[CP2K-user] A question about a formula in CP2K source code

[Junting]

Dear CP2K users and developers,

I am reading the source code of CP2K, and I have got lots of insights about 
formulas and algorithms from it. However, I have encountered a question 
when reading the aobasis/ai_verfc.F file in src directory.

The ai_verfc module is used for calculating the matrix <a|erfc(r)/r|b> in 
all-electron calculation. This term is divided into a nuclear term 
<a|1/r|b> and an error function term <a|erf(r)/r|b>. In the code prefactors 
corresponding these two terms are calculated by 

fnuc = 2.0_dp*pi*zetp*f0
ferf = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq*f0

and the basic s-orbital integrals are calculated by prefactors and 
incomplete gamma function

t = rcp2/zetp
CALL fgamma(nmax-1, t, f)
vnuc(1, 1, n) = fnuc*f(n-1)

t = -f4*rcp2/zetp
CALL fgamma(nmax-1, t, f)
verf(1, 1, n) = ferf*f(n-1)

Then the integrals over p or higher orbitals are calculated by a recursion 
procedure as is discussed in the literature: S. Obara and A. Saika, J. 
Chem. Phys. 84, 3963 (1986).

My question is that according to my derivations, my result for fnuc is 
consistent with the code, but my ferf is 
ferf = 2.0_dp*pi*zetp*f0*(-f4)**(n-0.5)
which is not consistent with the code, and is dependent on the parameter n 
in the incomplete gamma function f(n).

And it is reasonable that on the limit of zetc going to infinity, the error 
function term <a|erf(r)/r|b> should be the same as the nuclear term 
<a|1/r|b>. But according to source code, it seems not to be the case.

Therefore I want to ask, is this a potential mistake in the source code, or 
did I miss anything else in the derivations? I really appreciate it if 
someone could help me solve this question.

Best wishes,
Junting
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20200913/7683bbde/attachment.htm>