1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * HF/DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2011
9 * Copyright (c) 2010-2011, 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
18 #include <cstdio>
19 #include <cmath>
20 #include "gaunt.h"
21 #include "lmgrid.h"
22
23 extern "C" {
24 // 3j symbols
25 #include <gsl/gsl_sf_coupling.h>
26 }
27
gaunt_coefficient(int L,int M,int l,int m,int lp,int mp)28 double gaunt_coefficient(int L, int M, int l, int m, int lp, int mp) {
29 // First, compute square root factor
30 double res=sqrt((2*L+1)*(2*l+1)*(2*lp+1)/(4.0*M_PI));
31 // Plug in (l1 l2 l3 | 0 0 0), GSL uses half units
32 res*=gsl_sf_coupling_3j(2*L,2*l,2*lp,0,0,0);
33 // Plug in (l1 l2 l3 | m1 m2 m3)
34 res*=gsl_sf_coupling_3j(2*L,2*l,2*lp,-2*M,2*m,2*mp);
35 // and the phase factor
36 res*=pow(-1.0,M);
37
38 return res;
39 }
40
Gaunt()41 Gaunt::Gaunt() {
42 }
43
Gaunt(int Lmax,int lmax,int lpmax)44 Gaunt::Gaunt(int Lmax, int lmax, int lpmax) {
45 // Allocate storage
46 table=arma::cube(lmind(Lmax,Lmax)+1,lmind(lmax,lmax)+1,lmind(lpmax,lpmax)+1);
47
48 // Compute coefficients
49 for(int L=0;L<=Lmax;L++)
50 for(int M=-L;M<=L;M++)
51
52 for(int l=0;l<=lmax;l++)
53 for(int m=-l;m<=l;m++)
54
55 for(int lp=0;lp<=lpmax;lp++)
56 for(int mp=-lp;mp<=lp;mp++)
57 table(lmind(L,M),lmind(l,m),lmind(lp,mp))=gaunt_coefficient(L,M,l,m,lp,mp);
58 }
59
~Gaunt()60 Gaunt::~Gaunt() {
61 }
62
coeff(int L,int M,int l,int m,int lp,int mp) const63 double Gaunt::coeff(int L, int M, int l, int m, int lp, int mp) const {
64 return table(lmind(L,M),lmind(l,m),lmind(lp,mp));
65 }
66