SPC Water Model and Pair Correlation Function

LAMMPS LAMMPS Beginners
1

Hello,

I am using LAMMPS for the first time, and I want to make sure that the SPC Water model is fine, and also I am trying to get the Pair correlation function, but my data seems weird…
LMP File:
units real
atom_style full
dimension 3
boundary p p p
region box block 0.0 24.639 0.0 24.639 0.0 24.639 # ID, Type, lengths
create_box 2 box bond/types 1 angle/types 1 &
extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2 # nb types, ID, nb types

mass 1 15.9994
mass 2 1.008

pair_style lj/cut/coul/cut 10.0
pair_coeff 1 1 0.1553 3.166
pair_coeff 1 2 0.0 1.0
pair_coeff 2 2 0.0 1.0

bond_style zero
bond_coeff 1 1.0

angle_style zero
angle_coeff 1 104.5

molecule water spce.mol.txt
create_atoms 0 random 500 3454 NULL mol water 2537 overlap 1.33 # Randomly placed 1500 particules with Seed, molecule with name

fix rigid all shake 0.0001 10 10000 b 1 a 1
reset_timestep 0
velocity all create 300.0 5463576
fix integrate all nvt temp 294.0 394.0 100

#Compute RDF and save to file
compute myRDF all rdf 50 1 2
fix 1 all ave/time 100 1 100 c_myRDF[*] file tmp1.rdf mode vector
fix 2 all ave/time 100 1 100 c_myRDF[1] c_myRDF[2] c_myRDF[3] file tmp2.rdf mode vector

thermo_style custom step temp press etotal density pe ke
thermo 1000
run 10000 upto

write_data tip4p.data nocoeff

The mol file:

Water molecule. SPC/E geometry

3 atoms
2 bonds
1 angles

Coords

1 0.00000 -0.06461 0.00000
2 0.81649 0.51275 0.00000
3 -0.81649 0.51275 0.00000

Types

1 1 # O
2 2 # H
3 2 # H

Charges

1 -0.8476
2 0.4238
3 0.4238

Bonds

1 1 1 2
2 1 1 3

Angles

1 1 2 1 3

Shake Flags

1 1
2 1
3 1

Shake Atoms

1 1 2 3
2 1 2 3
3 1 2 3

Shake Bond Types

1 1 1 1
2 1 1 1
3 1 1 1

Special Bond Counts

1 2 0 0
2 1 1 0
3 1 1 0

Special Bonds

1 2 3
2 1 3
3 1 2

2

Here is an SPCE example:

https://www.lammps.org/bench/inputs/in.spce.txt

3

There are two obvious issues with your input.

  • You are running a temperature ramp. You can only compare to reference data if you have the same thermodynamic state (temperature and density)
  • There is no equilibration.
4

Hi, thank you for your answer, I tried to do a minimization, however, there was an issue with the shake command…

5

Thank you !

6

A minimization is no equilibration. This sounds like you need to find a proper tutor/mentor/adviser that can teach you MD simulation basics. A forum is not a good place that. You need someone that can look over your shoulder and guide you. Learning MD is more like learning a craft than science.

7

Yes, you are right… However, my adviser is at a conference this week… Then, do you mean that I should let my system at 294 Kelvins for some time, then I increase the temperature, then I let my system at 394 Kelvins for a few time again ? Sorry if this question is quite basic…

8

The classical paper (https://pubs.acs.org/doi/10.1021/jp003020w) uses 0.5-1.0 nanoseconds simulation time (0.5-1 million timesteps at a 1fs integration timestep) to calculate structural and dynamic data for various water models, and this was over 20 years ago.