LIBINT  2.6.0
comp_deriv_gauss.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_compderivgauss_h_
22 #define _libint2_src_bin_libint_compderivgauss_h_
23 
24 #include <generic_rr.h>
25 #include <set>
26 
27 using namespace std;
28 
29 namespace libint2 {
30 
32  static CR_DerivGauss_GenericInstantiator instance_;
33 
34  CR_DerivGauss_GenericInstantiator(); // this is a singleton
36 
37  // pairs of L,vectorize specify the instances of GenericGaussDeriv template to be created
38  std::set<std::pair<unsigned int, bool> > template_instances_;
39 
40  public:
41  static CR_DerivGauss_GenericInstantiator& instance();
42  void add(unsigned int L, bool vectorize);
43  };
44 
61  template <class IntType,
62  int part,
63  FunctionPosition where,
64  int trans_inv_part = -1,
65  FunctionPosition trans_inv_where = InBra>
66  class CR_DerivGauss : public GenericRecurrenceRelation< CR_DerivGauss<IntType,part,where,trans_inv_part,trans_inv_where>,
67  typename IntType::BasisFunctionType,
68  IntType >
69  {
70  private:
71  static constexpr auto trans_inv_oper = not IntType::OperType::Properties::odep;
72  // can use translational invariance iff the operator is translationally-invariant
73  static constexpr auto using_trans_inv = trans_inv_oper &&
74  (part == trans_inv_part) &&
75  (where == trans_inv_where);
76 
77  public:
78  typedef CR_DerivGauss ThisType;
79  typedef typename IntType::BasisFunctionType BasisFunctionType;
80  typedef IntType TargetType;
82  friend class GenericRecurrenceRelation<ThisType,BasisFunctionType,TargetType>;
83  // # of children varies:
84  // 1. translational invariance makes num_functions - 1 children
85  // 2. direct differentiation always makes at most 2 Gaussians
86  static constexpr auto max_nchildren = using_trans_inv ? (IntType::num_bf - 1) : 2u;
87 
88  using ParentType::Instance;
89 
91  static bool directional() { return true; }
92 
93  private:
94  using ParentType::RecurrenceRelation::expr_;
95  using ParentType::RecurrenceRelation::nflops_;
96  using ParentType::target_;
97  using ParentType::is_simple;
98 
100  CR_DerivGauss(const SafePtr<TargetType>&, unsigned int dir);
101 
102  static std::string descr() { return "CR_DerivGauss"; }
106  std::string generate_label() const
107  {
108  typedef typename TargetType::AuxIndexType mType;
109  static SafePtr<mType> aux0(new mType(0u));
110  ostringstream os;
111  os << descr() << "P" << part << to_string(where)
112  << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
113  return os.str();
114  }
115 
116 #if LIBINT_ENABLE_GENERIC_CODE
117  bool has_generic(const SafePtr<CompilationParameters>& cparams) const;
120  std::string generic_header() const;
122  std::string generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const;
123 #endif
124  };
125 
126  template <class IntType, int part, FunctionPosition where,
127  int trans_inv_part, FunctionPosition trans_inv_where>
128  CR_DerivGauss<IntType,part,where,trans_inv_part,trans_inv_where>::CR_DerivGauss(
129  const SafePtr< TargetType >& Tint,
130  unsigned int dir) :
131  ParentType(Tint,dir)
132  {
133  using namespace libint2::algebra;
134  using namespace libint2::prefactor;
135  using namespace libint2::braket;
136  typedef BasisFunctionType F;
137  const F& _1 = unit<F>(is_simple() ? dir : 0); // for shell sets increment the first index
138 
139  const typename IntType::AuxQuantaType& aux = Tint->aux();
140  const typename IntType::OperType& oper = Tint->oper();
141 
142  // WARNING assuming one function per position
143 
144  // the Gaussian must be differentiated in direction dir
145  {
146  if (where == InBra && Tint->bra(part,0).deriv().d(dir) == 0)
147  return;
148  if (where == InKet && Tint->ket(part,0).deriv().d(dir) == 0)
149  return;
150  }
151 
152  // if not using translational invariance ...
153  // can expand derivatives of primitive Gaussians only
154  if (not using_trans_inv) {
155  if (where == InBra && Tint->bra(part,0).contracted())
156  return;
157  if (where == InKet && Tint->ket(part,0).contracted())
158  return;
159  }
160 
161  typedef typename IntType::BraType IBraType;
162  typedef typename IntType::KetType IKetType;
163  IBraType* bra = new IBraType(Tint->bra());
164  IKetType* ket = new IKetType(Tint->ket());
165 
166  if (not using_trans_inv) { // differentiate
167 
168  if (where == InBra) {
169  F a(bra->member(part,0));
170 
171  // add a+1
172  F ap1(bra->member(part,0) + _1);
173  ap1.deriv().dec(dir);
174  bra->set_member(ap1,part,0);
175  auto int_ap1 = this->add_child(IntType::Instance(*bra,*ket,aux,oper));
176  bra->set_member(a,part,0);
177  if (is_simple()) {
178  std::ostringstream oss;
179  oss << "two_alpha" << part << "_bra";
180  expr_ = Scalar(oss.str()) * int_ap1; nflops_+=1;
181  }
182 
183  // See if a-1 exists
184  F am1(bra->member(part,0) - _1);
185  if (exists(am1)) {
186  am1.deriv().dec(dir);
187  bra->set_member(am1,part,0);
188  auto int_am1 = this->add_child(IntType::Instance(*bra,*ket,aux,oper));
189  bra->set_member(a,part,0);
190  if (is_simple()) {
191  expr_ -= Scalar(a[dir]) * int_am1; nflops_+=2;
192  }
193  }
194  return;
195  }
196 
197  if (where == InKet) {
198  F a(ket->member(part,0));
199 
200  // add a+1
201  F ap1(ket->member(part,0) + _1);
202  ap1.deriv().dec(dir);
203  ket->set_member(ap1,part,0);
204  auto int_ap1 = this->add_child(IntType::Instance(*bra,*ket,aux,oper));
205  ket->set_member(a,part,0);
206  if (is_simple()) {
207  std::ostringstream oss;
208  oss << "two_alpha" << part << "_ket";
209  expr_ = Scalar(oss.str()) * int_ap1; nflops_+=1;
210  }
211 
212  // See if a-1 exists
213  F am1(ket->member(part,0) - _1);
214  if (exists(am1)) {
215  am1.deriv().dec(dir);
216  ket->set_member(am1,part,0);
217  auto int_am1 = this->add_child(IntType::Instance(*bra,*ket,aux,oper));
218  ket->set_member(a,part,0);
219  if (is_simple()) {
220  expr_ -= Scalar(a[dir]) * int_am1; nflops_+=2;
221  }
222  }
223  return;
224  }
225  } else { // use translational invariance
226 
227  // remove one deriv quantum from the target function
228  if (where == InBra) bra->member(part,0).deriv().dec(dir);
229  if (where == InKet) ket->member(part,0).deriv().dec(dir);
230 
231  int term_count = 0;
232  for (int p = 0; p != IntType::num_particles; ++p) {
233  if (p != trans_inv_part || trans_inv_where != InBra) {
234  F a(bra->member(p, 0));
235  if (not a.is_unit()) {
236  F da(a);
237  da.deriv().inc(dir);
238  bra->set_member(da, p, 0);
239  auto int_da =
240  this->add_child(IntType::Instance(*bra, *ket, aux, oper));
241  bra->set_member(a, p, 0);
242  if (is_simple()) {
243  std::ostringstream oss;
244  if (term_count == 0)
245  expr_ = Scalar(-1) * int_da;
246  else
247  expr_ -= int_da;
248  ++term_count;
249  nflops_ += 1;
250  }
251  }
252  }
253  if (p != trans_inv_part || trans_inv_where != InKet) {
254  F a(ket->member(p, 0));
255  if (not a.is_unit()) {
256  F da(a);
257  da.deriv().inc(dir);
258  ket->set_member(da, p, 0);
259  auto int_da =
260  this->add_child(IntType::Instance(*bra, *ket, aux, oper));
261  ket->set_member(a, p, 0);
262  if (is_simple()) {
263  std::ostringstream oss;
264  if (term_count == 0)
265  expr_ = Scalar(-1) * int_da;
266  else
267  expr_ -= int_da;
268  ++term_count;
269  nflops_ += 1;
270  }
271  }
272  }
273  }
274  }
275 
276  return;
277  }
278 
279 #if LIBINT_ENABLE_GENERIC_CODE
280  template <class IntType, int part, FunctionPosition where,
281  int trans_inv_part, FunctionPosition trans_inv_where>
282  bool
283  CR_DerivGauss<IntType,part,where,trans_inv_part,trans_inv_where>::has_generic(
284  const SafePtr<CompilationParameters>& cparams
285  ) const
286  {
287  // not implemented generic relation for translational invariance yet
288  if (using_trans_inv) return false;
289 
290  if (TrivialBFSet<BasisFunctionType>::result)
291  return false;
292 
293  // generate generic code if the average quantum number > max_l_opt
294  const unsigned int max_opt_am = cparams->max_am_opt();
295  unsigned int am_total = 0;
296  unsigned int nfunctions = 0;
297  const unsigned int np = IntType::OperType::Properties::np;
298  for(unsigned int p=0; p<np; p++) {
299  unsigned int nbra = target_->bra().num_members(p);
300  for(unsigned int i=0; i<nbra; i++) {
301  am_total += target_->bra(p,i).norm();
302  ++nfunctions;
303  }
304  unsigned int nket = target_->ket().num_members(p);
305  for(unsigned int i=0; i<nket; i++) {
306  am_total += target_->ket(p,i).norm();
307  ++nfunctions;
308  }
309  }
310  if (am_total > max_opt_am*nfunctions)
311  return true;
312 
313  // else generate explicit code
314  return false;
315  }
316 
317  template <class IntType, int part, FunctionPosition where,
318  int trans_inv_part, FunctionPosition trans_inv_where>
319  std::string
320  CR_DerivGauss<IntType,part,where,trans_inv_part,trans_inv_where>::generic_header() const
321  {
322  return std::string("GenericGaussDeriv.h");
323  }
324 
325  template <class IntType, int part, FunctionPosition where,
326  int trans_inv_part, FunctionPosition trans_inv_where>
327  std::string
328  CR_DerivGauss<IntType,part,where,trans_inv_part,trans_inv_where>::generic_instance(
329  const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args
330  ) const
331  {
332  std::ostringstream oss;
333 
334  oss << "using namespace libint2;" << endl;
335 
336  BasisFunctionType sh(where == InBra ? target_->bra(part,0) : target_->ket(part,0));
337 
338  const unsigned int L = sh.norm();
339  const bool vectorize = (context->cparams()->max_vector_length() == 1) ? false : true;
340  oss << "libint2::GenericGaussDeriv<" << L << ","
341  << (vectorize ? "true" : "false")
342  << ">::compute(inteval";
343 
344  oss << "," << args->symbol(0); // target
345  const unsigned int nargs = args->n();
346  for(unsigned int a=1; a<nargs; a++) {
347  oss << "," << args->symbol(a);
348  }
349  // L == 0 => second argument not needed
350  if (nargs == 2)
351  oss << ",0";
352 
353  // then dimensions of basis function sets not involved in the transfer
354  unsigned int hsr = 1;
355  unsigned int lsr = 1;
356  const unsigned int np = IntType::OperType::Properties::np;
357  // a cleaner way to count the number of function sets referring
358  // to some particles is to construct a dummy integral and
359  // use subiterator policy
360  // WARNING !!!
361  for(int p=0; p<static_cast<int>(np); p++) {
362  unsigned int nbra = target_->bra().num_members(p);
363  assert(nbra == 1);
364  for(unsigned int i=0; i<nbra; i++) {
365  SubIterator* iter = target_->bra().member_subiter(p,i);
366  if (p < part || (p == part && where == InKet))
367  hsr *= iter->num_iter();
368  // skip p == part && where == InBra
369  if (p > part)
370  lsr *= iter->num_iter();
371  delete iter;
372  }
373  unsigned int nket = target_->ket().num_members(p);
374  assert(nket == 1);
375  for(unsigned int i=0; i<nket; i++) {
376  SubIterator* iter = target_->ket().member_subiter(p,i);
377  if (p < part)
378  hsr *= iter->num_iter();
379  // skip p == part && where == InKet
380  if (p > part || (p == part && where == InBra))
381  lsr *= iter->num_iter();
382  delete iter;
383  }
384  }
385  oss << "," << hsr << "," << lsr;
386 
387  // direction
388  oss << "," << this->dir();
389 
390  // two_alpha prefactor
391  oss << ",inteval->two_alpha" << part << "_" << (where == InBra ? "bra" : "ket");
392 
393  oss << ");";
394 
395  CR_DerivGauss_GenericInstantiator::instance().add(L, vectorize);
396 
397  return oss.str();
398  }
399 #endif // #if !LIBINT_ENABLE_GENERIC_CODE
400 
401 }; // namespace libint2
402 
403 #endif
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
RRImpl must inherit GenericRecurrenceRelation<RRImpl>
Definition: generic_rr.h:49
Compute relation for (geometric) derivative Gaussian ints of generic type IntType .
Definition: comp_deriv_gauss.h:66
these objects help to construct BraketPairs
Definition: braket.h:270
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
const SafePtr< DGVertex > & add_child(const SafePtr< DGVertex > &child)
add child
Definition: generic_rr.h:110
static bool directional()
always directional! Cartesian derivatives are applied in a particular direction
Definition: comp_deriv_gauss.h:91
Definition: comp_deriv_gauss.h:31