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

[Krack Matthias (PSI)]

Dear Junting

The integral routine ai_verfc works with normalised Cartesian Gaussian functions. CP2K performs the normalisation already during the initialisation phase after the Gaussian basis set information was read.

Matthias

From: cp... at googlegroups.com <cp... at googlegroups.com> On Behalf Of Junting
Sent: Donnerstag, 17. September 2020 19:13
To: cp2k <cp... at googlegroups.com>
Subject: Re: [CP2K:13916] A question about a formula in CP2K source code

Dear Matthias,

Sorry to bother you again. I have tried both one-electron integral and two-electron integral ways to derive <a|erf(alpha_c*r)/r|b>, and they gave me the same results, which is encouraging. However, there is one remaining question about the normalization factor.

The eq.30 of this work<https://pubs.rsc.org/en/content/articlelanding/2000/CP/b001167n#!divAbstract> shows that there is a nomalization factor N = (alpha_c^2/pi)^(3/2) when rewriting <a|erf(alpha_c*r)/r|b> to a three-center two-electron repulsion integral <ab||c>. But in the work of McMurchie and Davidson<https://www.sciencedirect.com/science/article/pii/002199917890092X>, the expression for lambda in eq. (3.31) seems to be for unnormalized case. If the lambda is multiplied by N = (alpha_c^2/pi)^(3/2), the result will be perfectly in agreement with my result of ferf (i.e. 2*pi*zetp*f0*sqrt(zetc*zetw)).

Since the ferf in aobasis/ai_verfc.F matches the the form of lambda, is it supposed to be mutiplied by N? Or has it already been multiplied somewhere else in the code?

Best,

Junting

在2020年9月14日星期一 UTC+8 下午10:29:52<Matthias Krack> 写道：
Dear Junting

The one-electron integral <a_c|-Z_c*erf(alpha*|r - R_c|)/(|r - R_c|)|b> can be rewritten as a three-center two-electron repulsion integral <ab||c> (see eq. 30 of this work<https://doi.org/10.1039/B001167N>).

HTH

Matthias

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

Dear Matthias,

Thank you so much for providing me with this insightful work!

However, Eq (3.31) is in Two-electron Integral part, and seems not to involve error function. The <a|erf(r)/r|b> term is one-electron integral, describing the interaction between electron and core, instead of eletron repulsion. I still do not understand how these two are related.

And from intuition, if zetc goes to infinity, ferf should become the same as fnuc. But according to source code, they seem quite different, which also makes me confused.

Best,

Junting
在2020年9月14日星期一 UTC+8 下午4:15:43<Matthias Krack> 写道：
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: c... at googlegroups.com <c... at googlegroups.com> On Behalf Of Junting
Sent: Montag, 14. September 2020 05:13
To: cp2k <c... 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.
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>.
--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/e6dafc7a-0a1b-49ec-bd14-86583294e012n%40googlegroups.com<https://groups.google.com/d/msgid/cp2k/e6dafc7a-0a1b-49ec-bd14-86583294e012n%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
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/4c5353a6-dd2f-48b7-848c-9e8a71bc6f34n%40googlegroups.com<https://groups.google.com/d/msgid/cp2k/4c5353a6-dd2f-48b7-848c-9e8a71bc6f34n%40googlegroups.com?utm_medium=email&utm_source=footer>.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20200918/209f41b2/attachment.htm>