src/nonbonded/SlaterIntegrals.hpp

/***************************************************************************
 *   This program 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.                                   *
 *                                                                         *
 *   This program 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/>.  *
 ***************************************************************************/

/***
 * This file was imported from qtpie, found at http://code.google.com/p/qtpie
 * and was modified minimally for use in OpenMD.
 *
 * No author attribution was found in the code, but it presumably is
 * the work of J. Chen and Todd J. Martinez.
 *
 * QTPIE (charge transfer with polarization current equilibration) is
 * a new charge model, similar to other charge models like QEq,
 * fluc-q, EEM or ABEEM. Unlike other existing charge models, however,
 * it is capable of describing both charge transfer and polarization
 * phenomena. It is also unique for its ability to describe
 * intermolecular charge transfer at reasonable computational cost.
 *
 * Good references to cite when using this code are:
 *
 * J. Chen and T. J. Martinez, "QTPIE: Charge transfer with
 * polarization current equalization. A fluctuating charge model with
 * correct asymptotics", Chemical Physics Letters 438 (2007),
 * 315-320. DOI: 10.1016/j.cplett.2007.02.065. Erratum: ibid, 463
 * (2008), 288. DOI: 10.1016/j.cplett.2008.08.060
 *
 * J. Chen, D. Hundertmark and T. J. Martinez, "A unified theoretical
 * framework for fluctuating-charge models in atom-space and in
 * bond-space", Journal of Chemical Physics 129 (2008), 214113. DOI:
 * 10.1063/1.3021400.
 *
 * J. Chen and T. J. Martinez, "Charge conservation in
 * electronegativity equalization and its implications for the
 * electrostatic properties of fluctuating-charge models", Journal of
 * Chemical Physics 131 (2009), 044114. DOI: 10.1063/1.3183167
 *
 * J. Chen and T. J. Martinez, "The dissociation catastrophe in
 * fluctuating-charge models and its implications for the concept of
 * atomic electronegativity", Progress in Theoretical Chemistry and
 * Physics, to appear. arXiv:0812.1543
 *
 * J. Chen, "Theory and applications of fluctuating-charge models",
 * PhD (chemical physics) thesis, University of Illinois at
 * Urbana-Champaign, Department of Chemistry, 2009.
 *
 * J. Chen and T. J. Martinez, "Size-extensive polarizabilities with
 * intermolecular charge transfer in a fluctuating-charge model", in
 * preparation. arXiv:0812.1544
 */

#include "config.h"
#include <cmath>
#include <cstdlib>
#include <iostream>
#include "math/Factorials.hpp"
#include "utils/NumericConstant.hpp"

#ifndef NONBONDED_SLATERINTEGRALS_HPP
#define NONBONDED_SLATERINTEGRALS_HPP

template <typename T> inline T sqr(T t) { return t*t; }
template <typename T> inline T mod(T x, T m)
{ return x<0 ? m - 1 - ((-x) - 1)%m : x%m; }

// #include "parameters.h"

/**
 * @brief Computes Rosen's Guillimer-Zener function A
 *  Computes Rosen's A integral, an auxiliary quantity needed to
 *  compute integrals involving Slater-type orbitals of s symmetry.
 * \f[
 *   A_n(\alpha) = \int_1^\infty x^n e^{-\alpha x}dx
 *   = \frac{n! e^{-\alpha}}{\alpha^{n+1}}\sum_{\nu=0}^n
 *   \frac{\alpha^\nu}{\nu!}
 * \f]
 * @param n - principal quantum number
 * @param alpha - Slater exponent 
 * @return the value of Rosen's A integral
 * @note N. Rosen, Phys. Rev., 38 (1931), 255
 */
inline RealType RosenA(int n, RealType a)
{
  RealType RosenA_ = 0.;
  if (a != 0.)
    {
      RealType Term = 1.;
      RosenA_ = Term;
      for (int nu=1; nu<=n; nu++)
        {
          Term *= a / nu;
          RosenA_ += Term;
        }
      RosenA_ = (RosenA_/Term) * (exp(-a)/a);
    }
  return RosenA_;
}

/**
 * @brief Computes Rosen's Guillimer-Zener function B
 *  Computes Rosen's B integral, an auxiliary quantity needed to
 *  compute integrals involving Slater-type orbitals of s symmetry.
 * \f[
 *  B_n(\alpha) = \int_{-1}^1 x^n e^{-\alpha x} dx
 *              = \frac{n!}{\alpha^{n+1}} 
 *                \sum_{\nu=0}^n \frac{e^\alpha(-\alpha)^\nu
 *                  - e^{-\alpha} \alpha^\nu}{\nu!}
 * \f]
 * @param n - principal quantum number
 * @param alpha - Slater exponent 
 * @return the value of Rosen's B integral
 * @note N. Rosen, Phys. Rev., 38 (1931), 255
 */
inline RealType RosenB(int n, RealType alpha)
{
  RealType TheSum, Term;
  RealType RosenB_, PSinhRosenA, PCoshRosenA, PHyperRosenA;
  int nu;
  bool IsPositive;
  if (alpha != 0.)
    {
      Term = 1.;
      TheSum = 1.;
      IsPositive = true;
                
      // These two expressions are (up to constant factors) equivalent
      // to computing the hyperbolic sine and cosine of a respectively
      // The series consists of adding up these terms in an alternating fashion
      PSinhRosenA =  exp(alpha) - exp(-alpha);
      PCoshRosenA = -exp(alpha) - exp(-alpha);
      TheSum = PSinhRosenA;
      for (unsigned nu=1; nu<=n; nu++)
        {
          if (IsPositive)
            {
              PHyperRosenA = PCoshRosenA;
              IsPositive = false;
            }
          else // term to add should be negative
            {
              PHyperRosenA = PSinhRosenA;
              IsPositive = true;
            }
          Term *= alpha / nu;
          TheSum += Term * PHyperRosenA;
        }
      RosenB_ = TheSum / (alpha*Term);
    }
  else // pathological case of a=0
    {
      printf("WARNING, a = 0 in RosenB\n");
      RosenB_ = (1. - pow(-1., n)) / (n + 1.);
    }
  return RosenB_;
}

/** @brief Computes Rosen's D combinatorial factor
 *  Computes Rosen's D factor, an auxiliary quantity needed to
 *  compute integrals involving Slater-type orbitals of s symmetry.
 *  \f[
 *  RosenD^{mn}_p = \sum_k (-1)^k \frac{m! n!}
 *                  {(p-k)!(m-p+k)!(n-k)!k!}
 *  \f]
 * @return the value of Rosen's D factor
 * @note N. Rosen, Phys. Rev., 38 (1931), 255
 */
inline RealType RosenD(int m, int n, int p)
{
  if (m+n+p > maxFact)
    {
      printf("Error, arguments exceed maximum factorial computed %d > %d\n", m+n+p, maxFact);
      ::exit(0);
    }
        
  RealType RosenD_ = 0;
  for (int k=max(p-m,0); k<=min(n,p); k++)
    {
      if (mod(k,2) == 0)
        RosenD_ += (fact[m] / (fact[p-k] * fact[m-p+k])) * (fact[n] / (fact[n-k] * fact[k]));
      else
        RosenD_ -= (fact[m] / ( fact[p-k] * fact[m-p+k])) * (fact[n] / (fact[n-k] * fact[k]));
    }
  return RosenD_;
}

/** @brief Computes Coulomb integral analytically over s-type STOs
 * Computes the two-center Coulomb integral over Slater-type orbitals of s symmetry.
 * @param a : Slater zeta exponent of first atom in inverse Bohr (au)
 * @param b : Slater zeta exponent of second atom in inverse Bohr (au)
 * @param m : principal quantum number of first atom
 * @param n : principal quantum number of second atom
 * @param R : internuclear distance in atomic units (bohr)
 * @return value of the Coulomb potential energy integral
 * @note N. Rosen, Phys. Rev., 38 (1931), 255
 * @note In Rosen's paper, this integral is known as K2.
 */
inline RealType sSTOCoulInt(RealType a, RealType b, int m, int n, RealType R)
{
  int nu, p;
  RealType x, K2;
  RealType Factor1, Factor2, Term, OneElectronTerm;
  RealType eps, epsi;
        
  // To speed up calculation, we terminate loop once contributions
  // to integral fall below the bound, epsilon
  RealType epsilon = 0.;

  // x is the argument of the auxiliary RosenA and RosenB functions
  x = 2. * a * R;
        
  // First compute the two-electron component
  RealType sSTOCoulInt_ = 0.;
  if (std::fabs(x) < OpenMD::NumericConstant::epsilon) // Pathological case
    {

      // This solution for the one-center coulomb integrals comes from 
      // Yoshiyuki Hase, Computers & Chemistry 9(4), pp. 285-287 (1985).
      
      RealType Term1 = fact[2*m - 1] / pow(2*a, 2*m);
      RealType Term2 = 0.;
      for (int nu = 1; nu <= 2*n; nu++) {
        Term2 += nu * pow(2*b, 2*n - nu) * fact[2*(m+n)-nu-1] / (fact[2*n-nu]*2*n * pow(2*(a+b), 2*(m+n)-nu));
      }
      sSTOCoulInt_ = pow(2*a, 2*m+1) * (Term1 - Term2) / fact[2*m];

      // Original QTPIE code for the one-center case is below.  Doesn't
      // appear to generate the correct one-center results.
      //
      // if ((a==b) && (m==n))
      //   {
      //     for (int nu=0; nu<=2*n-1; nu++)
      //       {
      //         K2 = 0.;
      //         for (unsigned p=0; p<=2*n+m; p++) K2 += 1. / fact[p];
      //         sSTOCoulInt_ += K2 * fact[2*n+m] / fact[m];
      //       }
      //     sSTOCoulInt_ = 2 * a / (n * fact[2*n]) * sSTOCoulInt_;
      //   }
      // else
      //   {
      //     // Not implemented
      //     printf("ERROR, sSTOCoulInt cannot compute from arguments\n");
      //     printf("a = %lf b = %lf m = %d n = %d R = %lf\n",a, b, m, n, R);
      //     exit(0);
      //   }

    }
  else
    {
      OneElectronTerm = 1./R + pow(x, 2*m)/(fact[2*m]*R)*
        ((x-2*m)*RosenA(2*m-1,x)-exp(-x)) + sSTOCoulInt_;
      eps = epsilon / OneElectronTerm;
      if (a == b)
        {
          // Apply Rosen (48)
          Factor1 = -a*pow(a*R, 2*m)/(n*fact[2*m]);
          for (int nu=0; nu<=2*n-1; nu++)
            {
              Factor2 = (2.*n-nu)/fact[nu]*pow(a*R,nu);
              epsi = eps / fabs(Factor1 * Factor2);
              K2 = 0.;
              for (int p=0; p<=m+(nu-1)/2; p++)
                {
                  Term = RosenD(2*m-1, nu, 2*p)/(2.*p+1.) *RosenA(2*m+nu-1-2*p,x);
                  K2 += Term;
                  if ((Term > 0) && (Term < epsi)) goto label1;
                }
              sSTOCoulInt_ +=  K2 * Factor2;
            }
        label1:
          sSTOCoulInt_ *= Factor1;
        }
      else
        {
          Factor1 = -a*pow(a*R,2*m)/(2.*n*fact[2*m]);
          epsi = eps/fabs(Factor1);
          if (b == 0.)
            printf("WARNING: b = 0 in sSTOCoulInt\n");
          else
            {
              // Apply Rosen (54)
              for (int nu=0; nu<=2*n-1; nu++)
                {
                  K2 = 0;
                  for (int p=0; p<=2*m+nu-1; p++)
                    K2=K2+RosenD(2*m-1,nu,p)*RosenB(p,R*(a-b))
                      *RosenA(2*m+nu-1-p,R*(a+b));
                  Term = K2*(2*n-nu)/fact[nu]*pow(b*R, nu);
                  sSTOCoulInt_ += Term;
                  if (fabs(Term) < epsi) goto label2;
                }
            label2:
              sSTOCoulInt_ *= Factor1;
            }
        }
      // Now add the one-electron term from Rosen (47) = Rosen (53)
      sSTOCoulInt_ += OneElectronTerm;
    }
  return sSTOCoulInt_;
}

/**
 * @brief Computes overlap integral analytically over s-type STOs
 *   Computes the overlap integral over two
 *   Slater-type orbitals of s symmetry.
 * @param a : Slater zeta exponent of first atom in inverse Bohr (au)
 * @param b : Slater zeta exponent of second atom in inverse Bohr (au)
 * @param m : principal quantum number of first atom
 * @param n : principal quantum number of second atom
 * @param R : internuclear distance in atomic units (bohr)
 * @return the value of the sSTOOvInt integral
 * @note N. Rosen, Phys. Rev., 38 (1931), 255
 * @note In the Rosen paper, this integral is known as I.
 */
inline RealType sSTOOvInt(RealType a, RealType b, int m, int n, RealType R)
{
  RealType Factor, Term, eps;
        
  // To speed up calculation, we terminate loop once contributions
  // to integral fall below the bound, epsilon
  RealType epsilon = 0.;
  RealType sSTOOvInt_ = 0.;
        
  if (a == b)
    {
      Factor = pow(a*R, m+n+1)/sqrt(fact[2*m]*fact[2*n]);
      eps = epsilon / fabs(Factor);
      for (int q=0; q<=(m+n)/2; q++)
        {
          Term = RosenD(m,n,2*q)/(2.*q+1.)*RosenA(m+n-2*q,a*R);
          sSTOOvInt_ += Term;
          if (fabs(Term) < eps) exit(0);
        }
      sSTOOvInt_ *= Factor;
    }
  else
    {
      Factor = 0.5*pow(a*R, m+0.5)*pow(b*R,n+0.5)
        /sqrt(fact[2*m]*fact[2*n]);
      eps = epsilon / fabs(Factor);
      for (int q=0; q<=m+n; q++)
        {
          Term = RosenD(m,n,q)*RosenB(q, R/2.*(a-b))
            * RosenA(m+n-q,R/2.*(a+b));
          sSTOOvInt_ += Term;
          if (fabs(Term) < eps) exit(0);
        }
      sSTOOvInt_ *= Factor;
    }
  return sSTOOvInt_;
}

/**
 * @brief Computes kinetic energy integral analytically over s-type STOs
 *   Computes the overlap integral over two Slater-type orbitals of s symmetry.
 * @param a : Slater zeta exponent of first atom in inverse Bohr (au)
 * @param b : Slater zeta exponent of second atom in inverse Bohr (au)
 * @param m : principal quantum number of first atom
 * @param n : principal quantum number of second atom
 * @param R : internuclear distance in atomic units (bohr)
 * @return the value of the kinetic energy integral
 * @note N. Rosen, Phys. Rev., 38 (1931), 255
 * @note untested
 */
inline RealType KinInt(RealType a, RealType b, int m, int n,RealType R)
{
  RealType KinInt_ = -0.5*b*b*sSTOOvInt(a, b, m, n, R);
  if (n > 0)
    {
      KinInt_ += b*b*pow(2*b/(2*b-1),0.5) * sSTOOvInt(a, b, m, n-1, R);
      if (n > 1) KinInt_ += pow(n*(n-1)/((n-0.5)*(n-1.5)), 0.5)
                   * sSTOOvInt(a, b, m, n-2, R);
    }
  return KinInt_;
}

/**
 * @brief Computes derivative of Coulomb integral with respect to the interatomic distance
 *   Computes the two-center Coulomb integral over Slater-type orbitals of s symmetry.
 * @param a: Slater zeta exponent of first atom in inverse Bohr (au)
 * @param b: Slater zeta exponent of second atom in inverse Bohr (au)
 * @param m: principal quantum number of first atom
 * @param n: principal quantum number of second atom
 * @param R: internuclear distance in atomic units (bohr)
 * @return the derivative of the Coulomb potential energy integral
 * @note Derived in QTPIE research notes, May 15 2007
 */
inline RealType sSTOCoulIntGrad(RealType a, RealType b, int m, int n, RealType R)
{
  RealType x, y, z, K2, TheSum;
  // x is the argument of the auxiliary RosenA and RosenB functions
  x = 2. * a * R;
        
  // First compute the two-electron component
  RealType sSTOCoulIntGrad_ = 0.;
  if (x==0) // Pathological case
    {
      printf("WARNING: argument given to sSTOCoulIntGrad is 0\n");
      printf("a = %lf R= %lf\n", a, R);
    }
  else
    {
      if (a == b)
        {
          TheSum = 0.;
          for (int nu=0; nu<=2*(n-1); nu++)
            {
              K2 = 0.;
              for (int p=0; p<=(m+nu)/2; p++)
                K2 += RosenD(2*m-1, nu+1, 2*p)/(2*p + 1.) * RosenA(2*m+nu-1-2*p, x);
              TheSum += (2*n-nu-1)/fact[nu]*pow(a*R, nu) * K2;
            }
          sSTOCoulIntGrad_ = -pow(a, 2*m+2)*pow(R, 2*m) /(n*fact[2*m])*TheSum;
          TheSum = 0.;
          for (int nu=0; nu<=2*n-1; nu++)
            {
              K2 = 0.;
              for (int p=0; p<=(m+nu-1)/2; p++)
                K2 += RosenD(2*m-1, nu, 2*p)/(2*p + 1.) * RosenA(2*m+nu-2*p, x);
              TheSum += (2*n-nu)/fact[nu]*pow(a*R,nu) * K2;
            }
          sSTOCoulIntGrad_ += 2*pow(a, 2*m+2)*pow(R, 2*m) /(n*fact[2*m])*TheSum;
        }
      else
        {
          // Slater exponents are different
          // First calculate some useful arguments
          y = R*(a+b);
          z = R*(a-b);
          TheSum = 0.;
          for (int nu=0; nu<=2*n-1; nu++)
            {
              K2 = 0.;
              for (int p=0; p<=2*m+nu; p++)
                K2 += RosenD(2*m-1, nu+1, p)
                  * RosenB(p,z)*RosenA(2*m+nu-p, y);
              TheSum += (2*n-nu-1)/fact[nu]*pow(b*R,nu) * K2;
            }
          sSTOCoulIntGrad_ = -b*pow(a,2*m+1)*pow(R,2*m)/
            (2*n*fact[2*m])*TheSum;
          TheSum = 0.;
          for (int nu=0; nu<=2*n; nu++)
            {
              K2 = 0.;
              for (int p=0; p<=2*m-1+nu; p++)
                K2 += RosenD(2*m-1, nu, p)
                  * ((a-b)*RosenB(p+1,z)*RosenA(2*m+nu-p-1, y)
                     +(a+b)*RosenB(p  ,z)*RosenA(2*m+nu-p  , y));
              TheSum += (2*n-nu)/fact[nu]*pow(b*R,nu) * K2;
            }
          sSTOCoulIntGrad_ += pow(a,2*m+1)*pow(R,2*m)/(2*n*fact[2*m])*TheSum;
        }
      // Now add one-electron terms and common term
      sSTOCoulIntGrad_ = sSTOCoulIntGrad_ - (2.*m+1.)/sqr(R)
        + 2.*m/R * sSTOCoulInt(a,b,m,n,R)
        + pow(x,2*m)/(fact[2*m]*sqr(R)) * ((2.*m+1.)*exp(-x)
                                           + 2.*m*(1.+2.*m-x)*RosenA(2*m-1,x));
    }
  return sSTOCoulIntGrad_;
}

/** 
 * @brief Computes gradient of overlap integral with respect to the interatomic diatance
 *   Computes the derivative of the overlap integral over two Slater-type orbitals of s symmetry.
 * @param a: Slater zeta exponent of first atom in inverse Bohr (au)
 * @param b: Slater zeta exponent of second atom in inverse Bohr (au)
 * @param m: principal quantum number of first atom
 * @param n: principal quantum number of second atom
 * @param R: internuclear distance in atomic units (bohr)
 * @return the derivative of the sSTOOvInt integral
 * @note Derived in QTPIE research notes, May 15 2007
 */
inline RealType sSTOOvIntGrad(RealType a, RealType b, int m, int n, RealType R)
{
  RealType w, x, y, z, TheSum;
        
  // Calculate first term
  RealType sSTOOvIntGrad_ = (m+n+1.)/R * sSTOOvInt(a, b, m, n, R);
        
  // Calculate remaining terms; answers depend on exponents 
  TheSum = 0.;
  x = a * R;
  if (a == b)
    {
      for (int q=0; q<=(m+n)/2; q++)
        TheSum += RosenD(m,n,2*q) / (2*q + 1.) * RosenA(m+n-2*q+1, x);
      sSTOOvIntGrad_ -= a*pow(x,m+n+1)/ sqrt(fact[2*m]*fact[2*n])*TheSum;
    }
  else
    {
      w = b*R;
      y = 0.5*R*(a+b);
      z = 0.5*R*(a-b);
      for (int q=0; q<m+n; q++)
        TheSum = TheSum + RosenD(m,n,q) *
          ((a-b)*RosenB(q+1,z)*RosenA(m+n-q  ,y)
           +(a+b)*RosenB(q  ,z)*RosenA(m+n-q+1,y));
      sSTOOvIntGrad_ -= 0.25*sqrt((pow(x, 2*m+1)*pow(w, 2*n+1))/(fact[2*m]*fact[2*n]))*TheSum;
    }
  return sSTOOvIntGrad_;
}

/**
 * @brief Calculates a Slater-type orbital exponent based on the hardness parameters
 * @param Hardness: chemical hardness in atomic units
 * @param        n: principal quantum number
 * @note Modified for use with OpenMD by Gezelter and Michalka.
 */
inline RealType getSTOZeta(int n, RealType hardness)
{
  //  Approximate the exact value of the constant of proportionality
  //  by its value at a very small distance epsilon
  //  since the exact R = 0 case has not be programmed
  RealType epsilon = 1.0e-8;
  
  // Assign orbital exponent  
  return pow(sSTOCoulInt(1., 1., n, n, epsilon) / hardness, -1./(3. + 2.*n));
}

#endif
Revision:	1766
Committed:	Thu Jul 5 17:08:25 2012 UTC (13 years ago) by gezelter
File size:	19245 byte(s)
Log Message:	Added Fluctuating Charge Langevin propagator, and made it the default fixed some errors on the one-center slater coulomb integrals, and some parameters in PhysicalConstants.
#	User	Rev	Content
1	gezelter	1718	/***************************************************************************
2			* This program is free software; you can redistribute it and/or modify *
3			* it under the terms of the GNU General Public License as published by *
4			* the Free Software Foundation; either version 3 of the License, or *
5			* (at your option) any later version. *
6			* *
7			* This program is distributed in the hope that it will be useful, *
8			* but WITHOUT ANY WARRANTY; without even the implied warranty of *
9			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
10			* GNU General Public License for more details. *
11			* *
12			* You should have received a copy of the GNU General Public License *
13			* along with this program; if not, see <http://www.gnu.org/licenses/>. *
14			***************************************************************************/
15
16			/***
17			* This file was imported from qtpie, found at http://code.google.com/p/qtpie
18			* and was modified minimally for use in OpenMD.
19			*
20			* No author attribution was found in the code, but it presumably is
21			* the work of J. Chen and Todd J. Martinez.
22			*
23			* QTPIE (charge transfer with polarization current equilibration) is
24			* a new charge model, similar to other charge models like QEq,
25			* fluc-q, EEM or ABEEM. Unlike other existing charge models, however,
26			* it is capable of describing both charge transfer and polarization
27			* phenomena. It is also unique for its ability to describe
28			* intermolecular charge transfer at reasonable computational cost.
29			*
30			* Good references to cite when using this code are:
31			*
32			* J. Chen and T. J. Martinez, "QTPIE: Charge transfer with
33			* polarization current equalization. A fluctuating charge model with
34			* correct asymptotics", Chemical Physics Letters 438 (2007),
35			* 315-320. DOI: 10.1016/j.cplett.2007.02.065. Erratum: ibid, 463
36			* (2008), 288. DOI: 10.1016/j.cplett.2008.08.060
37			*
38			* J. Chen, D. Hundertmark and T. J. Martinez, "A unified theoretical
39			* framework for fluctuating-charge models in atom-space and in
40			* bond-space", Journal of Chemical Physics 129 (2008), 214113. DOI:
41			* 10.1063/1.3021400.
42			*
43			* J. Chen and T. J. Martinez, "Charge conservation in
44			* electronegativity equalization and its implications for the
45			* electrostatic properties of fluctuating-charge models", Journal of
46			* Chemical Physics 131 (2009), 044114. DOI: 10.1063/1.3183167
47			*
48			* J. Chen and T. J. Martinez, "The dissociation catastrophe in
49			* fluctuating-charge models and its implications for the concept of
50			* atomic electronegativity", Progress in Theoretical Chemistry and
51			* Physics, to appear. arXiv:0812.1543
52			*
53			* J. Chen, "Theory and applications of fluctuating-charge models",
54			* PhD (chemical physics) thesis, University of Illinois at
55			* Urbana-Champaign, Department of Chemistry, 2009.
56			*
57			* J. Chen and T. J. Martinez, "Size-extensive polarizabilities with
58			* intermolecular charge transfer in a fluctuating-charge model", in
59			* preparation. arXiv:0812.1544
60			*/
61
62			#include "config.h"
63			#include <cmath>
64	gezelter	1749	#include <cstdlib>
65	jmichalk	1734	#include <iostream>
66	gezelter	1718	#include "math/Factorials.hpp"
67	gezelter	1766	#include "utils/NumericConstant.hpp"
68	gezelter	1718
69			#ifndef NONBONDED_SLATERINTEGRALS_HPP
70			#define NONBONDED_SLATERINTEGRALS_HPP
71
72			template <typename T> inline T sqr(T t) { return t*t; }
73			template <typename T> inline T mod(T x, T m)
74			{ return x<0 ? m - 1 - ((-x) - 1)%m : x%m; }
75
76			// #include "parameters.h"
77
78			/**
79			* @brief Computes Rosen's Guillimer-Zener function A
80			* Computes Rosen's A integral, an auxiliary quantity needed to
81			* compute integrals involving Slater-type orbitals of s symmetry.
82			* \f[
83			* A_n(\alpha) = \int_1^\infty x^n e^{-\alpha x}dx
84			* = \frac{n! e^{-\alpha}}{\alpha^{n+1}}\sum_{\nu=0}^n
85			* \frac{\alpha^\nu}{\nu!}
86			* \f]
87			* @param n - principal quantum number
88			* @param alpha - Slater exponent
89			* @return the value of Rosen's A integral
90			* @note N. Rosen, Phys. Rev., 38 (1931), 255
91			*/
92			inline RealType RosenA(int n, RealType a)
93			{
94			RealType RosenA_ = 0.;
95			if (a != 0.)
96			{
97			RealType Term = 1.;
98			RosenA_ = Term;
99	jmichalk	1734	for (int nu=1; nu<=n; nu++)
100	gezelter	1718	{
101			Term *= a / nu;
102			RosenA_ += Term;
103			}
104			RosenA_ = (RosenA_/Term) * (exp(-a)/a);
105			}
106			return RosenA_;
107			}
108
109			/**
110			* @brief Computes Rosen's Guillimer-Zener function B
111			* Computes Rosen's B integral, an auxiliary quantity needed to
112			* compute integrals involving Slater-type orbitals of s symmetry.
113			* \f[
114			* B_n(\alpha) = \int_{-1}^1 x^n e^{-\alpha x} dx
115			* = \frac{n!}{\alpha^{n+1}}
116			* \sum_{\nu=0}^n \frac{e^\alpha(-\alpha)^\nu
117			* - e^{-\alpha} \alpha^\nu}{\nu!}
118			* \f]
119			* @param n - principal quantum number
120			* @param alpha - Slater exponent
121			* @return the value of Rosen's B integral
122			* @note N. Rosen, Phys. Rev., 38 (1931), 255
123			*/
124			inline RealType RosenB(int n, RealType alpha)
125			{
126			RealType TheSum, Term;
127			RealType RosenB_, PSinhRosenA, PCoshRosenA, PHyperRosenA;
128			int nu;
129			bool IsPositive;
130			if (alpha != 0.)
131			{
132			Term = 1.;
133			TheSum = 1.;
134			IsPositive = true;
135
136			// These two expressions are (up to constant factors) equivalent
137			// to computing the hyperbolic sine and cosine of a respectively
138			// The series consists of adding up these terms in an alternating fashion
139			PSinhRosenA = exp(alpha) - exp(-alpha);
140			PCoshRosenA = -exp(alpha) - exp(-alpha);
141			TheSum = PSinhRosenA;
142			for (unsigned nu=1; nu<=n; nu++)
143			{
144			if (IsPositive)
145			{
146			PHyperRosenA = PCoshRosenA;
147			IsPositive = false;
148			}
149			else // term to add should be negative
150			{
151			PHyperRosenA = PSinhRosenA;
152			IsPositive = true;
153			}
154			Term *= alpha / nu;
155			TheSum += Term * PHyperRosenA;
156			}
157			RosenB_ = TheSum / (alpha*Term);
158			}
159			else // pathological case of a=0
160			{
161			printf("WARNING, a = 0 in RosenB\n");
162			RosenB_ = (1. - pow(-1., n)) / (n + 1.);
163			}
164			return RosenB_;
165			}
166
167			/** @brief Computes Rosen's D combinatorial factor
168			* Computes Rosen's D factor, an auxiliary quantity needed to
169			* compute integrals involving Slater-type orbitals of s symmetry.
170			* \f[
171			* RosenD^{mn}_p = \sum_k (-1)^k \frac{m! n!}
172			* {(p-k)!(m-p+k)!(n-k)!k!}
173			* \f]
174			* @return the value of Rosen's D factor
175			* @note N. Rosen, Phys. Rev., 38 (1931), 255
176			*/
177			inline RealType RosenD(int m, int n, int p)
178			{
179			if (m+n+p > maxFact)
180			{
181			printf("Error, arguments exceed maximum factorial computed %d > %d\n", m+n+p, maxFact);
182			::exit(0);
183			}
184
185			RealType RosenD_ = 0;
186			for (int k=max(p-m,0); k<=min(n,p); k++)
187			{
188			if (mod(k,2) == 0)
189			RosenD_ += (fact[m] / (fact[p-k] * fact[m-p+k])) * (fact[n] / (fact[n-k] * fact[k]));
190			else
191			RosenD_ -= (fact[m] / ( fact[p-k] * fact[m-p+k])) * (fact[n] / (fact[n-k] * fact[k]));
192			}
193			return RosenD_;
194			}
195
196			/** @brief Computes Coulomb integral analytically over s-type STOs
197			* Computes the two-center Coulomb integral over Slater-type orbitals of s symmetry.
198			* @param a : Slater zeta exponent of first atom in inverse Bohr (au)
199			* @param b : Slater zeta exponent of second atom in inverse Bohr (au)
200			* @param m : principal quantum number of first atom
201			* @param n : principal quantum number of second atom
202			* @param R : internuclear distance in atomic units (bohr)
203			* @return value of the Coulomb potential energy integral
204			* @note N. Rosen, Phys. Rev., 38 (1931), 255
205			* @note In Rosen's paper, this integral is known as K2.
206			*/
207			inline RealType sSTOCoulInt(RealType a, RealType b, int m, int n, RealType R)
208			{
209			int nu, p;
210			RealType x, K2;
211			RealType Factor1, Factor2, Term, OneElectronTerm;
212			RealType eps, epsi;
213
214			// To speed up calculation, we terminate loop once contributions
215			// to integral fall below the bound, epsilon
216			RealType epsilon = 0.;
217
218			// x is the argument of the auxiliary RosenA and RosenB functions
219			x = 2. * a * R;
220
221			// First compute the two-electron component
222			RealType sSTOCoulInt_ = 0.;
223	gezelter	1766	if (std::fabs(x) < OpenMD::NumericConstant::epsilon) // Pathological case
224	gezelter	1718	{
225	gezelter	1766
226			// This solution for the one-center coulomb integrals comes from
227			// Yoshiyuki Hase, Computers & Chemistry 9(4), pp. 285-287 (1985).
228
229			RealType Term1 = fact[2m - 1] / pow(2a, 2*m);
230			RealType Term2 = 0.;
231			for (int nu = 1; nu <= 2*n; nu++) {
232			Term2 += nu * pow(2b, 2n - nu) * fact[2(m+n)-nu-1] / (fact[2n-nu]2n * pow(2(a+b), 2(m+n)-nu));
233			}
234			sSTOCoulInt_ = pow(2a, 2m+1) * (Term1 - Term2) / fact[2*m];
235
236			// Original QTPIE code for the one-center case is below. Doesn't
237			// appear to generate the correct one-center results.
238			//
239			// if ((a==b) && (m==n))
240			// {
241			// for (int nu=0; nu<=2*n-1; nu++)
242			// {
243			// K2 = 0.;
244			// for (unsigned p=0; p<=2*n+m; p++) K2 += 1. / fact[p];
245			// sSTOCoulInt_ += K2 * fact[2*n+m] / fact[m];
246			// }
247			// sSTOCoulInt_ = 2 * a / (n * fact[2n]) sSTOCoulInt_;
248			// }
249			// else
250			// {
251			// // Not implemented
252			// printf("ERROR, sSTOCoulInt cannot compute from arguments\n");
253			// printf("a = %lf b = %lf m = %d n = %d R = %lf\n",a, b, m, n, R);
254			// exit(0);
255			// }
256
257	gezelter	1718	}
258			else
259			{
260			OneElectronTerm = 1./R + pow(x, 2m)/(fact[2m]R)
261			((x-2m)RosenA(2*m-1,x)-exp(-x)) + sSTOCoulInt_;
262			eps = epsilon / OneElectronTerm;
263			if (a == b)
264			{
265			// Apply Rosen (48)
266			Factor1 = -apow(aR, 2m)/(nfact[2*m]);
267			for (int nu=0; nu<=2*n-1; nu++)
268			{
269			Factor2 = (2.n-nu)/fact[nu]pow(a*R,nu);
270			epsi = eps / fabs(Factor1 * Factor2);
271			K2 = 0.;
272			for (int p=0; p<=m+(nu-1)/2; p++)
273			{
274			Term = RosenD(2m-1, nu, 2p)/(2.p+1.) RosenA(2m+nu-1-2p,x);
275			K2 += Term;
276			if ((Term > 0) && (Term < epsi)) goto label1;
277			}
278			sSTOCoulInt_ += K2 * Factor2;
279			}
280			label1:
281			sSTOCoulInt_ *= Factor1;
282			}
283			else
284			{
285			Factor1 = -apow(aR,2m)/(2.nfact[2m]);
286			epsi = eps/fabs(Factor1);
287			if (b == 0.)
288			printf("WARNING: b = 0 in sSTOCoulInt\n");
289			else
290			{
291			// Apply Rosen (54)
292			for (int nu=0; nu<=2*n-1; nu++)
293			{
294			K2 = 0;
295			for (int p=0; p<=2*m+nu-1; p++)
296			K2=K2+RosenD(2m-1,nu,p)RosenB(p,R*(a-b))
297			RosenA(2m+nu-1-p,R*(a+b));
298			Term = K2(2n-nu)/fact[nu]pow(bR, nu);
299			sSTOCoulInt_ += Term;
300			if (fabs(Term) < epsi) goto label2;
301			}
302			label2:
303			sSTOCoulInt_ *= Factor1;
304			}
305			}
306			// Now add the one-electron term from Rosen (47) = Rosen (53)
307			sSTOCoulInt_ += OneElectronTerm;
308			}
309			return sSTOCoulInt_;
310			}
311
312			/**
313			* @brief Computes overlap integral analytically over s-type STOs
314			* Computes the overlap integral over two
315			* Slater-type orbitals of s symmetry.
316			* @param a : Slater zeta exponent of first atom in inverse Bohr (au)
317			* @param b : Slater zeta exponent of second atom in inverse Bohr (au)
318			* @param m : principal quantum number of first atom
319			* @param n : principal quantum number of second atom
320			* @param R : internuclear distance in atomic units (bohr)
321			* @return the value of the sSTOOvInt integral
322			* @note N. Rosen, Phys. Rev., 38 (1931), 255
323			* @note In the Rosen paper, this integral is known as I.
324			*/
325			inline RealType sSTOOvInt(RealType a, RealType b, int m, int n, RealType R)
326			{
327			RealType Factor, Term, eps;
328
329			// To speed up calculation, we terminate loop once contributions
330			// to integral fall below the bound, epsilon
331			RealType epsilon = 0.;
332			RealType sSTOOvInt_ = 0.;
333
334			if (a == b)
335			{
336			Factor = pow(aR, m+n+1)/sqrt(fact[2m]fact[2n]);
337			eps = epsilon / fabs(Factor);
338			for (int q=0; q<=(m+n)/2; q++)
339			{
340			Term = RosenD(m,n,2q)/(2.q+1.)RosenA(m+n-2q,a*R);
341			sSTOOvInt_ += Term;
342			if (fabs(Term) < eps) exit(0);
343			}
344			sSTOOvInt_ *= Factor;
345			}
346			else
347			{
348			Factor = 0.5pow(aR, m+0.5)pow(bR,n+0.5)
349			/sqrt(fact[2m]fact[2*n]);
350			eps = epsilon / fabs(Factor);
351			for (int q=0; q<=m+n; q++)
352			{
353			Term = RosenD(m,n,q)RosenB(q, R/2.(a-b))
354			* RosenA(m+n-q,R/2.*(a+b));
355			sSTOOvInt_ += Term;
356			if (fabs(Term) < eps) exit(0);
357			}
358			sSTOOvInt_ *= Factor;
359			}
360			return sSTOOvInt_;
361			}
362
363			/**
364			* @brief Computes kinetic energy integral analytically over s-type STOs
365			* Computes the overlap integral over two Slater-type orbitals of s symmetry.
366			* @param a : Slater zeta exponent of first atom in inverse Bohr (au)
367			* @param b : Slater zeta exponent of second atom in inverse Bohr (au)
368			* @param m : principal quantum number of first atom
369			* @param n : principal quantum number of second atom
370			* @param R : internuclear distance in atomic units (bohr)
371			* @return the value of the kinetic energy integral
372			* @note N. Rosen, Phys. Rev., 38 (1931), 255
373			* @note untested
374			*/
375			inline RealType KinInt(RealType a, RealType b, int m, int n,RealType R)
376			{
377			RealType KinInt_ = -0.5bb*sSTOOvInt(a, b, m, n, R);
378			if (n > 0)
379			{
380			KinInt_ += bbpow(2b/(2b-1),0.5) * sSTOOvInt(a, b, m, n-1, R);
381			if (n > 1) KinInt_ += pow(n(n-1)/((n-0.5)(n-1.5)), 0.5)
382			* sSTOOvInt(a, b, m, n-2, R);
383			}
384			return KinInt_;
385			}
386
387			/**
388			* @brief Computes derivative of Coulomb integral with respect to the interatomic distance
389			* Computes the two-center Coulomb integral over Slater-type orbitals of s symmetry.
390			* @param a: Slater zeta exponent of first atom in inverse Bohr (au)
391			* @param b: Slater zeta exponent of second atom in inverse Bohr (au)
392			* @param m: principal quantum number of first atom
393			* @param n: principal quantum number of second atom
394			* @param R: internuclear distance in atomic units (bohr)
395			* @return the derivative of the Coulomb potential energy integral
396			* @note Derived in QTPIE research notes, May 15 2007
397			*/
398			inline RealType sSTOCoulIntGrad(RealType a, RealType b, int m, int n, RealType R)
399			{
400			RealType x, y, z, K2, TheSum;
401			// x is the argument of the auxiliary RosenA and RosenB functions
402			x = 2. * a * R;
403
404			// First compute the two-electron component
405			RealType sSTOCoulIntGrad_ = 0.;
406			if (x==0) // Pathological case
407			{
408			printf("WARNING: argument given to sSTOCoulIntGrad is 0\n");
409			printf("a = %lf R= %lf\n", a, R);
410			}
411			else
412			{
413			if (a == b)
414			{
415			TheSum = 0.;
416			for (int nu=0; nu<=2*(n-1); nu++)
417			{
418			K2 = 0.;
419			for (int p=0; p<=(m+nu)/2; p++)
420			K2 += RosenD(2m-1, nu+1, 2p)/(2p + 1.) RosenA(2m+nu-1-2p, x);
421			TheSum += (2n-nu-1)/fact[nu]pow(aR, nu) K2;
422			}
423			sSTOCoulIntGrad_ = -pow(a, 2m+2)pow(R, 2m) /(nfact[2m])TheSum;
424			TheSum = 0.;
425			for (int nu=0; nu<=2*n-1; nu++)
426			{
427			K2 = 0.;
428			for (int p=0; p<=(m+nu-1)/2; p++)
429			K2 += RosenD(2m-1, nu, 2p)/(2p + 1.) RosenA(2m+nu-2p, x);
430			TheSum += (2n-nu)/fact[nu]pow(aR,nu) K2;
431			}
432			sSTOCoulIntGrad_ += 2pow(a, 2m+2)pow(R, 2m) /(nfact[2m])*TheSum;
433			}
434			else
435			{
436			// Slater exponents are different
437			// First calculate some useful arguments
438			y = R*(a+b);
439			z = R*(a-b);
440			TheSum = 0.;
441			for (int nu=0; nu<=2*n-1; nu++)
442			{
443			K2 = 0.;
444			for (int p=0; p<=2*m+nu; p++)
445			K2 += RosenD(2*m-1, nu+1, p)
446			* RosenB(p,z)RosenA(2m+nu-p, y);
447			TheSum += (2n-nu-1)/fact[nu]pow(bR,nu) K2;
448			}
449			sSTOCoulIntGrad_ = -bpow(a,2m+1)pow(R,2m)/
450			(2nfact[2m])TheSum;
451			TheSum = 0.;
452			for (int nu=0; nu<=2*n; nu++)
453			{
454			K2 = 0.;
455			for (int p=0; p<=2*m-1+nu; p++)
456			K2 += RosenD(2*m-1, nu, p)
457			* ((a-b)RosenB(p+1,z)RosenA(2*m+nu-p-1, y)
458			+(a+b)RosenB(p ,z)RosenA(2*m+nu-p , y));
459			TheSum += (2n-nu)/fact[nu]pow(bR,nu) K2;
460			}
461			sSTOCoulIntGrad_ += pow(a,2m+1)pow(R,2m)/(2nfact[2m])*TheSum;
462			}
463			// Now add one-electron terms and common term
464			sSTOCoulIntGrad_ = sSTOCoulIntGrad_ - (2.*m+1.)/sqr(R)
465			+ 2.m/R sSTOCoulInt(a,b,m,n,R)
466			+ pow(x,2m)/(fact[2m]sqr(R)) ((2.m+1.)exp(-x)
467			+ 2.m(1.+2.m-x)RosenA(2*m-1,x));
468			}
469			return sSTOCoulIntGrad_;
470			}
471
472			/**
473			* @brief Computes gradient of overlap integral with respect to the interatomic diatance
474			* Computes the derivative of the overlap integral over two Slater-type orbitals of s symmetry.
475			* @param a: Slater zeta exponent of first atom in inverse Bohr (au)
476			* @param b: Slater zeta exponent of second atom in inverse Bohr (au)
477			* @param m: principal quantum number of first atom
478			* @param n: principal quantum number of second atom
479			* @param R: internuclear distance in atomic units (bohr)
480			* @return the derivative of the sSTOOvInt integral
481			* @note Derived in QTPIE research notes, May 15 2007
482			*/
483			inline RealType sSTOOvIntGrad(RealType a, RealType b, int m, int n, RealType R)
484			{
485			RealType w, x, y, z, TheSum;
486
487			// Calculate first term
488			RealType sSTOOvIntGrad_ = (m+n+1.)/R * sSTOOvInt(a, b, m, n, R);
489
490			// Calculate remaining terms; answers depend on exponents
491			TheSum = 0.;
492			x = a * R;
493			if (a == b)
494			{
495			for (int q=0; q<=(m+n)/2; q++)
496			TheSum += RosenD(m,n,2q) / (2q + 1.) * RosenA(m+n-2*q+1, x);
497			sSTOOvIntGrad_ -= apow(x,m+n+1)/ sqrt(fact[2m]fact[2n])*TheSum;
498			}
499			else
500			{
501			w = b*R;
502			y = 0.5R(a+b);
503			z = 0.5R(a-b);
504			for (int q=0; q<m+n; q++)
505			TheSum = TheSum + RosenD(m,n,q) *
506			((a-b)RosenB(q+1,z)RosenA(m+n-q ,y)
507			+(a+b)RosenB(q ,z)RosenA(m+n-q+1,y));
508			sSTOOvIntGrad_ -= 0.25sqrt((pow(x, 2m+1)pow(w, 2n+1))/(fact[2m]fact[2n]))TheSum;
509			}
510			return sSTOOvIntGrad_;
511			}
512
513			/**
514			* @brief Calculates a Slater-type orbital exponent based on the hardness parameters
515			* @param Hardness: chemical hardness in atomic units
516			* @param n: principal quantum number
517			* @note Modified for use with OpenMD by Gezelter and Michalka.
518			*/
519			inline RealType getSTOZeta(int n, RealType hardness)
520			{
521			// Approximate the exact value of the constant of proportionality
522			// by its value at a very small distance epsilon
523			// since the exact R = 0 case has not be programmed
524			RealType epsilon = 1.0e-8;
525
526			// Assign orbital exponent
527			return pow(sSTOCoulInt(1., 1., n, n, epsilon) / hardness, -1./(3. + 2.*n));
528			}
529
530			#endif