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 "basis.h"
17 #include "quadrature.h"
18 #include "polynomial_basis.h"
19 #include "chebyshev.h"
20 #include "../general/spherical_harmonics.h"
21 #include "../general/gaunt.h"
22 #include "utils.h"
23 #include "../general/scf_helpers.h"
24 #include <cassert>
25 #include <cfloat>
26 #include <helfem.h>
27 
28 #ifdef _OPENMP
29 #include <omp.h>
30 #endif
31 
32 namespace helfem {
33   namespace atomic {
34     namespace basis {
concatenate_grid(const arma::vec & left,const arma::vec & right)35       static arma::vec concatenate_grid(const arma::vec & left, const arma::vec & right) {
36         if(!left.n_elem)
37           return right;
38         if(!right.n_elem)
39           return left;
40 
41         if(left(0) != 0.0)
42           throw std::logic_error("left vector doesn't start from zero");
43         if(right(0) != 0.0)
44           throw std::logic_error("right vector doesn't start from zero");
45 
46         // Concatenated vector
47         arma::vec ret(left.n_elem + right.n_elem - 1);
48         ret.subvec(0,left.n_elem-1)=left;
49         ret.subvec(left.n_elem,ret.n_elem-1)=right.subvec(1,right.n_elem-1) + left(left.n_elem-1)*arma::ones<arma::vec>(right.n_elem-1);
50         return ret;
51       }
52 
normal_grid(int num_el,double rmax,int igrid,double zexp)53       arma::vec normal_grid(int num_el, double rmax, int igrid, double zexp) {
54         return utils::get_grid(rmax,num_el,igrid,zexp);
55       }
56 
finite_nuclear_grid(int num_el,double rmax,int igrid,double zexp,int num_el_nuc,double rnuc,int igrid_nuc,double zexp_nuc)57       arma::vec finite_nuclear_grid(int num_el, double rmax, int igrid, double zexp, int num_el_nuc, double rnuc, int igrid_nuc, double zexp_nuc) {
58         if(num_el_nuc) {
59           // Grid for the finite nucleus
60           arma::vec bnuc(utils::get_grid(rnuc,num_el_nuc,igrid_nuc,zexp_nuc));
61           // and the one for the electrons
62           arma::vec belec(utils::get_grid(rmax-rnuc,num_el,igrid,zexp));
63 
64           arma::vec bnucel(concatenate_grid(bnuc,bnuc));
65           return concatenate_grid(bnucel,belec);
66         } else {
67           return utils::get_grid(rmax,num_el,igrid,zexp);
68         }
69       }
70 
offcenter_nuclear_grid(int num_el0,int Zm,int Zlr,double Rhalf,int num_el,double rmax,int igrid,double zexp)71       arma::vec offcenter_nuclear_grid(int num_el0, int Zm, int Zlr, double Rhalf, int num_el, double rmax, int igrid, double zexp) {
72         // First boundary at
73         int b0used = (Zm != 0);
74         double b0=Zm*Rhalf/(Zm+Zlr);
75         // Second boundary at
76         int b1used = (Zlr != 0);
77         double b1=Rhalf;
78         // Last boundary at
79         double b2=rmax;
80 
81         printf("b0 = %e, b0used = %i\n",b0,b0used);
82         printf("b1 = %e, b1used = %i\n",b1,b1used);
83         printf("b2 = %e\n",b2);
84 
85         // Get grids
86         arma::vec bval0, bval1;
87         if(b0used) {
88           // 0 to b0
89           bval0=utils::get_grid(b0,num_el0,igrid,zexp);
90         }
91         if(b1used) {
92           // b0 to b1
93 
94           // Reverse grid to get tighter spacing around nucleus
95           bval1=-arma::reverse(utils::get_grid(b1-b0,num_el0,igrid,zexp));
96           bval1+=arma::ones<arma::vec>(bval1.n_elem)*(b1-b0);
97           // Assert numerical exactness
98           bval1(0)=0.0;
99           bval1(bval1.n_elem-1)=b1-b0;
100         }
101         arma::vec bval2=utils::get_grid(b2-b1,num_el,igrid,zexp);
102 
103         arma::vec bval;
104         if(b0used && b1used) {
105           bval=concatenate_grid(bval0,bval1);
106         } else if(b0used) {
107           bval=bval0;
108         } else if(b1used) {
109           bval=bval1;
110         }
111         if(b0used || b1used) {
112           bval=concatenate_grid(bval,bval2);
113         } else {
114           bval=bval2;
115         }
116 
117         return bval;
118       }
119 
form_grid(modelpotential::nuclear_model_t model,double Rrms,int Nelem,double Rmax,int igrid,double zexp,int Nelem0,int igrid0,double zexp0,int Z,int Zl,int Zr,double Rhalf)120       arma::vec form_grid(modelpotential::nuclear_model_t model, double Rrms, int Nelem, double Rmax, int igrid, double zexp, int Nelem0, int igrid0, double zexp0, int Z, int Zl, int Zr, double Rhalf) {
121         // Construct the radial basis
122         arma::vec bval;
123         if(model != modelpotential::POINT_NUCLEUS) {
124           printf("Finite-nucleus grid\n");
125 
126           if(Zl != 0 || Zr != 0)
127             throw std::logic_error("Off-center nuclei not supported in finite nucleus mode!\n");
128 
129           double rnuc;
130           if(model == modelpotential::HOLLOW_NUCLEUS)
131             rnuc = Rrms;
132           else if(model == modelpotential::SPHERICAL_NUCLEUS)
133             rnuc = sqrt(5.0/3.0)*Rrms;
134           else if(model == modelpotential::GAUSSIAN_NUCLEUS)
135             rnuc = 3*Rrms;
136           else
137             throw std::logic_error("Nuclear grid not handled!\n");
138 
139           bval=atomic::basis::finite_nuclear_grid(Nelem,Rmax,igrid,zexp,Nelem0,rnuc,igrid0,zexp0);
140 
141         } else if(Zl != 0 || Zr != 0) {
142           printf("Off-center grid\n");
143           bval=atomic::basis::offcenter_nuclear_grid(Nelem0,Z,std::max(Zl,Zr),Rhalf,Nelem,Rmax,igrid,zexp);
144         } else {
145           printf("Normal grid\n");
146           bval=atomic::basis::normal_grid(Nelem,Rmax,igrid,zexp);
147         }
148 
149         bval.print("Grid");
150 
151         return bval;
152       }
153 
angular_basis(int lmax,int mmax,arma::ivec & lval,arma::ivec & mval)154       void angular_basis(int lmax, int mmax, arma::ivec & lval, arma::ivec & mval) {
155         size_t nang=0;
156         for(int l=0;l<=lmax;l++)
157           nang+=2*std::min(mmax,l)+1;
158 
159         // Allocate memory
160         lval.zeros(nang);
161         mval.zeros(nang);
162 
163         // Store values
164         size_t iang=0;
165         for(int mabs=0;mabs<=mmax;mabs++)
166           for(int l=mabs;l<=lmax;l++) {
167             lval(iang)=l;
168             mval(iang)=mabs;
169             iang++;
170             if(mabs>0) {
171               lval(iang)=l;
172               mval(iang)=-mabs;
173               iang++;
174             }
175           }
176         if(iang!=lval.n_elem)
177           throw std::logic_error("Error.\n");
178       }
179     }
180   }
181 }
182