Difference in single point energy between GULP and Gromacs

GULP
1

Dear all,

I have computed the single point energy of a system consisting of 400 dimethyl ether molecules in a cubic box with unit cell length of 30 Å. I have used the General Amber Force Field (GAFF)[1] for this task. In GAFF, the 1-4 electrostatic and van der Waals interactions are scaled by 0.83333 and 0.5, respectively.

I computed the single point energy in both GULP and Gromacs, but the energy differs by 4791 kJ/mol between the two programs. I would like to understand the cause of this discrepancy.
In the following, I provide more details.

In the output of GULP, the following energy components are printed. GULP prints the energy in eV, but in the table below, I have added the energy in kJ/mol for easier comparison with Gromacs.

GULP eV kJ/mol
Interatomic potentials = 2261.656581 2.18217E+05
Three-body potentials = 25.19992868 2.43142E+03
Four-body potentials = 0.00007929 7.65032E-03
Improper torsions = 0 0.00000E+00
Out of plane potentials = 0 0.00000E+00
Monopole - monopole (real) = 87.23686723 8.41708E+03
Monopole - monopole (recip)= -78.72620441 -7.59592E+03
Monopole - monopole (total)= 8.51066283 8.21154E+02
Total lattice energy 2295.36725188 2.21469E+05

In the output of Gromacs, the energy components shown in the table below are printed.

Gromacs kJ/mol
Bond 7.5063E+02
Angle 2.4314E+03
Proper Dih. 7.6469E-03
LJ-14 3.2569E+03
Coulomb-14 3.8169E+03
LJ (SR) 2.1421E+05
Coulomb (SR) 1.4235E+03
Coul. recip. 3.7168E+02
Potential 2.2626E+05

The sum of all Coulomb components in Gromacs is 5.6120E+03 kJ/mol. The difference between this value and the “Monopole - monopole (total)” value in GULP is 4791 kJ/mol. The difference between the value of “Potential” in the Gromacs output and the value of “Total lattice energy” in the GULP output is also 4791 kJ/mol. I therefore conclude that the discrepancy between GULP and Gromacs is caused by the difference in how the Coulomb interactions are computed.

I wanted to attach the GULP input file, but I was not allowed. For this reason, I paste below the keyword line and the coulomb_subtract section. I believe these to parts of the input file should be the most relevant.

single molmec

coulomb_subtract intra x13 0.83333333
C4 C4 1.0 0.0 14.0
C4 H4 1.0 0.0 14.0
C4 O4 1.0 0.0 14.0
H4 H4 1.0 0.0 14.0
H4 O4 1.0 0.0 14.0
O4 O4 1.0 0.0 14.0

Is anyone able to spot any mistake in these two parts of the GULP input file?

Thanks for your help.

Best regards,

Torstein

[1] J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J Comput Chem 2004, 25, 1157-1174.

2

Hi Torstein,

As you deduced the issue seems to be how you are handling the Coulomb interactions. I’m guessing from your input setup that you created this yourself, rather than using the GAFF libraries that came with GULP? If you look in the libraries you should see that the Coulomb terms are handled differently:

coul o14 1/6
X X 7.0

What this means is the Coulomb subtract is only applied to 1-4 interactions and 1/6 of the interaction is removed. In your form of the Coulomb subtract, it will be applied to all atoms that are not 1-2 or 1-3 bonded & extended out to 14 Angstroms within a molecule (so you are subtracting a lot more than 1-4). In addition, you’ve scale the subtract by 5/6, rather than the interaction by 5/6. Hence this will give the wrong value. You can also see that the GULP way of doing this is more compact by using the atom wildcard (X) rather than having to specify each pair of interactions one by one.

Hope that helps,
Julian

PS I think the forum has a rule that you can’t upload files on your first post

3

Dear Julian,

Thanks a lot for your help.
I changed the coulomb_subtract line to

coulomb_subtract intra o14 1/6

and now the comparison between GULP and Gromacs is very good.
The difference between “Total lattice energy” in GULP and “Potential” in Gromacs is 0.16 kJ/mol.
The difference between “Monopole - monopole (total)” in GULP and the sum of all Coulomb components in Gromacs is 0.04 kJ/mol.

I think it is very good that Gromacs and GULP now give the same potential energy for the same system.
However, the real and reciprocal Coulomb energy in the two programs are different. In GULP, the reciprocal Coulomb energy is -7.59592E+03 kJ/mol while in Gromacs, the reciprocal Coulomb energy is 3.7168E+02 kJ/mol. Similarly, in GULP, the real Coulomb energy is 1.32080E+04 kJ/mol while in Gromacs, it is 5.24032E+03 kJ/mol.
I would like to understand why these two energy contributions are different in the two programs.

How are the real and reciprocal Coulomb energy computed in GULP?
Is it possible to control how much of the total Coulomb energy is real and how much is reciprocal?

Thanks again for your help.

Best regards,
Torstein

4

Dear Torstein,
The answer to your question as to how electrostatics is computed in periodic systems is via an Ewald sum. If you read about this method then you’ll understand why the split between real and reciprocal space is arbitrary and a matter of computational performance/convenience.
Best regards,
Julian

5

Dear Julian,

Thanks for the insight. I will read about the Ewald sum.

Best regards,
Torstein