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

[Igor Leontyev]

Ok, there is some progress in this subject.

1) The THOLE energy estimation in my initial post utilized very raff
approximation that contributions from C1-C2 and drude1-drude2 interactions
can be neglected comparably to (C1-drude2) and (C2-drude1)terms. A precise
calculation of the THOLE energy gives twice smaller absolute value:
Ethole=-6040 Kj/mol.

2) Calculation of the same structure in polarizable Drude model by CHARMM MD 
package supports this
value. Thanks Igor Vorobyov for the simulation.

3) The precise value deffers from the Gromacs output (see original post
Ethole=-4.22822e+04) even more significantly than one initially
estimated by Eq.(3). However, mdrun with "-debug" flag has shown that
the procedure "thole_pol" returns correct values of all THOLE energy
components, but the procedure was called 7 times instead of 1. You can check
the output energy -4.22822e+04 differs from the correct value "-6040" in 7
times. Running MD with different parallelization parameters revealeded that
the number of calls of "thole_pol" procedure and THOLE energy depends on the
number of used cpu. Thus, for np=4, ncalls=7, Ethole=-4.22822e+04; while,
for np=1, ncalls=12, Ethole=-7.24838e+04. The last THOLE energy value
differs from the correct number in 12 times precisely.

4) Thus, it looks like there is a problem in a code, possibly, in parallel
distribution of the THOLE energy calculation. Perhaps, this is the reason
why you haven't been able to reproduce literature data for Ethanol. I
localized the problem, but can not proceed deeper since I am not familiar
with gromacs code and parallel programming. I use Gromacs 3.3. The problem
may be fixed in version 4.0, but I can not check it since I don't have the
version installed. Hopefully, above will be helphul for someone who will be
interested in polarizable simulations.


>
> 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