Dear LAMMPS community,
I am unsure if this question relates to LAMMPS or general Physics. I have two thermostatted groups, one of each is gold (“hot region”) and the other one contains water (“cold region”). The “hot region” is maintained at 350 K while the “cold region” is maintained at 300K with Langevin thermostat. I have no overall energy drift in the system, but the absolute values of energies deposited by the thermostats into corresponding groups differ by about 60%. If this energy goes into, let’s say, rotation of the system as a whole, shouldn’t it be reflected in the total energy change? Right now it looks like the thermostats should be pumping energy into the system based on their net cumulative energy change, but this doesn’t seem to be reflected in the total energy. I appreciate any suggestions!
The computational setup is a spherical shell of SPC/E water surrounded by a 12-6 Lennard-Jones wall (wall/region with lj126). The “cold region” is ~ 2.3 nm thick and is shown by the orange color in the snapshot below. The snapshot is a slice of the whole system. The gold atoms are shown in yellow and the water atoms are shown in red and white. The graph demonstrates the total energy change of the system (Etot) relative to t = 0, the cumulative energy change for the hot thermostat (f_fHi), and the cumulative energy change for the cold thermostat (f_fLo).
It would have been helpful to see your input and the specific commands you were using (e.g. as a thermostat and which settings).
It should not. Please see the following statement from the fix nvt documentation:
The cumulative energy change in the system imposed by these fixes, via either thermostatting and/or barostatting, is included in the thermodynamic output keywords ecouple and econserve. See the thermo_style page for details.
So etotal does not account for the energy that was added to or removed from the system by the thermostat, but econserve does. While the cold part has to remove energy, the hot part has to add it and the sum of potential and kinetic energy of the whole system should be constant when the system is in (dynamic) equilibrium.
Dear Axel,
Thank you for your response!
My input script is below. I’ve also attached everything necessary to run the script in the .zip file, which can be accessed with the link:
#atom type 1 = O
#atom type 2 = H
#atom type 3 = Au
log log.out
variable Nprod equal 60000
variable Thi equal 350.0
variable Tlo equal 300.0
units metal
variable epsOAu equal 0.025576102438521528 #eV
variable epsOO equal 0.006737376507851235#eV
variable sigmaOAu equal 3.6 # A
variable epsHAu equal 0.0
variable sigmaHAu equal 0.0
variable cutLJ equal 9.0 #A
variable cutLJ equal 9.0 #A
variable kB equal 8.617333262e-5 #eV/K
variable N_dof_H equal 1.5947
variable N_dof_O equal 2.8106
atom_style full
boundary p p p
timestep 0.003
pair_style lj/cut/coul/long 9.0 9.0
special_bonds lj/coul 0.0 0.0 0.5
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
read_data spherical_cp_350K_thermalized.dat
pair_style hybrid lj/cut/coul/long 9.0 lj/cut 9.0 eam/alloy
pair_coeff 1 1 lj/cut/coul/long ${epsOO} 3.16555789 #O-O
pair_coeff 1 2 lj/cut/coul/long 0.0000 0.0000 # O-H
pair_coeff 2 2 lj/cut/coul/long 0.0000 0.0000 #H-H
pair_coeff 1 3 lj/cut ${epsOAu} 3.6# Au-O
pair_coeff 2 3 lj/cut ${epsHAu} 0.0# Au-H
pair_coeff * * eam/alloy Au-W1266.v1.eam.alloy NULL NULL Au #Au - Au
pair_modify tail no
kspace_style pppm 1.0e-6
bond_coeff 1 500000.00 1.000
angle_coeff 1 500000.0 109.47
neighbor 4.5 bin
group water type 1 2
group O type 1
group H type 2
group gold type 3
compute temp_water water temp
compute temp_gold gold temp
compute tKe all ke
compute tPe all pe
compute ke_H2O water ke/atom
compute ke_Au gold ke/atom
compute pe_H2O water pe/atom
compute pe_Au gold pe/atom
compute tKeAu gold reduce sum c_ke_Au
compute tPeAu gold reduce sum c_pe_Au
compute tKeH2O water reduce sum c_ke_H2O
compute tPeH2O water reduce sum c_pe_H2O
variable vE equal c_tKe+c_tPe
variable vEAu equal c_tKeAu+c_tPeAu
variable vEH2O equal c_tKeH2O+c_tPeH2O
variable time_cur equal time
variable temp_cur equal temp
variable temp_H2O equal c_temp_water
variable temp_gold equal c_temp_gold
variable Et equal etotal/atoms
variable Ep equal epair
variable Nevery equal 10
variable Nrepeat equal 25
variable Nfreq equal ${Nevery}*${Nrepeat}
variable sys_centx equal (xhi+xlo)/2.0
variable sys_centy equal (yhi+ylo)/2.0
variable sys_centz equal (zhi+zlo)/2.0+0.15
region r_Lo sphere ${sys_centx} ${sys_centy} ${sys_centz} 100.0 side out
run 0
group Tlo_group dynamic water region r_Lo every 10
run 0
fix fLo Tlo_group langevin ${Tlo} ${Tlo} 0.5 4838283 tally yes
fix fHi gold langevin ${Thi} ${Thi} 0.1 3272022 tally yes
thermo_style custom temp ke pe etotal v_temp_H2O v_temp_gold
thermo 500
restart 50000 system.restart
fix prt all ave/time ${Nevery} ${Nrepeat} ${Nfreq} v_temp_cur v_temp_H2O v_temp_gold v_Ep v_Et file ther_params.out
region r_wall sphere ${sys_centx} ${sys_centy} ${sys_centz} 123.475 side in
fix fnve all nve
fix 1 water shake 0.0001 20 0 b 1 a 1
fix fwall all wall/region r_wall lj126 ${epsOO} 3.1655 7.3
fix fE all ave/time 1 1 250 v_vE v_vEAu v_vEH2O f_fLo f_fHi f_fwall file out.energies
run ${Nprod}
write_data spherical_cp_350K_continued.dat nofix
