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

[David van der Spoel]

Igor Leontyev wrote:
> 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.
It's always good to be warned.

> 
> 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.]
> 
First, why are they sitting in the wrong place?
Second, why should there be a numerical divergence if the shell are on 
top of the nuclei?

> 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
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-developers-request at gromacs.org.


-- 
David van der Spoel, Ph.D., Professor of Biology
Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205. Fax: +4618511755.
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se