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

[Delmotte, Antoine]

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 
thanks in advance...

Antoine