diff --git a/src/basis.cpp b/src/basis.cpp index 9a0f65c8..01fb12df 100644 --- a/src/basis.cpp +++ b/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 GaussianShell::gradient_integral(const GaussianShell & rhs) const { + // Gradient matrix + std::vector 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 BasisSet::gradient_integral() const { + // Form gradient_integrals energy matrix + + // Size of basis set + size_t N=get_Nbf(); + + // Initialize matrix + std::vector 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 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> nuclear_data; for(size_t inuc=0;inuc 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 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; diff --git a/src/contrib/hneoci.cpp b/src/contrib/hneoci.cpp index 3b51eff9..7c2aaa2f 100644 --- a/src/contrib/hneoci.cpp +++ b/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 grade(basis.gradient_integral()); + std::vector 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 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 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 diff --git a/src/integrals.h b/src/integrals.h index ac7856fd..275a5376 100644 --- a/src/integrals.h +++ b/src/integrals.h @@ -17,6 +17,7 @@ #include "global.h" +#include #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 +std::tuple 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); diff --git a/src/obara-saika.cpp b/src/obara-saika.cpp index c24fbb70..a397c1c8 100644 --- a/src/obara-saika.cpp +++ b/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 gradient_int_os(double xa, double ya, double za, double zetaa, const std::vector & carta, double xb, double yb, double zb, double zetab, const std::vector & 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 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 huz(3); + for(size_t i=0;i<3;i++) + huz[i].zeros(carta.size(),cartb.size()); + + for(size_t i=0;i(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;i10*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;i1000*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 kinetic_int_pulay_os(double xa, double ya, double za, double zetaa, const std::vector & carta, double xb, double yb, double zb, double zetab, const std::vector & cartb) { // Compute shell of kinetic energy integrals diff --git a/src/obara-saika.h b/src/obara-saika.h index bd7dd66b..d854bef8 100644 --- a/src/obara-saika.h +++ b/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 kinetic_int_pulay_os(double xa, double ya, double za, double zetaa, const std::vector & carta, double xb, double yb, double zb, double zetab, const std::vector & cartb); +/// Compute shell of gradient integrals +std::vector gradient_int_os(double xa, double ya, double za, double zetaa, const std::vector & carta, double xb, double yb, double zb, double zetab, const std::vector & 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);