-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathintegrals.h
139 lines (110 loc) · 5.06 KB
/
integrals.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
//
// Created by Kshitij Surjuse on 5/9/23.
//
#ifndef HYLLERAAS_INTEGRALS_H
#define HYLLERAAS_INTEGRALS_H
#include "basis.h"
#include "k.h"
#include <Eigen/Dense>
namespace hylleraas {
typedef long double numeric_type;
typedef Eigen::Matrix<numeric_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Matrix;
/// computes overlap integral between basis function bf1 anf bf2
numeric_type S_ij(const basis &bf1, const basis &bf2) {
return K(bf1.n_ + bf2.n_, bf1.l_ + bf2.l_, bf1.m_ + bf2.m_,
bf1.alpha_, bf1.beta_, bf1.gamma_);
}
/// computes nu-ele attraction integral between basis function bf1 anf bf2
numeric_type Vne_ij(const basis &bf1, const basis &bf2) {
numeric_type one = 1.0;
auto k1 = K(bf1.n_ + bf2.n_ - one, bf1.l_ + bf2.l_, bf1.m_ + bf2.m_,
bf1.alpha_, bf1.beta_, bf1.gamma_);
auto k2 = K(bf1.n_ + bf2.n_, bf1.l_ + bf2.l_ - one, bf1.m_ + bf2.m_,
bf1.alpha_, bf1.beta_, bf1.gamma_);
return -one * bf1.Z_ * (k1 + k2);
}
/// computes ele-ele repulsion integral between basis function bf1 anf bf2
numeric_type Vee_ij(const basis &bf1, const basis &bf2) {
numeric_type one = 1.0;
return K(bf1.n_ + bf2.n_, bf1.l_ + bf2.l_, bf1.m_ + bf2.m_ - one,
bf1.alpha_, bf1.beta_, bf1.gamma_);
}
/// computes kinetic energy repulsion integral between basis function bf1 anf bf2
numeric_type T_ij(const basis &bf1, const basis &bf2) {
auto alpha = bf1.alpha_;
auto beta = bf1.beta_;
auto gamma = bf1.gamma_;
auto nj = bf2.n_;
auto lj = bf2.l_;
auto mj = bf2.m_;
auto screen_prefactor = [](numeric_type prefactor) {
return prefactor !=0.0;
};
auto Knlm = [&](numeric_type prefactor, numeric_type n, numeric_type l, numeric_type m) {
if(screen_prefactor(prefactor)){
return prefactor * K(bf1.n_ + bf2.n_ + n, bf1.l_ + bf2.l_ + l, bf1.m_ + bf2.m_ + m, alpha,
beta, gamma);
}
else{
numeric_type zero = 0.0;
return zero;
}
};
numeric_type Tij = 0.;
numeric_type zero, one, two, half,one_fourth,one_eighth ;
zero = 0.;
one = 1.;
two = 2.;
half = 0.5;
one_fourth = 0.25;
one_eighth = 0.125;
Tij = -one_eighth* (pow(alpha, two) + pow(beta, two) + pow(gamma, two)) * S_ij(bf1, bf2);
Tij += Knlm((half * nj * alpha + half * alpha), -one, zero, zero);
Tij += Knlm((half * lj * beta + half * alpha) ,zero, -one, zero);
Tij += Knlm((mj * gamma + gamma ),zero, zero, -one);
Tij -= Knlm((half * nj * (nj - one) + nj),-two, zero, zero);
Tij -= Knlm((half * lj * (lj - one) + lj),zero, -two, zero);
Tij -= Knlm((mj * (mj - one) + two * mj),zero, zero, -two);
Tij -= Knlm((one_eighth * alpha * gamma),-one, zero, one)
+ Knlm((one_eighth * alpha * gamma),one, zero, -one)
- Knlm((one_eighth * alpha * gamma),-one, two, -one);
Tij -= Knlm((one_eighth * beta * gamma),zero, -one, one)
+ Knlm((one_eighth * beta * gamma),zero, one, -one)
- Knlm((one_eighth * beta * gamma),two, -one, -one);
Tij += Knlm((one_fourth * nj * gamma),-two, zero, one)
+ Knlm((one_fourth * nj * gamma),zero, zero, -one)
- Knlm((one_fourth * nj * gamma),-two, two, -one);
Tij += Knlm((one_fourth * mj * alpha),-one, zero, zero)
+ Knlm((one_fourth * mj * alpha),one, zero, -two)
- Knlm((one_fourth * mj * alpha),-one, two, -two);
Tij -= Knlm((half * nj * mj),-two, zero, zero)
+ Knlm((half * nj * mj),zero, zero, -two)
- Knlm((half * nj * mj),-two, two, -two);
Tij += Knlm((one_fourth * lj * gamma),zero, -two, one)
+ Knlm((one_fourth * lj * gamma),zero, zero, -one)
- Knlm((one_fourth * lj * gamma),two, -two, -one);
Tij += Knlm((one_fourth * mj * beta),zero, -one , zero)
+ Knlm((one_fourth * mj * beta),zero, one, -two)
- Knlm((one_fourth * mj * beta),two,-one, -two);
Tij -= Knlm((half * lj * mj),zero, -two, zero)
+ Knlm((half * lj * mj),zero, zero, -two)
- Knlm((half * lj * mj),two, -two, -two);
return Tij;
}
/// computes integral matrix given a vector of basis functions
/// @param2 a reference to the function that computes the elements of the matrix
template<typename integral_type>
Matrix compute_integral(const std::vector<basis> &bfn, const integral_type &integral_ij) {
auto n = bfn.size();
Matrix matrix(n, n);
numeric_type zero = 0.0;
matrix.fill(zero);
for (auto i = 0; i < n; ++i) {
for (auto j = 0; j < n; ++j) {
matrix(i, j) = integral_ij(bfn[i], bfn[j]);
}
}
return matrix;
}
}
#endif //HYLLERAAS_INTEGRALS_H