LIBINT  2.6.0
multipole.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_multipole_h_
22 #define _libint2_src_bin_libint_multipole_h_
23 
24 #include <iostream>
25 #include <string>
26 #include <cassert>
27 #include <numeric>
28 #include <sstream>
29 #include <stdarray.h>
30 #include <smart_ptr.h>
31 #include <polyconstr.h>
32 #include <hashable.h>
33 #include <global_macros.h>
34 #include <util_types.h>
35 
36 namespace libint2 {
37 
40  template <unsigned NDIM = 3>
41  class CartesianMultipoleQuanta : public Hashable<LIBINT2_UINT_LEAST64,ReferToKey> {
42 
43  static_assert(NDIM == 1 || NDIM == 3, "CartesianMultipoleQuanta<NDIM>: only NDIM=1 or NDIM=3 are supported");
44  public:
45  CartesianMultipoleQuanta() : valid_(true) {
46  for(auto d=0u; d!=NDIM; ++d) n_[d] = 0u;
47  }
48  CartesianMultipoleQuanta(const CartesianMultipoleQuanta& other) : valid_(true) {
49  std::copy(other.n_, other.n_ + NDIM, n_);
50  }
51  CartesianMultipoleQuanta& operator=(const CartesianMultipoleQuanta& other) {
52  valid_ = other.valid_;
53  std::copy(other.n_, other.n_ + NDIM, n_);
54  return *this;
55  }
56  CartesianMultipoleQuanta& operator+=(const CartesianMultipoleQuanta& other) {
57  assert(valid_);
58  for(auto d=0u; d!=NDIM; ++d) n_[d] += other.n_[d];
59  return *this;
60  }
61  CartesianMultipoleQuanta& operator-=(const CartesianMultipoleQuanta& other) {
62  assert(valid_);
63  for(auto d=0u; d!=NDIM; ++d) n_[d] -= other.n_[d];
64  return *this;
65  }
66  ~CartesianMultipoleQuanta() = default;
67 
69  unsigned int operator[](unsigned int xyz) const {
70  assert(xyz < NDIM);
71  return n_[xyz];
72  }
74  void inc(unsigned int xyz, unsigned int c = 1u) {
75  assert(xyz < NDIM);
76  assert(valid_);
77  n_[xyz] += c;
78  }
80  void dec(unsigned int xyz, unsigned int c = 1u) {
81  assert(xyz < NDIM);
82  //assert(valid_);
83  if (n_[xyz] >= c)
84  n_[xyz] -= c;
85  else
86  valid_ = false;
87  }
89  unsigned int norm() const {
90  return std::accumulate(n_, n_+NDIM, 0u);
91  }
93  bool zero() const { return norm() == 0; }
95  bool valid() const { return valid_; }
97  LIBINT2_UINT_LEAST64 key() const {
98  if (NDIM == 3u) {
99  unsigned nxy = n_[1] + n_[2];
100  unsigned l = nxy + n_[0];
101  LIBINT2_UINT_LEAST64 key = nxy*(nxy+1)/2 + n_[2];
102  const auto result = key + key_l_offset.at(l);
103  assert(result < max_key);
104  return result;
105  }
106  if (NDIM == 1u) {
107  const auto result = n_[0];
108  assert(result < max_key);
109  return result;
110  }
111  assert(false);
112  }
114  std::string label() const {
115  char result[NDIM+1];
116  for(auto xyz=0u; xyz<NDIM; ++xyz) result[xyz] = '0' + n_[xyz];
117  result[NDIM] = '\0';
118  return std::string(result);
119  }
120 
121  /* ---------------
122  * key computation
123  * --------------- */
124  const constexpr static unsigned max_qn = LIBINT_CARTGAUSS_MAX_AM;
125 
126  // nkeys[L] is the number of possible CartesianMultipoleQuanta with \c norm() == L
127  // \note for NDIM=1 nkeys[L] = 1
128  // \note for NDIM=3 nkeys[L] = (L+1)(L+2)/2
129 
132  const static unsigned max_key = NDIM == 3 ? (1 + max_qn)*(2 + max_qn)*(3 + max_qn)/6 : (1+max_qn);
133 
135  void print(std::ostream& os = std::cout) const;
136 
137  protected:
138 
140  void invalidate() { valid_ = false; }
141 
142  private:
143  unsigned int n_[NDIM];
144  bool valid_; // indicates valid/invalid state
148 
149  };
150 
151  // explicit instantiation declaration, definitions are multipole.cc
152  template<>
153  std::array<LIBINT2_UINT_LEAST64, CartesianMultipoleQuanta<1u>::max_qn+1> CartesianMultipoleQuanta<1u>::key_l_offset;
154  template<>
155  std::array<LIBINT2_UINT_LEAST64, CartesianMultipoleQuanta<3u>::max_qn+1> CartesianMultipoleQuanta<3u>::key_l_offset;
156 
157  template <unsigned NDIM>
158  CartesianMultipoleQuanta<NDIM> operator-(const CartesianMultipoleQuanta<NDIM>& A, const CartesianMultipoleQuanta<NDIM>& B) {
159  CartesianMultipoleQuanta<NDIM> Diff(A);
160  for(unsigned int xyz=0; xyz<3; ++xyz)
161  Diff.dec(xyz,B[xyz]);
162  return Diff;
163  }
164 
165  template <unsigned NDIM>
166  bool operator==(const CartesianMultipoleQuanta<NDIM>& A, const CartesianMultipoleQuanta<NDIM>& B) {
167  for(unsigned d=0; d!=NDIM; ++d)
168  if (A[d] != B[d])
169  return false;
170  return true;
171  }
172 
174  template <unsigned NDIM> inline bool exists(const CartesianMultipoleQuanta<NDIM>& A) { return A.valid(); }
175 
176 
178 
185  class SphericalMultipoleQuanta : public Hashable<LIBINT2_UINT_LEAST64,ReferToKey> {
186  public:
187  enum Sign {plus, minus};
188  const constexpr static unsigned max_qn = LIBINT_CARTGAUSS_MAX_AM;
189 
194  : SphericalMultipoleQuanta(l, m, m < 0 ? Sign::minus : Sign::plus) {}
196  SphericalMultipoleQuanta(int l, int m, Sign sign)
197  : l_(l),
198  m_(std::abs(m)),
199  sign_(sign),
200  valid_(true),
201  phase_(sign == Sign::plus && m < 0 ? -1 : 1) {
202  if (l < 0) valid_ = false;
203  if (m_ > l_) valid_ = false;
204  // N^-_{0,0} = 0
205  if (sign_ == Sign::minus && m_ == 0) valid_ = false;
206  }
207 
208  int l() const { assert(valid_); return static_cast<int>(l_); }
209  int m() const { assert(valid_); return static_cast<int>(m_); }
210  Sign sign() const { assert(valid_); return sign_; }
211  bool valid() const { return valid_; }
212  int phase() const { assert(valid_); return phase_; }
214  bool is_precomputed() const { assert(valid_); return sign_ == Sign::plus && l_ == 0 && m_ == 0; }
215  int value() const { assert(is_precomputed()); return 1; }
216 
217  const static unsigned max_key = (1 + max_qn) * (1 + max_qn);
218 
220  LIBINT2_UINT_LEAST64 key() const {
221  assert(valid_);
222  const auto result = l_*l_ + (sign_ == Sign::plus ? (l_ + m_) : (l_ - m_));
223  assert(result < max_key);
224  return result;
225  }
226 
227  private:
228  // N^-_{l,m} is symmetric is stored as N^-_{l,|m|}
229  int l_;
230  int m_;
231  Sign sign_;
232  bool valid_;
233  int phase_;
234  };
235 
237  inline bool exists(const SphericalMultipoleQuanta& A) { return A.valid(); }
238 
240  inline SphericalMultipoleQuanta::Sign flip(const SphericalMultipoleQuanta::Sign& s) {
241  return s == SphericalMultipoleQuanta::Sign::minus ? SphericalMultipoleQuanta::Sign::plus : SphericalMultipoleQuanta::Sign::minus;
242  }
243 
244 } // namespace libint2
245 
246 #endif
247 
void print(std::ostream &os=std::cout) const
Print out the content.
Represents quantum numbers of cartesian multipole operator.
Definition: multipole.h:41
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
Objects of Hashable<T> class provide hashing function key() which computes keys of type KeyType.
Definition: hashable.h:73
SphericalMultipoleQuanta(int l, int m)
constructs if , otherwise constructs
Definition: multipole.h:193
Represents quantum numbers of real spherical multipole operator defined in Eqs.
Definition: multipole.h:185
LIBINT2_UINT_LEAST64 key() const
Implements Hashable<unsigned>::key()
Definition: multipole.h:220
SphericalMultipoleQuanta(int l, int m, Sign sign)
constructs
Definition: multipole.h:196
bool zero() const
norm() == 0
Definition: multipole.h:93
static const unsigned max_key
The range of keys is [0,max_key).
Definition: multipole.h:132
LIBINT2_UINT_LEAST64 key() const
Implements Hashable<unsigned>::key()
Definition: multipole.h:97
void inc(unsigned int xyz, unsigned int c=1u)
Add c quanta along xyz.
Definition: multipole.h:74
bool valid() const
Return false if this object is invalid.
Definition: multipole.h:95
unsigned int operator[](unsigned int xyz) const
returns the number of quanta along xyz
Definition: multipole.h:69
void invalidate()
make this object invalid
Definition: multipole.h:140
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
unsigned int norm() const
Returns the sum of quantum numbers.
Definition: multipole.h:89
bool is_precomputed() const
Definition: multipole.h:214
SphericalMultipoleQuanta::Sign flip(const SphericalMultipoleQuanta::Sign &s)
plus <-> minus
Definition: multipole.h:240
SphericalMultipoleQuanta()
constructs an object in default (unusable) state
Definition: multipole.h:191
std::string label() const
Return a compact label.
Definition: multipole.h:114
void dec(unsigned int xyz, unsigned int c=1u)
Subtract c quanta along xyz. If impossible, invalidate the object, but do not change its quanta!
Definition: multipole.h:80