# Modelling water using the TIP3P flexible system [closed]

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_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
$$$$
`
• Possible duplicate regarding cutoffs: mattermodeling.stackexchange.com/questions/1859/… – Godzilla Aug 16 '20 at 6:04
• I would suggest that, yes. People usually use about 1.0-1.4 nm cutoffs with TIP3P so 0.5 is definitely outside of what's commonly used. That being said, I don't think that's your only problem, so you will have to reproduce the problem on a simple system, otherwise it will be very difficult to debug. – Godzilla Aug 16 '20 at 21:47
• Hi, I would recommend to see Carlos Vega's benchmark papers on water models and try to replicate their results. A 30**3 angstrom box seems too large for a first calculation, if this your first time trying MD simulations. I would recommend starting from standard 216-water box; then energy minimization with NVT equilibration. Then finally NPT production run. Keep in mind, density in NPT equilibrations will take a long time to converge, so you will possibly need to run several ns long simulations. – mykd Aug 17 '20 at 11:02
• About the pressure: when simulating an NPT ensemble, the average pressure is the one you fix, but the instantaneous one can take any value. The same holds true for the density. In general, you should first equilibrate your system and then do production runs that you need to average your observables on. – lr1985 Aug 18 '20 at 13:40
• @Godzilla I suggested before that we close such questions so that people don't waste time on them later: mattermodeling.meta.stackexchange.com/q/162/5. It seems my last comment here was a waste of time, now that I see that the asker has not been online since August 30th, and didn't respond to the last comment here on August 18th. If you support my idea, I suggest writing an answer there, about the desire to have "abandoned question" listed as a close reason. – Nike Dattani Sep 13 '20 at 19:26