nonconservation of energy with fix ehex and fix nve command in n-heptane/benzene mixture

LAMMPS LAMMPS Mailing List Mirror
1

Dear all，
When I tried to simulate the NEMD processes with fix ehex and fix nve command, the temperature and energy of the whole system are improving with simulation time. In my mixture, n-heptane is set as non-rigid molecule and benzene is set as rigid molecule, each of them has 500 molecules and the mixture has 1000 molecules totally. Both of the n-heptane and benzene use Trappe-UA force field, and the particles has no charge. Therefore, I did not use kspace command in the simulation. When the system become stable in 303.15 K and 1 atm, I tried to used fix ehex command to make a hot region in the middle of the box and two cold regions in the boundaries. The benzene use fix rigid/nve/small command and n-heptane use fix nve command to make the whole system perform nve dynamics. However, I found that the temperature and energy of the whole system are improving. I have tried to test the energy conservation in nve ensemble without ehex algorithm when the simulation timestep set as 2 fs and 1 fs, both simulations are in good energy conservation. But when I add the fix ehex command, the energy is improving both when the timestep set as 2 fs and 1 fs. Could you give me some advice to avoid this problem? here is my input file and log file. I also plot the drawing of energy vs timestep.

########################################################################## input file

#restart file for benzene-n-heptane (0.5) model with number of 1000

log 0.5-nc7-3-nve.log
echo both
read_restart restart.1667000.nvt

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes exclude molecule/intra benzene
reset_timestep 0

################################################ test NVE (timestep=2fs)
timestep 2
fix 1 benzene rigid/nve/small molecule
fix 2 heptane nve

thermo 1000
thermo_style custom step temp ke pe etotal press

run 5000000
############################################### test NVE (timestep=1fs)
timestep 1
run 5000000

unfix 1
unfix 2

############################################### fix ehex process
reset_timestep 0

restart 2000000 restart.*.ehex1

variable delta equal (zhi-zlo)/20
variable zhi_T1 equal zhi-(zhi-zlo)/2+{delta} variable zlo_T1 equal zhi-(zhi-zlo)/2-{delta}
variable zhi_T21 equal zlo+{delta} variable zlo_T21 equal zlo variable zhi_T22 equal zhi variable zlo_T22 equal zhi-{delta}

region T1_region block INF INF INF INF {zlo_T1} {zhi_T1}
region T21_region block INF INF INF INF {zlo_T21} {zhi_T21}
region T22_region block INF INF INF INF {zlo_T22} {zhi_T22}

fix 1 benzene rigid/nve/small molecule
fix 2 heptane nve

fix T1_fix all ehex 1 0.008 region T1_region
fix T21_fix all ehex 1 -0.004 region T21_region
fix T22_fix all ehex 1 -0.004 region T22_region

thermo 1000
thermo_style custom step temp ke pe etotal c_T1 c_T21 c_T22
run 1000000

截图1.png

截图2.png

2

Hi Xiaoyu,

The current implementation of fix ehex assumes that you use it in combination with either fix shake or fix rattle when applying it to rigid bodies (see doc pages). If using fix rattle is an option for the system you study, then I’d suggest you try that in combination with fix ehex.

Best wishes,
Peter

截图1.png

截图2.png

3

Dear Peter,

Thank you for your reply! In my NEMD simulation with “fix ehex”, I only used “fix rigid/nve/small” command to make the benzene molecules being in rigid bodies. Do you mean I should use “fix rattle” to constrain the bonds and angles of the Benzene molecules with “fix ehex” command? If I use “fix rattle” and “fix ehex”, should I use “fix nve” to replace the “fix rigid/nve/small” to the all particles in the mixture to make the energy conserved? by the way, the force field I used is TraPPE-UA, there are only 6 sites in a benzene molecule and each site is (CH2). Is “fix rattle” used to constrain the C-C bonds and C-C-C angles in the benzene molecule?

Yours sincerely,
Xiaoyu Chen

截图1.png

截图2.png

4

Hi Xiaoyu,

The operator splitting analysis underlying fix ehex is very specific to velocity Verlet integration (standard fix nve). I have never tested it in combination with other integration schemes, such as fix rigid/nve/small. If it’s possible to constrain molecules with 6 sites using fix rattle or fix shake, then I would try that next. The doc pages probably talk about the limitations. If it’s not possible, then I’d explore other thermostatting options.

Best wishes,
Peter

截图1.png

截图2.png

5

Dear Peter,

Thanks so much for your reply. Last time after your email, I have tried to use “fix ehex” to the n-heptane molecules only and the simulation works well. I think it may have the same pattern of the NPT dynamics of the mixture with rigid bodies and non-rigid particles in lammps. I hope it can provide a new idea about how to deal with the mixture with some rigid bodies in a NEMD process using ehex algorithm.

Yours sincerely,

Xiaoyu Chen