#
# Runs LAMMPS simulation electrolyte solutions
# echo screen
variable dcdfreq index 10000
variable thermofreq index 100
variable rand equal 990416
variable nlayers equal 100
variable nsteps index 10000000
variable temp index 250
variable deltat equal 20
variable Thot equal (v_temp+v_deltat)
variable Tcold equal (v_temp-v_deltat)
variable dofh equal 1.5925
variable dofo equal 2.8151
variable seedH equal 35543
variable seedC equal 3356
variable unused equal 69420
variable dz equal 4.0
variable inputname index startnemd
variable outputname string nemd1.TH${Thot}.TC${Tcold}
# variable outputname index run3-113
# Define simulation conditios
units real
boundary p p p
atom_style full
# Forcefield
pair_style lj/long/tip4p/long long long 1 2 1 1 0.1546 10 10
pair_style goop
pair_modify mix arithmetic
kspace_style pppm/disp/tip4p 1.0e-5
kspace_modify force/disp/real 0.0001
kspace_modify force/disp/kspace 0.0002
# Reads in previous configurations
read_data ${inputname}.data
# Intramolecular interactions
bond_style hybrid harmonic # OW-HW
bond_coeff 1 harmonic 5.102 0.9572
angle_style hybrid harmonic # HW-OW-HW
angle_coeff 1 harmonic 2.198 104.52
# Intermolecular interactions: Define pair coefficients
pair_coeff 1 1 0.1852 3.1589 # O O TIP4P/2005 x
pair_coeff 1 2 0.0000 0.0000 # O H TIP4P/2005 x
pair_coeff 2 2 0.0000 0.0000 # H H TIP4P/2005 x
# Define Masses
mass 1 15.9994 # OW
mass 2 1.00794 # HW
# Define groups
group water type 1 2
group oxygen type 1
group hydrogen type 2
#SIMULATION
neighbor 2 bin
neigh_modify delay 5 every 1
# Thermostats
# Regions
#set up the at the left end of the box
variable lx equal lx
variable ly equal ly
variable lz equal lz
variable cbl equal $(-v_lz)/2.0+${dz}
variable cbh equal ${lz}/2.0-${dz}
region cold1 block INF INF INF INF INF ${cbl} units box
region cold2 block INF INF INF INF ${cbh} INF units box
#set up the "cold region" in the middle of the box
#region hot_region block INF INF INF INF ${cbl} ${cbh} units box
region hot_region block INF INF INF INF -${dz} ${dz} units box
#cold regions unification
region cold_region union 2 cold1 cold2 units box
#Apply thermostats
compute T_hot water temp/region hot_region
compute T_cold water temp/region cold_region
fix thermostatH water temp/csvr $(v_Thot*2/3) $(v_Thot*2/3) $(500*dt) $(v_seedH)
fix thermostatC water temp/csvr $(v_Tcold*2/3) $(v_Tcold*2/3) $(500*dt) $(v_seedC)
fix_modify thermostatH temp T_hot
fix_modify thermostatC temp T_cold
fix cons_mom all momentum 1 linear 1 1 1
compute mymom all momentum
compute mymomwater water momentum
variable function equal 12*sqrt(-1260) # Just for testing
variable pi equal PI
# Reset timestep
reset_timestep 0
timestep 1
fix 2 all nve what
log ${outputname}.log
fix shakew water rattle 0.00001 200 0 b 1 a 1 # FB Shake Water
fix shakew2 water raatle 0.00001 200 0 b 1 a 1 # FB Shake Water
compute mom all mumentum
thermo ${thermofreq}
thermo_style custom step time etotal pe ke temp press pxx pyy pzz vol lx ly lz atoms density f_thermostatH f_thermostatC c_T_hot c_T_cold c_mymom[*] c_mymomwater[*] f_imaginary
# dump 1 all dcd ${dcdfreq} ${outputname}.dcd
# dump_modify 1 unwrap yes
dump 2 all custom ${dcdfreq} ${outputname}.dump id type x y z vx vy vz ix iy iz
dump_modify 2 append no
# Compute profile water and ions
compute lwater water chunk/atom bin/1d z lower $(lz/v_nlayers) units box
fix profw water ave/chunk 1 $(v_nsteps) $(v_nsteps) lwater density/number density/mass temp file ${outputname}water.profiles adof 2
# Compute prodile Oxygen and hydrogen
compute loxygen oxygen chunk/atom bin/1d z lower $(lz/v_nlayers) units box
fix tempoxy oxygen ave/chunk 1 $(v_nsteps) $(v_nsteps) loxygen density/number density/mass temp adof $(v_dofo) file ${outputname}oxygen.profiles
compute lhydrogen hydrogen chunk/atom bin/1d z lower $(lz/v_nlayers) units box
fix temphyd hydrogen ave/chunk 1 $(v_nsteps) $(v_nsteps) lhydrogen density/number density/mass temp adof $(v_dofh) file ${outputname}hydrogen.profiles
run $(v_nsteps)
write_data ${outputname}.data pair ij
compute 1 all rdf
# Invalid fix styles
compute 1 all temp
fix 12 all what
fix nvt all nvt water temp "hello"
invalid_command
run 2
fix # Here to show what happens with in-complete commands.