# [gmx-users] Electrostatics and van Der Waals in Gromacs

Delmotte, Antoine antoine.delmotte09 at imperial.ac.uk
Fri Aug 24 17:47:57 CEST 2012

```Dear Gromacs users,

I am trying to get a full understanding of how the energy of
electrostatics and VDW interactions are calculated in Gromacs, but I am
not sure I have the right picture in mind. Could someone let me know if
the procedure below is correct?

So, let's say I have a PDB file, and I want to calculate the
electrostatic interaction between a pair of atoms. Can I do this:

1) Take the charges q_1 and q_2 of the two atoms from those in the
"aminoacids.rtp" file of the force field.
Example:
From aminoacids.rtp:
[ ALA ]
[ atoms ]
N    opls_238   -0.500     1
H    opls_241    0.300     1
CA    opls_224B   0.140     1
HA    opls_140    0.060     1
CB    opls_135   -0.180     2
If atom 1 is CA from alanine, I would choose q_1 = 0.140?

2) Choose a value for the dielectric constant
=> I read that epsilon_r = 4 is a reasonable value for the inside
of a protein, is that correct?

3) Calculate the distance r_12 between the two atoms

4) Simply evaluate:  E_electrostatics = 138.935 * q1 * q2 / (epsilon_r
* r_12);

Is that the way Gromacs is calculating the actual energy or am I missing
something? Is this at least a good approximation, or not at all? Is the
aminoacids.rtp the right place to look for the two charges? Is epsilon_r
= 4 reasonable?

Similarly, if I want to calculate the van der Waals energy between my
two atoms, can I do this:

1) Take the values of epsilon_1, epsilon_2, sigma_1  and sigma_2 from
the last two columns of the ffbonded.itp for the line corresponding to
my two atoms/

Example:
From ffnonbonded.itp:
; name  bond_type    mass    charge   ptype sigma      epsilon
opls_001   C   6      12.01100     0.500       A 3.75000e-01
4.39320e-01 ; SIG
opls_002   O   8      15.99940    -0.500       A 2.96000e-01
8.78640e-01 ; SIG
(...)
if atom 1 is opls_001 and atom 2 opls_002, epsilon_1 = 4.39320e-01,
sigma_1 = 3.75000e-01, etc...

2) calculate epsilon = sqrt(epsilon_1*epsilon_2) and sigma =
sqrt(sigma_1*sigma_2)

3) Calculate the distance r_12 between the two atoms

4) Simply evaluate E_vdw = 4*epsilon * [ (sigma/r_12)^12 -
(sigma/r_12)^6   ]

(I am aware that there are other ways to do step 2) and other potentials
than L-J, but let's assume they are those my forcefield is using).

Same question: Is that how gromacs is calculating E_vdw or am I missing
something important here?

One last detail: is there a specific term for hydrogen bonds or their
energy is simply included in the sum of electrostatics and VDW?

My apologies for this long and probably very basic question and many