[lammps-users] TIP4P not running in parallel

I am attempting to generate and equilibrate a box of TIP4P water. I consistently get the error: "ERROR on proc 0: TIP4P hydrogen is missing" while running in parallel, but not in serial. I don't think the problem is due to an incorrectly formatted file. If so, how could it run in serial? I think there could be a problem in the initial geometry obtained from "Materials Studio".

Does anyone know of software that can generate boxes of TIP4P water? Ultimately, I need to use this software to solvate metal ions and metal surfaces for corrosion studies. VMD has a nice solvation tool, but it generates TIP3P geometries, and is not amenable to my needs.

Thanks in advance.

I am attempting to generate and equilibrate a box of TIP4P water. I consistently get the error: "ERROR on proc 0: TIP4P hydrogen is missing" while running in parallel, but not in serial. I don't think the problem is due to an incorrectly formatted file. If so, how could it run in serial? I think there could be a problem in the initial geometry obtained from "Materials Studio".

Does anyone know of software that can generate boxes of TIP4P water? Ultimately, I need to use this software to solvate metal ions and metal surfaces for corrosion studies. VMD has a nice solvation tool, but it generates TIP3P geometries, and is not amenable to my needs.

christopher,

i disagree, the tip3p geometries in VMD are perfectly good for what
you want to do. the H-O-H geometry in TIP3P and TIP4P is identical
and that is all you need to input to LAMMPS. after using the solvate
tool, you will need to equilibrate your water in any case and then the
difference between TIP3P and TIP4P is small compared to the disruption
that you generate through the solvation. i strongly recommend to start
the equilibration process with a langevin thermostat to quickly dissipate
any soundwaves travelling through your system from the collapse around
the "hole" that the solvate tool will cut around your solute.

cheers,
   axel.

Axel,
I have tried to implement your recommendation of using the Langevin thermostat and nve integration. However, using the output of VMD's solvate plugin, I find that the system cannot even run one step without the "missing TIP4P Hydrogen" error. Running in serial, with langevin and nve I can run one step but the pressure is infinite.

I think that the source of this problem may be due to the initial geometry. My reasoning is as follows: Running the TIP3P model with the LJ parameters and charges in the lammps manual, the VMD solvate geometry behaves well. However, after changing the LJ parameters and charges to TIP4P, the pressure is infinite. Is there any way for me to generate TIP4P systems?

Is there any problem in my input file?
units real
dimension 3
boundary p p p
atom_style full
read_data solvate.dat
pair_style lj/cut/coul/long/tip4p 2 1 1 1 0.125 8.0
pair_coeff 1 1*2 0.000 0.000
pair_coeff 2 2 0.16275 3.16435
bond_style harmonic
bond_coeff 1 450 0.9572
angle_style harmonic
angle_coeff 1 55 104.52
kspace_style pppm/tip4p 1.0e-5

velocity all create 298.0 12345689 dist uniform
fix 1 all shake 1e-6 500 50000 m 1 a 1
neighbor 0.3 bin
neigh_modify delay 0 every 2 check yes

#output
thermo_style custom step temp pe ke etotal press
thermo_modify norm no
restart 50000 restart.tip4p

#run
timestep 0.5
fix 3 all langevin 298.0 298.0 300 123456
fix 4 all nve
run 600000

Christopher:

It sounds like the issue is with the numbering of atoms in the LAMMPS data file. If the format isn’t O-H-H for each molecule (in that order), the system will throw the error you mentioned.

If that isn’t the problem, could you post a (truncated) sample of your LAMMPS data file?

–AEI

Ahmed,
The data file is a short version of that produced by VMD's TOPOTOOLS plugin. As far as I can tell, it matches with what the LAMMPS manual specified for TIP4P input.

LAMMPS data file. CGCMM style. generated by VMD/TopoTools v1.1 on Sun Oct 17 14:56:14 EDT 2010
20535 atoms
13690 bonds
6845 angles
0 dihedrals
0 impropers
2 atom types
1 bond types
1 angle types
0 dihedral types
0 improper types
-0.048999 59.951001 xlo xhi
-0.047999 59.952001 ylo yhi
0.015501 60.015501 zlo zhi

Masses

1 1.007940 # HT
2 15.999400 # OT

Atoms

1 1 2 -1.048400 12.756000 53.344002 3.856000 # OT TIP4
2 1 1 0.524200 12.527000 53.366001 4.833000 # HT TIP4
3 1 1 0.524200 12.039000 53.747002 3.454000 # HT TIP4
4 2 2 -1.048400 51.402000 7.122000 13.392000 # OT TIP4
5 2 1 0.524200 51.609001 7.299000 12.429000 # HT TIP4
6 2 1 0.524200 50.471001 7.406000 13.589000 # HT TIP4

Bonds

1 1 1 2
2 1 1 3
3 1 4 5
4 1 4 6

Angles

1 1 2 1 3
2 1 5 4 6

christopher,

you problem is the extremely small "skin" for the neighborlists of 0.3 angstrom.
a number between 2.0 and 4.0 would be much more reasonable.

i just generated a data file from the VMD water box and
with that change to your input, the run continues nicely in parallel.

since you are using shake, you can also crank up the time step
to 2.0 fs for more throughput.

cheers,
     axel.

One addendum - a small skin distance should only be a problem
if you don't reneighbor frequently enough, i.e. neigh_modify delay 0
every 1 check yes
is the safest setting.

Steve

One addendum - a small skin distance should only be a problem
if you don't reneighbor frequently enough, i.e. neigh_modify delay 0
every 1 check yes
is the safest setting.

but isn't TIP4P different in that respect?

you have to have all three atoms for each
molecule available so you can compute
the location of point M from them. so with
a too short skin, the pair style cannot find
all three coordinates that it needs. i did
test with "delay 0 and every 1 check yes"
and it (still) fails in parallel. for TIP4P (in real
units!) you would need to have a skin of at
least 1.0.

cheers,
     axel.

I also find that running with a skin distance of <= 1.0 fails (parallel and serial) even if "delay 0 and every 1 check yes" is set with a time step as small as 0.05fs. But it will run with a skin of 2.0.

I would like to make it clear that I am using the waterbox generated by VMD (which is TIP3P) with the long-range lammps TIP4P parameters. I had to change the charges in the VMD output to reflect the TIP4P model. Consequently, the box is not at equilibrium and (with the input file attached) the initial pressure is ~40 Mbar.

Thank you all for your continuing interest in this matter.

Regards,
Chris

#### TIP4P INPUT #####
units real
dimension 3
boundary p p p
atom_style full
read_data vmdsolvate.dat

pair_style lj/cut/coul/long/tip4p 2 1 1 1 0.125 10.0
pair_coeff 1 1*2 0.000 0.000
pair_coeff 2 2 0.16275 3.16435
bond_style harmonic
bond_coeff 1 450 0.9572
angle_style harmonic
angle_coeff 1 55 104.52
kspace_style pppm/tip4p 1.0e-6

velocity all create 298.0 12345689 dist uniform
fix 1 all shake 1e-6 500 0 m 1 a 1
fix 3 all langevin 298.0 298.0 1000 123458
fix 2 all nve

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes
thermo 100
thermo_style custom step temp pe ke etotal press ecoul
thermo_modify norm no flush yes

#run variables
timestep 0.05
run 10000

I also find that running with a skin distance of <= 1.0 fails (parallel and serial) even if "delay 0 and every 1 check yes" is set with a time step as small as 0.05fs. But it will run with a skin of 2.0.

I would like to make it clear that I am using the waterbox generated by VMD (which is TIP3P) with the long-range lammps TIP4P parameters. I had to change the charges in the VMD output to reflect the TIP4P model. Consequently, the box is not at equilibrium and (with the input file attached) the initial pressure is ~40 Mbar.

you can also easily reset the charges from within lammps.
eg. with:

read_data data.water

group ox type 2
group hy type 1

set group ox charge -1.04
set group hy charge 0.52

of course the initial values for pressure and temperature
will be off due to the difference between the water potentials
and using shake. but with a proper use of the langevin
thermostat, that should dissipate quickly.

btw: i would choose a shorter time constant
during the initial equilibration for more efficient
"rattling the system into place". during that period,
you could not only crank up the time step, but also
reduce the pair cutoff and kspace convergence to
massively speed up the equilibration process
before going back to more conservative values
and an nvt/npt nose-hoover thermostat.

btw: if you are interested in squeezing out the
maximum performance, please note that the
tip4p pair style benefits more than others from
the OpenMP parallelization in my lammps-icms
branch (as it needs more computations per iteration).
so if you are running on a multi-core cluster, that
code with pair_style lj/cut/coul/long/tip4p/omp
might give you an addition performance boost.

cheers,
    axel.

Hi:

I'm running a simulation with pairwise, bonded, and angle interactions. I need to use the fix bond/break command to remove bonds in the system, but this does not affect the angle interactions.

Just wondering if there is a workaround (however inefficient) that can turn off the angle interactions between atoms that include the broken bond, but only when the bond breaks.

Thanks.

Gopinath.

The O atom only needs to find the 2 H atoms
within the force cutoff plus the skin. I presume the
cutoff is plenty long enough to include the H atoms
in a water. Am I missing something here?

Steve