[molpro-user] Calculation of C6 and C9 dispersion coefficients in propag_util.F and new_propag_util.F

[Paul Jerabek]

Dear Molpro community,

I'm interested in calculation non-additive dispersion coefficients for
homo-nuclear, closed-shell three-body systems. They can be calculated
from the multipole polarizabilities:

Z_{l_1, l_2, l_3} = (1 / \pi) * int_0^{\infty} (alpha^{l_1}(iw)
alpha^{l_1}(iw) alpha^{l_1}(iw))

[ Equation (13) in https://arxiv.org/abs/physics/0505131 ]

To my understanding, one requires the dynamic multipole polarizabilities
for that and needs to integrate over them.

Molpro provides such an functionality via

{CCSD;  CPROP,PROPAGATOR=1,EOMPROP=1, PROP_ORDER=3, HIGHW=0, DISPCOEF=12}

which I basically understand: Calculation of polarizabilities for
different imaginary frequencies and then Laguerre-Gauss summation with
certain weights.

Molpro even calculates the C_6 and C_9 coefficients (which are the
two-body and three-body dipole-dipole and dipole-dipole-dipole
interactions, respectively) and I thought it would be a good idea to try
to reproduce those numbers and then go to higher interactions (like e.g.
dipole-dipole-quadrupole).

For my attached example of Argon, the Molpro results

C_6 = 66.0
C_9 = 535.9

are in the right ball-park compared to literature:

C_6 =  61.833-64.543  [Ref: Table IV in
http://aip.scitation.org/doi/abs/10.1063]/1.463012]
C_9 =  484.2-521.7  [Ref: Table V (Z_{111} * 3 because of the
definition) http://aip.scitation.org/doi/abs/10.1063/1.463012]


But here the problems start: I don't really understand why the C_6 and
C_9 coefficients are calculated the way they are!

By looking into the source of the modules

src/eom/propag_util.F
src/eom/new_propag_util.F

I could deduce that

C_6 = 2/3 * w0 * 1/PI * Sum[ (alpha_xx+alpha_yy+alpha_zz)^2 *
weighting_factor * 1/((1+x)^2) ]

C_9 = 2/9 * w0 * 1/PI * Sum[ (alpha_xx+alpha_yy+alpha_zz)^3 *
weighting_factor * 1/((1+x)^2) ]

(with w0=0.3 <- took my a while to find this!)

and then I could reproduce the numbers.


But how does this relate to the usual definitions? E.g., for C_6 it
should be something like

C_6 = 3/PI * int_0^{\infty} (alpha^{a}(iw) alpha^{b}(iw))

[Eq. (26) in http://pubs.acs.org/doi/abs/10.1021/acs.chemrev.5b00533]


I understand that there is some averaging over the components of the
polarizabilities involved, but I cannot figure out what's going on in
these equation used in the Molpro module.

Help would be really appreciated!

Thanks in advance and best wishes,

	Paul

-- 
Dr. Paul Jerabek
Schwerdtfeger Group, Centre for Theoretical Chemistry and Physics
Bob Tindall Bldg., NZ Institute for Advanced Study
Massey University (Albany Campus)
Private Bag 102904
North Shore MSC, Auckland
New Zealand
Phone +64 9 414 0800 ext. 41698


-- 
Dr. Paul Jerabek
Fachbereich Chemie, AK Frenking, Theoretische Chemie
Philipps-Universität Marburg
Hans-Meerwein-Strasse
D-35032 Marburg
Germany

Phone: +49-(0)6421-28-27001
-------------- next part --------------
***,Ar finite field calculations

geometry={Ar};                       !z-matrix input
basis=aug-cc-pCVQZ                         !define default basis

hf

{CCSD; core, 0;  CPROP,PROPAGATOR=1,EOMPROP=1, PROP_ORDER=3, HIGHW=0, DISPCOEF=12}