From e688afbf24a21d3edac8e1f3c2fd6f48f87212f3 Mon Sep 17 00:00:00 2001 From: Ivo Filot Date: Tue, 5 Dec 2023 21:18:49 +0100 Subject: [PATCH] Further pruning Boost and GLM --- pyqint/cgf.cpp | 46 ++++---- pyqint/cgf.h | 52 +++++----- pyqint/factorials.h | 28 +++++ pyqint/integrals.cpp | 242 +++++++++++++++++++------------------------ pyqint/integrals.h | 115 ++++++++++---------- pyqint/vec3.h | 157 ++++++++++++++++++++++++++++ tests/test_cgf.py | 2 +- 7 files changed, 396 insertions(+), 246 deletions(-) create mode 100644 pyqint/factorials.h create mode 100644 pyqint/vec3.h diff --git a/pyqint/cgf.cpp b/pyqint/cgf.cpp index fd41564..a95fa81 100644 --- a/pyqint/cgf.cpp +++ b/pyqint/cgf.cpp @@ -22,7 +22,7 @@ #include "cgf.h" GTO::GTO(double _c, - const vec3& _position, // position (unit = Bohr) + const Vec3& _position, // position (unit = Bohr) double _alpha, unsigned int _l, unsigned int _m, @@ -51,7 +51,7 @@ GTO::GTO(double _c, l(_l), m(_m), n(_n), - position(vec3(_x,_y,_z)) { + position(Vec3(_x,_y,_z)) { // calculate the normalization constant this->calculate_normalization_constant(); @@ -61,12 +61,12 @@ GTO::GTO(double _c, * @fn get_amp * @brief Gets the amplitude of the GTO * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return const double amplitude */ -const double GTO::get_amp(const vec3& r) const { - double r2 = (r - this->position).squaredNorm(); +const double GTO::get_amp(const Vec3& r) const { + double r2 = (r - this->position).norm2(); return this->norm * std::pow(r[0]-this->position[0], l) * @@ -79,11 +79,11 @@ const double GTO::get_amp(const vec3& r) const { * @fn get_gradient * @brief Gets the gradient of the GTO * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return gradient */ -vec3 GTO::get_grad(const vec3& r) const { +Vec3 GTO::get_grad(const Vec3& r) const { const double ex = std::exp(-this->alpha * std::pow(r[0]-this->position[0],2)); const double fx = std::pow(r[0] - this->position[0], this->l) * ex; @@ -107,7 +107,7 @@ vec3 GTO::get_grad(const vec3& r) const { gz += std::pow(r[2] - this->position[2], this->n-1) * ez; } - return vec3(this->norm * gx * fy * fz, + return Vec3(this->norm * gx * fy * fz, this->norm * fx * gy * fz, this->norm * fx * fy * gz); } @@ -119,15 +119,13 @@ vec3 GTO::get_grad(const vec3& r) const { * @return void */ void GTO::calculate_normalization_constant() { - static const double pi = boost::math::constants::pi(); - double nom = std::pow(2.0, 2.0 * (l + m + n) + 3.0 / 2.0) * std::pow(alpha, (l + m + n) + 3.0 / 2.0); - double denom = (l < 1 ? 1 : boost::math::double_factorial(2 * l - 1) )* - (m < 1 ? 1 : boost::math::double_factorial(2 * m - 1) )* - (n < 1 ? 1 : boost::math::double_factorial(2 * n - 1) )* - std::pow(pi, 3.0 / 2.0); + double denom = (l < 1 ? 1 : double_factorial(2 * l - 1) )* + (m < 1 ? 1 : double_factorial(2 * m - 1) )* + (n < 1 ? 1 : double_factorial(2 * n - 1) )* + std::pow(M_PI, 3.0 / 2.0); this->norm = std::sqrt(nom / denom); } @@ -139,7 +137,7 @@ void GTO::calculate_normalization_constant() { * @return CGF */ CGF::CGF(): - r(vec3(0,0,0)) { + r(Vec3(0,0,0)) { // do nothing } @@ -150,7 +148,7 @@ CGF::CGF(): * @return CGF */ CGF::CGF(double x, double y, double z) : - r(vec3(x,y,z)) {} + r(Vec3(x,y,z)) {} /* * @fn CGF @@ -158,7 +156,7 @@ CGF::CGF(double x, double y, double z) : * * @return CGF */ -CGF::CGF(const vec3& _r): +CGF::CGF(const Vec3& _r): r(_r) { // do nothing } @@ -167,11 +165,11 @@ CGF::CGF(const vec3& _r): * @fn get_amp * @brief Gets the amplitude of the CGF * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return const double amplitude */ -const double CGF::get_amp(const vec3& r) const { +const double CGF::get_amp(const Vec3& r) const { double sum = 0.0; for(const auto& gto : this->gtos) { @@ -185,12 +183,12 @@ const double CGF::get_amp(const vec3& r) const { * @fn get_amp * @brief Gets the gradient of the CGF * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return gradient */ -std::vector CGF::get_grad(const vec3& r) const { - vec3 sum = vec3(0,0,0); +std::vector CGF::get_grad(const Vec3& r) const { + Vec3 sum = Vec3(0,0,0); for(const auto& gto : this->gtos) { sum += gto.get_coefficient() * gto.get_grad(r); @@ -206,7 +204,7 @@ std::vector CGF::get_grad(const vec3& r) const { * @param unsigned int type type of the orbital (see above for the list) * @param double alpha alpha value * @param double c coefficient - * @param const vec3& vec3 position + * @param const Vec3& Vec3 position * * @return void */ @@ -286,7 +284,7 @@ void CGF::add_gto(double c, * * @return void */ -void CGF::set_position(const vec3 &pos) { +void CGF::set_position(const Vec3 &pos) { this->r = pos; for(unsigned int i=0; igtos.size(); i++) { diff --git a/pyqint/cgf.h b/pyqint/cgf.h index c16dfd7..2a725fb 100644 --- a/pyqint/cgf.h +++ b/pyqint/cgf.h @@ -22,10 +22,8 @@ #pragma once #include -#include -#include - -typedef Eigen::Vector3d vec3; +#include "factorials.h" +#include "vec3.h" /* * Gaussian Type Orbital @@ -42,7 +40,7 @@ class GTO { // Gaussian Type Orbital double c; // coefficient double alpha; // alpha value in the exponent unsigned int l,m,n; // powers of the polynomial - vec3 position; // position vector (unit = Bohr) + Vec3 position; // position vector (unit = Bohr) double norm; // normalization constant public: @@ -87,7 +85,7 @@ class GTO { // Gaussian Type Orbital * returns double value of the incomplete Gamma Function */ GTO(double _c, - const vec3& _position, + const Vec3& _position, double _alpha, unsigned int _l, unsigned int _m, @@ -161,9 +159,9 @@ class GTO { // Gaussian Type Orbital * @fn get_position * @brief Get the center of the Gaussian * - * @return const double vec3 + * @return const double Vec3 */ - inline const vec3& get_position() const { + inline const Vec3& get_position() const { return this->position; } @@ -171,33 +169,33 @@ class GTO { // Gaussian Type Orbital * @fn get_amp * @brief Gets the amplitude of the GTO * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return const double amplitude */ - const double get_amp(const vec3& r) const; + const double get_amp(const Vec3& r) const; /* * @fn get_amp * @brief Gets the amplitude of the GTO * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return const double amplitude */ inline double get_amp(double x, double y, double z) const { - return this->get_amp(vec3(x,y,z)); + return this->get_amp(Vec3(x,y,z)); } /* * @fn get_gradient * @brief Gets the gradient of the GTO * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return gradient */ - vec3 get_grad(const vec3& r) const; + Vec3 get_grad(const Vec3& r) const; /* * @fn set_position @@ -205,7 +203,7 @@ class GTO { // Gaussian Type Orbital * * @return void */ - inline void set_position(const vec3& _position) { + inline void set_position(const Vec3& _position) { this->position = _position; } @@ -223,7 +221,7 @@ class GTO { // Gaussian Type Orbital class CGF { // Contracted Gaussian Function private: std::vector gtos; // vector holding all gtos - vec3 r; // position of the CGF + Vec3 r; // position of the CGF public: /* @@ -248,7 +246,7 @@ class CGF { // Contracted Gaussian Function * * @return CGF */ - CGF(const vec3& _r); + CGF(const Vec3& _r); // type of GTOs to add enum{ @@ -271,7 +269,7 @@ class CGF { // Contracted Gaussian Function * * @return Position */ - inline const vec3& get_r() const { + inline const Vec3& get_r() const { return r; } @@ -325,44 +323,44 @@ class CGF { // Contracted Gaussian Function * @fn get_amp * @brief Gets the amplitude of the CGF * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return const double amplitude */ - const double get_amp(const vec3& r) const; + const double get_amp(const Vec3& r) const; /* * @fn get_amp * @brief Gets the amplitude of the GTO * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return const double amplitude */ inline double get_amp(double x, double y, double z) const { - return this->get_amp(vec3(x,y,z)); + return this->get_amp(Vec3(x,y,z)); } /* * @fn get_grad * @brief Gets the gradient of the CGF * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return gradient */ - std::vector get_grad(const vec3& r) const; + std::vector get_grad(const Vec3& r) const; /* * @fn get_grad * @brief Gets the gradient of the CGF * - * @param vec3 r coordinates + * @param Vec3 r coordinates * * @return gradient */ inline std::vector get_grad(double x, double y, double z) const { - return this->get_grad(vec3(x,y,z)); + return this->get_grad(Vec3(x,y,z)); } /* @@ -405,5 +403,5 @@ class CGF { // Contracted Gaussian Function * * @return void */ - void set_position(const vec3 &pos); + void set_position(const Vec3 &pos); }; diff --git a/pyqint/factorials.h b/pyqint/factorials.h new file mode 100644 index 0000000..2fba5b8 --- /dev/null +++ b/pyqint/factorials.h @@ -0,0 +1,28 @@ +#ifndef _FACTORIALS_H +#define _FACTORIALS_H + +static double factorial(unsigned int n) { + static const double ans[] = { + 1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,479001600,6227020800,87178291200,1307674368000 + }; + + if(n > 15) { + return n * factorial(n-1); + } else { + return ans[n]; + } +} + +static double double_factorial(unsigned int n) { + static const double ans[] = { + 1,1,2,3,8,15,48,105,384,945,3840,10395,46080,135135,645120,2027025 + }; + + if(n > 15) { + return n * double_factorial(n-2); + } else { + return ans[n]; + } +} + +#endif \ No newline at end of file diff --git a/pyqint/integrals.cpp b/pyqint/integrals.cpp index 4f136cf..e6b8780 100644 --- a/pyqint/integrals.cpp +++ b/pyqint/integrals.cpp @@ -42,22 +42,20 @@ std::vector Integrator::evaluate_cgfs(const std::vector& cgfs, size_t sz = cgfs.size(); - // Construct 2x2 matrices to hold values for the overlap, - // kinetic and two nuclear integral values, respectively. - Eigen::MatrixXd S = Eigen::MatrixXd::Zero(sz, sz); - Eigen::MatrixXd T = Eigen::MatrixXd::Zero(sz, sz); - Eigen::MatrixXd V = Eigen::MatrixXd::Zero(sz, sz); - std::vector Vn(charges.size(), Eigen::MatrixXd::Zero(sz, sz)); + std::vector S(sz*sz, 0.0); + std::vector T(sz*sz, 0.0); + std::vector V(sz*sz, 0.0); + std::vector Vn(charges.size() * sz * sz, 0.0); // calculate the integral values using the integrator class #pragma omp parallel for schedule(dynamic) for(int i=0; i<(int)sz; i++) { // have to use signed int for MSVC OpenMP here for(unsigned int j=i; joverlap(cgfs[i], cgfs[j]); - T(i,j) = T(j,i) = this->kinetic(cgfs[i], cgfs[j]); + S[i*sz+j] = S[j*sz+i] = this->overlap(cgfs[i], cgfs[j]); + T[i*sz+j] = T[j*sz+i] = this->kinetic(cgfs[i], cgfs[j]); for(unsigned int k=0; knuclear(cgfs[i], cgfs[j], vec3(px[k], py[k], pz[k]), charges[k]); + Vn[k*sz*sz + i*sz+j] = Vn[k*sz*sz + j*sz+i] = this->nuclear(cgfs[i], cgfs[j], Vec3(px[k], py[k], pz[k]), charges[k]); } } } @@ -66,7 +64,7 @@ std::vector Integrator::evaluate_cgfs(const std::vector& cgfs, for(unsigned int i=0; i Integrator::evaluate_geometric_derivatives(const std::vector // Construct 2x2 matrices to hold values for the overlap, // kinetic and two nuclear integral values, respectively. - auto S = std::vector(charges.size() * 3, Eigen::MatrixXd::Zero(sz, sz)); - auto T = std::vector(charges.size() * 3, Eigen::MatrixXd::Zero(sz, sz)); - auto V = std::vector(charges.size() * 3, Eigen::MatrixXd::Zero(sz, sz)); + std::vector S(charges.size() * 3 * sz * sz, 0.0); + std::vector T(charges.size() * 3 * sz * sz, 0.0); + std::vector V(charges.size() * 3 * sz * sz, 0.0); // calculate the integral values using the integrator class for(unsigned int n=0; n Vn(charges.size(), Eigen::MatrixXd::Zero(sz, sz)); + std::vector Vn(charges.size() * sz * sz, 0.0); #pragma omp parallel for schedule(dynamic) for(int i=0; i<(int)sz; i++) { // have to use signed int for MSVC OpenMP here for(int j=0; j<(int)sz; j++) { - S[n*3+k](i,j) = this->overlap_deriv(cgfs[i], cgfs[j], vec3(px[n], py[n], pz[n]), k); - T[n*3+k](i,j) = this->kinetic_deriv(cgfs[i], cgfs[j], vec3(px[n], py[n], pz[n]), k); + S[(n*3+k) * sz * sz + i * sz + j] = this->overlap_deriv(cgfs[i], cgfs[j], Vec3(px[n], py[n], pz[n]), k); + T[(n*3+k) * sz * sz + i * sz + j] = this->kinetic_deriv(cgfs[i], cgfs[j], Vec3(px[n], py[n], pz[n]), k); for(unsigned int l=0; lnuclear_deriv(cgfs[i], cgfs[j], - vec3(px[l], py[l], pz[l]), + Vn[l*sz*sz + i*sz + j] = this->nuclear_deriv(cgfs[i], cgfs[j], + Vec3(px[l], py[l], pz[l]), charges[l], - vec3(px[n], py[n], pz[n]), + Vec3(px[n], py[n], pz[n]), k); } } @@ -187,7 +185,7 @@ std::vector Integrator::evaluate_geometric_derivatives(const std::vector for(unsigned int i=0; i Integrator::evaluate_geometric_derivatives(const std::vector const size_t l = jobs[s][4]; const size_t n = jobs[s][5]; const size_t d = jobs[s][6]; - tedouble[idx] = this->repulsion_deriv(cgfs[i], cgfs[j], cgfs[k], cgfs[l], vec3(px[n], py[n], pz[n]), d); + tedouble[idx] = this->repulsion_deriv(cgfs[i], cgfs[j], cgfs[k], cgfs[l], Vec3(px[n], py[n], pz[n]), d); } // package everything into results vector, will be unpacked in @@ -247,13 +245,13 @@ std::vector Integrator::evaluate_geometric_derivatives(const std::vector for(size_t n=0; n @@ -792,14 +790,14 @@ double Integrator::repulsion(const CGF &cgf1,const CGF &cgf2,const CGF &cgf3,con * @return double value of the repulsion integral */ double Integrator::repulsion_deriv(const CGF &cgf1, const CGF &cgf2, const CGF &cgf3, const CGF &cgf4, - const vec3& nucleus, unsigned int coord) const { + const Vec3& nucleus, unsigned int coord) const { double sum = 0; // check if cgf originates from nucleus - bool cgf1_nuc = (cgf1.get_r() - nucleus).squaredNorm() < 0.0001; - bool cgf2_nuc = (cgf2.get_r() - nucleus).squaredNorm() < 0.0001; - bool cgf3_nuc = (cgf3.get_r() - nucleus).squaredNorm() < 0.0001; - bool cgf4_nuc = (cgf4.get_r() - nucleus).squaredNorm() < 0.0001; + bool cgf1_nuc = (cgf1.get_r() - nucleus).norm2() < 0.0001; + bool cgf2_nuc = (cgf2.get_r() - nucleus).norm2() < 0.0001; + bool cgf3_nuc = (cgf3.get_r() - nucleus).norm2() < 0.0001; + bool cgf4_nuc = (cgf4.get_r() - nucleus).norm2() < 0.0001; // early exit if(cgf1_nuc == cgf2_nuc && cgf2_nuc == cgf3_nuc && cgf3_nuc == cgf4_nuc) { @@ -904,27 +902,25 @@ double Integrator::repulsion_deriv(const GTO& gto1, const GTO& gto2, const GTO & * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param double alpha2 Gaussian exponent of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * * Calculates the value of < gto1 | gto2 > * * @return double value of the overlap integral */ -double Integrator::overlap(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const vec3 &a, - double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const vec3 &b) const { - - static const double pi = boost::math::constants::pi(); +double Integrator::overlap(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const Vec3 &a, + double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const Vec3 &b) const { - double rab2 = (a-b).squaredNorm(); + double rab2 = (a-b).norm2(); double gamma = alpha1 + alpha2; - vec3 p = this->gaussian_product_center(alpha1, a, alpha2, b); + Vec3 p = this->gaussian_product_center(alpha1, a, alpha2, b); - double pre = std::pow(pi / gamma, 1.5) * std::exp(-alpha1 * alpha2 * rab2 / gamma); + double pre = std::pow(M_PI / gamma, 1.5) * std::exp(-alpha1 * alpha2 * rab2 / gamma); double wx = this->overlap_1D(l1, l2, p[0]-a[0], p[0]-b[0], gamma); double wy = this->overlap_1D(m1, m2, p[1]-a[1], p[1]-b[1], gamma); double wz = this->overlap_1D(n1, n2, p[2]-a[2], p[2]-b[2], gamma); @@ -939,31 +935,30 @@ double Integrator::overlap(double alpha1, unsigned int l1, unsigned int m1, unsi * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param double alpha2 Gaussian exponent of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * * @return double value of the overlap integral */ -double Integrator::dipole(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const vec3 &a, - double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const vec3 &b, +double Integrator::dipole(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const Vec3 &a, + double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const Vec3 &b, unsigned int cc, double cref) const { - static const double pi = boost::math::constants::pi(); - double rab2 = (a-b).squaredNorm(); + double rab2 = (a-b).norm2(); double gamma = alpha1 + alpha2; // determine new product center - vec3 p = this->gaussian_product_center(alpha1, a, alpha2, b); + Vec3 p = this->gaussian_product_center(alpha1, a, alpha2, b); // determine correcting pre-factor - double pre = std::pow(pi / gamma, 1.5) * std::exp(-alpha1 * alpha2 * rab2 / gamma); + double pre = std::pow(M_PI / gamma, 1.5) * std::exp(-alpha1 * alpha2 * rab2 / gamma); // construct adjusted triple product - vec3 w; + Vec3 w; const std::array o1({l1,m1,n1}); const std::array o2({l2,m2,n2}); @@ -972,7 +967,7 @@ double Integrator::dipole(double alpha1, unsigned int l1, unsigned int m1, unsig } // create copy and adjust the "cc" value - vec3 wd = w; + Vec3 wd = w; wd[cc] = this->overlap_1D(o1[cc], o2[cc]+1, p[cc]-a[cc], p[cc]-b[cc], gamma); @@ -999,7 +994,7 @@ double Integrator::overlap_1D(int l1, int l2, double x1, double x2, double gamma for(int i=0; i < (1 + std::floor(0.5 * (l1 + l2))); i++) { sum += this->binomial_prefactor(2*i, l1, l2, x1, x2) * - (i == 0 ? 1 : this->double_factorial(2 * i - 1) ) / + (i == 0 ? 1 : double_factorial(2 * i - 1) ) / std::pow(2 * gamma, i); } @@ -1011,14 +1006,14 @@ double Integrator::overlap_1D(int l1, int l2, double x1, double x2, double gamma * * @param double alpha1 Gaussian exponent of the first GTO * @param double alpha2 Gaussian exponent of the second GTO - * @param const vec3 a Center of the first GTO - * @param const vec3 b Center of the second GTO + * @param const Vec3 a Center of the first GTO + * @param const Vec3 b Center of the second GTO * * * @return new gaussian product center */ -vec3 Integrator::gaussian_product_center(double alpha1, const vec3& a, - double alpha2, const vec3& b) const { +Vec3 Integrator::gaussian_product_center(double alpha1, const Vec3& a, + double alpha2, const Vec3& b) const { return (alpha1 * a + alpha2 * b) / (alpha1 + alpha2); } @@ -1042,20 +1037,19 @@ double Integrator::binomial(int a, int b) const { if( (a < 0) | (b < 0) | (a-b < 0) ) { return 1.0; } - return this->factorial(a) / (this->factorial(b) * this->factorial(a-b)); + return factorial(a) / (factorial(b) * factorial(a-b)); } -double Integrator::nuclear(const vec3& a, int l1, int m1, int n1, double alpha1, - const vec3& b, int l2, int m2, int n2, - double alpha2, const vec3& c) const { +double Integrator::nuclear(const Vec3& a, int l1, int m1, int n1, double alpha1, + const Vec3& b, int l2, int m2, int n2, + double alpha2, const Vec3& c) const { - static const double pi = boost::math::constants::pi(); double gamma = alpha1 + alpha2; - vec3 p = gaussian_product_center(alpha1, a, alpha2, b); - double rab2 = (a-b).squaredNorm(); - double rcp2 = (c-p).squaredNorm(); + Vec3 p = gaussian_product_center(alpha1, a, alpha2, b); + double rab2 = (a-b).norm2(); + double rcp2 = (c-p).norm2(); std::vector ax = A_array(l1, l2, p[0]-a[0], p[0]-b[0], p[0]-c[0], gamma); std::vector ay = A_array(m1, m2, p[1]-a[1], p[1]-b[1], p[1]-c[1], gamma); @@ -1071,7 +1065,7 @@ double Integrator::nuclear(const vec3& a, int l1, int m1, int n1, double alpha1, } } - return -2.0 * pi / gamma * std::exp(-alpha1*alpha2*rab2/gamma) * sum; + return -2.0 * M_PI / gamma * std::exp(-alpha1*alpha2*rab2/gamma) * sum; } std::vector Integrator::A_array(const int l1, const int l2, const double pa, const double pb, const double cp, const double g) const { @@ -1092,40 +1086,39 @@ std::vector Integrator::A_array(const int l1, const int l2, const double double Integrator::A_term(const int i, const int r, const int u, const int l1, const int l2, const double pax, const double pbx, const double cpx, const double gamma) const { return std::pow(-1,i) * this->binomial_prefactor(i,l1,l2,pax,pbx)* - std::pow(-1,u) * this->factorial(i)*std::pow(cpx,i-2*r-2*u)* - std::pow(0.25/gamma,r+u)/this->factorial(r)/this->factorial(u)/this->factorial(i-2*r-2*u); + std::pow(-1,u) * factorial(i)*std::pow(cpx,i-2*r-2*u)* + std::pow(0.25/gamma,r+u)/factorial(r)/factorial(u)/factorial(i-2*r-2*u); } /** * @brief Performs nuclear integral evaluation * - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO * @param double alpha1 Gaussian exponent of the first GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO * @param double alpha2 Gaussian exponent of the second GTO - * @param vec3 c + * @param Vec3 c * * @return double value of the nuclear integral */ -double Integrator::repulsion(const vec3 &a, const int la, const int ma, const int na, const double alphaa, - const vec3 &b, const int lb, const int mb, const int nb, const double alphab, - const vec3 &c, const int lc, const int mc, const int nc, const double alphac, - const vec3 &d, const int ld, const int md, const int nd, const double alphad) const { +double Integrator::repulsion(const Vec3 &a, const int la, const int ma, const int na, const double alphaa, + const Vec3 &b, const int lb, const int mb, const int nb, const double alphab, + const Vec3 &c, const int lc, const int mc, const int nc, const double alphac, + const Vec3 &d, const int ld, const int md, const int nd, const double alphad) const { - static const double pi = boost::math::constants::pi(); - double rab2 = (a-b).squaredNorm(); - double rcd2 = (c-d).squaredNorm(); + double rab2 = (a-b).norm2(); + double rcd2 = (c-d).norm2(); - vec3 p = gaussian_product_center(alphaa, a, alphab, b); - vec3 q = gaussian_product_center(alphac, c, alphad, d); + Vec3 p = gaussian_product_center(alphaa, a, alphab, b); + Vec3 q = gaussian_product_center(alphac, c, alphad, d); - double rpq2 = (p-q).squaredNorm(); + double rpq2 = (p-q).norm2(); double gamma1 = alphaa + alphab; double gamma2 = alphac + alphad; @@ -1144,7 +1137,7 @@ double Integrator::repulsion(const vec3 &a, const int la, const int ma, const in } } - return 2.0 * std::pow(pi,2.5)/(gamma1*gamma2*std::sqrt(gamma1+gamma2))* + return 2.0 * std::pow(M_PI,2.5)/(gamma1*gamma2*std::sqrt(gamma1+gamma2))* std::exp(-alphaa*alphab*rab2/gamma1)* std::exp(-alphac*alphad*rcd2/gamma2)*sum; } @@ -1154,33 +1147,32 @@ double Integrator::repulsion(const vec3 &a, const int la, const int ma, const in * * This function uses function-level caching of the Fgamma function * - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO * @param double alpha1 Gaussian exponent of the first GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO * @param double alpha2 Gaussian exponent of the second GTO - * @param vec3 c + * @param Vec3 c * * @return double value of the nuclear integral */ -double Integrator::repulsion_fgamma_cached(const vec3 &a, const int la, const int ma, const int na, const double alphaa, - const vec3 &b, const int lb, const int mb, const int nb, const double alphab, - const vec3 &c, const int lc, const int mc, const int nc, const double alphac, - const vec3 &d, const int ld, const int md, const int nd, const double alphad) const { +double Integrator::repulsion_fgamma_cached(const Vec3 &a, const int la, const int ma, const int na, const double alphaa, + const Vec3 &b, const int lb, const int mb, const int nb, const double alphab, + const Vec3 &c, const int lc, const int mc, const int nc, const double alphac, + const Vec3 &d, const int ld, const int md, const int nd, const double alphad) const { - static const double pi = boost::math::constants::pi(); - double rab2 = (a-b).squaredNorm(); - double rcd2 = (c-d).squaredNorm(); + double rab2 = (a-b).norm2(); + double rcd2 = (c-d).norm2(); - vec3 p = gaussian_product_center(alphaa, a, alphab, b); - vec3 q = gaussian_product_center(alphac, c, alphad, d); + Vec3 p = gaussian_product_center(alphaa, a, alphab, b); + Vec3 q = gaussian_product_center(alphac, c, alphad, d); - double rpq2 = (p-q).squaredNorm(); + double rpq2 = (p-q).norm2(); double gamma1 = alphaa + alphab; double gamma2 = alphac + alphad; @@ -1208,7 +1200,7 @@ double Integrator::repulsion_fgamma_cached(const vec3 &a, const int la, const in } } - return 2.0 * std::pow(pi,2.5)/(gamma1*gamma2*std::sqrt(gamma1+gamma2))* + return 2.0 * std::pow(M_PI,2.5)/(gamma1*gamma2*std::sqrt(gamma1+gamma2))* std::exp(-alphaa*alphab*rab2/gamma1)* std::exp(-alphac*alphad*rcd2/gamma2)*sum; } @@ -1255,7 +1247,7 @@ double Integrator::B0(int i, int r, double g) const { } double Integrator::fact_ratio2(unsigned int a, unsigned int b) const { - return this->factorial(a) / this->factorial(b) / this->factorial(a - 2*b); + return factorial(a) / factorial(b) / factorial(a - 2*b); } size_t Integrator::teindex(size_t i, size_t j, size_t k, size_t l) const { @@ -1276,30 +1268,6 @@ size_t Integrator::teindex(size_t i, size_t j, size_t k, size_t l) const { return ij * (ij + 1) / 2 + kl; } -double Integrator::factorial(unsigned int n) const { - static const double ans[] = { - 1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,479001600,6227020800,87178291200,1307674368000 - }; - - if(n > 15) { - return n * this->factorial(n-1); - } else { - return ans[n]; - } -} - -double Integrator::double_factorial(unsigned int n) const { - static const double ans[] = { - 1,1,2,3,8,15,48,105,384,945,3840,10395,46080,135135,645120,2027025 - }; - - if(n > 15) { - return n * this->double_factorial(n-2); - } else { - return ans[n]; - } -} - void Integrator::init() { std::unordered_map map{ {200505,"2.5"}, diff --git a/pyqint/integrals.h b/pyqint/integrals.h index 85a8cdf..af5f671 100644 --- a/pyqint/integrals.h +++ b/pyqint/integrals.h @@ -23,9 +23,14 @@ #include #include +#include +#include #include "gamma.h" #include "cgf.h" +#include "factorials.h" + +typedef std::vector> MATDOUBLE; class Integrator { private: @@ -163,7 +168,7 @@ class Integrator { */ double overlap_deriv(const CGF& cgf1, const CGF& cgf2, - const vec3& nucleus, + const Vec3& nucleus, unsigned int coord) const; /** @@ -222,7 +227,7 @@ class Integrator { // expanded interface for Cython inline double overlap_deriv(const CGF& cgf1, const CGF& cgf2, double cx, double cy, double cz, unsigned int coord) const { - return this->overlap_deriv(cgf1, cgf2, vec3(cx, cy, cz), coord); + return this->overlap_deriv(cgf1, cgf2, Vec3(cx, cy, cz), coord); } /************************************************************************** @@ -267,7 +272,7 @@ class Integrator { */ double kinetic_deriv(const CGF& cgf1, const CGF& cgf2, - const vec3& nucleus, + const Vec3& nucleus, unsigned int coord) const; @@ -289,7 +294,7 @@ class Integrator { const CGF& cgf2, double cx, double cy, double cz, unsigned int coord) const { - return this->kinetic_deriv(cgf1, cgf2, vec3(cx, cy, cz), coord); + return this->kinetic_deriv(cgf1, cgf2, Vec3(cx, cy, cz), coord); } /** @@ -325,7 +330,7 @@ class Integrator { */ double nuclear(const CGF &cgf1, const CGF &cgf2, - const vec3& nucleus, + const Vec3& nucleus, unsigned int charge) const; /** @@ -341,7 +346,7 @@ class Integrator { */ double nuclear_gto(const GTO >o1, const GTO >o2, - const vec3& nucleus) const; + const Vec3& nucleus) const; /** @@ -364,7 +369,7 @@ class Integrator { const CGF &cgf2, double cx, double cy, double cz, unsigned int charge) const { - return this->nuclear(cgf1, cgf2, vec3(cx, cy, cz), charge); + return this->nuclear(cgf1, cgf2, Vec3(cx, cy, cz), charge); } /** @@ -372,20 +377,20 @@ class Integrator { * * @param const CGF& cgf1 Contracted Gaussian Function * @param const CGF& cgf2 Contracted Gaussian Function - * @param const vec3 nucleus Position of the nucleus + * @param const Vec3 nucleus Position of the nucleus * @param unsigned int charge charge of the nucleus in a.u. * * Calculates the value of < cgf1 | V | cgf2 > * * @return double value of the nuclear integral */ - double nuclear_deriv(const CGF &cgf1, const CGF &cgf2, const vec3& nucleus, unsigned int charge, - const vec3& nucderiv, unsigned int coord) const; + double nuclear_deriv(const CGF &cgf1, const CGF &cgf2, const Vec3& nucleus, unsigned int charge, + const Vec3& nucderiv, unsigned int coord) const; // expanded notation for Cython interface inline double nuclear_deriv(const CGF &cgf1, const CGF &cgf2, double cx, double cy, double cz, unsigned int charge, double dx, double dy, double dz, unsigned int coord) const { - return this->nuclear_deriv(cgf1, cgf2, vec3(cx, cy, cz), charge, vec3(dx, dy, dz), coord); + return this->nuclear_deriv(cgf1, cgf2, Vec3(cx, cy, cz), charge, Vec3(dx, dy, dz), coord); } /** @@ -400,7 +405,7 @@ class Integrator { * @return double value of the nuclear integral */ inline double nuclear_gto(const GTO >o1, const GTO >o2, double cx, double cy, double cz) const { - return this->nuclear_gto(gto1, gto2, vec3(cx, cy, cz)); + return this->nuclear_gto(gto1, gto2, Vec3(cx, cy, cz)); } /************************************************************************** @@ -442,7 +447,7 @@ class Integrator { * @param const CGF& cgf2 Contracted Gaussian Function * @param const CGF& cgf3 Contracted Gaussian Function * @param const CGF& cgf4 Contracted Gaussian Function - * @param const vec3& nucleus Nucleus coordinates + * @param const Vec3& nucleus Nucleus coordinates * @param unsigned int coord Derivative direction * * Calculates the value of d/dcx < cgf1 | cgf2 | cgf3 | cgf4 > @@ -450,11 +455,11 @@ class Integrator { * @return double value of the repulsion integral */ double repulsion_deriv(const CGF &cgf1,const CGF &cgf2,const CGF &cgf3,const CGF &cgf4, - const vec3& nucleus, unsigned int coord) const; + const Vec3& nucleus, unsigned int coord) const; inline double repulsion_deriv(const CGF &cgf1,const CGF &cgf2,const CGF &cgf3,const CGF &cgf4, double cx, double cy, double cz, unsigned int coord) const { - return this->repulsion_deriv(cgf1, cgf2, cgf3, cgf4, vec3(cx, cy, cz), coord); + return this->repulsion_deriv(cgf1, cgf2, cgf3, cgf4, Vec3(cx, cy, cz), coord); } /** @@ -488,17 +493,17 @@ class Integrator { * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param double alpha2 Gaussian exponent of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * * @return double value of the overlap integral */ - double overlap(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const vec3 &a, - double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const vec3 &b) const; + double overlap(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const Vec3 &a, + double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const Vec3 &b) const; /** * @brief Performs overlap integral evaluation @@ -507,87 +512,87 @@ class Integrator { * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param double alpha2 Gaussian exponent of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * * @return double value of the overlap integral */ - double dipole(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const vec3 &a, - double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const vec3 &b, + double dipole(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const Vec3 &a, + double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const Vec3 &b, unsigned int cc, double cref = 0.0) const; /** * @brief Performs nuclear integral evaluation * - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO * @param double alpha1 Gaussian exponent of the first GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO * @param double alpha2 Gaussian exponent of the second GTO - * @param vec3 c + * @param Vec3 c * * @return double value of the nuclear integral */ - double nuclear(const vec3& a, + double nuclear(const Vec3& a, int l1, int m1, int n1, double alpha1, - const vec3& b, + const Vec3& b, int l2, int m2, int n2, double alpha2, - const vec3& c) const; + const Vec3& c) const; /** * @brief Performs nuclear integral evaluation * - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO * @param double alpha1 Gaussian exponent of the first GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO * @param double alpha2 Gaussian exponent of the second GTO - * @param vec3 c Nuclear position + * @param Vec3 c Nuclear position * @param coord Cartesian direction to derive nuclear coordinate towards * * @return double value of the nuclear integral derived towards nuclear coordinate */ - double nuclear_deriv_op(const vec3& a, int l1, int m1, int n1, double alpha1, - const vec3& b, int l2, int m2, int n2, - double alpha2, const vec3& c, unsigned int coord) const; + double nuclear_deriv_op(const Vec3& a, int l1, int m1, int n1, double alpha1, + const Vec3& b, int l2, int m2, int n2, + double alpha2, const Vec3& c, unsigned int coord) const; /** * @brief Performs nuclear integral evaluation * - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO * @param double alpha1 Gaussian exponent of the first GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO * @param double alpha2 Gaussian exponent of the second GTO - * @param vec3 c + * @param Vec3 c * * @return double value of the nuclear integral */ - double repulsion(const vec3 &a, const int la, const int ma, const int na, const double alphaa, - const vec3 &b, const int lb, const int mb, const int nb, const double alphab, - const vec3 &c, const int lc, const int mc, const int nc, const double alphac, - const vec3 &d, const int ld, const int md, const int nd, const double alphad) const; + double repulsion(const Vec3 &a, const int la, const int ma, const int na, const double alphaa, + const Vec3 &b, const int lb, const int mb, const int nb, const double alphab, + const Vec3 &c, const int lc, const int mc, const int nc, const double alphac, + const Vec3 &d, const int ld, const int md, const int nd, const double alphad) const; /** * @brief Performs nuclear integral evaluation including caching of Fgamma @@ -596,24 +601,24 @@ class Integrator { * was suggested in https://github.com/ifilot/hfcxx/issues/8, but explicit unit testing * actually shows not appreciable difference in speed. * - * @param vec3 a Center of the Gaussian orbital of the first GTO + * @param Vec3 a Center of the Gaussian orbital of the first GTO * @param unsigned int l1 Power of x component of the polynomial of the first GTO * @param unsigned int m1 Power of y component of the polynomial of the first GTO * @param unsigned int n1 Power of z component of the polynomial of the first GTO * @param double alpha1 Gaussian exponent of the first GTO - * @param vec3 b Center of the Gaussian orbital of the second GTO + * @param Vec3 b Center of the Gaussian orbital of the second GTO * @param unsigned int l2 Power of x component of the polynomial of the second GTO * @param unsigned int m2 Power of y component of the polynomial of the second GTO * @param unsigned int n2 Power of z component of the polynomial of the second GTO * @param double alpha2 Gaussian exponent of the second GTO - * @param vec3 c + * @param Vec3 c * * @return double value of the nuclear integral */ - double repulsion_fgamma_cached(const vec3 &a, const int la, const int ma, const int na, const double alphaa, - const vec3 &b, const int lb, const int mb, const int nb, const double alphab, - const vec3 &c, const int lc, const int mc, const int nc, const double alphac, - const vec3 &d, const int ld, const int md, const int nd, const double alphad) const; + double repulsion_fgamma_cached(const Vec3 &a, const int la, const int ma, const int na, const double alphaa, + const Vec3 &b, const int lb, const int mb, const int nb, const double alphab, + const Vec3 &c, const int lc, const int mc, const int nc, const double alphac, + const Vec3 &d, const int ld, const int md, const int nd, const double alphad) const; /** * @brief Calculates one dimensional overlap integral @@ -639,14 +644,14 @@ class Integrator { * * @param double alpha1 Gaussian exponent of the first GTO * @param double alpha2 Gaussian exponent of the second GTO - * @param const vec3 a Center of the first GTO - * @param const vec3 b Center of the second GTO + * @param const Vec3 a Center of the first GTO + * @param const Vec3 b Center of the second GTO * * * @return new gaussian product center */ - vec3 gaussian_product_center(double alpha1, const vec3 &a, - double alpha2, const vec3 &b) const; + Vec3 gaussian_product_center(double alpha1, const Vec3 &a, + double alpha2, const Vec3 &b) const; double binomial_prefactor(int s, int ia, int ib, double xpa, double xpb) const; @@ -679,9 +684,5 @@ class Integrator { double B0(int i, int r, double q) const; double fact_ratio2(unsigned int a, unsigned int b) const; - double factorial(unsigned int n) const; - - double double_factorial(unsigned int n) const; - void init(); }; diff --git a/pyqint/vec3.h b/pyqint/vec3.h new file mode 100644 index 0000000..58c0265 --- /dev/null +++ b/pyqint/vec3.h @@ -0,0 +1,157 @@ +#ifndef _VEC3_H +#define _VEC3_H + +typedef float mat33[3][3]; + +/** + * Custom 3-vector class +*/ +class Vec3 { +public: + float x = 0.0; + float y = 0.0; + float z = 0.0; + + /** + * Default constructor + */ + Vec3() {} + + /** + * Initialization constructor + */ + Vec3(float _x, float _y, float _z) : x(_x), y(_y), z(_z) {} + + float& operator[](int n) { + switch(n) { + case 0: + return this->x; + case 1: + return this->y; + case 2: + return this->z; + default: // should never reach this + return this->z; + } + } + + const float& operator[](int n) const { + switch(n) { + case 0: + return this->x; + case 1: + return this->y; + case 2: + return this->z; + default: // should never reach this + return this->z; + } + } + + /** + * Multiplication operation between matrix and vector + */ + friend Vec3 operator*(const mat33& lhs, const Vec3& rhs) { + Vec3 r; + r.x = lhs[0][0] * rhs.x + lhs[1][0] * rhs.y + lhs[2][0] * rhs.z; + r.y = lhs[0][1] * rhs.x + lhs[1][1] * rhs.y + lhs[2][1] * rhs.z; + r.z = lhs[0][2] * rhs.x + lhs[1][2] * rhs.y + lhs[2][2] * rhs.z; + + return r; + } + + /** + * Multiplication operation between scalar and vector + */ + friend Vec3 operator*(float v, const Vec3& rhs) { + return Vec3(v * rhs.x, v * rhs.y, v * rhs.z); + } + + /** + * Multiplication operation between vector and scalar + */ + friend Vec3 operator*(const Vec3& rhs, float v) { + return Vec3(v * rhs.x, v * rhs.y, v * rhs.z); + } + + /** + * Return normalized vector + */ + Vec3 normalized() const { + float l = std::sqrt(this->x * this->x + this->y * this->y + this->z * this->z); + return Vec3(this->x / l, this->y / l, this->z / l); + } + + /** + * Return dot product between two vectors + */ + float dot(const Vec3& rhs) const { + return this->x * rhs.x + this->y * rhs.y + this->z * rhs.z; + } + + /** + * Addition operation between two vectors + */ + friend Vec3 operator+(const Vec3& lhs, const Vec3& rhs) { + return Vec3(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z); + } + + /** + * Subtraction operation between two vectors + */ + friend Vec3 operator-(const Vec3& lhs, const Vec3& rhs) { + return Vec3(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z); + } + + /** + * Addition assignment operation + */ + void operator+=(const Vec3& rhs) { + this->x += rhs.x; + this->y += rhs.y; + this->z += rhs.z; + } + + /** + * Subtraction assignment operation + */ + void operator-=(const Vec3& rhs) { + this->x -= rhs.x; + this->y -= rhs.y; + this->z -= rhs.z; + } + + /** + * Divide vector by a scalar operation + */ + friend Vec3 operator/(const Vec3& rhs, float v) { + return Vec3(rhs.x / v, rhs.y / v, rhs.z / v); + } + + /** + * Calculate cross product between two vectors + */ + Vec3 cross(const Vec3& rhs) const { + return Vec3( + this->y * rhs.z - this->z * rhs.y, + this->z * rhs.x - this->x * rhs.z, + this->x * rhs.y - this->y * rhs.x + ); + } + + /** + * Calculate squared sum of coefficients + */ + double norm2() const { + return this->x * this->x + this->y * this->y + this->z * this->z; + } + + /** + * Calculate product of coefficients + */ + double prod() const { + return this->x * this->y * this->z; + } +}; + +#endif // _VEC3_H diff --git a/tests/test_cgf.py b/tests/test_cgf.py index de3b980..951e55c 100644 --- a/tests/test_cgf.py +++ b/tests/test_cgf.py @@ -66,7 +66,7 @@ def testPlotGrid(self): res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x))) ans = np.load(os.path.join(os.path.dirname(__file__), 'results', 'h2o_orb_1b2.npy')) - np.testing.assert_almost_equal(res, ans, 8) + np.testing.assert_almost_equal(res, ans, 6) if __name__ == '__main__': unittest.main()