I am simulating TIP4P/2005 water model. Everything seems to run correctly as I reproduce the liquid density from the original TIP4P/2005 reference. However what puzzles me is how big is the kspace energy (elong) relative to short-range Coulomb (ecoul). I am getting elong = ca -285000 kcal/mol and ecoul = ca 231000 kcal/mol. When simulating similar systems in GROMACS I was used to reciprocal sum values being two orders of magnitude smaller than the direct sum. I tried decreasing the value of gewald manually to 0.1 A^{-1} which should put more weight on direct sum, but I get elong = ca -139000 kcal/mol and ecoul = only about 85000 kcal/mol. Does anybody understand, why kspace sum is so big? I wanted to attach my input files but apparently new users are not allowed to do that.
@Jan_J
Some points to consider:
- absolute energy values have no meaning in classical MD, only energy differences between different states or their gradients (forces) have a meaning
- You should not look at the real-space versus kspace values anyway and only consider the sum of the two (with the caveat of the previous statement).
- Gromacs uses a smooth PME variant while LAMMPS has either PPPM or a plain explicit Ewald sum. Each of these may have a different reference energy, which in turn may also depend on the details of your simulation setup
- It is usually a bad idea to change the individual long-range settings. If you want to change the balance between real space and reciprocal space, you should just change the coulomb cutoff, everything else will be computed from that cutoff and the convergence parameter.
- In either code, if you want to compare the accuracy, you should look at the limit of the forces, compare codes or e.g. with coul/cut using a very, very large cutoff (which requires a system with a small number of atoms).
Be aware that GROMACS has Gnarly Routines Optimised Mighty Aggressively Chasing Speed. In particular, its default behaviour is to tune the short-range / long-range electrostatics balance for the first few thousand steps of the simulation, increasing the cutoff as necessary where the cost of larger neighbour lists is offset by the reduced parallelism required for K-space FFTs.
LAMMPS is more modular* – the KSpace class does not have permission to modify the neighbourlist cutoff willy-nilly, from what I remember. Nonetheless, you can emulate GROMACS’s “warm-up” by setting the Coulombic cutoff much larger than the non-Coulombic cutoff and then using fix tune/kspace (fix tune/kspace command — LAMMPS documentation) to adjust the KSpace parameters (see caveats from this old forum post).
*LAMMPS has more “compartmentalized” objects in its object-oriented programming style. In my opinion this makes it easier to maintain and to add new or different features to, while making more esoteric tricks more tricky. Not better or worse than GROMACS, but a different tool for different (worlds of) jobs.
To get into the specifics: an inverse length of 0.1 per angstrom is one over a length of 10 angstroms, so that (essentially) each point charge is now modelled, short-range, as a Gaussian charge of width \sigma = 10A (with the long-range calculation compensating for the difference).
But a water model usually has a 12A cutoff. And a Gaussian curve has about a third of its distribution outside \pm 1\sigma! So the short-range Coulombic calculation with those settings would be expected to miss 1-10% of the Gaussian charges’ interactions.
Of course, pure water has no free charges, so in practice and on average most of those interactions would cancel out (being half repulsive and half attractive) and you would still get somewhat usable results. But you should expect such a simulation to show poor energy conservation and, with enough statistics, meaningfully different behaviour from the standard model.