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

[David van der Spoel]

David van der Spoel wrote:
> 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:

Maybe you can post your topology.

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