Tip4p water molecule

LAMMPS LAMMPS General Discussion
1

Dear LAMMPS users,

I want to use TIP4P/2005 water molecules in LAMMPS. I used Moltemplate with OPLS force field to create water.lt. I have tried the following two methods.

I would appreciate if I have your advice on the following questions/problems:

• Using the first method, I can run on GPU. However, I can run only on the CPU if I apply the second method. Is there any way to run it on GPU?
• In the first approach, the system initially has high energy and pressure, and minimization does not work. How to fix very high energy when minimization does not work?
• In the first approach, If I replace the fix rigid with a fix (not rigid) and use fix shake, the bonds and atoms are lost, and the simulation ends with an error. Does the result of fix shake + fix differ from fix rigid, or are they the same?
• Which of these methods is correct and should be implemented?

First Method

I created water.pdb, including one oxygen, two hydrogens, and one massless atom. I set the charges, pair coefficients, pair style, and kspace style as follows. I then used fix rigid/npt in the simulation.

water inherits OPLSAA{
write(“Data Atoms”) {
$atom:O1 $mol @atom:65 0.00 0.0 0.00000 0.21451
$atom:H1 $mol @atom:66 0.00 0.0 0.75695 -0.37137
$atom:H2 $mol @atom:66 0.00 0.0 -0.75695 -0.37137
$atom:N1 $mol @atom:67 0.00 0.0 0.00000 0.05991
}

write(“Data Bond List”) {
$bond:O1H1 $atom:O1 $atom:H1
$bond:O1H2 $atom:O1 $atom:H2
$bond:O1N1 $atom:O1 $atom:N1
}

}
set type 65 charge 0 # oxygen
set type 66 charge 0.5564 # hydrogen
set type 67 charge -1.1128 # massless
bond_coeff 231 0.0 0.9572
pair_coeff 65 65 0.1852 3.1589
pair_coeff 66 66 0.0 0.0
pair_coeff 67 67 0.0 0.0
angle_coeff 475 0.0 104.52
angle_coeff 476 0.0 52.26
………………….
pair_style lj/cut/coul/long 12 14
……………………
kspace_style pppm 1e-5.0000

Second Method
I created water.pdb, in which the position of three atoms, including one oxygen and two hydrogens. I applied both of these pair styles: lj/cut/tip4p/long/ and lj/cut/tip4p/long/gpu

COMPND UNNAMED
AUTHOR GENERATED BY OPEN BABEL 2.3.90
HETATM 1 O HOH 1 0.000 0.000 0.215 1.00 0.00 O
HETATM 2 H HOH 1 0.000 0.757 -0.371 1.00 0.00 H
HETATM 3 H HOH 1 0.000 -0.757 -0.371 1.00 0.00 H
CONECT 2 1
CONECT 3 1
END

2

That is not correct. You can run on the GPU with the second method, but only offload the pair style to the GPU. That is actually in many cases the faster choice (the more MPI processes the more likely). Even though a GPU version of pppm (but not pppm/tip4p) exists, only part of the process is GPU accelerated and not very efficiently so. Since the GPU package can run the pair style on the GPU and all other force computations on the CPU concurrently, this is quite efficient. For that reason we have the option “pair/only” for the package gpu command implemented, so that the gpu suffix is only applied to the pair style.

While the minimize command may not work, you can still do simulated annealing, i.e. run MD with a sufficiently short time step and fix viscous to remove energy until the system has settled and you can run normal MD to further equilibrate it. Mind you, you cannot use the minimize command with fix shake either unless you use a very recent version of LAMMPS. For older versions you have to temporarily increase the force constant for rigid bonds and angles by multiple orders of magnitude during minimization. In current versions of LAMMPS this is done automatically as part of fix shake. For details, see its documentation.

Fix shake cannot work on 4-site clusters as required by the TIP4P geometry, so your input must be incorrect anyway.

If applied correctly, the trajectories obtained via fix rigid or fix shake should be statistically equivalent although not identical. If you can use fix shake, you can usually use a much larger time step and thus require less time to accumulate a fixed amount of trajectory data.

In principle, either approach is correct, if applied correctly. However, there are many subtle details that matter, so it is easy to make mistakes and then the results would be bogus.