multipole.h

/*
 *  Copyright (C) 2004-2021 Edward F. Valeev
 *
 *  This file is part of Libint.
 *
 *  Libint is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  Libint is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with Libint.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef _libint2_src_bin_libint_multipole_h_
#define _libint2_src_bin_libint_multipole_h_
 
#include <iostream>
#include <string>
#include <cassert>
#include <numeric>
#include <sstream>
#include <stdarray.h>
#include <smart_ptr.h>
#include <polyconstr.h>
#include <hashable.h>
#include <global_macros.h>
#include <util_types.h>
 
namespace libint2 {
 
  template <unsigned NDIM = 3>
  class CartesianMultipoleQuanta : public Hashable<LIBINT2_UINT_LEAST64,ReferToKey> {
 
    static_assert(NDIM == 1 || NDIM == 3, "CartesianMultipoleQuanta<NDIM>: only NDIM=1 or NDIM=3 are supported");
  public:
    CartesianMultipoleQuanta() : valid_(true) {
      for(auto d=0u; d!=NDIM; ++d) n_[d] = 0u;
    }
    CartesianMultipoleQuanta(const CartesianMultipoleQuanta& other) : valid_(true) {
      std::copy(other.n_, other.n_ + NDIM, n_);
    }
    CartesianMultipoleQuanta& operator=(const CartesianMultipoleQuanta& other) {
      valid_ = other.valid_;
      std::copy(other.n_, other.n_ + NDIM, n_);
      return *this;
    }
    CartesianMultipoleQuanta& operator+=(const CartesianMultipoleQuanta& other) {
      assert(valid_);
      for(auto d=0u; d!=NDIM; ++d) n_[d] += other.n_[d];
      return *this;
    }
    CartesianMultipoleQuanta& operator-=(const CartesianMultipoleQuanta& other) {
      assert(valid_);
      for(auto d=0u; d!=NDIM; ++d) n_[d] -= other.n_[d];
      return *this;
    }
    ~CartesianMultipoleQuanta() = default;
 
    unsigned int operator[](unsigned int xyz) const {
      assert(xyz < NDIM);
      return n_[xyz];
    }
    void inc(unsigned int xyz, unsigned int c = 1u) {
      assert(xyz < NDIM);
      assert(valid_);
      n_[xyz] += c;
    }
    void dec(unsigned int xyz, unsigned int c = 1u) {
      assert(xyz < NDIM);
      //assert(valid_);
      if (n_[xyz] >= c)
        n_[xyz] -= c;
      else
        valid_ = false;
    }
    unsigned int norm() const {
      return std::accumulate(n_, n_+NDIM, 0u);
    }
    bool zero() const { return norm() == 0; }
    bool valid() const { return valid_; }
    LIBINT2_UINT_LEAST64 key() const {
      if (NDIM == 3u) {
        unsigned nxy = n_[1] + n_[2];
        unsigned l = nxy + n_[0];
        LIBINT2_UINT_LEAST64 key = nxy*(nxy+1)/2 + n_[2];
        const auto result = key + key_l_offset.at(l);
        assert(result < max_key);
        return result;
      }
      if (NDIM == 1u) {
        const auto result = n_[0];
        assert(result < max_key);
        return result;
      }
      assert(false);
    }
    std::string label() const {
      char result[NDIM+1];
      for(auto xyz=0u; xyz<NDIM; ++xyz) result[xyz] = '0' + n_[xyz];
      result[NDIM] = '\0';
      return std::string(result);
    }
 
    /* ---------------
     * key computation
     * --------------- */
    constexpr static unsigned max_qn = LIBINT_CARTGAUSS_MAX_AM;
 
    // nkeys[L] is the number of possible CartesianMultipoleQuanta with \c norm() == L
    // \note for NDIM=1 nkeys[L] = 1
    // \note for NDIM=3 nkeys[L] = (L+1)(L+2)/2
 
    const static unsigned max_key = NDIM == 3 ? (1 + max_qn)*(2 + max_qn)*(3 + max_qn)/6 : (1+max_qn);
 
    void print(std::ostream& os = std::cout) const;
 
  protected:
 
    void invalidate() { valid_ = false; }
 
  private:
    unsigned int n_[NDIM];
    bool valid_;  // indicates valid/invalid state
    static std::array<LIBINT2_UINT_LEAST64, CartesianMultipoleQuanta::max_qn+1> key_l_offset;
 
  };
 
  // explicit instantiation declaration, definitions are multipole.cc
  template<>
  std::array<LIBINT2_UINT_LEAST64, CartesianMultipoleQuanta<1u>::max_qn+1> CartesianMultipoleQuanta<1u>::key_l_offset;
  template<>
  std::array<LIBINT2_UINT_LEAST64, CartesianMultipoleQuanta<3u>::max_qn+1> CartesianMultipoleQuanta<3u>::key_l_offset;
 
  template <unsigned NDIM>
  CartesianMultipoleQuanta<NDIM> operator-(const CartesianMultipoleQuanta<NDIM>& A, const CartesianMultipoleQuanta<NDIM>& B) {
    CartesianMultipoleQuanta<NDIM> Diff(A);
    for(unsigned int xyz=0; xyz<3; ++xyz)
      Diff.dec(xyz,B[xyz]);
    return Diff;
  }
 
  template <unsigned NDIM>
  bool operator==(const CartesianMultipoleQuanta<NDIM>& A, const CartesianMultipoleQuanta<NDIM>& B) {
    for(unsigned d=0; d!=NDIM; ++d)
      if (A[d] != B[d])
        return false;
    return true;
  }
 
  template <unsigned NDIM> inline bool exists(const CartesianMultipoleQuanta<NDIM>& A) { return A.valid(); }
 
 
 
    class SphericalMultipoleQuanta : public Hashable<LIBINT2_UINT_LEAST64,ReferToKey> {
     public:
      enum Sign {plus, minus};
      constexpr static unsigned max_qn = LIBINT_CARTGAUSS_MAX_AM;
 
      SphericalMultipoleQuanta() : SphericalMultipoleQuanta(0,1) {}
      SphericalMultipoleQuanta(int l, int m)
          : SphericalMultipoleQuanta(l, m, m < 0 ? Sign::minus : Sign::plus) {}
      SphericalMultipoleQuanta(int l, int m, Sign sign)
          : l_(l),
            m_(std::abs(m)),
            sign_(sign),
            valid_(true),
            phase_(sign == Sign::plus && m < 0 ? -1 : 1) {
        if (l < 0) valid_ = false;
        if (m_ > l_) valid_ = false;
        // N^-_{0,0} = 0
        if (sign_ == Sign::minus && m_ == 0) valid_ = false;
      }
 
      int l() const { assert(valid_); return static_cast<int>(l_); }
      int m() const { assert(valid_); return static_cast<int>(m_); }
      Sign sign() const { assert(valid_); return sign_; }
      bool valid() const { return valid_; }
      int phase() const { assert(valid_); return phase_; }
      bool is_precomputed() const { assert(valid_); return sign_ == Sign::plus && l_ == 0 && m_ == 0; }
      int value() const { assert(is_precomputed()); return 1; }
 
      const static unsigned max_key = (1 + max_qn) * (1 + max_qn);
 
      LIBINT2_UINT_LEAST64 key() const {
        assert(valid_);
        const auto result = l_*l_ + (sign_ == Sign::plus ? (l_ + m_) : (l_ - m_));
        assert(result < max_key);
        return result;
      }
 
     private:
      // N^-_{l,m} is symmetric is stored as N^-_{l,|m|}
      int l_;
      int m_;
      Sign sign_;
      bool valid_;
      int phase_;
    };
 
    inline bool exists(const SphericalMultipoleQuanta& A) { return A.valid(); }
 
    inline SphericalMultipoleQuanta::Sign flip(const SphericalMultipoleQuanta::Sign& s) {
      return s == SphericalMultipoleQuanta::Sign::minus ? SphericalMultipoleQuanta::Sign::plus : SphericalMultipoleQuanta::Sign::minus;
    }
 
}  // namespace libint2
 
#endif