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