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:

Reference: M. Orsi, Comparative assessment of the ELBA coarse-grained

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

units real
atom_style full
read_data data.singleSPC
include forcefield.SPC
pair_modify tail no

replicate 10 10 10

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

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
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[1]
variable yPress equal c_thermo_press[2]
variable zPress equal c_thermo_press[3]

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