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