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 #include <cstdio>
17 #include <cmath>
18 #include "gaunt.h"
19 
20 extern "C" {
21   // 3j symbols
22 #include <gsl/gsl_sf_coupling.h>
23 }
24 
25 /// Index of (l,m) in tables: l^2 + l + m
26 #define genind(l,m) ( ((size_t) (l))*(size_t (l)) + (size_t) (l) + (size_t) (m))
27 /// Index of (l,m) in m limited table
28 #define LMind(L,M) ( ((size_t) (L))*(size_t (2*Mmax+1)) + (size_t) (Mmax) + (size_t) (M))
29 #define lmind(L,M) ( ((size_t) (L))*(size_t (2*mmax+1)) + (size_t) (mmax) + (size_t) (M))
30 #define lpmpind(L,M) ( ((size_t) (L))*(size_t (2*mpmax+1)) + (size_t) (mpmax) + (size_t) (M))
31 
32 namespace helfem {
33   namespace gaunt {
34 
gaunt_coefficient(int L,int M,int l,int m,int lp,int mp)35     double gaunt_coefficient(int L, int M, int l, int m, int lp, int mp) {
36       if(M != m+mp)
37         // Coefficient vanishes if l values don't match
38         return 0.0;
39       if(L < abs(l-lp) || L > l+lp)
40         // Can't couple
41         return 0.0;
42 
43       // First, compute square root factor
44       double res=sqrt((2*L+1)*(2*l+1)*(2*lp+1)/(4.0*M_PI));
45       // Plug in (l1 l2 l3 | 0 0 0), GSL uses half units
46       res*=gsl_sf_coupling_3j(2*L,2*l,2*lp,0,0,0);
47       // Plug in (l1 l2 l3 | m1 m2 m3)
48       res*=gsl_sf_coupling_3j(2*L,2*l,2*lp,-2*M,2*m,2*mp);
49       // and the phase factor
50       res*=pow(-1.0,M);
51 
52       return res;
53     }
54 
modified_gaunt_coefficient(int lj,int mj,int L,int M,int li,int mi)55     double modified_gaunt_coefficient(int lj, int mj, int L, int M, int li, int mi) {
56       static const double const0(2.0/3.0*sqrt(M_PI));
57       static const double const2(4.0/15.0*sqrt(5.0*M_PI));
58 
59       // Coupling through Y_0^0
60       double cpl0(gaunt_coefficient(L,M,0,0,L,M)*gaunt_coefficient(lj,mj,li,mi,L,M));
61 
62       // Coupling through Y_2^0
63       double cpl2=0.0;
64       for(int Lp=std::max(std::max(L-2,0),std::abs(M));Lp<=L+2;Lp++)
65         cpl2+=gaunt_coefficient(Lp,M,2,0,L,M)*gaunt_coefficient(lj,mj,li,mi,Lp,M);
66 
67       return const0*cpl0+const2*cpl2;
68     }
69 
Gaunt()70     Gaunt::Gaunt() {
71     }
72 
Gaunt(int Lmax,int lmax,int lpmax)73     Gaunt::Gaunt(int Lmax, int lmax, int lpmax) {
74       // Allocate storage
75       mlimit=false;
76       table=arma::cube(genind(Lmax,Lmax)+1,genind(lmax,lmax)+1,genind(lpmax,lpmax)+1);
77 
78       // Compute coefficients
79 #ifdef _OPENMP
80 #pragma omp parallel for collapse(3)
81 #endif
82       for(int L=0;L<=Lmax;L++)
83 	  for(int l=0;l<=lmax;l++)
84             for(int lp=0;lp<=lpmax;lp++)
85               for(int M=-L;M<=L;M++)
86                 for(int m=-l;m<=l;m++)
87                   for(int mp=-lp;mp<=lp;mp++)
88                     table(genind(L,M),genind(l,m),genind(lp,mp))=gaunt_coefficient(L,M,l,m,lp,mp);
89     }
90 
Gaunt(int Lmax,int Mmax_,int lmax,int mmax_,int lpmax,int mpmax_)91     Gaunt::Gaunt(int Lmax, int Mmax_, int lmax, int mmax_, int lpmax, int mpmax_) : Mmax(Mmax_), mmax(mmax_), mpmax(mpmax_) {
92       // Allocate storage
93       mlimit=true;
94       table=arma::cube(LMind(Lmax,Mmax)+1,lmind(lmax,mmax)+1,lpmpind(lpmax,mpmax)+1);
95 
96       // Compute coefficients
97 #ifdef _OPENMP
98 #pragma omp parallel for collapse(3)
99 #endif
100       for(int L=0;L<=Lmax;L++)
101 	  for(int l=0;l<=lmax;l++)
102             for(int lp=0;lp<=lpmax;lp++)
103               for(int M=-Mmax;M<=Mmax;M++)
104                 for(int m=-mmax;m<=mmax;m++)
105                   for(int mp=-mpmax;mp<=mpmax;mp++)
106                     table(LMind(L,M),lmind(l,m),lpmpind(lp,mp))=gaunt_coefficient(L,M,l,m,lp,mp);
107     }
108 
~Gaunt()109     Gaunt::~Gaunt() {
110     }
111 
coeff(int L,int M,int l,int m,int lp,int mp) const112     double Gaunt::coeff(int L, int M, int l, int m, int lp, int mp) const {
113       if(std::abs(M)>L) return 0.0;
114       if(std::abs(m)>l) return 0.0;
115       if(std::abs(mp)>lp) return 0.0;
116 
117       size_t irow, icol, islice;
118       if(mlimit) {
119         irow=LMind(L,M);
120         icol=lmind(l,m);
121         islice=lpmpind(lp,mp);
122       } else {
123         irow=genind(L,M);
124         icol=genind(l,m);
125         islice=genind(lp,mp);
126       }
127 
128 #ifndef ARMA_NO_DEBUG
129       if(irow>table.n_rows) {
130         std::ostringstream oss;
131         oss << "Row index overflow for coeff(" << L << "," << M << "," << l << "," << m << "," << lp << "," << mp << ")!\n";
132         oss << "Wanted element at (" << irow << "," << icol << "," << islice << ") but table is " << table.n_rows << " x " << table.n_cols << " x " << table.n_slices << "\n";
133         throw std::logic_error(oss.str());
134       }
135 
136       if(icol>table.n_cols) {
137         std::ostringstream oss;
138         oss << "Column index overflow for coeff(" << L << "," << M << "," << l << "," << m << "," << lp << "," << mp << ")!\n";
139         oss << "Wanted element at (" << irow << "," << icol << "," << islice << ") but table is " << table.n_rows << " x " << table.n_cols << " x " << table.n_slices << "\n";
140         throw std::logic_error(oss.str());
141       }
142 
143       if(islice>table.n_slices) {
144         std::ostringstream oss;
145         oss << "Slice index overflow for coeff(" << L << "," << M << "," << l << "," << m << "," << lp << "," << mp << ")!\n";
146         oss << "Wanted element at (" << irow << "," << icol << "," << islice << ") but table is " << table.n_rows << " x " << table.n_cols << " x " << table.n_slices << "\n";
147         throw std::logic_error(oss.str());
148       }
149 #endif
150 
151       return table(irow,icol,islice);
152     }
153 
cosine_coupling(int lj,int mj,int li,int mi) const154     double Gaunt::cosine_coupling(int lj, int mj, int li, int mi) const {
155       // cos th = const1 * Y_1^0
156       static const double const1(2.0*sqrt(M_PI/3.0));
157       return const1*coeff(lj,mj,1,0,li,mi);
158     }
159 
cosine2_coupling(int lj,int mj,int li,int mi) const160     double Gaunt::cosine2_coupling(int lj, int mj, int li, int mi) const {
161       // cos^2 th = const0 * Y_0^0 + const2 * Y_2^0
162       static const double const0(2.0/3.0*sqrt(M_PI));
163       static const double const2(4.0/15.0*sqrt(5.0*M_PI));
164       return const0*coeff(lj,mj,0,0,li,mi) + const2*coeff(lj,mj,2,0,li,mi);
165     }
166 
mod_coeff(int lj,int mj,int L,int M,int li,int mi) const167     double Gaunt::mod_coeff(int lj, int mj, int L, int M, int li, int mi) const {
168       static const double const0(2.0/3.0*sqrt(M_PI));
169       static const double const2(4.0/15.0*sqrt(5.0*M_PI));
170 
171       // Coupling through Y_0^0
172       double cpl0(coeff(L,M,0,0,L,M)*coeff(lj,mj,li,mi,L,M));
173 
174       // Coupling through Y_2^0
175       double cpl2=0.0;
176       for(int Lp=std::max(std::max(L-2,0),std::abs(M));Lp<=L+2;Lp++)
177         cpl2+=coeff(Lp,M,2,0,L,M)*coeff(lj,mj,li,mi,Lp,M);
178 
179       return const0*cpl0+const2*cpl2;
180     }
181 
cosine3_coupling(int lj,int mj,int li,int mi) const182     double Gaunt::cosine3_coupling(int lj, int mj, int li, int mi) const {
183       // cos^3 th = const1 * Y_1^0 + const3 * Y_3^0
184       static const double const1(2.0/5.0*sqrt(3.0*M_PI));
185       static const double const3(4.0/35.0*sqrt(7.0*M_PI));
186       return const1*coeff(lj,mj,1,0,li,mi) + const3*coeff(lj,mj,3,0,li,mi);
187     }
188 
cosine4_coupling(int lj,int mj,int li,int mi) const189     double Gaunt::cosine4_coupling(int lj, int mj, int li, int mi) const {
190       // cos^4 th = const0 * Y_0^0 + const2 * Y_2^0 + const4 * Y_4^0
191       static const double const0(2.0/5.0*sqrt(M_PI));
192       static const double const2(8.0/35.0*sqrt(5.0*M_PI));
193       static const double const4(16.0/105.0*sqrt(M_PI));
194       return const0*coeff(lj,mj,0,0,li,mi) + const2*coeff(lj,mj,2,0,li,mi) + const4*coeff(lj,mj,4,0,li,mi);
195     }
196 
cosine5_coupling(int lj,int mj,int li,int mi) const197     double Gaunt::cosine5_coupling(int lj, int mj, int li, int mi) const {
198       // cos^5 th = const1 * Y_1^0 + const3 * Y_3^0 + const5 * Y_5^0
199       static const double const1(2.0/7.0*sqrt(3.0*M_PI));
200       static const double const3(8.0/63.0*sqrt(7.0*M_PI));
201       static const double const5(16.0/693.0*sqrt(11.0*M_PI));
202       return const1*coeff(lj,mj,1,0,li,mi) + const3*coeff(lj,mj,3,0,li,mi) + const5*coeff(lj,mj,5,0,li,mi);
203     }
204 
sine2_coupling(int lj,int mj,int li,int mi) const205     double Gaunt::sine2_coupling(int lj, int mj, int li, int mi) const {
206       // sin^2 th = const0 * Y_0^0 + const2 * Y_2^0
207       static const double const0(4.0/3.0*sqrt(M_PI));
208       static const double const2(-4.0/15.0*sqrt(5.0*M_PI));
209       return const0*coeff(lj,mj,0,0,li,mi) + const2*coeff(lj,mj,2,0,li,mi);
210     }
211 
cosine2_sine2_coupling(int lj,int mj,int li,int mi) const212     double Gaunt::cosine2_sine2_coupling(int lj, int mj, int li, int mi) const {
213       // cos^2 th sin^2 th = const0 * Y_0^0 + const2 * Y_2^0 + const4 * Y_4^0
214       static const double const0(4.0/15.0*sqrt(M_PI));
215       static const double const2(4.0/105.0*sqrt(5.0*M_PI));
216       static const double const4(-16.0/105.0*sqrt(M_PI));
217       return const0*coeff(lj,mj,0,0,li,mi) + const2*coeff(lj,mj,2,0,li,mi) + const4*coeff(lj,mj,4,0,li,mi);
218     }
219   }
220 }
221