LIBINT  2.1.0-stable
comp_deriv_gauss.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_compderivgauss_h_
21 #define _libint2_src_bin_libint_compderivgauss_h_
22 
23 #include <generic_rr.h>
24 #include <set>
25 
26 using namespace std;
27 
28 namespace libint2 {
29 
31  static CR_DerivGauss_GenericInstantiator instance_;
32 
33  CR_DerivGauss_GenericInstantiator(); // this is a singleton
35 
36  // pairs of L,vectorize specify the instances of GenericGaussDeriv template to be created
37  std::set<std::pair<unsigned int, bool> > template_instances_;
38 
39  public:
40  static CR_DerivGauss_GenericInstantiator& instance();
41  void add(unsigned int L, bool vectorize);
42  };
43 
47  template <class IntType, int part, FunctionPosition where>
48  class CR_DerivGauss : public GenericRecurrenceRelation< CR_DerivGauss<IntType,part,where>,
49  typename IntType::BasisFunctionType,
50  IntType >
51  {
52  public:
53  typedef CR_DerivGauss ThisType;
54  typedef typename IntType::BasisFunctionType BasisFunctionType;
55  typedef IntType TargetType;
57  friend class GenericRecurrenceRelation<ThisType,BasisFunctionType,TargetType>;
58  static const unsigned int max_nchildren = 2;
59 
60  using ParentType::Instance;
61 
63  static bool directional() { return true; }
64 
65  private:
66  using ParentType::RecurrenceRelation::expr_;
67  using ParentType::RecurrenceRelation::nflops_;
68  using ParentType::target_;
69  using ParentType::is_simple;
70 
72  CR_DerivGauss(const SafePtr<TargetType>&, unsigned int dir);
73 
74  static std::string descr() { return "CR_DerivGauss"; }
78  std::string generate_label() const
79  {
80  typedef typename TargetType::AuxIndexType mType;
81  static SafePtr<mType> aux0(new mType(0u));
82  ostringstream os;
83  os << descr() << "P" << part << to_string(where)
84  << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
85  return os.str();
86  }
87 
88 #if LIBINT_ENABLE_GENERIC_CODE
89  bool has_generic(const SafePtr<CompilationParameters>& cparams) const;
92  std::string generic_header() const;
94  std::string generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const;
95 #endif
96  };
97 
98  template <class IntType, int part, FunctionPosition where>
99  CR_DerivGauss<IntType,part,where>::CR_DerivGauss(const SafePtr< TargetType >& Tint,
100  unsigned int dir) :
101  ParentType(Tint,dir)
102  {
103  using namespace libint2::algebra;
104  using namespace libint2::prefactor;
105  using namespace libint2::braket;
106  typedef BasisFunctionType F;
107  const F& _1 = unit<F>(is_simple() ? dir : 0); // for shell sets increment the first index
108 
109  const typename IntType::AuxQuantaType& aux = Tint->aux();
110  const typename IntType::OperType& oper = Tint->oper();
111 
112  // WARNING assuming one function per position
113  { // can't apply to contracted basis functions
114  if (where == InBra && Tint->bra(part,0).contracted())
115  return;
116  if (where == InKet && Tint->ket(part,0).contracted())
117  return;
118  }
119 
120  // the Gaussian must be differentiated in direction dir
121  {
122  if (where == InBra && Tint->bra(part,0).deriv().d(dir) == 0)
123  return;
124  if (where == InKet && Tint->ket(part,0).deriv().d(dir) == 0)
125  return;
126  }
127 
128  typedef typename IntType::BraType IBraType;
129  typedef typename IntType::KetType IKetType;
130  IBraType* bra = new IBraType(Tint->bra());
131  IKetType* ket = new IKetType(Tint->ket());
132 
133  if (where == InBra) {
134  F a(bra->member(part,0));
135 
136  // add a+1
137  F ap1(bra->member(part,0) + _1);
138  ap1.deriv().dec(dir);
139  bra->set_member(ap1,part,0);
140  auto int_ap1 = this->add_child(IntType::Instance(*bra,*ket,aux,oper));
141  bra->set_member(a,part,0);
142  if (is_simple()) {
143  std::ostringstream oss;
144  oss << "two_alpha" << part << "_bra";
145  expr_ = Scalar(oss.str()) * int_ap1; nflops_+=1;
146  }
147 
148  // See if a-1 exists
149  F am1(bra->member(part,0) - _1);
150  if (exists(am1)) {
151  am1.deriv().dec(dir);
152  bra->set_member(am1,part,0);
153  auto int_am1 = this->add_child(IntType::Instance(*bra,*ket,aux,oper));
154  bra->set_member(a,part,0);
155  if (is_simple()) {
156  expr_ -= Scalar(a[dir]) * int_am1; nflops_+=2;
157  }
158  }
159  return;
160  }
161 
162  if (where == InKet) {
163  F a(ket->member(part,0));
164 
165  // add a+1
166  F ap1(ket->member(part,0) + _1);
167  ap1.deriv().dec(dir);
168  ket->set_member(ap1,part,0);
169  auto int_ap1 = this->add_child(IntType::Instance(*bra,*ket,aux,oper));
170  ket->set_member(a,part,0);
171  if (is_simple()) {
172  std::ostringstream oss;
173  oss << "two_alpha" << part << "_ket";
174  expr_ = Scalar(oss.str()) * int_ap1; nflops_+=1;
175  }
176 
177  // See if a-1 exists
178  F am1(ket->member(part,0) - _1);
179  if (exists(am1)) {
180  am1.deriv().dec(dir);
181  ket->set_member(am1,part,0);
182  auto int_am1 = this->add_child(IntType::Instance(*bra,*ket,aux,oper));
183  ket->set_member(a,part,0);
184  if (is_simple()) {
185  expr_ -= Scalar(a[dir]) * int_am1; nflops_+=2;
186  }
187  }
188  return;
189  }
190 
191  return;
192  }
193 
194 #if LIBINT_ENABLE_GENERIC_CODE
195  template <class IntType, int part, FunctionPosition where>
196  bool
197  CR_DerivGauss<IntType,part,where>::has_generic(const SafePtr<CompilationParameters>& cparams) const
198  {
200  return false;
201 
202  // generate generic code if the average quantum number > max_l_opt
203  const unsigned int max_opt_am = cparams->max_am_opt();
204  unsigned int am_total = 0;
205  unsigned int nfunctions = 0;
206  const unsigned int np = IntType::OperType::Properties::np;
207  for(unsigned int p=0; p<np; p++) {
208  unsigned int nbra = target_->bra().num_members(p);
209  for(unsigned int i=0; i<nbra; i++) {
210  am_total += target_->bra(p,i).norm();
211  ++nfunctions;
212  }
213  unsigned int nket = target_->ket().num_members(p);
214  for(unsigned int i=0; i<nket; i++) {
215  am_total += target_->ket(p,i).norm();
216  ++nfunctions;
217  }
218  }
219  if (am_total > max_opt_am*nfunctions)
220  return true;
221 
222  // else generate explicit code
223  return false;
224  }
225 
226  template <class IntType, int part, FunctionPosition where>
227  std::string
229  {
230  return std::string("GenericGaussDeriv.h");
231  }
232 
233  template <class IntType, int part, FunctionPosition where>
234  std::string
235  CR_DerivGauss<IntType,part,where>::generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const
236  {
237  std::ostringstream oss;
238 
239  oss << "using namespace libint2;" << endl;
240 
241  BasisFunctionType sh(where == InBra ? target_->bra(part,0) : target_->ket(part,0));
242 
243  const unsigned int L = sh.norm();
244  const bool vectorize = (context->cparams()->max_vector_length() == 1) ? false : true;
245  oss << "libint2::GenericGaussDeriv<" << L << ","
246  << (vectorize ? "true" : "false")
247  << ">::compute(inteval";
248 
249  oss << "," << args->symbol(0); // target
250  const unsigned int nargs = args->n();
251  for(unsigned int a=1; a<nargs; a++) {
252  oss << "," << args->symbol(a);
253  }
254  // L == 0 => second argument not needed
255  if (nargs == 2)
256  oss << ",0";
257 
258  // then dimensions of basis function sets not involved in the transfer
259  unsigned int hsr = 1;
260  unsigned int lsr = 1;
261  const unsigned int np = IntType::OperType::Properties::np;
262  // a cleaner way to count the number of function sets referring
263  // to some particles is to construct a dummy integral and
264  // use subiterator policy
265  // WARNING !!!
266  for(int p=0; p<static_cast<int>(np); p++) {
267  unsigned int nbra = target_->bra().num_members(p);
268  assert(nbra == 1);
269  for(unsigned int i=0; i<nbra; i++) {
270  SubIterator* iter = target_->bra().member_subiter(p,i);
271  if (p < part || (p == part && where == InKet))
272  hsr *= iter->num_iter();
273  // skip p == part && where == InBra
274  if (p > part)
275  lsr *= iter->num_iter();
276  delete iter;
277  }
278  unsigned int nket = target_->ket().num_members(p);
279  assert(nket == 1);
280  for(unsigned int i=0; i<nket; i++) {
281  SubIterator* iter = target_->ket().member_subiter(p,i);
282  if (p < part)
283  hsr *= iter->num_iter();
284  // skip p == part && where == InKet
285  if (p > part || (p == part && where == InBra))
286  lsr *= iter->num_iter();
287  delete iter;
288  }
289  }
290  oss << "," << hsr << "," << lsr;
291 
292  // direction
293  oss << "," << this->dir();
294 
295  // two_alpha prefactor
296  oss << ",inteval->two_alpha" << part << "_" << (where == InBra ? "bra" : "ket");
297 
298  oss << ");";
299 
300  CR_DerivGauss_GenericInstantiator::instance().add(L, vectorize);
301 
302  return oss.str();
303  }
304 #endif // #if !LIBINT_ENABLE_GENERIC_CODE
305 
306 }; // namespace libint2
307 
308 #endif
TrivialBFSet<T> defines static member result, which is true if T is a basis function set consisting o...
Definition: bfset.h:879
static bool directional()
always directional! Cartesian derivatives are applied in a particular direction
Definition: comp_deriv_gauss.h:63
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
Definition: stdarray.h:18
Definition: prefactors.h:159
RRImpl must inherit GenericRecurrenceRelation<RRImpl>
Definition: generic_rr.h:48
QuantumNumbersA<T,N> is a set of N quantum numbers of type T implemented in terms of a C-style array...
Definition: quanta.h:198
compute relation for derivative 2-e ERI.
Definition: comp_deriv_gauss.h:48
std::string to_string(const T &x)
Converts x to its string representation.
Definition: entity.h:71
Definition: algebra.h:32
these objects help to construct BraketPairs
Definition: braket.h:442
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:91
Iterator provides a base class for all object iterator classes.
Definition: iter.h:44
virtual unsigned int num_iter() const =0
Returns a number of iterations (number of elements in a set over which to iterate).
Definition: comp_deriv_gauss.h:30