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

[Krack Matthias (PSI)]

Dear Junting

Did you compare your results with the work of McMurchie and Davidson<https://www.sciencedirect.com/science/article/pii/002199917890092X> from 1978? The expression for lambda in eq. (3.31) seems to match the implemented ferf for the basic integral.

Best

Matthias

From: cp... at googlegroups.com <cp... at googlegroups.com> On Behalf Of Junting
Sent: Montag, 14. September 2020 05:13
To: cp2k <cp... at googlegroups.com>
Subject: [CP2K:13866] A question about a formula in CP2K source code

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
--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp... at googlegroups.com<mailto:cp... at googlegroups.com>.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/89673088-6dba-4f9f-9cf8-464569253964o%40googlegroups.com<https://groups.google.com/d/msgid/cp2k/89673088-6dba-4f9f-9cf8-464569253964o%40googlegroups.com?utm_medium=email&utm_source=footer>.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20200914/1d3b469f/attachment.htm>