Dear all,

I am exploring the capabilities of LAMMPS for water simulations. For typical water models, the CPU requirements should scale linearly with the number of distances that need to be calculated for each pair of molecules, at least approximately. In SPC like models, there are 9 such distances (all pairs O-O, O-H and H-H), while for TIP4P models there are 10 distances (the previous 9 pairs + the distance between LJ interaction sites on each molecule). However, when I run short MD simulations using both models, TIP4P takes almost twice as long as SPC. Anybody knows why this is so? Apparently, Gromacs gets the right scaling.

The final lines of the log files are:

For SPC:

Dear all,

I am exploring the capabilities of LAMMPS for water simulations. For typical

water models, the CPU requirements should scale linearly with the number of

distances that need to be calculated for each pair of molecules, at least

approximately. In SPC like models, there are 9 such distances (all pairs

O-O, O-H and H-H), while for TIP4P models there are 10 distances (the

previous 9 pairs + the distance between LJ interaction sites on each

molecule). However, when I run short MD simulations using both models, TIP4P

takes almost twice as long as SPC. Anybody knows why this is so? Apparently,

Gromacs gets the right scaling.

the main difference is that gromacs supports dummy atoms on a general

level while LAMMPS doesn't.

there is not only the number of distances that matter but also the

effort to project the M-point forces back on O and H and do it

correctly so for flexible geometries and different M offsets. i don't

know the details of the gromacs implementation of tip4p, but i would

not be surprised if there are some corners that are being cut for

better performance at the expense of some generality and accuracy.

axel.

also, i cannot reproduce your timing differences when doing a real

apples to apples comparison. there is about 50% overhead in my tests.

please see the attached log files. this is the same density and the

same cutoff settings. TIP4P has to have a few more neighbors, because

it needs to have the coordinates of the hydrogen atoms within reach to

compute the position of M where with SPC/E they don't matter. this is

managed automatically by the tip4p pair style and you can see the

larger number of neighbors in the log.

log.tip4p-shake.gz (1.68 KB)

log.spce.gz (1.7 KB)

Dear Axel,

Thanks a lot for your detailed answer.

I was a little puzzled by the difference between your results and mine. It

took me some time to find out where the difference was. If I run in the NVT

ensemble I get:

SPC

Files.tgz (3.92 KB)

Dear Axel,

[...]

There is a huge difference (100% in pair time). I am attaching the full log

files.

Do you have any idea why this is happening? Are pair forces being calculated

yes, i have.

a second time for the TIP4P model in order to get the pressure?

no. what happens is that for SPC/E LAMMPS can use a much more

efficient way to compute the stress tensor since it has the

coordinates for all particles and then the virial stress can be

computed from f_i * r_i with a single loop over all atoms. however,

for TIP4P, this is not possible since the location of the M point is

computed on the fly and then the forces and virial stress

contributions have to be tallied for each pair of atoms.

i've implemented some optimizations to the process in the styles in

USER-OMP and it would be interesting to see how this works out. you

don't need to compile with OpenMP, just install the package and use

-suffix omp when you run. those optimizations are independent of using

OpenMP.

axel.