Hex: src/radial.h Source File

| Hex 2.2 Hydrogen-electron collision solver

radial.h
Go to the documentation of this file.
1 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2 // //
3 // / / / / __ \ \ / / //
4 // / /__ / / / _ \ \ \/ / //
5 // / ___ / | |/_/ / /\ \ //
6 // / / / / \_\ / / \ \ //
7 // //
8 // //
9 // Copyright (c) 2017, Jakub Benda, Charles University in Prague //
10 // //
11 // MIT License: //
12 // //
13 // Permission is hereby granted, free of charge, to any person obtaining a //
14 // copy of this software and associated documentation files (the "Software"), //
15 // to deal in the Software without restriction, including without limitation //
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense, //
17 // and/or sell copies of the Software, and to permit persons to whom the //
18 // Software is furnished to do so, subject to the following conditions: //
19 // //
20 // The above copyright notice and this permission notice shall be included //
21 // in all copies or substantial portions of the Software. //
22 // //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS //
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, //
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE //
26 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, //
27 // WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF //
28 // OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. //
29 // //
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31
32 #ifndef HEX_ECS_RADIAL_H
33 #define HEX_ECS_RADIAL_H
34
35 // --------------------------------------------------------------------------------- //
36
37 #include <cmath>
38 #include <functional>
39 #include <vector>
40
41 // --------------------------------------------------------------------------------- //
42
43 #include "hex-arrays.h"
44 #include "hex-symbandmatrix.h"
45 #include "hex-special.h"
46
47 // --------------------------------------------------------------------------------- //
48
49 #include "bspline.h"
50 #include "gauss.h"
51 #include "io.h"
52 #include "parallel.h"
53
54 // --------------------------------------------------------------------------------- //
55
56 #define EXPANSION_QUADRATURE_POINTS 20
57
58 // --------------------------------------------------------------------------------- //
59
60 // Complex -> Complex function.
61 typedef std::function<Complex(Complex)> CCFunction;
62
63 // Complex x Complex -> Complex function.
64 typedef std::function<Complex(Complex,Complex)> C2CFunction;
65
66 // --------------------------------------------------------------------------------- //
67
89 class RadialIntegrals
90 {
91 public:
92
99 RadialIntegrals
100 (
101 Bspline const & bspline_x,
102 Bspline const & bspline_y,
103 int Nlambdas
104 );
105
114 void setupOneElectronIntegrals (bool sharedscratch = false, bool iammaster = true);
115
122 void setupOneElectronIntegrals (Parallel const & par, CommandLine const & cmd);
123
130 void setupTwoElectronIntegrals (Parallel const & par, CommandLine const & cmd);
131
137 void verbose (bool v) { verbose_ = v; }
138
145 int maxlambda () const { return Nlambdas_ - 1; }
146
147 // compute overlap matrix of radial function
148 Complex computeOverlapMatrixElement_iknot (Bspline const & bspline, GaussLegendre const & g, int i, int j, CCFunction func, int iknot) const;
149 Complex computeOverlapMatrixElement (Bspline const & bspline, GaussLegendre const & g, int i, int j, CCFunction func) const;
150 SymBandMatrix<Complex> computeOverlapMatrix (Bspline const & bspline, CCFunction func) const;
151
160 static CsrMatrix<LU_int_t,Complex> computeOverlapMatrix
161 (
162 Bspline const & bspline1,
163 Bspline const & bspline2,
164 Real R_left,
165 Real R_right
166 );
167
177 Complex computeD_iknot
178 (
179 Bspline const & bspline,
180 GaussLegendre const & g,
181 int i, int j, int iknot
182 ) const;
183
191 Complex computeD
192 (
193 Bspline const & bspline,
194 GaussLegendre const & g,
195 int i, int j
196 ) const;
197
210 Complex computeM_iknot
211 (
212 Bspline const & bspline,
213 GaussLegendre const & g,
214 int a, int i, int j,
215 int iknot,
216 Real rmin, Real rmax,
217 Real scale
218 ) const;
219
229 Complex computeM
230 (
231 Bspline const & bspline,
232 GaussLegendre const & g,
233 int a, int i, int j,
234 Real rmin, Real rmax,
235 bool scale = false
236 ) const;
237
247 Complex computeM
248 (
249 Bspline const & bspline,
250 GaussLegendre const & g,
251 int a, int i, int j,
252 int begin_knot, int end_knot,
253 Real rmin, Real rmax,
254 bool scale = false
255 ) const;
256
268 cArray computeMi
269 (
270 Bspline const & bspline,
271 GaussLegendre const & g,
272 int a,
273 Real rmin, Real rmax
274 ) const;
275
283 Complex computeR (int lambda, int i, int j, int k, int l) const;
284
285 void diagonalR (int lambda);
286 void coupledR (int lambda);
287
321 Complex computeRtri
322 (
323 int lambda,
324 Bspline const & xspline, GaussLegendre const & xg,
325 int k, int l, int iknotx, Real tx1, Real tx2,
326 Bspline const & yspline, GaussLegendre const & yg,
327 int m, int n, int iknoty, Real ty1, Real ty2
328 ) const;
329
336 SymBandMatrix<Complex> calc_R_tr_dia_block
337 (
338 unsigned lambda,
339 int i, int k
340 ) const;
341
352 void apply_R_matrix
353 (
354 unsigned int lambda,
355 Complex a, const cArrayView src,
356 Complex b, cArrayView dst,
357 int x_row_limit = -1,
358 int y_row_limit = -1
359 ) const;
360
368 cArray overlap
369 (
370 Bspline const & bspline,
371 GaussLegendre const & g,
372 CCFunction funct
373 ) const;
374
384 cArray overlap
385 (
386 Bspline const & bsplinex,
387 GaussLegendre const & gx,
388 Bspline const & bspliney,
389 GaussLegendre const & gy,
390 C2CFunction funct
391 ) const;
392
403 static cArray overlapP
404 (
405 Bspline const & bspline,
406 GaussLegendre const & g,
407 Real Z, int n, int l
408 );
409
419 cArray overlapj
420 (
421 Bspline const & bspline,
422 GaussLegendre const & g,
423 int maxell,
424 const rArrayView vk,
425 bool fast_bessel = false
426 ) const;
427
428 // Return reference to the B-spline object.
429 Bspline const & bspline () const { return bspline_x_; }
430 Bspline const & bspline_x () const { return bspline_x_; }
431 Bspline const & bspline_y () const { return bspline_y_; }
432
433 // Return the Gauss-Legendre integrator object.
434 GaussLegendre const & gaussleg () const { return g_x_; }
435 GaussLegendre const & gaussleg_x () const { return g_x_; }
436 GaussLegendre const & gaussleg_y () const { return g_y_; }
437
438 #define OneElectronMatrixAccessors(M) \
439 SymBandMatrix<Complex> const & M () const { return M##_x_; } \
440 SymBandMatrix<Complex> const & M##_x () const { return M##_x_; } \
441 SymBandMatrix<Complex> const & M##_y () const { return M##_y_; } \
442 Complex M##_x (std::size_t i, std::size_t j) const { return M##_x_(i,j); } \
443 Complex M##_y (std::size_t i, std::size_t j) const { return M##_y_(i,j); } \
444
445 // Access the precomputed one-electron overlap matrices.
446 OneElectronMatrixAccessors(D)
447 OneElectronMatrixAccessors(S)
448 OneElectronMatrixAccessors(Mm1)
449 OneElectronMatrixAccessors(Mm2)
450
451 #define OneElectronMatrixArrayAccessors(M) \
452 SymBandMatrix<Complex> const & M (int L) const { return M##_x_[L]; } \
453 SymBandMatrix<Complex> const & M##_x (int L) const { return M##_x_[L]; } \
454 SymBandMatrix<Complex> const & M##_y (int L) const { return M##_y_[L]; } \
455 Complex M##_x (int L, std::size_t i, std::size_t j) const { return M##_x_[L](i,j); } \
456 Complex M##_y (int L, std::size_t i, std::size_t j) const { return M##_y_[L](i,j); } \
457
458 // Access the precomputed scaled full integral moments of order L / -L-1.
459 OneElectronMatrixArrayAccessors(Mtr_L)
460 OneElectronMatrixArrayAccessors(Mtr_mLm1)
461
462 #define OneElectronPartialMatrixAccessors(M) \
463 cArrayView M##_x (int L = -1) const \
464 { \
465 if (L < 0) return M##_x_; \
466 std::size_t mi_size = bspline_x_.Nspline() * (2 * bspline_x_.order() + 1) * (bspline_x_.order() + 1); \
467 return cArrayView (M##_x_, L * mi_size, mi_size); \
468 } \
469 cArrayView M##_y (int L = -1) const \
470 { \
471 if (L < 0) return M##_y_; \
472 std::size_t mi_size = bspline_y_.Nspline() * (2 * bspline_y_.order() + 1) * (bspline_y_.order() + 1); \
473 return cArrayView (M##_y_, L * mi_size, mi_size); \
474 }
475
476 // Access the precomputed scaled partial overlap matrices of order L / -L-1.
477 OneElectronPartialMatrixAccessors(Mitr_L)
478 OneElectronPartialMatrixAccessors(Mitr_mLm1)
479
480 // Return reference to the precomputed matrix of two-electron integrals for given multipole.
481 BlockSymBandMatrix<Complex> const & R_tr_dia (unsigned i) const { return R_tr_dia_[i]; }
482 CsrMatrix<LU_int_t,Complex> const & R_coupled (unsigned i) const { return R_coupled_[i]; }
483
484 private:
485
486 // B-spline environment
487 Bspline bspline_x_;
488 Bspline bspline_y_;
489
490 // Asymptotic radiuses
491 // - These are somewhat arbitrary radiuses that restrict position dependence of the Coulombic
492 // interaction. When either of the coordinates exceeds the limit, it is assumed to be exactly equal
493 // to the bound instead, making the potential constant with respect to further change of that
494 // coordinate. These bounds are sufficiently far from the real grid to avoid corruption of the
495 // physical solution. In single-domain problem it just stops position dependence of the potential
496 // for large distances, which is irrelevant. The importance lies in the multi-domain solution,
497 // where the potential may diverge in prepended complex grid.
498 Real rxmin_, rymin_, rxmax_, rymax_;
499
500 // Gauss-Legendre integrator
501 GaussLegendre g_x_;
502 GaussLegendre g_y_;
503
504 // one-electron moment and overlap matrices
505 SymBandMatrix<Complex> D_x_, S_x_, Mm1_x_, Mm2_x_;
506 SymBandMatrix<Complex> D_y_, S_y_, Mm1_y_, Mm2_y_;
507
508 // one-electron full integral moments for various orders (used to calculate R-integrals)
509 std::vector<SymBandMatrix<Complex>> Mtr_L_x_, Mtr_mLm1_x_;
510 std::vector<SymBandMatrix<Complex>> Mtr_L_y_, Mtr_mLm1_y_;
511
512 // partial one-electron integral moments for various orders (used to calculate R-integrals)
513 cArray Mitr_L_x_, Mitr_mLm1_x_;
514 cArray Mitr_L_y_, Mitr_mLm1_y_;
515
516 // two-electron integral matrices
517 std::vector<BlockSymBandMatrix<Complex>> R_tr_dia_;
518 std::vector<CsrMatrix<LU_int_t,Complex>> R_coupled_;
519
520 // verbose output
521 bool verbose_;
522
523 // number of multipole matrices
524 int Nlambdas_;
525 };
526
527 #endif // HEX_ECS_RADIAL_H
RadialIntegrals::computeM_iknot
Complex computeM_iknot(Bspline const &bspline, GaussLegendre const &g, int a, int i, int j, int iknot, Real rmin, Real rmax, Real scale) const
Partial integral moment.
Definition: radial.cpp:409
RadialIntegrals
Class that calculates and manages the radial integrals for given B-spline bases.
Definition: radial.h:89
RadialIntegrals::gaussleg
GaussLegendre const & gaussleg() const
Definition: radial.h:434
RadialIntegrals::computeOverlapMatrixElement
Complex computeOverlapMatrixElement(Bspline const &bspline, GaussLegendre const &g, int i, int j, CCFunction func) const
Definition: radial.cpp:284
RadialIntegrals::computeMi
cArray computeMi(Bspline const &bspline, GaussLegendre const &g, int a, Real rmin, Real rmax) const
Partial integral moments.
Definition: radial.cpp:86
RadialIntegrals::R_coupled
CsrMatrix< LU_int_t, Complex > const & R_coupled(unsigned i) const
Definition: radial.h:482
RadialIntegrals::overlapP
static cArray overlapP(Bspline const &bspline, GaussLegendre const &g, Real Z, int n, int l)
Compute P-overlaps.
Definition: radial.cpp:962
bspline.h
RadialIntegrals::setupTwoElectronIntegrals
void setupTwoElectronIntegrals(Parallel const &par, CommandLine const &cmd)
Calculate the two-electron integral matrix.
Definition: radial.cpp:628
RadialIntegrals::computeD
Complex computeD(Bspline const &bspline, GaussLegendre const &g, int i, int j) const
Derivative overlap.
Definition: radial.cpp:226
Bspline
B-spline environment.
Definition: bspline.h:57
RadialIntegrals::overlap
cArray overlap(Bspline const &bspline, GaussLegendre const &g, CCFunction funct) const
Compute B-spline overlaps of arbitrary one-dimensional function.
Definition: radial.cpp:840
OneElectronPartialMatrixAccessors
#define OneElectronPartialMatrixAccessors(M)
CsrMatrix
Definition: luft.h:42
GaussLegendre
Gauss-Legendre quadrature.
Definition: gauss.h:79
CCFunction
std::function< Complex(Complex)> CCFunction
Definition: radial.h:61
RadialIntegrals::bspline_y
Bspline const & bspline_y() const
Definition: radial.h:431
RadialIntegrals::computeOverlapMatrixElement_iknot
Complex computeOverlapMatrixElement_iknot(Bspline const &bspline, GaussLegendre const &g, int i, int j, CCFunction func, int iknot) const
Definition: radial.cpp:247
io.h
RadialIntegrals::computeD_iknot
Complex computeD_iknot(Bspline const &bspline, GaussLegendre const &g, int i, int j, int iknot) const
Partial derivative overlap.
Definition: radial.cpp:187
RadialIntegrals::computeOverlapMatrix
SymBandMatrix< Complex > computeOverlapMatrix(Bspline const &bspline, CCFunction func) const
Definition: radial.cpp:305
RadialIntegrals::computeM
Complex computeM(Bspline const &bspline, GaussLegendre const &g, int a, int i, int j, Real rmin, Real rmax, bool scale=false) const
Integral moments.
Definition: radial.cpp:449
CommandLine
Command line parameters.
Definition: io.h:57
RadialIntegrals::setupOneElectronIntegrals
void setupOneElectronIntegrals(bool sharedscratch=false, bool iammaster=true)
Calculate one-electron integral matrices.
Definition: radial.cpp:551
RadialIntegrals::gaussleg_x
GaussLegendre const & gaussleg_x() const
Definition: radial.h:435
RadialIntegrals::maxlambda
int maxlambda() const
Maximal multipole moment.
Definition: radial.h:145
RadialIntegrals::computeRtri
Complex computeRtri(int lambda, Bspline const &xspline, GaussLegendre const &xg, int k, int l, int iknotx, Real tx1, Real tx2, Bspline const &yspline, GaussLegendre const &yg, int m, int n, int iknoty, Real ty1, Real ty2) const
Triangular R-integral.
Definition: slater.cpp:44
gauss.h
RadialIntegrals::overlapj
cArray overlapj(Bspline const &bspline, GaussLegendre const &g, int maxell, const rArrayView vk, bool fast_bessel=false) const
Compute j-overlaps.
Definition: radial.cpp:1021
C2CFunction
std::function< Complex(Complex, Complex)> C2CFunction
Definition: radial.h:64
parallel.h
RadialIntegrals::coupledR
void coupledR(int lambda)
Definition: slater.cpp:111
OneElectronMatrixArrayAccessors
#define OneElectronMatrixArrayAccessors(M)
RadialIntegrals::bspline_x
Bspline const & bspline_x() const
Definition: radial.h:430
Parallel
MPI info.
Definition: parallel.h:50
RadialIntegrals::bspline
Bspline const & bspline() const
Definition: radial.h:429
RadialIntegrals::RadialIntegrals
RadialIntegrals(Bspline const &bspline_x, Bspline const &bspline_y, int Nlambdas)
Constructor.
Definition: radial.cpp:56
RadialIntegrals::diagonalR
void diagonalR(int lambda)
RadialIntegrals::apply_R_matrix
void apply_R_matrix(unsigned int lambda, Complex a, const cArrayView src, Complex b, cArrayView dst, int x_row_limit=-1, int y_row_limit=-1) const
Multiply vector by matrix of two-electron integrals.
Definition: radial.cpp:768
RadialIntegrals::verbose
void verbose(bool v)
Verbosity control.
Definition: radial.h:137
RadialIntegrals::OneElectronMatrixAccessors
OneElectronMatrixAccessors(D) OneElectronMatrixAccessors(S) OneElectronMatrixAccessors(Mm1) OneElectronMatrixAccessors(Mm2) OneElectronMatrixArrayAccessors(Mtr_L) OneElectronMatrixArrayAccessors(Mtr_mLm1) OneElectronPartialMatrixAccessors(Mitr_L) OneElectronPartialMatrixAccessors(Mitr_mLm1) BlockSymBandMatrix< Complex > const &R_tr_dia(unsigned i) const
Definition: radial.h:446
RadialIntegrals::calc_R_tr_dia_block
SymBandMatrix< Complex > calc_R_tr_dia_block(unsigned lambda, int i, int k) const
Calculate particular sub-matrix of the radial integrals matrix.
Definition: radial.cpp:750
RadialIntegrals::computeR
Complex computeR(int lambda, int i, int j, int k, int l) const
Two-electron integral for multipole .
Definition: slater.cpp:213
RadialIntegrals::gaussleg_y
GaussLegendre const & gaussleg_y() const
Definition: radial.h:436