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