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