There is an error in the GROMACS documentation on this page:
It is stated incorrectly that the Lennard Jones long range correction energy for the SPC water system is -0.75 KJ/mole, when in fact it is 3 times less at -0.25KJ/mole. The virial correction is correctly stated as being 0.75KJ/mole.
Can someone correct the documentation, and also confirm that the actual GROMACS program calculates Lennard Jones long range correction energies in a correct manner?
The values in the manual are indeed incorrect. But how did you get to 0.25 kJ/mol per molecule? I get 0.5 kJ/mol per molecule.
This code is covered in our integration and regression tests and values have been compared against theoretical values and other codes, so I’m pretty sure the code is correct.
Doing the calculation step-by-step:
The formula for the longe range energy (assuming g(r)=1 for r > r_c) is:
V_lr = (-2/3) * pi * N * rho * C_6 * r_c^-3 where * denotes multiplication.
Taking each term:
For SPC water, the Lennard Jones oxygen site parameters are:
epsilon = 78.2 Kelvin, and sigma = 3.1656 Angstroms
The corresponding C_6 parameter is therefore 4 * epsilon * sigma^6
= 4 * (78.2 Kelvin) * (3.1656 Angstroms)^6
= 314776.98 Kelvin*Angstroms^6
= 2.6171849e-57 (KJ/mol)metres^6
(There are no Lennard Jones parameters for the hydrogen sites, so there is just one LJ site per water molecule.)
r_c^-3 = (9 Angstroms)^-3 = 1.3717421e27 metres^-3
rho = N/V, so if N = Avagado’s number, N_A = 6.02214076e23 then the mass of N_A water molecules is 18.015 grams.
Since the density of water is 1g/cm^3 then we have a volume of 18.015 cm^3 = 1.8015e-5 metres^3.
Therefore N/V = 55509.298 metres^-3
Therefore we have finally:
V_lr = (-2/3) * pi * N * rho * C_6 * r_c^-3
= (-2/3) * pi * N_A * ( 55509.298 metres^-3 ) * ( 2.6171849e-57 (KJ/mol)metres^6 ) * ( 1.3717421e27 metres^-3 )
= -0.25135187 KJ/mol
To double check, the factor to get from V_lr to the pressure correction P_lr is just 2/V so we have:
P_lr = -0.25135187 KJ/mol * (2/V)
= -0.25135187 KJ/mol * 2 / ( 1.8015e-5 metres^3 )
= -27940.732 (KJ/mol)metres^-3
= -27940.732 * 1e3 (J/mol)metres^-3 / 1e5 pascals = -279.40732 bar
I cannot see how you get V_lr = 0.5 KJ/mol as the long range correction energy.