I am modelling water using a TIP3P flexible system. I am using the pair, bond and angle parameters given on the LAMMPS webpage. I am trying to run an $NPT$ simulation at $298\,K$ and $1 \,\mathrm{atm}$ - so I will be using real units in my LAMMPS simulation.
I have created a 30-by-30-by-30 angstrom box and I am placing my water molecules in a lattice in this box. I am imposing periodic boundary conditions on all sides of the box. I am using LJ and Coulomb interactions to model my system, but I am not sure what to use as a cut-off distance for the TIP3P simulation.
I expect to see a density of 0.982 g/cc, but I am not seeing that. Also, I expect pressure to be 1 atm, but my pressure in my dump output is huge - and is also negative sometimes.
What am I doing wrong? Any advice would be appreciated to debug my code.
variable data_name index sys.data
variable settings_name index sys.settings
variable nsteps index 10000
variable avg_freq index 500
variable coords_freq index 2000
variable thermo_freq index 500
variable dump4avg index 50
variable Tinit index 300
variable T0 index 298
variable Tf index 298
variable vseed1 index 8453
variable vseed2 index 8892
#===========================================================
# GENERAL PROCEDURES
#===========================================================
units real
dimension 3
newton on
boundary p p p
atom_style full
#===========================================================
# FORCE FIELD DEFINITION
#===========================================================
pair_style lj/cut/coul/cut 5 5
bond_style harmonic
angle_style harmonic
dihedral_style none
kspace_style none
improper_style none
pair_modify mix geometric shift yes
#===========================================================
# SETUP SIMULATIONS
#===========================================================
# READ IN COEFFICIENTS/COORDINATES/TOPOLOGY
read_data ${data_name} #read sys.data
#include ${settings_name} #read sys.settings
# SET RUN PARAMETERS
timestep 0.01 #real-units, in femtoseconds
run_style verlet # Velocity-Verlet integrator
# SET OUTPUTS
thermo_style custom step temp vol density etotal pe ke enthalpy press #lammps predefined keywords
thermo_modify format float %14.6f #format in log file
thermo ${thermo_freq} #store all of the above information every ${thermo_freq} time step. ${x} refers to the variable x
# DECLARE RELEVANT OUTPUT VARIABLES
variable my_step equal step #these are LAMMPS predefined words for timestep, temperature, density, potential energy, and so on
variable my_temp equal temp
variable my_rho equal density
variable my_pe equal pe
variable my_ke equal ke
variable my_etot equal etotal
variable my_ent equal enthalpy
variable my_P equal press
variable my_vol equal vol
fix averages all ave/time ${dump4avg} $(v_avg_freq/v_dump4avg) ${avg_freq} v_my_temp v_my_etot v_my_pe v_my_ke v_my_ent v_my_P file thermo.avg
#write a file called thermo.avg, compute averages of
# CREATE COORDINATE DUMPS
dump crds all custom ${coords_freq} coords.lammpstrj id type x y z vx vy vz
dump_modify crds sort id
bond_coeff * 450 0.9572 #not sure why I have to do this here
angle_coeff * 55 104.52
#===========================================================
# RUN DYNAMICS
#===========================================================
velocity all create ${Tinit} ${vseed1} mom yes rot yes #zero out momentum and any rotation
fix dynamics all npt temp ${T0} ${Tf} 298 iso 1.0 1.0 100.0
run ${nsteps}
write_data ${data_name}.end pair ii
unfix dynamics
```
