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 "lcao.h"
17 #include <cmath>
18 
19 // For factorials
20 extern "C" {
21 #include <gsl/gsl_sf_gamma.h>
22 }
23 
24 namespace helfem {
25   namespace lcao {
double_factorial(unsigned int n)26     static double double_factorial(unsigned int n) {
27       return gsl_sf_doublefact(n);
28     }
factorial(unsigned int n)29     static double factorial(unsigned int n) {
30       return gsl_sf_fact(n);
31     }
32 
33     /// Evaluate radial GTO
radial_GTO(double r,int l,double alpha)34     double radial_GTO(double r, int l, double alpha) {
35       return std::pow(2,l+2) * std::pow(alpha,(2*l+3)/4.0) * std::pow(r,l) * exp(-alpha*r*r) / ( std::pow(2.0*M_PI,0.25) * sqrt(double_factorial(2*l+1)));
36     }
37 
radial_GTO(const arma::vec & r,int l,const arma::vec & alpha)38     arma::mat radial_GTO(const arma::vec & r, int l, const arma::vec & alpha) {
39       arma::mat gto(r.n_elem, alpha.n_elem);
40       for(size_t j=0;j<alpha.n_elem;j++)
41         for(size_t i=0;i<r.n_elem;i++)
42           gto(i,j)=radial_GTO(r(i),l,alpha(j));
43       return gto;
44     }
45 
46     /// Evaluate radial STO
radial_STO(double r,int l,double zeta)47     double radial_STO(double r, int l, double zeta) {
48       return std::pow(2*zeta,l+1.5)/sqrt(factorial(2*l+2)) * std::pow(r,l) * exp(-zeta*r);
49     }
50 
radial_STO(const arma::vec & r,int l,const arma::vec & alpha)51     arma::mat radial_STO(const arma::vec & r, int l, const arma::vec & alpha) {
52       arma::mat gto(r.n_elem, alpha.n_elem);
53       for(size_t j=0;j<alpha.n_elem;j++)
54         for(size_t i=0;i<r.n_elem;i++)
55           gto(i,j)=radial_STO(r(i),l,alpha(j));
56       return gto;
57     }
58   }
59 }
60