Dear all,
I am simulating a NPzT system. Water in confined in two solid walls, of which one is fixed and the other one acts as a piston to control the system pressure with the “fix aveforce” command. In my previous post, I followed Axel’s advices. I realized the temperature control of the piston wall. But a new problem turns up now. In the equilibration period, if I thermostat the solid atoms and water simultaneously, it works well, However, if I changed water into NVE ensemble while maintain langevin thermostat for two solid walls, water goes up to a extremely high temperature and expand quickly. The temperature of water is supposed to fluatute around 300K since the piston and substrate walls are kept at 300K. I feel frustrated for being unable to find the reason. It would be highly appreaciated if anyone could kindly point out my problem. The input script and the system configuration are attached below for your information.
Best regards,
Huaqiang
variable Ts equal 300.0
variable Ts_end equal 400.0
variable dt equal 2.0
variable spring_factor equal 67.6425
variable Tdamp equal 100dt
variable Pdamp equal 1000dt
---------- Initialize Simulation ---------------------
clear
echo screen
units real
dimension 3
boundary p p s
atom_style full
#atom_modify map array
bond_style harmonic
angle_style harmonic
#processors * * 1
---------- Create Atoms ---------------------
read_data system.data
---------- forcefield ---------------
variable epl_CuCu equal 4.72
variable sigma_CuCu equal 1.9297
variable epl_OO equal 0.1553
variable sigma_OO equal 3.166
variable epl_CuO equal sqrt({epl_CuCu}*{epl_OO})
variable sigma_CuO equal ({sigma_OO}+{sigma_CuCu})/2
variable epl_CuO2 equal 0.1*${epl_CuO}
mass 1 63.546 # Cu
mass 2 63.546 # Cu
mass 3 15.9994 # O
mass 4 1.00794 # H
pair_style hybrid lj/cut 10.0 lj/cut/coul/cut 9.0 9.0
pair_modify shift yes
#------------Piston-substrate interation--------------
pair_coeff 1 1 lj/cut {epl_CuCu} {sigma_CuCu}
pair_coeff 2 2 lj/cut {epl_CuCu} {sigma_CuCu}
pair_coeff 1 2 lj/cut {epl_CuCu} {sigma_CuCu}
#------------Water-water interation--------------
pair_coeff 3 3 lj/cut/coul/cut {epl_OO} {sigma_OO}
pair_coeff 3 4 lj/cut/coul/cut 0.0 0.0
pair_coeff 4 4 lj/cut/coul/cut 0.0 0.0
#------------Substrate-water interation--------------
pair_coeff 1 3 lj/cut {epl_CuO} {sigma_CuO} # Cu-O
pair_coeff 1 4 lj/cut 0.0 0.0 # Cu-H
#------------Piston-water interation--------------
pair_coeff 2 3 lj/cut {epl_CuO} {sigma_CuO} # Cu-O
pair_coeff 2 4 lj/cut 0.0 0.0 # Cu-H
bond_style harmonic
bond_coeff 1 450.0 1.0
angle_style harmonic
angle_coeff 1 55.0 109.47
---------- Group and Region Definition ---------------------
group substrate type 1
group piston type 2
group solid union substrate piston
group water type 3 4
region fixed block INF INF INF INF INF 1 units box
region phan block INF INF INF INF 1 6 units box
group fixed region fixed
group phan region phan
group gthermo subtract all fixed
group gnvtpre subtract all water fixed
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes exclude group fixed fixed
---------- Run Minimization ---------------------
fix spring solid spring/self ${spring_factor}
fix freeze fixed setforce 0 0 0 #fix the atoms at bottom
fix pistonkeep piston setforce 0.0 0.0 0.0
thermo 1000
thermo_style custom step pe press pxx pyy pzz
min_style cg
minimize 1e-6 1e-8 5000 10000
write_data system_min.data
reset_timestep 0
---------- Run 1st Equilibration ---------------------
unfix spring
unfix pistonkeep
velocity all create ${Ts} 13425 mom yes rot yes dist gaussian loop local
velocity fixed set 0.0 0.0 0.0 sum no units box
velocity piston set NULL NULL 0.0 sum no units box
fix springsub substrate spring/self {spring_factor} fix springpis piston spring/self {spring_factor} xy
fix fShakeSPCE water shake 0.0001 10 0 b 1 a 1
dump mydump1 all custom 10000 equil.*.lammpstrj id type x y z vx vy vz
dump_modify mydump1 sort id
fix loadbalance all balance 10000 1.0 shift xyz 10 1.1 weight time 0.8 out balance.output
fix fNVEall gthermo nve
group gequil union substrate water
fix fNVTlan1 gthermo langevin {Ts} {Ts} ${Tdamp} 54687 tally yes
variable pressure equal 1.0 # atmosphere pressure, 1.01325 bar
variable num_piston equal “count(piston)”
variable force equal -{pressure}*lx*ly*1.4579e-5 #1 kcal/(mol A3) = 6.95e9 Pa = 68591.167 atmospheres, 1 atmospheres = 1.4579e-5 kcal/(mol A3)
#fix pistonkeep piston setforce 0.0 0.0 NULL
fix conspress1 piston aveforce NULL NULL {force}
timestep ${dt}
run 100000
---------- Run 2nd Equilibration ---------------------
unfix fNVTlan1
unfix conspress1
fix fNVTlan2 substrate langevin {Ts} {Ts} {Tdamp} 5734 tally yes
fix fNVTlan3 piston langevin {Ts} {Ts} {Tdamp} 5734 tally yes
fix conspress2 piston aveforce NULL NULL ${force}
run 100000
write_data system_equil.data
write_restart equil.restart
#reset_timestep 0
----------------- Main NEMD simulation: Temperature ramp in the phantom layers -----------------
undump mydump1
dump mydump all custom 10000 run.*.lammpstrj id type x y z vx vy vz
dump_modify mydump sort id
unfix fNVTlan2
unfix fNVTlan3
unfix conspress2
fix tempramp phan langevin {Ts_end} {Ts_end} {Tdamp} 4464 tally yes
fix fNVTlan4 piston langevin {Ts} {Ts} {Tdamp} 5734 tally yes
fix conspress3 piston aveforce NULL NULL ${force}
