Implement removal of center-of-mass terms, not sure if this is 100% c…

…orrect

src/basis.cpp

@@ -974,6 +974,43 @@ arma::mat GaussianShell::kinetic(const GaussianShell & rhs) const {
  return T;
 }
 
+// Calculate gradient matrix element between basis functions
+std::vector<arma::mat> GaussianShell::gradient_integral(const GaussianShell & rhs) const {
+ // Gradient matrix
+ std::vector<arma::mat> T(3);
+ for(size_t i=0;i<3;i++)
+ T[i].zeros(cart.size(),rhs.cart.size());
+
+ // Coordinates
+ double xa=cen.x;
+ double ya=cen.y;
+ double za=cen.z;
+
+ double xb=rhs.cen.x;
+ double yb=rhs.cen.y;
+ double zb=rhs.cen.z;
+
+ for(size_t ixl=0;ixl<c.size();ixl++)
+ for(size_t ixr=0;ixr<rhs.c.size();ixr++) {
+ auto itg = gradient_int_os(xa,ya,za,c[ixl].z,cart,xb,yb,zb,rhs.c[ixr].z,rhs.cart);
+ for(size_t ic=0;ic<3;ic++)
+ T[ic]+=c[ixl].c*rhs.c[ixr].c*itg[ic];
+ }
+
+ // Transformation to spherical harmonics. Left side:
+ if(uselm) {
+ for(size_t ic=0;ic<3;ic++)
+ T[ic]=transmat*T[ic];
+ }
+ // Right side
+ if(rhs.uselm) {
+ for(size_t ic=0;ic<3;ic++)
+ T[ic]=T[ic]*arma::trans(rhs.transmat);
+ }
+
+ return T;
+}
+
 
 // Calculate nuclear attraction matrix element between basis functions
 arma::mat GaussianShell::nuclear(double cx, double cy, double cz, const GaussianShell & rhs) const {
@@ -2316,6 +2353,39 @@ arma::mat BasisSet::kinetic() const {
  return T;
 }
 
+std::vector<arma::mat> BasisSet::gradient_integral() const {
+ // Form gradient_integrals energy matrix
+
+ // Size of basis set
+ size_t N=get_Nbf();
+
+ // Initialize matrix
+ std::vector<arma::mat> T(3);
+ for(size_t ic=0; ic<3;ic++)
+ T[ic].zeros(N,N);
+
+ // Loop over shells
+#ifdef _OPENMP
+#pragma omp parallel for schedule(dynamic)
+#endif
+ for(size_t ip=0;ip<shellpairs.size();ip++) {
+ // Shells in pair
+ size_t i=shellpairs[ip].is;
+ size_t j=shellpairs[ip].js;
+
+ // Get partial gradient_integral energy matrix
+ std::vector<arma::mat> tmp=shells[i].gradient_integral(shells[j]);
+
+ // Store result
+ for(size_t ic=0;ic<3;ic++) {
+ T[ic].submat(shells[i].get_first_ind(),shells[j].get_first_ind(),shells[i].get_last_ind(),shells[j].get_last_ind())=tmp[ic];
+ T[ic].submat(shells[j].get_first_ind(),shells[i].get_first_ind(),shells[j].get_last_ind(),shells[i].get_last_ind())=arma::trans(tmp[ic]);
+ }
+ }
+
+ return T;
+}
+
 arma::mat BasisSet::nuclear() const {
  std::vector<std::tuple<int,double,double,double>> nuclear_data;
  for(size_t inuc=0;inuc<nuclei.size();inuc++) {

src/basis.h

@@ -445,6 +445,8 @@ class BasisSet {
  arma::mat coulomb_overlap(const BasisSet & rhs) const;
  /// Calculate kinetic energy matrix
  arma::mat kinetic() const;
+ /// Calculate gradient integral matrix
+ std::vector<arma::mat> gradient_integral() const;
  /// Calculate nuclear repulsion matrix
  arma::mat nuclear() const;
  /// Calculate nuclear repulsion matrix with given nuclei
@@ -644,6 +646,8 @@ class GaussianShell {
  arma::mat coulomb_overlap(const GaussianShell & rhs) const;
  /// Calculate kinetic energy matrix between shells
  arma::mat kinetic(const GaussianShell & rhs) const;
+ /// Calculate gradient integral between shells
+ std::vector<arma::mat> gradient_integral(const GaussianShell & rhs) const;
  /// Calculate nuclear repulsion matrix between shells
  arma::mat nuclear(double cx, double cy, double cz, const GaussianShell & rhs) const;
 

src/contrib/hneoci.cpp

@@ -114,6 +114,7 @@ int main_guarded(int argc, char **argv) {
  settings.add_int("Verbosity", "Verboseness level", 5);
  settings.add_bool("H2", "Run H2+ instead of H atom?", false);
  settings.add_double("H2BondLength", "Bond length for H2+ in a.u.", 2.0);
+ settings.add_bool("RemoveCOM", "Remove COM terms", false);
 
  // Parse settings
  settings.parse(std::string(argv[1]),true);
@@ -125,6 +126,10 @@ int main_guarded(int argc, char **argv) {
  double proton_mass = settings.get_double("ProtonMass");
  bool dimer = settings.get_bool("H2");
  double R = settings.get_double("H2BondLength");
+ bool removecom = settings.get_bool("RemoveCOM");
+
+ // Total mass of the system
+ double MT = 1+proton_mass;
 
  // Read in basis sets
  BasisSetLibrary baslib;
@@ -186,9 +191,19 @@ int main_guarded(int argc, char **argv) {
  arma::mat Se(basis.overlap());
  arma::mat Sp(pbasis.overlap());
 
+ // Get gradient integrals
+ std::vector<arma::mat> grade(basis.gradient_integral());
+ std::vector<arma::mat> gradp(pbasis.gradient_integral());
+
  // Kinetic energy matrices
- arma::mat T = basis.kinetic();
- arma::mat Tp = pbasis.kinetic()/proton_mass;
+ arma::mat T, Tp;
+ if(removecom) {
+ T = (1.0-1.0/MT) * basis.kinetic();
+ Tp = (1.0/proton_mass - 1.0/MT)* pbasis.kinetic();
+ } else {
+ T = basis.kinetic();
+ Tp = pbasis.kinetic()/proton_mass;
+ }
 
  // Nuclear attraction
  arma::mat V(T.n_rows, T.n_cols, arma::fill::zeros);
@@ -215,7 +230,6 @@ int main_guarded(int argc, char **argv) {
  arma::vec E_bo;
  arma::mat C_bo;
  arma::eig_sym(E_bo, C_bo, H_bo);
-
  arma::mat Xe_bo = Xe * C_bo;
  E_bo.print("Electronic spectrum without quantum proton");
 
@@ -225,7 +239,7 @@ int main_guarded(int argc, char **argv) {
  arma::mat Cp_bo;
  arma::eig_sym(Ep_bo, Cp_bo, Hp_bo);
  arma::mat Xp_bo = Xp * Cp_bo;
- Ep_bo.print("Nucleon spectrum");
+ Ep_bo.print("Quantum proton spectrum");
 
  // Compute the two-electron integrals
  size_t e_nbf = basis.get_Nbf();
@@ -295,7 +309,6 @@ int main_guarded(int argc, char **argv) {
  size_t I = I_start+II;
  size_t J = J_start+JJ;
 
- //double element = -eris[((ii*N_j+jj)*N_I+II)*N_J+JJ];
  double element = -eris[((ii*N_j+jj)*N_I+II)*N_J+JJ];
  // (ij|IJ)
  V_ao(BFIDX(i,I),BFIDX(j,J)) = element;
@@ -311,6 +324,20 @@ int main_guarded(int argc, char **argv) {
  printf("Finished e-p integrals\n");
  fflush(stdout);
 
+ if(removecom) {
+ for(size_t i=0; i<e_nbf; i++)
+ for(size_t j=0; j<e_nbf; j++)
+ for(size_t I=0; I<p_nbf; I++)
+ for(size_t J=0; J<p_nbf; J++) {
+ double el = 0.0;
+ for(size_t ic=0; ic<3; ic++)
+ el += grade[ic](i,j) * gradp[ic](I,J);
+ V_ao(BFIDX(j,J),BFIDX(i,I)) += 1.0/MT * el;
+ }
+ printf("Finished nabla.nabla terms\n");
+ fflush(stdout);
+ }
+
  // Compute transformation matrix to go to normalized AO basis
  arma::mat ao_to_orth(e_nbf*p_nbf,e_nmo*p_nmo,arma::fill::zeros);
 #pragma omp parallel for collapse(4)

src/integrals.cpp

@@ -176,6 +176,24 @@ double kinetic_int(double xa, double ya, double za, double zetaa, int la, int ma
 }
 
 
+std::tuple<double,double,double> gradient_int(double xa, double ya, double za, double zetaa, int la, int ma, int na, double xb, double yb, double zb, double zetab, int lb, int mb, int nb) {
+
+ // Acting on exponent gives -2 x alpha
+ double xint=2.0*zetab*overlap_int(xa,ya,za,zetaa,la,ma,na,xb,yb,zb,zetab,lb+1,mb,nb);
+ double yint=2.0*zetab*overlap_int(xa,ya,za,zetaa,la,ma,na,xb,yb,zb,zetab,lb,mb+1,nb);
+ double zint=2.0*zetab*overlap_int(xa,ya,za,zetaa,la,ma,na,xb,yb,zb,zetab,lb,mb,nb+1);
+ // Acting on polynomial gives a factor
+ if(lb>=1)
+ xint-=lb*overlap_int(xa,ya,za,zetaa,la,ma,na,xb,yb,zb,zetab,lb-1,mb,nb);
+ if(mb>=1)
+ yint-=mb*overlap_int(xa,ya,za,zetaa,la,ma,na,xb,yb,zb,zetab,lb,mb-1,nb);
+ if(nb>=1)
+ zint-=nb*overlap_int(xa,ya,za,zetaa,la,ma,na,xb,yb,zb,zetab,lb,mb,nb-1);
+
+ return std::make_tuple(xint,yint,zint);
+}
+
+
 // Compute A array for nucler attraction integral (THO 2.18)
 std::vector<double> A_array(int la, int lb, double PAx, double PBx, double PCx, double zeta) {
  // The idea behind the array is that although in principle it has three

src/integrals.h

@@ -17,6 +17,7 @@
 
 
 #include "global.h"
+#include <tuple>
 
 #ifndef ERKALE_INTEGRALS
 #define ERKALE_INTEGRALS
@@ -40,6 +41,8 @@ double overlap_int(double xa, double ya, double za, double zetaa, int la, int ma
 
 /// Calculate kinetic energy integral
 double kinetic_int(double xa, double ya, double za, double zetaa, int la, int ma, int na, double xb, double yb, double zb, double zetab, int lb, int mb, int nb);
+/// Calculate gradient integral <u|nabla|v>
+std::tuple<double,double,double> gradient_int(double xa, double ya, double za, double zetaa, int la, int ma, int na, double xb, double yb, double zb, double zetab, int lb, int mb, int nb);
 
 /// Calculate nuclear attraction integral
 double nuclear_int(double xa, double ya, double za, double zetaa, int la, int ma, int na, double xnuc, double ynuc, double znuc, double xb, double yb, double zb, double zetab, int lb, int mb, int nb);

src/obara-saika.cpp

@@ -329,6 +329,113 @@ arma::mat kinetic_int_os(double xa, double ya, double za, double zetaa, const st
  return T;
 }
 
+std::vector<arma::mat> gradient_int_os(double xa, double ya, double za, double zetaa, const std::vector<shellf_t> & carta, double xb, double yb, double zb, double zetab, const std::vector<shellf_t> & cartb) {
+ // Compute shell of kinetic energy integrals
+
+ // Angular momenta of shells
+ int am_a=carta[0].l+carta[0].m+carta[0].n;
+ int am_b=cartb[0].l+cartb[0].m+cartb[0].n;
+
+ // Returned matrix
+ std::vector<arma::mat> T(3);
+ for(size_t i=0;i<3;i++)
+ T[i].zeros(carta.size(),cartb.size());
+
+ // Get 1d overlap integrals
+ arma::mat ox_arr=overlap_ints_1d(xa,xb,zetaa,zetab,am_a,am_b);
+ arma::mat oy_arr=overlap_ints_1d(ya,yb,zetaa,zetab,am_a,am_b);
+ arma::mat oz_arr=overlap_ints_1d(za,zb,zetaa,zetab,am_a,am_b);
+
+ // Get kinetic energy integrals
+ arma::mat kx_arr=derivative_ints_1d(xa,xb,zetaa,zetab,am_a,am_b,1);
+ arma::mat ky_arr=derivative_ints_1d(ya,yb,zetaa,zetab,am_a,am_b,1);
+ arma::mat kz_arr=derivative_ints_1d(za,zb,zetaa,zetab,am_a,am_b,1);
+
+ double ox, oy, oz;
+ double kx, ky, kz;
+
+ int la, ma, na;
+ int lb, mb, nb;
+
+ double anorm, bnorm;
+
+ for(size_t i=0;i<carta.size();i++) {
+ anorm=carta[i].relnorm;
+
+ la=carta[i].l;
+ ma=carta[i].m;
+ na=carta[i].n;
+
+ for(size_t j=0;j<cartb.size();j++) {
+ lb=cartb[j].l;
+ mb=cartb[j].m;
+ nb=cartb[j].n;
+
+ bnorm=cartb[j].relnorm;
+
+ ox=ox_arr(la,lb);
+ oy=oy_arr(ma,mb);
+ oz=oz_arr(na,nb);
+
+ kx=kx_arr(la,lb);
+ ky=ky_arr(ma,mb);
+ kz=kz_arr(na,nb);
+
+ T[0](i,j)=anorm*bnorm*kx*oy*oz;
+ T[1](i,j)=anorm*bnorm*ox*ky*oz;
+ T[2](i,j)=anorm*bnorm*ox*oy*kz;
+ }
+ }
+
+#ifdef DEBUG
+
+ std::vector<arma::mat> huz(3);
+ for(size_t i=0;i<3;i++)
+ huz[i].zeros(carta.size(),cartb.size());
+
+ for(size_t i=0;i<carta.size();i++) {
+ la=carta[i].l;
+ ma=carta[i].m;
+ na=carta[i].n;
+ anorm=carta[i].relnorm;
+
+ for(size_t j=0;j<cartb.size();j++) {
+ lb=cartb[j].l;
+ mb=cartb[j].m;
+ nb=cartb[j].n;
+ bnorm=cartb[j].relnorm;
+
+ auto res = gradient_int(xa,ya,za,zetaa,la,ma,na,xb,yb,zb,zetab,lb,mb,nb);
+ huz[0](i,j)=anorm*bnorm*std::get<0>(res);
+ huz[1](i,j)=anorm*bnorm*std::get<1>(res);
+ huz[2](i,j)=anorm*bnorm*std::get<2>(res);
+ }
+ }
+
+ int diff=0;
+ for(size_t ic=0;ic<3;ic++)
+ for(size_t i=0;i<carta.size();i++)
+ for(size_t j=0;j<cartb.size();j++)
+ if(fabs(T[ic](i,j)-huz[ic](i,j))>10*DBL_EPSILON*fabs(huz[ic](i,j)))
+ diff++;
+
+ if(diff==0)
+ //printf("Computed shell of KE (%e,%e,%e)-(%e,%e,%e) with zeta=(%e,%e) and am=(%i,%i), the results match.\n",xa,ya,za,xb,yb,zb,zetaa,zetab,am_a,am_b);
+ ;
+ else
+ for(size_t ic=0;ic<3;ic++)
+ for(size_t i=0;i<carta.size();i++)
+ for(size_t j=0;j<cartb.size();j++)
+ if(fabs(T[ic](i,j)-huz[ic](i,j))>1000*DBL_EPSILON*fabs(huz[ic](i,j))) {
+ printf("Computed gradient ic=%i (%e,%e,%e)-(%e,%e,%e) with zeta=(%e,%e) and am=(%i,%i,%i)-(%i,%i,%i)\n",ic,xa,ya,za,xb,yb,zb,zetaa,zetab,la,ma,na,lb,mb,nb);
+ printf("Huzinaga gives %e, OS gives %e.\n",huz[ic](i,j),T[ic](i,j));
+ }
+
+#endif
+
+ return T;
+}
+
 std::vector<arma::mat> kinetic_int_pulay_os(double xa, double ya, double za, double zetaa, const std::vector<shellf_t> & carta, double xb, double yb, double zb, double zetab, const std::vector<shellf_t> & cartb) {
  // Compute shell of kinetic energy integrals
 

src/obara-saika.h

@@ -50,6 +50,8 @@ double kinetic_int_1d(double xa, double xb, double zetaa, double zetab, int la,
 /// Compute shell of kinetic energy integral derivatives (for Pulay force). Order same as in overlap_int_pulay_os
 std::vector<arma::mat> kinetic_int_pulay_os(double xa, double ya, double za, double zetaa, const std::vector<shellf_t> & carta, double xb, double yb, double zb, double zetab, const std::vector<shellf_t> & cartb);
 
+/// Compute shell of gradient integrals
+std::vector<arma::mat> gradient_int_os(double xa, double ya, double za, double zetaa, const std::vector<shellf_t> & carta, double xb, double yb, double zb, double zetab, const std::vector<shellf_t> & cartb);
 
 /// Compute matrix element of derivative operator
 double derivative_int_1d(double xa, double xb, double zetaa, double zetab, int la, int lb, int eval);