[gmx-developers] Polarizable simulations in Gromacs (THOLE interaction)

[Igor Leontyev]

I would like to carry out simulation of benzene solvent in the polarizable
Drude model reported in the paper from Benoit Roux and MacKerell J. Phys.
Chem. B, 2007, 111 (11), pp 2873-2885
(http://pubs.acs.org/doi/full/10.1021/jp0663614).
Although, it is announced that the polarizable Drude (shell particles) model
is implemented in Gromacs, the polarizable  simulation protocol is poorly
described in Gromacs Manual. A description of key topology elements for
definition of polarizable dipole-dipole interactions is omitted, e.g. fields
[ polarization ] and [ thole_polarization ]. COULD SOMEONE GIVE A HINT WHERE
IS A CLEAR DESCRIPTION OF THE GROMACS POLARIZABLE MD PROTOCOL ? The most
informative source regarding the subject is an example of the topology file
given here (
http://www.gromacs.org/pipermail/gmx-developers/2006-March/001571.html )
However, it is emphasized that authors "haven't been able to reproduce
literature Data" with this topology file.

I have prepared topology benzene.itp file according to the example above but
it also produces undesirable result. The THOLE interaction energy seems to
be incorrect. Please, comment my test below (the topology benzene.itp file
is also attached):

System:
Single benzene molecule, drudes (shell particles) are placed only on Carbon
atoms

Configuration:
Minimized geometry; all drudes seat in the positions of own atoms, except
only 2 neighboring drude particles. Drude-1 is seating in the position of C2
and drude-2 is seating in the position of C1. [To avoid a numerical
divergence a location of the position-exchanged drude particles actually are
slightly (0.001nm) different from the position of the neighboring heavy
atom.]

THOLE interaction:
For shuch configuration we have non-zero THOLE interaction between only
charges C1, Dride-1, C2 and Drude-2 since all other dipoles are ZERO (drudes
are seating in positions of heavy atoms). Moreover, C1-C2 and Drude1-Drude2
THOLE interaction is negligible comparably to the terms induced by
interaction of atoms C1-Drude2 and C2-Drude-1 because the distance
r(C1-Drude2)=r(C2-Drude-1)=0.001nm << r(C1-C2)=r(Drude1-Drude2)=0.137nm.
Thus, for this case THOLE energy is determined only by 2 interactions:
C1-Drude2 and C2-Drude-1, which are equal to each other due to symmetry.

In the limit of small interatomic distance rij the THOLE functional:
Vthole=qi*qj/rij*{ 1 - (1+uij/2)*exp(-u12)}, where  uij = rij *
a/(alphai*alphaj)^(1/6), is reduced to the relation:

Vthole=qi*qj *1/2*a/(alphai*alphaj)^(1/6)              (1)

In our case: alphai=alphaj=alpha, qi=-qj=q(D) and there are 2 equal
interaction pairs C1-Drude2 and C2-Drude-1 that gives:

Vthole=-q(D)*q(D)*a/(alpha)^(1/3)                      (3)

In our case q(D)= -2.0356282 and alpha= 1.376e-03,
that produces, Vthole=-13457.9 Kj/mol ,

while Gromacs gives few times larger -42282 Kj/mol:


  Energies (kJ/mol)
           Bond            U-B    Proper Dih.        LJ (SR)   Coulomb (SR)
    3.39546e+00    1.10876e-02    8.07184e-02    4.44678e+01   -4.84596e+01
   Polarization     Thole Pol.      Potential    Kinetic En.   Total Energy
    8.15375e+03   -4.22822e+04   -3.41289e+04    0.00000e+00   -3.41289e+04
    Temperature Pressure (bar)
    0.00000e+00    0.00000e+00
------------------------------------------

Just for the comparison, in the case when all drudes are seating in the
positions of atoms, THOLE enegy is zero and output is

   Energies (kJ/mol)
           Bond            U-B    Proper Dih.        LJ (SR)   Coulomb (SR)
    3.39546e+00    1.10876e-02    8.07184e-02    4.44678e+01    3.73386e+00
   Polarization     Thole Pol.      Potential    Kinetic En.   Total Energy
    0.00000e+00    0.00000e+00    5.16889e+01    0.00000e+00    5.16889e+01
    Temperature Pressure (bar)
    0.00000e+00    0.00000e+00
------------------------------------------


Igor Leontyev
-------------- next part --------------
A non-text attachment was scrubbed...
Name: benzene.itp
Type: application/octet-stream
Size: 9800 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20090112/5a8e72bb/attachment.obj>