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

[Igor Leontyev]

>> However, it is emphasized that authors "haven't been able to reproduce
>> literature Data" with this topology file.
>David van der Spoel wrote:
>It's always good to be warned.
No doubts. I appreciate this.

>>Please, comment my test below
>David van der Spoel wrote:
>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?

1) It's just a test to figure out how THOLE interaction is implemented in
Gromacs. The interaction should be correctly calculated for any
configuration of Drudes (shell particles), isn't it? If it's true then I
have a freedom to chose the most appropriate configuration to simplify
analytical calculations. The test configuration which I used allows to
obtain simple expression for calculation of the THOLE energy, see my Eq(3):
> Vthole=-q(D)*q(D)*a/(alpha)^(1/3)                      (3)
BTW, for a given force constant, the test configuration is energetically
more favorable than that where all drudes are located on top of own nuclei.

2) There are no divergence if shell are on a top of own nuclei but in my
test configuration 2 shell particles are on top of the neighboring nuclei.
In this case we have zero distances r(C1-Drude2) and r(C2-Drude-1) that lead
to numerical divergence of THOLE interaction.


>David van der Spoel wrote:
>Maybe you can post your topology.
The benzene molecular topology is defined in the benzene.itp file which has
been attached to my initial post. Here I attached the structure file and
bellow is the system topology:

[ defaults ]
; nbfunc    comb-rule  gen-pairs   fudgeLJ     fudgeQQ
  1    2                  no    0.5         0.5
; Include forcefield parameters
#include "benzene.itp"
[ system ]
; Name
One polarizable benzene molecule

[ molecules ]
; Compound        #mols
   BNZ             1

---------------------------------------------------------------------------------------
Thanks David for your response.
Best regards,
Igor Leontyev









>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.pdb
Type: application/octet-stream
Size: 1367 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20090112/e3d0e665/attachment.obj>