LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] About simulating Water-vapor system.
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] About simulating Water-vapor system.

From: Rubayat Bin Shahadat <rubayat37@...24...>
Date: Thu, 12 Apr 2018 20:38:17 +0600

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.5*lx
variable xHi equal 0.5*lx
variable yLo equal -0.5*ly
variable yHi equal 0.5*ly
variable zLo equal -0.5*lz
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.0*lz
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[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.5*lz*(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

Muhammad Rubayat Bin Shahadat
Dept. of Mechanical Engineering
Hajee Mohammad Danesh Science and Technology University, Dinajpur.