Program Listing for File ints.h

↰ Return to documentation for file (gbasis/ints.h)

// HORTON: Helpful Open-source Research TOol for N-fermion systems.
// Copyright (C) 2011-2017 The HORTON Development Team
//
// This file is part of HORTON.
//
// HORTON 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.
//
// HORTON 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 this program; if not, see <http://www.gnu.org/licenses/>
//
//--

#ifndef GBASIS_INTS_H_
#define GBASIS_INTS_H_

#include "libint2.h"
#include "calc.h"
#include "iter_pow.h"

class GB2Integral : public GBCalculator {
 protected:
  long shell_type0, shell_type1;
  const double *r0, *r1;
  IterPow2 i2p;
 public:
  explicit GB2Integral(long max_shell_type)
      : GBCalculator(max_shell_type, 1, 2), shell_type0(0), shell_type1(0),
        r0(NULL), r1(NULL), i2p() {}
  void reset(long shell_type0, long shell_type1, const double *r0, const double *r1);

  virtual void add(double coeff, double alpha0, double alpha1, const double *scales0,
                   const double *scales1) = 0;

  void cart_to_pure();

  const long get_shell_type0() const { return shell_type0; }

  const long get_shell_type1() const { return shell_type1; }
};

class GB2OverlapIntegral : public GB2Integral {
 public:
  explicit GB2OverlapIntegral(long max_shell_type) : GB2Integral(max_shell_type) {}

  virtual void
  add(double coeff, double alpha0, double alpha1, const double *scales0, const double *scales1);
};

class GB2KineticIntegral : public GB2Integral {
 public:
  explicit GB2KineticIntegral(long max_shell_type) : GB2Integral(max_shell_type) {}

  virtual void
  add(double coeff, double alpha0, double alpha1, const double *scales0, const double *scales1);
};

class GB2AttractionIntegral : public GB2Integral {
 private:
  double *charges;    // !< Array with values of the nuclear charges.
  double *centers;    // !< The centers where the charges are located.
  long ncharge;       // !< Number of nuclear charges.

  double *work_g0;    // !< Temporary array to store intermediate results.
  double *work_g1;    // !< Temporary array to store intermediate results.
  double *work_g2;    // !< Temporary array to store intermediate results.
  double *work_boys;  // !< Temporary array to store the laplace of the interaction potential.

 public:
  GB2AttractionIntegral(long max_shell_type, double *charges, double *centers, long ncharge);

  ~GB2AttractionIntegral();

  virtual void add(double coeff, double alpha0, double alpha1, const double *scales0,
                   const double *scales1);

  virtual void laplace_of_potential(double gamma, double arg, long mmax, double *output) = 0;
};

class GB2NuclearAttractionIntegral : public GB2AttractionIntegral {
 public:
  explicit GB2NuclearAttractionIntegral(long max_shell_type, double *charges,
                                        double *centers, long ncharge)
      : GB2AttractionIntegral(max_shell_type, charges, centers, ncharge) {}

  virtual void laplace_of_potential(double gamma, double arg, long mmax, double *output);
};

class GB2ErfAttractionIntegral : public GB2AttractionIntegral {
 public:
  GB2ErfAttractionIntegral(long max_shell_type, double *charges, double *centers,
                           long ncharge, double mu)
      : GB2AttractionIntegral(max_shell_type, charges, centers, ncharge), mu(mu) {}

  virtual void laplace_of_potential(double gamma, double arg, long mmax, double *output);

  const double get_mu() const { return mu; }  // !< The range-separation parameter.

 private:
  double mu;  // !< The range-separation parameter.
};

class GB2GaussAttractionIntegral : public GB2AttractionIntegral {
 public:
  GB2GaussAttractionIntegral(long max_shell_type, double *charges, double *centers, long ncharge,
                             double c, double alpha)
      : GB2AttractionIntegral(max_shell_type, charges, centers, ncharge), c(c),
        alpha(alpha) {}

  virtual void laplace_of_potential(double gamma, double arg, long mmax, double *output);

  const double get_c() const { return c; }  // !< Coefficient of the gaussian.
  const double get_alpha() const { return alpha; }  // !< Exponential parameter of the gaussian.

 private:
  double c;  // !< Coefficient of the gaussian.
  double alpha;  // !< Exponential parameter of the gaussian.
};

class GB2MomentIntegral : public GB2Integral {
 private:
  long *xyz;          // !< Powers for x, y and z of the multipole moment.
  double *center;     // !< The origin w.r.t. to which the multipole moment is computed.

 public:
  GB2MomentIntegral(long max_shell_type, long *xyz, double *center);

  virtual void add(double coeff, double alpha0, double alpha1,
                   const double *scales0, const double *scales1);
};

class GB4Integral : public GBCalculator {
 public:
  explicit GB4Integral(long max_shell_type)
    : GBCalculator(max_shell_type, 1, 4), shell_type0(0), shell_type1(0),
      shell_type2(0), shell_type3(0), r0(NULL), r1(NULL), r2(NULL), r3(NULL) {}

  virtual void reset(long shell_type0, long shell_type1, long shell_type2,
                     long shell_type3, const double *r0, const double *r1,
                     const double *r2, const double *r3);

  virtual void add(double coeff, double alpha0, double alpha1, double alpha2,
                   double alpha3, const double *scales0, const double *scales1,
                   const double *scales2, const double *scales3) = 0;

  // ! Transform the results in the work array from Cartesian to pure functions where needed.
  void cart_to_pure();

  const long get_shell_type0() const { return shell_type0; }  // !< Shell type of contraction 0
  const long get_shell_type1() const { return shell_type1; }  // !< Shell type of contraction 1
  const long get_shell_type2() const { return shell_type2; }  // !< Shell type of contraction 2
  const long get_shell_type3() const { return shell_type3; }  // !< Shell type of contraction 3

 protected:
  long shell_type0;  // !< Shell type of contraction 0
  long shell_type1;  // !< Shell type of contraction 1
  long shell_type2;  // !< Shell type of contraction 2
  long shell_type3;  // !< Shell type of contraction 3
  const double *r0;  // !< Center of contraction 0
  const double *r1;  // !< Center of contraction 1
  const double *r2;  // !< Center of contraction 2
  const double *r3;  // !< Center of contraction 3
};

typedef struct {
  unsigned int am;  // !< Shell type in LibInt conventions.
  const double *r;  // !< Center of a primitive shell.
  double alpha;     // !< Exponent of a primitive shell.
} libint_arg_t;

class GB4IntegralLibInt : public GB4Integral {
 public:
  explicit GB4IntegralLibInt(long max_shell_type);

  ~GB4IntegralLibInt();

  virtual void reset(long shell_type0, long shell_type1, long shell_type2, long shell_type3,
                     const double *r0, const double *r1, const double *r2, const double *r3);

  virtual void add(double coeff, double alpha0, double alpha1, double alpha2, double alpha3,
                   const double *scales0, const double *scales1, const double *scales2,
                   const double *scales3);

  virtual void laplace_of_potential(double prefac, double rho, double t, long mmax,
                                    double *output) = 0;

 private:
  Libint_eri_t erieval;         // !< LibInt runtime object.
  libint_arg_t libint_args[4];  // !< Arguments (shell info) for libint.
  long order[4];                // !< Re-ordering of shells for compatibility with LibInt.
  double ab[3];                 // !< Relative vector from shell 2 to 0 (LibInt order).
  double cd[3];                 // !< Relative vector from shell 3 to 1 (LibInt order).
  double ab2;                   // !< Norm squared of ab.
  double cd2;                   // !< Norm squared of cd.
};

class GB4ElectronRepulsionIntegralLibInt : public GB4IntegralLibInt {
 public:
  explicit GB4ElectronRepulsionIntegralLibInt(long max_shell_type)
      : GB4IntegralLibInt(max_shell_type) {}

  virtual void laplace_of_potential(double prefac, double rho, double t, long mmax,
                                    double *output);
};

class GB4ErfIntegralLibInt : public GB4IntegralLibInt {
 public:
  GB4ErfIntegralLibInt(long max_shell_type, double mu)
      : GB4IntegralLibInt(max_shell_type), mu(mu) {}

  virtual void laplace_of_potential(double prefac, double rho, double t, long mmax,
                                    double *output);

  const double get_mu() const { return mu; }  // !< The range-separation parameter.

 private:
  double mu;  // !< The range-separation parameter.
};

class GB4GaussIntegralLibInt : public GB4IntegralLibInt {
 public:
  GB4GaussIntegralLibInt(long max_shell_type, double c, double alpha)
      : GB4IntegralLibInt(max_shell_type), c(c), alpha(alpha) {}

  virtual void laplace_of_potential(double prefac, double rho, double t, long mmax,
                                    double *output);

  const double get_c() const { return c; }  // !< Coefficient of the gaussian.
  const double get_alpha() const { return alpha; }  // !< Exponential parameter of the gaussian.

 private:
  double c;  // !< Coefficient of the gaussian.
  double alpha;  // !< Exponential parameter of the gaussian.
};

class GB4RAlphaIntegralLibInt : public GB4IntegralLibInt {
 public:
  GB4RAlphaIntegralLibInt(long max_shell_type, double alpha)
      : GB4IntegralLibInt(max_shell_type), alpha(alpha) {}

  virtual void laplace_of_potential(double prefac, double rho, double t, long mmax,
                                    double *output);

  const double get_alpha() const { return alpha; }  // !< The power of r.

 private:
  double alpha;  // !< The power of r.
};

class GB4DeltaIntegralLibInt : public GB4IntegralLibInt {
 public:
  explicit GB4DeltaIntegralLibInt(long max_shell_type)
      : GB4IntegralLibInt(max_shell_type) {}

  virtual void laplace_of_potential(double prefac, double rho, double t, long mmax,
                                    double* output);
};


class GB4DIntegralLibInt : public GB4Integral {
 public:
  explicit GB4DIntegralLibInt(long max_shell_type);
  ~GB4DIntegralLibInt();

  virtual void reset(long shell_type0, long shell_type1, long shell_type2, long shell_type3,
                     const double* r0, const double* r1, const double* r2, const double* r3);
  virtual void add(double coeff, double alpha0, double alpha1, double alpha2, double alpha3,
                   const double* scales0, const double* scales1, const double* scales2,
                   const double* scales3);

  virtual void laplace_of_potential(double prefac, double rho, double t, double* p,
                                    double* q, long mmax, double* output) = 0;

 private:
  Libint_eri_t erieval;
  libint_arg_t libint_args[4];
  long order[4];
  double ab[3];
  double cd[3];
  double ab2;
  double cd2;
};

class GB4IntraDensIntegralLibInt : public GB4DIntegralLibInt {
 public:
  GB4IntraDensIntegralLibInt(long max_shell_type, double* point)
      : GB4DIntegralLibInt(max_shell_type), point(point) {}

  virtual void laplace_of_potential(double prefac, double rho, double t, double* p,
                                    double* q, long mmax, double* output);

 private:
  double* point;
};

#endif  // GBASIS_INTS_H_