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 

units       real
dimension   3   
newton      on  
boundary    p p p   
atom_style  full    
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 

read_data ${data_name} #read sys.data
#include ${settings_name} #read sys.settings

timestep    0.01        #real-units, in femtoseconds
run_style   verlet      # Velocity-Verlet integrator

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

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

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

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
  • 2
    $\begingroup$ Possible duplicate regarding cutoffs: mattermodeling.stackexchange.com/questions/1859/… $\endgroup$ – Godzilla Aug 16 '20 at 6:04
  • 1
    $\begingroup$ 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. $\endgroup$ – Godzilla Aug 16 '20 at 21:47
  • 1
    $\begingroup$ 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. $\endgroup$ – mykd Aug 17 '20 at 11:02
  • 1
    $\begingroup$ 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. $\endgroup$ – lr1985 Aug 18 '20 at 13:40
  • 1
    $\begingroup$ @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. $\endgroup$ – Nike Dattani Sep 13 '20 at 19:26

The author seems to have abandoned this question but I will summarize what the comments are saying.

It seems this issue was likely related to the cutoff applied to the TIP3P model. This resulted in pressure and density values that were incorrect. Since the author was doing an NPT simulation and reported an incorrect pressure was seen (the average pressure is fixed), it is also likely that they were taking instantaneous values rather than the average value which may also contribute to the error seen.


Not the answer you're looking for? Browse other questions tagged or ask your own question.