LIBINT  2.6.0
integral.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_integral_h_
22 #define _libint2_src_bin_libint_integral_h_
23 
24 #include <cassert>
25 #include <smart_ptr.h>
26 #include <dgvertex.h>
27 #include <oper.h>
28 #include <iter.h>
29 #include <policy_spec.h>
30 #include <quanta.h>
31 #include <equiv.h>
32 #include <singl_stack.h>
33 #include <class_registry.h>
34 #include <global_macros.h>
35 
36 #if USE_BRAKET_H
37 # include <braket.h>
38 #endif
39 
40 using namespace std;
41 
42 extern long living_count;
43 
44 namespace libint2 {
45 
50  template <class BasisFunctionSet> class IntegralSet {
51 
52  public:
53  virtual ~IntegralSet() {};
54 
55 #if USE_BRAKET_H
56  virtual unsigned int num_part() const =0;
59  virtual unsigned int num_func_bra(unsigned int p) const =0;
61  virtual unsigned int num_func_ket(unsigned int p) const =0;
63  virtual const BasisFunctionSet& bra(unsigned int p, unsigned int i) const =0;
65  virtual const BasisFunctionSet& ket(unsigned int p, unsigned int i) const =0;
67  virtual BasisFunctionSet& bra(unsigned int p, unsigned int i) =0;
69  virtual BasisFunctionSet& ket(unsigned int p, unsigned int i) =0;
70 #else
71  virtual unsigned int np() const =0;
74  virtual const SafePtr<BasisFunctionSet> bra(unsigned int p, unsigned int i) const =0;
76  virtual const SafePtr<BasisFunctionSet> ket(unsigned int p, unsigned int i) const =0;
77 #endif
78  };
79 
90  template <class Oper, class BFS, class BraSetType, class KetSetType, class AuxQuanta = EmptySet>
92  public IntegralSet<BFS>, public DGVertex,
93  public EnableSafePtrFromThis< GenIntegralSet<Oper,BFS,BraSetType,KetSetType,AuxQuanta> >
94  {
95  public:
96  typedef GenIntegralSet this_type;
103 #if USE_INT_KEY_TO_HASH
104  typedef mpz_class key_type;
105 #else
106  typedef std::string key_type;
107 #endif
111  typedef Oper OperType;
113  typedef typename BraSetType::bfs_type BasisFunctionType;
115  static constexpr auto num_particles = BraSetType::num_particles;
116  static_assert(BraSetType::num_particles == KetSetType::num_particles,
117  "number of particles in bra and ket must be same");
119  static constexpr auto num_bf = BraSetType::num_bf + KetSetType::num_bf;
120 
127  virtual ~GenIntegralSet();
129 
131  static const SafePtr<GenIntegralSet> Instance(const BraSetType& bra, const KetSetType& ket, const AuxQuanta& aux, const Oper& oper=Oper());
132 
134  virtual bool operator==(const GenIntegralSet&) const;
136  bool equiv(const SafePtr<DGVertex>& v) const
137  {
138  return PtrComp::equiv(this,v);
139  }
141  virtual unsigned int size() const;
143  virtual const std::string& label() const;
145  virtual const std::string& id() const;
147  virtual std::string description() const;
148 
150  unsigned int num_part() const { return OperType::Properties::np; }
152  virtual unsigned int num_func_bra(unsigned int p) const { return bra_.num_members(p); }
154  virtual unsigned int num_func_ket(unsigned int p) const { return ket_.num_members(p); }
156  typename BraSetType::bfs_cref bra(unsigned int p, unsigned int i) const;
158  typename KetSetType::bfs_cref ket(unsigned int p, unsigned int i) const;
160  typename BraSetType::bfs_ref bra(unsigned int p, unsigned int i);
162  typename KetSetType::bfs_ref ket(unsigned int p, unsigned int i);
163 
164  typedef BraSetType BraType;
165  typedef KetSetType KetType;
166  typedef Oper OperatorType;
167  typedef AuxQuanta AuxQuantaType;
168 
170  const SafePtr<Oper> oper() const;
172  const BraType& bra() const;
174  const KetType& ket() const;
176  const SafePtr<AuxQuanta> aux() const;
177 
179  DGVertex::KeyReturnType key() const {
180  if (key_ == 0) compute_key();
181  return key_;
182  }
183 
185  void unregister() const;
186 
188  virtual bool auto_unroll() const { return false; }
189 
190  protected:
191  // Basic Integral constructor. It is protected so that derived classes don't have to behave like singletons
192  GenIntegralSet(const Oper& oper, const BraSetType& bra, const KetSetType& ket, const AuxQuanta& aux);
194  static key_type compute_key(const Oper& O, const BraType& bra, const KetType& ket, const AuxQuanta& aux) {
195 #define TEST_KEYTYPE_SAFETY 0
196 #if TEST_KEYTYPE_SAFETY
197  key_type remainder = UINT64_MAX;
198  remainder /= (key_type)aux.max_key(); assert(remainder != 0);
199  remainder /= (key_type)ket.max_key(); assert(remainder != 0);
200  remainder /= (key_type)bra.max_key(); assert(remainder != 0);
201  remainder /= (key_type)O.max_key; assert(remainder != 0);
202 #endif
203 
204  key_type key;
205 
206  key = ( (key_type(O.key()) * KeyTypes::cast(bra.max_key()) + KeyTypes::cast(bra.key()) ) * KeyTypes::cast(ket.max_key()) +
207  KeyTypes::cast(ket.key()) ) * KeyTypes::cast(aux.max_key()) + KeyTypes::cast(aux.key());
208 
209  return key;
210 
211  }
212 
213  BraSetType bra_;
214  KetSetType ket_;
215 
217  void set_size(unsigned int sz);
219  virtual bool this_precomputed() const { return false; }
221  void reset_cache() { key_ = 0; size_ = 0; }
222 
223  private:
224  //
225  // All integrals are Singletons by nature, therefore they must be treated as such
226  // 1) No public constructors are provided
227  // 2) protected members are provided to implement Singleton-type functionality
228  //
229  GenIntegralSet();
231  // Copy is not permitted
232  GenIntegralSet& operator=(const GenIntegralSet& source);
233 
234  // This is used to manage GenIntegralSet objects as singletons
235  static SingletonManagerType singl_manager_;
236 
237  // The operator needs to be a real object rather than real type to be able to construct a SubIterator, etc.
238  SafePtr<Oper> O_;
239  // Same for AuxQuanta
240  SafePtr<AuxQuanta> aux_;
241 
242  // size of the integral set
243  mutable unsigned int size_;
244 
245  // unique label -- no 2 GenIntegralSet instances can have the same label!
246  mutable std::string label_;
247  // generates label_
248  std::string generate_label() const;
249  // key
250  mutable key_type key_;
251 
253  void compute_key() const {
254  key_ = compute_key(*(O_.get()),bra_,ket_,*(aux_.get()));
255  }
256 
257  };
258 
259 #if USE_INT_KEY_TO_HASH
260  // I use key() to hash GenIntegralSet. Therefore compute_key() must return unique keys!
261  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
262  typename GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::SingletonManagerType
263  GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::singl_manager_(&GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::key);
264 #else
265  // I use label() to hash GenIntegralSet. Therefore labels must be unique!
266  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
267  typename GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::SingletonManagerType
268  GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::singl_manager_(&GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::label);
269 #endif
270 
271  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
272  GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::GenIntegralSet(const Op& oper, const BraSetType& bra, const KetSetType& ket, const AuxQuanta& aux) :
273  DGVertex(ClassInfo<GenIntegralSet>::Instance().id()), bra_(bra), ket_(ket), O_(SafePtr<Op>(new Op(oper))), aux_(SafePtr<AuxQuanta>(new AuxQuanta(aux))),
274  size_(0), label_()
275  {
276  if (Op::Properties::np != bra.num_part())
277  throw std::runtime_error("GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::GenIntegralSet(bra,ket) -- number of particles in bra doesn't match that in the operator");
278  if (Op::Properties::np != ket.num_part())
279  throw std::runtime_error("GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::GenIntegralSet(bra,ket) -- number of particles in ket doesn't match that in the operator");
280  compute_key();
281 #if DEBUG
282  std::cout << "GenIntegralSet: constructed " << label() << std::endl;
283  std::cout << "GenIntegralSet: living_count = " << ++living_count << std::endl;
284 #endif
285  }
286 
287  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
289  {
290 #if DEBUG
291  std::cout << "GenIntegralSet: destructed " << label() << std::endl;
292  std::cout << "GenIntegralSet: living_count = " << --living_count << std::endl;
293 #endif
294  }
295 
296  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
297  const SafePtr< GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta> >
298  GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::Instance(const BraSetType& bra, const KetSetType& ket, const AuxQuanta& aux, const Op& oper)
299  {
300  typedef typename SingletonManagerType::value_type map_value_type;
301  key_type key = compute_key(oper,bra,ket,aux);
302  const map_value_type& val = singl_manager_.find(key);
303  if (!val.second) {
304  SafePtr<this_type> this_int(new this_type(oper,bra,ket,aux));
305  // Use singl_manager_ to make sure this is a new object of this type
306  const map_value_type& val = singl_manager_.find(this_int);
307  val.second->instid_ = val.first;
308  this_int.reset();
309  return val.second;
310  }
311  return val.second;
312  }
313 
314  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
315  typename BraSetType::bfs_cref
317  {
318  return bra_.member(p,i);
319  }
320 
321  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
322  typename KetSetType::bfs_cref
324  {
325  return ket_.member(p,i);
326  }
327 
328  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
329  typename BraSetType::bfs_ref
331  {
332  reset_cache();
333  return bra_.member(p,i);
334  }
335 
336  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
337  typename KetSetType::bfs_ref
339  {
340  reset_cache();
341  return ket_.member(p,i);
342  }
343 
344  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
345  const typename GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::BraType&
347  {
348  return bra_;
349  }
350 
351  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
352  const typename GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::KetType&
354  {
355  return ket_;
356  }
357 
358  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
359  const SafePtr<Op>
361  {
362  return O_;
363  }
364 
365  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
366  const SafePtr<AuxQuanta>
368  {
369  return aux_;
370  }
371 
372  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
373  bool
375  {
376 #if USE_INT_KEY_TO_COMPARE
377  return key() == a.key();
378 #else
379  bool aux_equiv = PtrEquiv<AuxQuanta>::equiv(aux_,a.aux_);
380  if (!aux_equiv) return false;
381  bool oper_equiv = PtrEquiv<Op>::equiv(O_,a.O_);
382  bool bra_equiv = PtrEquiv<BraSetType>::equiv(bra_,a.bra_);
383  if (!bra_equiv) return false;
384  bool ket_equiv = PtrEquiv<KetSetType>::equiv(ket_,a.ket_);
385  if (!ket_equiv) return false;
386  return true;
387 #endif
388  }
389 
390  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
391  void
393  {
394  singl_manager_.remove(const_pointer_cast<this_type,const this_type>(EnableSafePtrFromThis<this_type>::SafePtr_from_this()));
395  }
396 
397  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
398  unsigned int
400  {
401  if (size_ > 0)
402  return size_;
403 
404 #if COMPUTE_SIZE_DIRECTLY
405  // WARNING: will not work if aux_ includes "sets" of numbers, i.e. when a quantum number can be a set of numbers
406  // but I don't at the moment see why this would be needed
407  size_ = bra_.size() * ket_.size() * O_->num_oper();
408 #else
409  // compute size
410  SafePtr<this_type> this_ptr = const_pointer_cast<this_type,const this_type>(EnableSafePtrFromThis<GenIntegralSet>::SafePtr_from_this());
411  SafePtr< SubIteratorBase<this_type> > siter(new SubIteratorBase<this_type>(this_ptr));
412  size_ = siter->num_iter();
413  if (size_ == 0)
414  throw std::runtime_error("GenIntegralSet<Op,BFS,BraSetType,KetSetType,AuxQuanta>::size() -- size is 0");
415 #endif
416 
417  return size_;
418  }
419 
420  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
421  void
423  {
424  size_ = sz;
425  }
426 
427  template <class BraSetType, class KetSetType, class AuxQuanta, class Op>
428  std::string
429  genintegralset_label(const BraSetType& bra, const KetSetType& ket, const SafePtr<AuxQuanta>& aux, const SafePtr<Op>& O)
430  {
431  ostringstream os;
432  os << "< ";
433  for(unsigned int p=0; p<Op::Properties::np; p++) {
434  unsigned int nbra = bra.num_members(p);
435  for(unsigned int i=0; i<nbra; i++)
436 #if USE_BRAKET_H
437  os << bra.member(p,i).label() << "(" << p << ") ";
438 #else
439  os << bra.member(p,i)->label() << "(" << p << ") ";
440 #endif
441  }
442  os << "| " << O->label() << " | ";
443  for(unsigned int p=0; p<Op::Properties::np; p++) {
444  unsigned int nket = ket.num_members(p);
445  for(unsigned int i=0; i<nket; i++)
446 #if USE_BRAKET_H
447  os << ket.member(p,i).label() << "(" << p << ") ";
448 #else
449  os << ket.member(p,i)->label() << "(" << p << ") ";
450 #endif
451  }
452  os << "> ^ { " << aux->label() << " }";
453  return os.str();
454  }
455 
456  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
457  const std::string&
459  {
460  if (label_.empty())
461  label_ = generate_label();
462  return label_;
463  }
464 
465  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
466  std::string
468  {
469  return genintegralset_label(bra_,ket_,aux_,O_);
470  }
471 
472  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
473  const std::string&
475  {
476  return label();
477  }
478 
479  template <class Op, class BFS, class BraSetType, class KetSetType, class AuxQuanta>
480  std::string
482  {
483  ostringstream os;
484  os << " GenIntegralSet: " << label();
485  const std::string descr = os.str();
486  return descr;
487  }
488 
489 };
490 
491 #endif
virtual unsigned int size() const
Specialization of DGVertex::size()
Definition: integral.h:399
GenIntegralSet< typename Oper::iter_type, BFS, typename BraSetType::iter_type, typename KetSetType::iter_type, typename AuxQuanta::iter_type > iter_type
GenIntegralSet is a set of these subobjects.
Definition: integral.h:98
IntegralSet< BFS > parent_type
GenIntegralSet is derived from IntegralSet.
Definition: integral.h:100
SingletonStack<T,KeyType> helps to implement Singleton-like objects of type T.
Definition: singl_stack.h:43
bool equiv(const SafePtr< DGVertex > &v) const
Specialization of DGVertex::equiv()
Definition: integral.h:136
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
GenIntegralSet is a set of integrals over functions derived from BFS.
Definition: integral.h:91
This is a vertex of a Directed Graph (DG)
Definition: dgvertex.h:43
SubIteratorBase<T> provides a base class for a sub-iterator class for T.
Definition: iter.h:73
virtual bool auto_unroll() const
If consists of precomputed elements, override this to return true.
Definition: integral.h:188
PtrEquiv<T> provides a set of comparison functions named 'equiv' which take as arguments a mix of ref...
Definition: equiv.h:36
Oper is OperSet characterized by properties Props.
Definition: oper.h:90
BraSetType::bfs_type BasisFunctionType
This is the real type of basis functions.
Definition: integral.h:113
PtrEquiv< GenIntegralSet > PtrComp
This type provides comparison operations on pointers to GenIntegralSet.
Definition: integral.h:102
DGVertex::KeyReturnType key() const
Implements Hashable::key()
Definition: integral.h:179
LIBINT2_UINT_LEAST64 key() const
Implements Hashable::key()
Definition: braket.h:177
virtual unsigned int num_func_bra(unsigned int p) const
Implementation of IntegralSet::num_func_bra.
Definition: integral.h:152
unsigned int num_part() const
Implementation of IntegralSet::num_part.
Definition: integral.h:150
LIBINT2_UINT_LEAST64 max_key() const
key lies in range [0,max_key())
Definition: braket.h:190
This is an abstract base for sets of all types of integrals.
Definition: integral.h:50
virtual bool this_precomputed() const
Specialization of DGVertex::this_precomputed()
Definition: integral.h:219
virtual unsigned int num_func_ket(unsigned int p) const
Implementation of IntegralSet::num_func_ket.
Definition: integral.h:154
void reset_cache()
Resets all cached values.
Definition: integral.h:221
ArrayBraket is a lightweight implementation of Braket concept.
Definition: braket.h:38
Oper OperType
This is the type of the operator.
Definition: integral.h:111
static key_type compute_key(const Oper &O, const BraType &bra, const KetType &ket, const AuxQuanta &aux)
computes a key. it's protected so that derived classes can use it to implement smart caching in const...
Definition: integral.h:194