Unphysical Temperature Rise in Water Flash Evaporation Simulation with LAMMPS

LAMMPS LAMMPS General Discussion
1

Dear LAMMPS Community,

I am writing to seek advice regarding an unphysical issue I have encountered in my recent molecular dynamics simulations of flash evaporation.

Over the past few months, I have been studying flash evaporation processes using LAMMPS. My current focus is on the flash evaporation of water molecules. The simulation setup is as follows:

I created a cubic box of liquid water.

The system was first relaxed in the NVT ensemble at 360 K until equilibrium was reached.

After equilibration, I switched to the NVE ensemble and suddenly expanded the simulation box in the z-direction to simulate rapid decompression, aiming to induce flash evaporation.

According to the principles of flash evaporation, the temperature of the system should decrease due to the latent heat absorbed during the phase transition. However, in my simulations, I observed a continuous increase in temperature, which is clearly unphysical.

I have double-checked my input script and parameters, including the intermolecular potentials (I am using the SPC/E water model) and the decompression rate, but the issue persists. I suspect there might be something overlooked in the ensemble transition, the box expansion method, or the energy conservation settings in NVE.

Has anyone encountered a similar problem or have insights into what might be causing this behavior? Any suggestions on how to resolve this issue would be greatly appreciated.

Thank you in advance for your time and assistance. I look forward to your valuable feedback.

Best regards,

Yu Xing

2

I think the problem here is not the behavior of the system, but your method of simulating the evaporation. By expanding the box (assuming that you use the change_box command applied to the group all) you also scale the atom positions. That changes your system from a well equilibrated system to a system, where all atoms instantly have a significant change in their environment and potential energy and are far from equilibrium. For a rigid water potential you have the additional complication that your constraints are suddenly violated by that scaling, too. Even if you created an empty group, you would still have an unphysical instant change in the environment of the atoms where you “cut” the box due to periodic boundaries.

Have you surveyed the published literature to see what others have done to simulate this process?

This is not my area of expertise, but common sense tells me that to see what you want to observe, you would have to use fix deform and deform the system slowly so that the molecules can stay in equilibrium. Even if that would require a “very long” simulation as far as atomic scale simulations go, it would still be quick compared to macroscopic experiments (i.e. nano- or microseconds).

3

Try this experiment out:

Suppose I have a LAMMPS data file for a fully equilibrated 1nm cube of fully periodic water (and all the atoms are inside the cube – assume we don’t have any issues with bonds crossing periodic boundaries).

Now suppose I edit the data file and manually change the z-length from 10 angstroms to 20 angstroms, and I still run a fully periodic simulation. What would this new system represent? It would represent 1nm-thick slabs of water alternating with 1nm-thick slabs of vacuum.

What is the total potential energy per simulation box? It must be much higher than before, because each simulation box now contains two 1-nm^2 water-vacuum interfaces under significant surface tension. (You can confirm this by outputting the potential energy with the thermo output from a run 0 command.)

If the potential energy supplied externally through a discontinuous configurational change is enough to “pay for” the latent heat of evaporation, then it makes sense that you would see the kinetic temperature increase instead of decrease. This is the same point @akohlmey was making but with a much more graphic example (I am not saying your script actually does this – but it is likely doing something similar, as he said.)

1 Like
4

Dear Professor,
Thank you so much for taking the time to respond to my questions and for your valuable advice. I am truly honored to receive your reply.
I must admit that my familiarity with the literature is still quite limited, and my understanding of molecular dynamics remains incomplete. When faced with a challenge like this, I realize my current ability is not yet sufficient to resolve it on my own.
Following your suggestions, I have actively tried to improve my approach. Specifically, I switched from using “change box” to “fix deform.” Due to the presence of long-range electrostatic interactions handled by PPPM, I also implemented a reflecting wall in the z-direction that adjusts dynamically along with the deformation.
Unfortunately, the issue of rising temperature continues to persist. After running the simulation for two full days—now reaching 6,727,000 steps—the length in the z-direction has expanded from 39.9 Å to 154.46 Å, and the temperature has increased from 360 K to 758 K. Given these results, I feel I must seek further guidance in finding a new solution.
Once again, thank you for your patience and support. I sincerely hope to learn from your expertise and would be deeply grateful for any further suggestions you might offer. The files I have prepared are attached below for your reference.

Respectfully yours,
Yu Xing

Basic settings

units real
dimension 3
boundary p p f
atom_style full
timestep 2.0
#-----

Neighbor list settings

neighbor 2 bin
neigh_modify delay 0 every 1 check yes
#-----

Simulation box

region whole block -18.95 20.95 -18.95 20.95 -18.95 20.95
region liquid block -17.95 19.95 -17.95 19.95 -17.95 19.95

create_box 2 whole bond/types 1 angle/types 1 extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
#-----

Atom types

mass 1 15.9994 # O
mass 2 1.008 # H
#-----

Water molecule definition

molecule H2O spce.mol
#-----

create atoms

lattice bcc 3.788
create_atoms 0 region liquid mol H2O 1000

write_data water.data
#-----

Group definitions

group liquid region liquid
#-----

Force field settings

SPC/E water model

pair_style lj/cut/coul/long 9.0
pair_coeff 1 1 0.1553 3.166 # O-O
pair_coeff 1 2 0.0 0.0 # O-H
pair_coeff 2 2 0.0 0.0 # H-H

Bond and angle parameters

bond_style harmonic
bond_coeff 1 1000.0 1.0 # O-H bond

angle_style harmonic
angle_coeff 1 100.0 109.47 #1=Hw-Ow-Hw

Long-range electrostatics

kspace_style pppm 1.0e-4
kspace_modify slab 3.0
#-----
fix rigid all shake 0.0001 20 0 b 1 a 1
fix zwalls all wall/reflect zhi EDGE zlo EDGE units box
#-----

Energy minimization

min_style cg
minimize 1.0e-1 1.0e-3 100 1000
reset_timestep 0

===

Initial velocities

velocity all create 360 12345 rot yes dist gaussian

NVT integration

fix nvt liquid nvt temp 360 360 100

Thermodynamic output for NVT

thermo 1000
thermo_style custom step press pe temp ke density lx ly lz

Trajectory output for NVT

dump nvt all custom 1000 nvt.lammpstrj id type x y z vx vy vz
dump_modify nvt sort id

Run NVT

run 1000000

write_restart restart.equil

variable target_lz equal 1000
variable rate equal 0.00001
variable interval equal 1000
variable loops equal 48000
#-----
fix deform all deform 1 z vel 0.00001 units box remap v
fix walls all wall/reflect zhi EDGE zlo EDGE units box
fix nve all nve
#-----
thermo 1000
thermo_style custom step time temp press density lz
dump flash all custom 1000 flash.lammpstrj id type x y z
dump_modify flash sort id
#-----
print “Starting main loop…”
variable i loop {loops} label main_loop run {interval}
#-----
unfix walls
fix walls all wall/reflect zhi EDGE zlo EDGE units box
next i
jump SELF main_loop
#-----
unfix deform
run 200000
write_data flash_final.data
write_restart restart.flash
print “Simulation completed successfully!”

5

Dear Professor,
Thank you very much for your thoughtful reply and for raising the important point about the potential energy increase exceeding the latent heat of vaporization.
I sincerely appreciate your insights, as I have been earnestly trying to achieve the physical phenomenon of flash evaporation through various simulation approaches, though without success thus far.
You are absolutely right to point out the consequences of artificially expanding the simulation box. In my original setup with a water box, when I abruptly change the box dimensions, the log files indeed show a sudden decrease in pressure and density, while the kinetic and potential energies remain relatively unchanged in the short term.
My goal has been to achieve flash evaporation through rapid depressurization, but I realize my current understanding may be insufficient. I will carefully reconsider the implications of this box expansion method and explore alternative approaches to better address this issue.
I am aware that there may be additional factors I need to consider, and I would be grateful for any further suggestions or references you might recommend on this topic.
Thank you again for your valuable guidance.

Best regards,
YuXing

6

I am sorry, but I don’t see a question about LAMMPS here. What you are asking about is advice on how to do your research. That is outside the scope of this forum and not something that I (or others) have the time to assist you with. That is a topic for a discussion with your adviser or supervisor, or a suitable collaborator that has the sufficient expertise (I don’t. I already told you that this kind of research is not my area of expertise and that my response was based simply on common sense).

If you want, we can re-categorize this discussion to Science Talk - Materials Science Community Discourse

7

Hi YuXing,

The only additional advice I can offer is to carefully consider what the thermodynamic conditions are for the process you intend to simulate, and see how other molecular dynamics studies have accounted for them. In particular, it sounds like you are trying to simulate an isenthalpic process, and there are references online to that kind of simulation.

As @akohlmey has said, this would not really be a discussion about the specific use of LAMMPS in this simulation and so it might be a bit difficult for you to get additional help on this thread.