# About simulating Water-vapor system.

Dear all,

I am trying to simulate the water vapor system by SPC/E model using lammps 15. But the fix-shake command is not working. I have installed rigid package in lammps 15 as shake command is a part of rigid package but still the system is not working. It would be great if anyone explain the error and suggest me a way to solve it. Thank you so much.

Following is the input script:

# model for water, Molecular Physics (2014), 112, 1566-1576

units real
atom_style full
include forcefield.SPC
pair_modify tail no

replicate 10 10 10

variable xLo equal -0.5lx
variable xHi equal 0.5
lx
variable yLo equal -0.5ly
variable yHi equal 0.5
ly
variable zLo equal -0.5lz
variable zHi equal 0.5
lz

# Recenter system coords about the origin (0,0,0):

displace_atoms all move {xLo} {yLo} \${zLo} units box

# Adjust box dimensions:

change_box all x final {xLo} {xHi} units box
change_box all y final {yLo} {yHi} units box
change_box all z final {zLo} {zHi} units box

# Increase z edge to yield “Vacuum|Water|Vacuum” system:

variable zLoNew equal -2.0lz
variable zHiNew equal 2.0
lz
change_box all z final {zLoNew} {zHiNew} units box

variable Nblock equal 10000
variable Nrun equal 5.0*{Nblock} variable Ndump equal {Nrun}/100
variable Ne equal 10
variable Nr equal {Nblock}/{Ne}
variable Nf equal \${Nrun}
variable Dz equal 0.1

variable A_in_m equal 1e-10 # Angstrom in meter
variable atm_in_Pa equal 101325 # note: 1 Pa = 1 N/m^2
variable N_in_mN equal 1e3 # Newton in milliNewton
variable convFac equal {A_in_m}*{atm_in_Pa}*\${N_in_mN}

variable Text equal 300

group hydrogen type 1
group oxygen type 2

velocity all create \${Text} 1234

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

timestep 0.01

#fix loadBalance all balance 20 z 5 1.05

fix constrain all 0.01 20 100 b 1 a 1
fix integrate all nvt temp {Text} {Text} 100
fix removeMomentum all momentum 1 linear 1 1 1

compute T all temp
fix TAve all ave/time {Ne} {Nr} \${Nf} c_T

variable P equal press
fix PAve all ave/time {Ne} {Nr} \${Nf} v_P

variable xPress equal c_thermo_press
variable yPress equal c_thermo_press
variable zPress equal c_thermo_press

# Evaluate and average surface tension in mN/m:

variable st equal 0.5lz(v_zPress-0.5*(v_xPress+v_yPress))*{convFac} fix st all ave/time {Ne} {Nr} {Nf} v_st

variable xyArea equal lx*ly

thermo_style custom step temp f_TAve press f_PAve v_xyArea lz f_st
thermo_modify flush yes
thermo \${Nf}

compute cO oxygen chunk/atom bin/1d z center {Dz} units box fix numDensO oxygen ave/chunk {Ne} {Nr} {Nf} cO density/number &
file numDensO.zProfile

compute cH hydrogen chunk/atom bin/1d z center {Dz} units box fix numDensH hydrogen ave/chunk {Ne} {Nr} {Nf} cH density/number &
file numDensH.zProfile

compute A all chunk/atom bin/1d z center \${Dz} units box

compute charges all property/atom q
fix qDens all ave/chunk {Ne} {Nr} \${Nf} A c_charges file qDens.zProfile

dump trj all atom \${Ndump} wat.lammpstrj
dump_modify trj scale no sort id

run \${Nrun}

#write_restart restart.wat