############################################################################### # Description: # ------------ # # This file is a LAMMPS input script for equilibrating TIP4P/ice water # with an initial lattice structure at a temperature of # 300 K. The output of the simulation are two files. # The file "avg_nve_post_heat.dat" contains the average temperature and # and total system energy, respectively, of the last NVE simulation. # The binary restart file "nve.res" contains the last configuration # of that run which serves as an input for the NEMD simulation. # # Dependencies: # ------------- # # The input file waterandgraphite which contains more than 4000 TIP4P/ice molecules above the graphite. # ############################################################################### ################################################# # Variable declaracations ################################################# # number of timesteps for the two NVE runs: # - NNVEpre: before correcting the energy # - NNVEpost: after correcting the energy variable NNVEpre equal 1000 variable NNVEpost equal 200000 # number of timesteps for the short velocity rescaling run variable Nrescale equal 10000 # number of timesteps for the NVT simulations variable NNVT equal 1000000 # number of timesteps for correcting the energy variable Nadd equal 50 # sampling frequency for temperature and energy variable Nsamp equal 10 # seed for generating random initial velocities variable seed equal 1001 # timestep variable dt equal 1. # reference temperature variable Tref equal 300. variable Tice equal 220. # damping parameter for Nosé-Hoover thermostat variable Tdamp equal 1000 # parameters for the TIP4P/ice model variable epsOO equal 0.21084 variable sigOO equal 3.1668 variable theta equal 104.52 # cutoff radii for Coulomb and LJ interactions, respectively variable cutC equal 8.5 variable cutLJ equal 11 # specification of units, spatial dimensions, boundary conditions and atom-style units real dimension 3 boundary p p p neighbor 3.0 bin neigh_modify delay 0 every 1 check yes atom_style full read_data water_on_graphite.data # group atoms to molecules group O type 1 group H type 2 group water type 1 2 group graphite type 3 # specify pairwise interactions with long-range Ewald solver pair_style hybrid lj/cut 12 lj/cut/tip4p/long 1 2 1 1 0.1577 ${cutLJ} ${cutC} pair_coeff * * lj/cut/tip4p/long 0 0 pair_coeff 1 1 lj/cut/tip4p/long ${epsOO} ${sigOO} pair_coeff 3 3 lj/cut 0.0663 3.5812 14 pair_coeff 1 3 lj/cut 0.13 3.2 12 kspace_style pppm/tip4p 1.e-5 # specify harmonic intra-molecular interactions # NOTE: this will have no effect as we use SHAKE bond_style harmonic angle_style harmonic dihedral_style opls improper_style harmonic bond_coeff 1 0 0.9572 bond_coeff 2 610 1.42028 angle_coeff 1 0 ${theta} angle_coeff 2 80.0 120.0 dihedral_coeff 1 0.0 7.25 0.0 0.0 improper_coeff 1 5.0 180.0 # use SHAKE with a precision of 1e10. # NOTE: We use SHAKE instead of RATTLE for the equilibration run, # because the LAMMPS implementation of RATTLE is not fully # compatible with fix nvt. # (The constraints are only satisfied exactly with fix nve.) fix 1 water shake 1e-10 400 0 b 1 a 1 # use standard correction terms for the trucated tail of the LJ potential pair_modify tail yes # compute the temperature, kinetic energy, potential energy and # total energy, respectively compute cT water temp compute cKe water ke compute cPe all pe variable E equal c_cKe+c_cPe dump dp1 all atom 1000 dump.lammpstrj ################################################# # Short run with velocity rescaling # to drive the system close to the target # temperature. ################################################# reset_timestep 0 fix fNVE1 all nve # rescale velocities every 10 timesteps fix 2 graphite setforce 0 0 0 fix fHEATING water temp/rescale 10 ${Tref} ${Tref} 0 1 timestep ${dt} # output timestep, temperature, energies and pressure to console every 100 steps thermo 100 thermo_style custom step temp ke pe etotal press thermo_modify temp cT run ${Nrescale} unfix fHEATING unfix fNVE1 ################################################# # NVT simulation with Nosé-Hover thermostat ################################################# reset_timestep 0 fix fNVT water nvt temp ${Tref} ${Tref} ${Tdamp} timestep ${dt} # average temperature and total energy during that run variable repeatNVT equal ${NNVT}/${Nsamp} fix fAvgNVT water ave/time ${Nsamp} ${repeatNVT} ${NNVT} c_cT v_E run ${NNVT} unfix fNVT # Note: use $() here fore immediate evaluation as the fix will be deleted afterwards variable Tavg_NVT equal $(f_fAvgNVT[1]) variable Eavg_NVT equal $(f_fAvgNVT[2]) unfix fAvgNVT ################################################# # NVE simulation ################################################# reset_timestep 0 variable repeatNVEpre equal ${NNVEpre}/${Nsamp} fix fAvgNVEpre water ave/time ${Nsamp} ${repeatNVEpre} ${NNVEpre} c_cT v_E fix fNVEpre water nve run ${NNVEpre} # Note: use $ here fore immediate evaluation as the fix will be deleted afterwards variable Tavg_NVE equal $(f_fAvgNVEpre[1]) variable Eavg_NVE equal $(f_fAvgNVEpre[2]) unfix fAvgNVEpre unfix fNVEpre # output averages to console print "AVERAGES:" print "_NVT = ${Tavg_NVT}" print "_NVT = ${Eavg_NVT}" print "_NVE = ${Tavg_NVE}" print "_NVE = ${Eavg_NVE}" ################################################# # Adjust the energy of the last configuration ################################################# variable dQ equal ${Eavg_NVT}-${Eavg_NVE} print "dQ = ${dQ}" # define box as region variable xlo equal xlo variable ylo equal ylo variable xhi equal xhi variable yhi equal yhi variable zlo equal zlo variable zhi equal zhi region boxreg block ${xlo} ${xhi} ${ylo} ${yhi} ${zlo} ${zhi} # calculate rate variable dQdt equal ${dQ}/${Nadd}/${dt} # add heat fix fHEAT water ehex 1 ${dQdt} region boxreg fix fNVEADD water nve run ${Nadd} unfix fNVEADD unfix fHEAT ################################################# # Second NVE simulation with correct temperature ################################################# reset_timestep 0 fix fNVEpost water nve variable repeatNVEpost equal ${NNVEpost}/${Nsamp} fix fAvgNVEpost water ave/time ${Nsamp} ${repeatNVEpost} ${NNVEpost} c_cT v_E fix fAvgNVEpostFile water ave/time ${Nsamp} 1 ${Nsamp} c_cT v_E file avg_nve_post_heat.dat run ${NNVEpost} write_restart "nve.res" variable Tavg_NVE3 equal $(f_fAvgNVEpost[1]) variable Eavg_NVE3 equal $(f_fAvgNVEpost[2]) unfix fAvgNVEpost unfix fAvgNVEpostFile unfix fNVEpost print "AVERAGES:" print "_NVE = ${Tavg_NVE3}" print "_NVE = ${Eavg_NVE3}" ################################################# # Short run with velocity rescaling # to drive the system close to the target ice # temperature 220K. ################################################# reset_timestep 0 fix fNVE1 water nve # rescale velocities every 10 timesteps fix fHEATING water temp/rescale 10 220 220 0 1 timestep 2 # output timestep, temperature, energies and pressure to console every 100 steps thermo 100 thermo_style custom step temp ke pe etotal press thermo_modify temp cT run ${Nrescale} unfix fHEATING unfix fNVE1 unfix 2 ################################################# # NVT simulation with Nosé-Hover thermostat ################################################# reset_timestep 0 fix fNVT water nvt temp 220 220 ${Tdamp} timestep 2 # average temperature and total energy during that run variable repeatNVT equal ${NNVT}/${Nsamp} fix fAvgNVT all ave/time ${Nsamp} ${repeatNVT} ${NNVT} c_cT v_E fix 2 graphite setforce 0 0 0 run 600000000 unfix fNVT # Note: use $() here fore immediate evaluation as the fix will be deleted afterwards variable Tavg_NVT equal $(f_fAvgNVT[1]) variable Eavg_NVT equal $(f_fAvgNVT[2]) unfix fAvgNVT