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

[Junting]

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>
> .
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20200917/855b030c/attachment.htm>