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