[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
thanks in advance...
Antoine
More information about the gromacs.org_gmx-users
mailing list