LIBINT 2.7.2
gaussoper.h
1/*
2 * Copyright (C) 2004-2021 Edward F. Valeev
3 *
4 * This file is part of Libint.
5 *
6 * Libint is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * Libint is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with Libint. If not, see <http://www.gnu.org/licenses/>.
18 *
19 */
20
21#ifndef _libint2_src_bin_libint_gaussoper_h_
22#define _libint2_src_bin_libint_gaussoper_h_
23
24#include <boost/type_traits/is_same.hpp>
25#include <braket.h>
26#include <prefactors.h>
27#include <global_macros.h>
28
29namespace libint2 {
30
32 template <class F, BraketType BKType>
33 LinearCombination<
34 SafePtr<DGVertex>,
35 BraketPair<F,BKType>
37 if (BKType == CBra || BKType == CKet)
38 throw std::logic_error("R12vec_dot_Nabla1 can only be applied to physicists brakets");
39
40 const char* zeta = (BKType == PBra) ? "zeta_A" : "zeta_B";
41 const char* XY = (BKType == PBra) ? "AC" : "BD";
42
43 const auto& f = bkt[0];
44 const auto& g = bkt[1];
46 ResultType result;
47
48 using namespace libint2::prefactor;
49 if (f.norm())
50 result += make_pair(Scalar((double)f.norm()),
52
53 const unsigned int nxyz = boost::is_same<F,CGF>::value ? 3 : 1;
54 for(unsigned int xyz=0; xyz<nxyz; ++xyz) {
55 using namespace libint2::algebra;
56 F _1 = unit<F>(xyz);
57 const F& fm1 = f - _1;
58 const F& gp1 = g + _1;
59 if (exists(fm1)) {
60 const double f_xyz = (double)(f.qn(xyz));
61 result += make_pair(Scalar(-1.0*f_xyz),
62 BraketPair<F,BKType>(fm1,gp1));
63 result += make_pair(Scalar(f_xyz)*Vector(XY)[xyz],
65 }
66 const F& fp1 = f + _1;
67 const F& fp2 = fp1 + _1;
68 result += make_pair(Scalar(-2.0) * Scalar(zeta),
70 result += make_pair(Scalar(2.0) * Scalar(zeta),
71 BraketPair<F,BKType>(fp1,gp1));
72 result += make_pair(Scalar(-2.0) * Scalar(zeta) * Vector(XY)[xyz],
74 }
75
76 return result;
77 }
78
80 template <class F, BraketType BKType>
81 LinearCombination<
82 SafePtr<DGVertex>,
83 BraketPair<F,BKType>
85 if (BKType == CBra || BKType == CKet)
86 throw std::logic_error("R12vec_dot_Nabla2 can only be applied to physicists brakets");
87
88 const char* zeta = (BKType == PBra) ? "zeta_C" : "zeta_D";
89 const char* XY = (BKType == PBra) ? "AC" : "BD";
90
91 const F& f = bkt[0];
92 const F& g = bkt[1];
94 ResultType result;
95
96 using namespace libint2::prefactor;
97 if (g.norm())
98 result += make_pair(Scalar(-1.0*g.norm()),
100
101 const unsigned int nxyz = boost::is_same<F,CGF>::value ? 3 : 1;
102 for(unsigned int xyz=0; xyz<nxyz; ++xyz) {
103 using namespace libint2::algebra;
104 F _1 = unit<F>(xyz);
105 const F& fp1 = f + _1;
106 const F& gm1 = g - _1;
107 if (exists(gm1)) {
108 const double g_xyz = (double)(g.qn(xyz));
109 result += make_pair(Scalar(g_xyz),
110 BraketPair<F,BKType>(fp1,gm1));
111 result += make_pair(Scalar(g_xyz)*Vector(XY)[xyz],
112 BraketPair<F,BKType>(f,gm1));
113 }
114 const F& gp1 = g + _1;
115 const F& gp2 = gp1 + _1;
116 result += make_pair(Scalar(-2.0) * Scalar(zeta),
117 BraketPair<F,BKType>(fp1,gp1));
118 result += make_pair(Scalar(2.0) * Scalar(zeta),
119 BraketPair<F,BKType>(f,gp2));
120 result += make_pair(Scalar(-2.0) * Scalar(zeta) * Vector(XY)[xyz],
121 BraketPair<F,BKType>(f,gp1));
122 }
123
124 return result;
125 }
126
128 template <class F, BraketType BKType>
129 LinearCombination<
130 SafePtr<DGVertex>,
131 BraketPair<F,BKType>
132 > Nabla1(const BraketPair<F,BKType>& bkt, int xyz) {
133 if (BKType == CBra || BKType == CKet)
134 throw std::logic_error("Nabla1 can only be applied to physicists brakets");
135
136 const char* zeta = (BKType == PBra) ? "zeta_A" : "zeta_B";
137
138 const F& f = bkt[0];
139 const F& g = bkt[1];
141 ResultType result;
142
143 using namespace libint2::prefactor;
144 using namespace libint2::algebra;
145 // if applying to shell pair, change angular momentum along 0
146 const unsigned int dir = boost::is_same<F,CGF>::value ? xyz : 0;
147 F _1 = unit<F>(dir);
148 const F& fm1 = f - _1;
149 if (exists(fm1)) {
150 const double f_xyz = (double)(f.qn(dir));
151 result += make_pair(Scalar(f_xyz),
152 BraketPair<F,BKType>(fm1,g));
153 }
154 const F& fp1 = f + _1;
155 result += make_pair(Scalar(-2.0) * Scalar(zeta),
156 BraketPair<F,BKType>(fp1,g));
157 return result;
158 }
159
161 template <class F, BraketType BKType>
162 LinearCombination<
163 SafePtr<DGVertex>,
164 BraketPair<F,BKType>
165 > Nabla2(const BraketPair<F,BKType>& bkt, int xyz) {
166 if (BKType == CBra || BKType == CKet)
167 throw std::logic_error("Nabla2 can only be applied to physicists brakets");
168
169 const char* zeta = (BKType == PBra) ? "zeta_C" : "zeta_D";
170
171 const F& f = bkt[0];
172 const F& g = bkt[1];
174 ResultType result;
175
176 using namespace libint2::prefactor;
177 using namespace libint2::algebra;
178 // if applying to shell pair, change angular momentum along 0
179 const unsigned int dir = boost::is_same<F,CGF>::value ? xyz : 0;
180 F _1 = unit<F>(dir);
181 const F& gm1 = g - _1;
182 if (exists(gm1)) {
183 const double g_xyz = (double)(g.qn(dir));
184 result += make_pair(Scalar(g_xyz),
185 BraketPair<F,BKType>(f,gm1));
186 }
187 const F& gp1 = g + _1;
188 result += make_pair(Scalar(-2.0) * Scalar(zeta),
189 BraketPair<F,BKType>(f,gp1));
190 return result;
191 }
192
194 template <class F, BraketType BKType>
195 LinearCombination<
196 SafePtr<DGVertex>,
197 BraketPair<F,BKType>
199 unsigned int xyz) {
200 if (BKType == CBra || BKType == CKet)
201 throw std::logic_error("R12v can only be applied to physicists brakets");
202
203 const char* XY = (BKType == PBra) ? "AC" : "BD";
204
205 const F& f = bkt[0];
206 const F& g = bkt[1];
208 ResultType result;
209
210 using namespace libint2::prefactor;
211 using namespace libint2::algebra;
212 // if applying to shell pair, change angular momentum along 0
213 const unsigned int dir = boost::is_same<F,CGF>::value ? xyz : 0;
214 F _1 = unit<F>(dir);
215 const F& fp1 = f + _1;
216 const F& gp1 = g + _1;
217 result += make_pair(Scalar(1.0),
218 BraketPair<F,BKType>(fp1,g));
219 result += make_pair(Scalar(-1.0),
220 BraketPair<F,BKType>(f,gp1));
221 result += make_pair(Vector(XY)[xyz],
223 return result;
224 }
225
226};
227
228#endif // header guard
BraketPair is a trimmed down version of ArrayBraket specialized for same-particle or different-partic...
Definition: src/bin/libint/braket.h:244
represents linear combination of objects of type T with coefficients of type C
Definition: algebra.h:222
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
LinearCombination< SafePtr< DGVertex >, BraketPair< F, BKType > > Nabla1(const BraketPair< F, BKType > &bkt, int xyz)
Applies Nabla1 to a physicists' braket.
Definition: gaussoper.h:132
LinearCombination< SafePtr< DGVertex >, BraketPair< F, BKType > > Nabla2(const BraketPair< F, BKType > &bkt, int xyz)
Applies Nabla2 to a physicists' braket.
Definition: gaussoper.h:165
LinearCombination< SafePtr< DGVertex >, BraketPair< F, BKType > > R12vec_dot_Nabla2(const BraketPair< F, BKType > &bkt)
Applies R12vec_dot_Nabla2 to a physicists' braket.
Definition: gaussoper.h:84
LinearCombination< SafePtr< DGVertex >, BraketPair< F, BKType > > R12vec_dot_Nabla1(const BraketPair< F, BKType > &bkt)
Applies R12vec_dot_Nabla1 to a physicists' braket.
Definition: gaussoper.h:36
LinearCombination< SafePtr< DGVertex >, BraketPair< F, BKType > > R12v(const BraketPair< F, BKType > &bkt, unsigned int xyz)
Applies R12v to a physicists' braket.
Definition: gaussoper.h:198