LIBINT  2.6.0
gaussoper.h
1 /*
2  * Copyright (C) 2004-2019 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 
29 namespace 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],
64  BraketPair<F,BKType>(fm1,g));
65  }
66  const F& fp1 = f + _1;
67  const F& fp2 = fp1 + _1;
68  result += make_pair(Scalar(-2.0) * Scalar(zeta),
69  BraketPair<F,BKType>(fp2,g));
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],
73  BraketPair<F,BKType>(fp1,g));
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],
222  BraketPair<F,BKType>(f,g));
223  return result;
224  }
225 
226 };
227 
228 #endif // header guard
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
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
LinearCombination< SafePtr< DGVertex >, BraketPair< F, BKType > > Nabla1(const BraketPair< F, BKType > &bkt, int xyz)
Applies Nabla1 to a physicists' braket.
Definition: gaussoper.h:132
represents linear combination of objects of type T with coefficients of type C
Definition: algebra.h:223
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
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
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
LinearCombination< SafePtr< DGVertex >, BraketPair< F, BKType > > Nabla2(const BraketPair< F, BKType > &bkt, int xyz)
Applies Nabla2 to a physicists' braket.
Definition: gaussoper.h:165
BraketPair is a trimmed down version of ArrayBraket specialized for same-particle or different-partic...
Definition: braket.h:244