Equilibrium of water

Dear LAMMPS users,

Hello! I’m a new user to LAMMPS and a beginner of MD. Recently, I’m assigned a job to simulate the alignment of water molecules under laser electric field. I start with the process of equilibrium of water molecules. I searched related input scripts for a while and found Dr. Mario Orsi has made a simulation of bulk water. The address is http://www.orsi.sems.qmul.ac.uk/downloads.html. I slightly modify his code with changing the amount of water molecules to 256 with fcc arrangement. My input file is attached at the end of the post. Water model is TIP4P-2005. I choose NPT ensemble with temperature 298k and pressure 1atm.

Currently, the problem is the temperature reaches 298k quickly and oscillates around 3K but the pressure oscillates violently. I have run the code for 4ns but the oscillation still exists. The time average pressure oscillates between -100 atm to 100 atm. I’m curious if it is needed to do barostate. How do we define pressure for liquid? Since the pe oscillates between -11.3 to -11.4, could I consider it as equilibrated?

Another problem is could LAMMPS calculate polarization of water molecules? I find two ways in journal paper. One is to change charges of atoms to compensate the induced dipole. The other way is to directly calculate the added energy induced by polarization.

Thanks in advance.

Ps. input scripts for equilibrium of water

----------------- Init Section -----------------

units real
atom_style full

pair_style hybrid lj/cut/tip4p/long 1 2 1 1 0.1546 13.0
pair_modify pair lj/cut/tip4p/long tail yes

bond_style hybrid harmonic
angle_style hybrid harmonic
kspace_style pppm/tip4p 1.0e-5

read_data “water.equil.data”

bond_coeff 1 harmonic 0.0 0.9572
angle_coeff 1 harmonic 0.0 104.52

O-O nonbonded interactions

pair_coeff 1 1 lj/cut/tip4p/long 0.18520 3.1589

O-H and H-H interactions

pair_coeff 2 2 lj/cut/tip4p/long 0.00 0.00
pair_coeff 1 2 lj/cut/tip4p/long 0.00 0.00

Define a group for the tip4p water molecules:

group tip4p_def type 1 2

Define variable

variable Teq equal 298.0
variable Peq equal 1.0
variable Tstep equal 2.0

variable Nrun equal 100000
variable Nf equal {Nrun}/100 variable Ne equal 10 variable Nr equal {Nf}/{Ne} variable Ndump equal {Nrun}/100

variable watMoleMass equal 18.0153 # /(g/mol)
variable nAvog equal 6.0221415e23 # Avogadro’s number
variable watMoleculeMass equal ({watMoleMass}/{nAvog}) # /(g/molecule)
variable A3_in_cm3 equal 1e-24 # Angstrom^3 in cm^3
variable nAtoms equal atoms
variable nMolecules equal v_nAtoms/3

Initialize molecule

velocity all create ${Teq} 1234

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

fix fShakeTIP4P tip4p_def shake 1e-4 100 0 b 1 a 1
fix fxnpt all npt temp {Teq} {Teq} (100.0*v_Tstep) iso {Peq} {Peq} (1000.0*v_Tstep)
fix removeMomentum all momentum 1 linear 1 1 1

timestep ${Tstep}

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

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

compute PE all pe pair kspace
variable PE_Mol equal c_PE/v_nMolecules
fix PEAve_Mol all ave/time {Ne} {Nr} ${Nf} v_PE_Mol

variable Dens equal v_nMolecules*{watMoleculeMass}/(vol*{A3_in_cm3})
fix DensAve all ave/time {Ne} {Nr} ${Nf} v_Dens file wat.dens

restart (v_Nrun/10) water.equil.*.restart thermo {Nf}
thermo_style custom step time temp f_TempAve press f_PressAve f_PEAve_Mol f_DensAve
thermo_modify flush yes

dump trj all atom {Ndump} wat.equil.trj dump_modify trj append yes run {Nrun}

write_restart water.equil.end