There is a contradiction in the system’s energy after `fix evaporate`

Since I raised the issue of an anomaly in the system’s energy after the fix evaporate in the simulation last time, I still have not figured out the issue. Therefore, I am giving a more detailed post here again, and thanks Mr. Axel for his patient guidance.

Run in LAMMPS (2 Aug 2023 - Update 1), on CentOS with 64 cores.

The simulation I am conducting is of water evaporating from a Cu surface at room temperature, as shown in the figure.

The water molecules that drift from the surface to the top are considered evaporated, and I use the fix evaporate command to remove them. Before I proceed with the evaporation, I confirm that the system has already reached thermal equilibrium.

I conducted a simulation as follows: The total number of steps in the simulation was 15 million. In the first phase, which consisted of the initial 5 million steps, I did not use fix evaporate and simply carried out a normal fix nve preheat. Then, in the second phase, I used fix evaporate for 5 million steps to simulate evaporation. In the final phase, the last 5 million steps, I turned off fix evaporate and the entire system was fixed with nve. Normally, during the fix evaporate phase, both the kinetic and potential energy of the system should gradually decrease (due to the deletion of water molecules). However, this was not the case:

This is the Pe-steps curve.

This is the Ke-steps curve.

In the second phase, after using fix evaporate, the potential energy of the system decreased, but the total kinetic energy of the system did not change, which is a very contradictory result. I checked the temperature and found that the temperature did not change significantly throughout the process.

here is the atom-steps figure.

I am not clear about the reason for this contradiction, I think it is a LAMMPS bug for fix evaporate and here is my input:

units			metal
dimension		3
boundary		p p p
timestep		0.001
atom_style		full
neighbor		3.0 bin
neigh_modify	delay 0 every 1 check yes


kspace_style		pppm 1.0e-5

pair_style 		hybrid eam/fs lj/cut/coul/long 10
bond_style		zero
angle_style		zero

pair_coeff      * * eam/fs Cu1.eam.fs NULL NULL Cu
pair_coeff      1 1  lj/cut/coul/long 0.00659558  3.1507 # o
pair_coeff      2 2  lj/cut/coul/long 0.0   0 # h
pair_coeff      1 2  lj/cut/coul/long 0.0   0 #o-h
pair_coeff      1 3  lj/cut/coul/long 0.052520151	2.720475282 #Cu-o
pair_coeff      2 3  lj/cut/coul/long 0 0 #Cu-h
bond_coeff	1    0.9574 # o*-h*
angle_coeff 	1    104.52 # h*-o*-h*

region		fixed block EDGE EDGE EDGE EDGE EDGE -10 units box
region		evap_space block EDGE EDGE EDGE EDGE 90 EDGE
region          elf_space block EDGE EDGE EDGE EDGE 1 80

group			Water type 1 2 
group			fixed region fixed

set type 1 charge -0.834
set type 2 charge 0.417

compute      peng   all pe/atom
compute      patoms Water reduce sum c_peng
compute      keng   all ke/atom
compute      katoms Water reduce sum c_keng

variable     wateng equal "c_patoms+c_katoms"

compute      wtemp Water temp
compute      wtempcom Water temp/com

thermo			1000
thermo_style    custom step temp atoms v_wateng c_wtemp c_wtempcom c_patoms c_katoms etotal ke pe

fix 		freeze fixed setforce 0 0 0

fix		001 Water shake 0.0001 100 0 b 1 a 1
#fix             evp Water evaporate 10000 30 evap_space 34255 molecule yes
fix		fnvt Water nvt temp 295.35 295.35 0.1
run		5000000

fix             evp Water evaporate 20000 30 evap_space 34255 molecule yes
dump            1 all atom 2000
run 		5000000

unfix		evp
run		5000000

Could you share also the data file?

I think the jump in the potential energy after the end of evaporation is due to the fact that the evaporation was conducted using the nvt time integrator, and not the nve as you stated.

My guess is that as the number of atoms changes, the thermostat adds energy to the system but without accounting for the decreasing number of molecules. In the third run, that excess energy is removed and hence you see a discontinuity in the kinetic energy.

Could you try running the simulation with the following modifications?

compute_modify wtemp dynamic/dof yes
fix		001 Water shake 0.0001 100 0 b 1 a 1
fix		fnvt Water nvt temp 295.35 295.35 0.1
fix_modify fnvt temp wtemp
run		5000000

And keep the rest the same.


OK, here is my data file and Cu EAM potential
Cu1.eam.fs (738.7 KB)

data.out (1.6 MB)

And I will try this code to run the simmulation.

Thank you, If you hadn’t reminded me, I would have thought I was doing NVE integration all along. You’re right, in fact, the entire simulation is NVT integration.

this data file is wrong, here is the data file (4.5 MB)

Thank you. Let us know how the energy behaves with the new input file.

I experimented with a smaller sample of 512 TIP3P water molecules and found the simulation stuck after removing the first molecule. There is no crash or error message reported: LAMMPS keeps running in the background, but no I/O after the first evaporation event (see the Atoms column):

   Step         TotEng         E_vdwl         E_coul         E_long         E_bond        E_angle        Atoms        KinEng         PotEng          Temp          Press          Volume        Density          CPU  
     89800  -2460.6624      852.34101     -4489.8286      0              0              0                   1536   1176.8252     -3637.4876      385.92422      66.708391      32381.915      0.47299395     167.96439    
     90000  -4040.9463      888.30118     -6113.4338      0              0              0                   1533   1184.1863     -5225.1326      389.0989      -744.43469      32381.915      0.47207013     168.80912 

Here are the files to reproduce this issue. Tested on 7Feb24 (Update 1), 21Nov23, and 2Jun22. (78.3 KB) (2.8 KB)
water02.psf (147.4 KB)

PS The same issue happens with pair_style lj/cut/coul/long and pair_style lj/cut/coul/cut.

Hello Otello, I tried the parameters you provided and conducted a molecular simulation.

The results showed that everything became consistent. Thank you very much.

Here is the Pe

Here is the Ke

I made slight adjustments to the step numbers in each stage, and from the graph, it seems that there are no apparent issues. This parameter scheme is based on the NVT simulation. I would like to ask, if I want to use an NVE integrator, what modifications should I make to my parameters?

Thank you for sharing your results. If you want to carry out the evaporation step without thermostatting, you need the following:

# 1. Thermalisation
compute_modify wtemp dynamic/dof yes
fix		001 Water shake 0.0001 100 0 b 1 a 1
fix		fnvt Water nvt temp 295.35 295.35 $(100*dt)
fix_modify fnvt temp wtemp
run		5000000

# 2. Adiabatic evaporation.
unfix fnvt
fix		fnve Water nve
fix             evp Water evaporate 20000 30 evap_space 34255 molecule yes
dump            1 all atom 2000
run 		5000000

# 3. Post-evaporation equilibration .
unfix		evp
run		5000000

I expect the temperature to drop in the evaporation phase as it mimics an adiabatic volume increase.