1 /*
2  *                This source code is part of
3  *
4  *                          HelFEM
5  *                             -
6  * Finite element methods for electronic structure calculations on small systems
7  *
8  * Written by Susi Lehtola, 2018-
9  * Copyright (c) 2018- Susi Lehtola
10  *
11  * This program is free software; you can redistribute it and/or
12  * modify it under the terms of the GNU General Public License
13  * as published by the Free Software Foundation; either version 2
14  * of the License, or (at your option) any later version.
15  */
16 
17 #include "erfc_expn.h"
18 #include <cmath>
19 #include <cfloat>
20 #include <stdexcept>
21 #include <cstdio>
22 
23 // For factorials
24 extern "C" {
25 #include <gsl/gsl_sf_gamma.h>
26 }
27 
28 namespace helfem {
29   namespace atomic {
30     namespace erfc_expn {
double_factorial(unsigned int n)31       inline static double double_factorial(unsigned int n) {
32         if(n==0)
33           return 1.0;
34 
35         return gsl_sf_doublefact(n);
36       }
37 
factorial(unsigned int n)38       inline static double factorial(unsigned int n) {
39         if(n==0)
40           return 1.0;
41 
42         return gsl_sf_fact(n);
43       }
44 
choose(int n,int m)45       inline static double choose(int n, int m) {
46         // Special cases
47         if(n==-1)
48           // choose(-1,m) = (-1)^m
49           return std::pow(-1.0,m);
50         if(n==0)
51           // choose(0,m) = 0 except for choose(0,0) = 1
52           return m==0 ? 1.0 : 0.0;
53         if(m==0)
54           // choose(n,0) = 1 for all n
55           return 1.0;
56         if(m==1)
57           // choose(n,1) = n for all n
58           return n;
59         if(n>0 && m>0 && m>n)
60           // choose(n,m) = 0 for m>n positive
61           return 0.0;
62 
63         // Negative binomials
64         if(n<0) {
65           return choose(n+m-1,m)*std::pow(-1,m);
66         } else {
67           return gsl_sf_choose(n,m);
68         }
69       }
70 
71       // Angyan et al, equation (22)
Fn(unsigned int n,double Xi,double xi)72       inline static double Fn(unsigned int n, double Xi, double xi) {
73         // Exponential factors
74         double explus(std::exp(-std::pow(Xi+xi,2)));
75         double exminus(std::exp(-std::pow(Xi-xi,2)));
76 
77         // Prefactor
78         double prefac(-1.0/(4.0*Xi*xi));
79 
80         double F=0.0;
81         // Looks like there's a typo in the equation: I can't make the
82         // function match the equations in the appendix unless the
83         // lower limit is 0 instead of 1.
84         for(unsigned int p=0;p<=n;p++) {
85           F += std::pow(prefac,p+1) * (factorial(n+p)/(factorial(p)*factorial(n-p))) * (std::pow(-1,n-p) * explus - exminus);
86         }
87         // Apply prefactor
88         return 2.0/sqrt(M_PI)*F;
89       }
90 
91       // Angyan et al, equation (24)
Hn(unsigned int n,double Xi,double xi)92       inline static double Hn(unsigned int n, double Xi, double xi) {
93         if(Xi<xi)
94           throw std::logic_error("Xi < xi");
95 
96         double Xi2np1=std::pow(Xi,2*n+1);
97         double xi2np1=std::pow(xi,2*n+1);
98 
99         double Hn = (Xi2np1+xi2np1)*erfc(Xi+xi) - (Xi2np1-xi2np1)*erfc(Xi-xi);
100         return Hn/(2.0*std::pow(xi*Xi,n+1));
101       }
102 
103       // Angyan et al, equation (21)
Phi_general(unsigned int n,double Xi,double xi)104       double Phi_general(unsigned int n, double Xi, double xi) {
105         // Make sure arguments are in the correct order
106         if(Xi < xi)
107           std::swap(Xi,xi);
108 
109         double Fnarr[n];
110         for(unsigned int i=0;i<n;i++)
111           Fnarr[i]=Fn(i,Xi,xi);
112 
113         double sum = 0.0;
114         for(unsigned int m=1;m<=n;m++) {
115           double Xim(std::pow(Xi,m));
116           double xim(std::pow(xi,m));
117           sum += Fnarr[n-m]*((Xim*Xim + xim*xim)/(Xim*xim));
118         }
119 
120         return Fn(n,Xi,xi) + sum + Hn(n,Xi,xi);
121       }
122 
123       // Angyan et al, equations 28 and 29
Dnk(int n,int k,double Xi)124       double Dnk(int n, int k, double Xi) {
125         // Prefactor
126         double prefac = std::exp(-std::pow(Xi,2))/sqrt(M_PI)*std::pow(2,n+1)*std::pow(Xi,2*n+1);
127 
128         double D = 0.0;
129         if(k==0) {
130           // Compute the sum
131           double sum = 0.0;
132           for(int m=1;m<=n;m++)
133             sum += 1.0/(double_factorial(2*(n-m)+1)*std::pow(2*Xi*Xi,m));
134 
135           D = erfc(Xi) + prefac*sum;
136         } else {
137           // Compute the sum
138           double sum = 0.0;
139           for(int m=1;m<=k;m++)
140             sum += choose(m-k-1,m-1)*std::pow(2*Xi*Xi,k-m)/double_factorial(2*(n+k-m)+1);
141 
142           D = prefac * (2.0*n+1.0)/(factorial(k)*(2.0*(n+k)+1.0)) * sum;
143         }
144 
145         return D;
146       }
147 
148       // Angyan et al, equation (30), evaluated to full numerical precision
Phi_short(unsigned int n,double Xi,double xi)149       double Phi_short(unsigned int n, double Xi, double xi) {
150         // Make sure arguments are in the correct order
151         if(Xi < xi)
152           std::swap(Xi,xi);
153 
154         double Phi = 0.0;
155         double dPhi;
156         unsigned int k;
157         double tol=DBL_EPSILON;
158         for(k=0; k<=30; k+=2) {
159           // Unroll odd values so that we don't truncate too soon by
160           // accident
161           dPhi = Dnk(n,k  ,Xi)*std::pow(xi,n+2*k)
162             +    Dnk(n,k+1,Xi)*std::pow(xi,n+2*(k+1));
163           Phi += dPhi;
164           if(std::abs(dPhi) < tol*std::abs(Phi)) break;
165         }
166         if(std::abs(dPhi) >= tol*std::abs(Phi))
167           fprintf(stderr,"Warning - short-range Phi not converged, ratio %e\n",dPhi/Phi);
168         //fprintf(stderr,"Phi%i(%e,%e) converged in %u iterations\n",n,Xi,xi,k);
169 
170         return Phi/std::pow(Xi,n+1);
171       }
172 
173       // Wrapper
Phi(unsigned int n,double Xi,double xi)174       double Phi(unsigned int n, double Xi, double xi) {
175         // Make sure arguments are in the correct order
176         if(Xi < xi)
177           std::swap(Xi,xi);
178 
179         // See text on top of page 8624 of Angyan et al
180         if(xi < 0.4 || (Xi < 0.5 && xi < 2*Xi)) {
181           // Short-range Taylor polynomial
182           return Phi_short(n,Xi,xi);
183         } else {
184           // General expansion, susceptible to numerical noise for
185           // small arguments
186           return Phi_general(n,Xi,xi);
187         }
188       }
189     }
190   }
191 }
192