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