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 17 #include "twodquadrature.h" 18 #include "chebyshev.h" 19 #include "../general/lcao.h" 20 #include "../general/model_potential.h" 21 #include "utils.h" 22 23 namespace helfem { 24 namespace diatomic { 25 namespace twodquad { TwoDGridWorker()26 TwoDGridWorker::TwoDGridWorker() { 27 } 28 TwoDGridWorker(const helfem::diatomic::basis::TwoDBasis * basp_,int lang)29 TwoDGridWorker::TwoDGridWorker(const helfem::diatomic::basis::TwoDBasis * basp_, int lang) : basp(basp_) { 30 // Get angular grid 31 chebyshev::chebyshev(lang,cth,wang); 32 } 33 ~TwoDGridWorker()34 TwoDGridWorker::~TwoDGridWorker() { 35 } 36 compute_bf(size_t iel,size_t irad,int m_)37 void TwoDGridWorker::compute_bf(size_t iel, size_t irad, int m_) { 38 // Store m 39 m=m_; 40 // Update function list 41 bf_ind=basp->bf_list(iel,m); 42 43 // Get radial weights. Only do one radial quadrature point at a 44 // time, since this is an easy way to save a lot of memory. 45 r.zeros(1); 46 r(0)=basp->get_r(iel)(irad); 47 wrad.zeros(1); 48 wrad(0)=basp->get_wrad(iel)(irad); 49 50 double Rhalf(basp->get_Rhalf()); 51 52 // Calculate helpers 53 arma::vec shmu(arma::sinh(r)); 54 55 arma::vec sth(cth.n_elem); 56 for(size_t ia=0;ia<cth.n_elem;ia++) 57 sth(ia)=sqrt(1.0 - cth(ia)*cth(ia)); 58 59 // Update total weights 60 wtot.zeros(wrad.n_elem*wang.n_elem); 61 for(size_t ia=0;ia<wang.n_elem;ia++) 62 for(size_t ir=0;ir<wrad.n_elem;ir++) { 63 size_t idx=ia*wrad.n_elem+ir; 64 // sin(th) is already contained within wang, but we don't want to divide by it since it may be zero. Phi integrals yield 2 pi 65 wtot(idx)=2.0*M_PI*wang(ia)*wrad(ir)*std::pow(Rhalf,3)*shmu(ir)*(std::pow(shmu(ir),2)+std::pow(sth(ia),2)); 66 } 67 68 // Compute basis function values 69 bf.zeros(bf_ind.n_elem,wtot.n_elem); 70 // Loop over angular grid 71 #ifdef _OPENMP 72 #pragma omp parallel for 73 #endif 74 for(size_t ia=0;ia<cth.n_elem;ia++) { 75 // Evaluate basis functions at angular point 76 arma::mat abf(basp->eval_bf(iel, irad, cth(ia), m)); 77 if(abf.n_cols != bf_ind.n_elem) { 78 std::ostringstream oss; 79 oss << "Mismatch! Have " << bf_ind.n_elem << " basis function indices but " << abf.n_cols << " basis functions!\n"; 80 throw std::logic_error(oss.str()); 81 } 82 // Store functions 83 bf.cols(ia*wrad.n_elem,(ia+1)*wrad.n_elem-1)=arma::trans(abf); 84 } 85 } 86 model_potential(const modelpotential::ModelPotential * p1,const modelpotential::ModelPotential * p2)87 void TwoDGridWorker::model_potential(const modelpotential::ModelPotential * p1, const modelpotential::ModelPotential * p2) { 88 double Rhalf(basp->get_Rhalf()); 89 arma::vec chmu(arma::cosh(r)); 90 91 itg.zeros(1,wtot.n_elem); 92 for(size_t ia=0;ia<wang.n_elem;ia++) 93 for(size_t ir=0;ir<wrad.n_elem;ir++) { 94 size_t idx=ia*wrad.n_elem+ir; 95 96 double r1=Rhalf*(chmu(ir) + cth(ia)); 97 double r2=Rhalf*(chmu(ir) - cth(ia)); 98 99 double V1(p1->V(r1)); 100 double V2(p2->V(r2)); 101 if(std::isnormal(V1)) 102 itg(idx)+=V1; 103 if(std::isnormal(V2)) 104 itg(idx)+=V2; 105 } 106 } 107 unit_pot()108 void TwoDGridWorker::unit_pot() { 109 itg.ones(1,wtot.n_elem); 110 } 111 gto(int l,const arma::vec & expn,probe_t p)112 void TwoDGridWorker::gto(int l, const arma::vec & expn, probe_t p) { 113 double Rhalf(basp->get_Rhalf()); 114 arma::vec chmu(arma::cosh(r)); 115 116 itg.zeros(expn.n_elem,wtot.n_elem); 117 if(p==PROBE_LEFT) { 118 for(size_t ia=0;ia<wang.n_elem;ia++) 119 for(size_t ir=0;ir<wrad.n_elem;ir++) { 120 size_t idx=ia*wrad.n_elem+ir; 121 122 double ra(Rhalf*(chmu(ir) + cth(ia))); 123 for(size_t ix=0;ix<expn.n_elem;ix++) 124 itg(ix,idx)=lcao::radial_GTO(ra,l,expn(ix)); 125 } 126 127 } else if(p==PROBE_RIGHT) { 128 for(size_t ia=0;ia<wang.n_elem;ia++) 129 for(size_t ir=0;ir<wrad.n_elem;ir++) { 130 size_t idx=ia*wrad.n_elem+ir; 131 132 double rb(Rhalf*(chmu(ir) - cth(ia))); 133 for(size_t ix=0;ix<expn.n_elem;ix++) 134 itg(ix,idx)=lcao::radial_GTO(rb,l,expn(ix)); 135 } 136 137 } else if(p==PROBE_MIDDLE) { 138 for(size_t ia=0;ia<wang.n_elem;ia++) 139 for(size_t ir=0;ir<wrad.n_elem;ir++) { 140 size_t idx=ia*wrad.n_elem+ir; 141 142 double rc(Rhalf*sqrt(std::pow(chmu(ir),2) + std::pow(cth(ia),2) -1.0)); 143 for(size_t ix=0;ix<expn.n_elem;ix++) 144 itg(ix,idx)=lcao::radial_GTO(rc,l,expn(ix)); 145 } 146 } 147 148 // Assure normalization 149 itg/=sqrt(4.0*M_PI); 150 } 151 sto(int l,const arma::vec & expn,probe_t p)152 void TwoDGridWorker::sto(int l, const arma::vec & expn, probe_t p) { 153 double Rhalf(basp->get_Rhalf()); 154 arma::vec chmu(arma::cosh(r)); 155 156 itg.zeros(expn.n_elem,wtot.n_elem); 157 if(p==PROBE_LEFT) { 158 for(size_t ia=0;ia<wang.n_elem;ia++) 159 for(size_t ir=0;ir<wrad.n_elem;ir++) { 160 size_t idx=ia*wrad.n_elem+ir; 161 double ra(Rhalf*(chmu(ir) + cth(ia))); 162 for(size_t ix=0;ix<expn.n_elem;ix++) 163 itg(ix,idx)=lcao::radial_STO(ra,l,expn(ix)); 164 } 165 166 } else if(p==PROBE_RIGHT) { 167 for(size_t ia=0;ia<wang.n_elem;ia++) 168 for(size_t ir=0;ir<wrad.n_elem;ir++) { 169 size_t idx=ia*wrad.n_elem+ir; 170 double rb(Rhalf*(chmu(ir) - cth(ia))); 171 for(size_t ix=0;ix<expn.n_elem;ix++) 172 itg(ix,idx)=lcao::radial_STO(rb,l,expn(ix)); 173 } 174 175 } else if(p==PROBE_MIDDLE) { 176 for(size_t ia=0;ia<wang.n_elem;ia++) 177 for(size_t ir=0;ir<wrad.n_elem;ir++) { 178 size_t idx=ia*wrad.n_elem+ir; 179 double rc(Rhalf*sqrt(std::pow(chmu(ir),2) + std::pow(cth(ia),2) -1.0)); 180 for(size_t ix=0;ix<expn.n_elem;ix++) 181 itg(ix,idx)=lcao::radial_STO(rc,l,expn(ix)); 182 } 183 } 184 185 // Assure normalization 186 itg/=sqrt(4.0*M_PI); 187 } 188 eval_pot(arma::mat & Vo) const189 void TwoDGridWorker::eval_pot(arma::mat & Vo) const { 190 if(itg.n_rows != 1) 191 throw std::logic_error("Should only have one column in integrand!\n"); 192 Vo.submat(bf_ind,bf_ind)+=bf*arma::diagmat(itg%wtot)*arma::trans(bf); 193 } 194 eval_proj(arma::mat & Vo) const195 void TwoDGridWorker::eval_proj(arma::mat & Vo) const { 196 Vo.cols(bf_ind)+=itg*arma::diagmat(wtot)*arma::trans(bf); 197 } 198 eval_proj_overlap(arma::mat & Vo) const199 void TwoDGridWorker::eval_proj_overlap(arma::mat & Vo) const { 200 Vo+=itg*arma::diagmat(wtot)*arma::trans(itg); 201 } 202 TwoDGrid()203 TwoDGrid::TwoDGrid() { 204 } 205 TwoDGrid(const helfem::diatomic::basis::TwoDBasis * basp_,int lang_)206 TwoDGrid::TwoDGrid(const helfem::diatomic::basis::TwoDBasis * basp_, int lang_) : basp(basp_), lang(lang_) { 207 } 208 ~TwoDGrid()209 TwoDGrid::~TwoDGrid() { 210 } 211 model_potential(const modelpotential::ModelPotential * p1,const modelpotential::ModelPotential * p2)212 arma::mat TwoDGrid::model_potential(const modelpotential::ModelPotential * p1, const modelpotential::ModelPotential * p2) { 213 arma::mat H; 214 H.zeros(basp->Ndummy(),basp->Ndummy()); 215 216 // Get unique m values in basis set 217 arma::ivec muni(basp->get_mval()); 218 muni=muni(arma::find_unique(muni)); 219 { 220 TwoDGridWorker grid(basp,lang); 221 222 for(size_t im=0;im<muni.n_elem;im++) { 223 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 224 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 225 grid.compute_bf(iel,irad,muni(im)); 226 grid.model_potential(p1, p2); 227 grid.eval_pot(H); 228 } 229 } 230 } 231 } 232 233 H=basp->remove_boundaries(H); 234 235 return H; 236 } 237 overlap()238 arma::mat TwoDGrid::overlap() { 239 arma::mat S; 240 S.zeros(basp->Ndummy(),basp->Ndummy()); 241 242 // Get unique m values in basis set 243 arma::ivec muni(basp->get_mval()); 244 muni=muni(arma::find_unique(muni)); 245 { 246 TwoDGridWorker grid(basp,lang); 247 248 for(size_t im=0;im<muni.n_elem;im++) { 249 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 250 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 251 grid.compute_bf(iel,irad,muni(im)); 252 grid.unit_pot(); 253 grid.eval_pot(S); 254 } 255 } 256 } 257 } 258 259 S=basp->remove_boundaries(S); 260 261 return S; 262 } 263 gto_projection(int l,int m,const arma::vec & expn,probe_t p)264 arma::mat TwoDGrid::gto_projection(int l, int m, const arma::vec & expn, probe_t p) { 265 arma::mat S; 266 S.zeros(expn.n_elem,basp->Ndummy()); 267 TwoDGridWorker grid(basp,lang); 268 269 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 270 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 271 grid.compute_bf(iel,irad,m); 272 grid.gto(l, expn, p); 273 grid.eval_proj(S); 274 } 275 } 276 277 S=S.cols(basp->pure_indices()); 278 279 return S; 280 } 281 gto_overlap(int l,int m,const arma::vec & expn,probe_t p)282 arma::mat TwoDGrid::gto_overlap(int l, int m, const arma::vec & expn, probe_t p) { 283 arma::mat S; 284 S.zeros(expn.n_elem,expn.n_elem); 285 TwoDGridWorker grid(basp,lang); 286 287 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 288 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 289 grid.compute_bf(iel,irad,m); 290 grid.gto(l, expn, p); 291 grid.eval_proj_overlap(S); 292 } 293 } 294 295 return S; 296 } 297 sto_projection(int l,int m,const arma::vec & expn,probe_t p)298 arma::mat TwoDGrid::sto_projection(int l, int m, const arma::vec & expn, probe_t p) { 299 arma::mat S; 300 S.zeros(expn.n_elem,basp->Ndummy()); 301 TwoDGridWorker grid(basp,lang); 302 303 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 304 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 305 grid.compute_bf(iel,irad,m); 306 grid.sto(l, expn, p); 307 grid.eval_proj(S); 308 } 309 } 310 311 S=S.cols(basp->pure_indices()); 312 313 return S; 314 } 315 sto_overlap(int l,int m,const arma::vec & expn,probe_t p)316 arma::mat TwoDGrid::sto_overlap(int l, int m, const arma::vec & expn, probe_t p) { 317 arma::mat S; 318 S.zeros(expn.n_elem,expn.n_elem); 319 TwoDGridWorker grid(basp,lang); 320 321 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 322 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 323 grid.compute_bf(iel,irad,m); 324 grid.sto(l, expn, p); 325 grid.eval_proj_overlap(S); 326 } 327 } 328 329 return S; 330 } 331 } 332 } 333 } 334