Skip to content

Commit

Permalink
Further pruning Boost and GLM
Browse files Browse the repository at this point in the history
  • Loading branch information
@ifilot
ifilot committed Dec 5, 2023
1 parent 0c35f19 commit e688afb
Show file tree
Hide file tree
Showing 7 changed files with 396 additions and 246 deletions.
46 changes: 22 additions & 24 deletions pyqint/cgf.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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();
Expand All @@ -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) *
Expand All @@ -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;

Expand All @@ -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);
}
Expand All @@ -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>();

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<double>(2 * l - 1) )*
(m < 1 ? 1 : boost::math::double_factorial<double>(2 * m - 1) )*
(n < 1 ? 1 : boost::math::double_factorial<double>(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);
}
Expand All @@ -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
}

Expand All @@ -150,15 +148,15 @@ CGF::CGF():
* @return CGF
*/
CGF::CGF(double x, double y, double z) :
r(vec3(x,y,z)) {}
r(Vec3(x,y,z)) {}

/*
* @fn CGF
* @brief Default constructor
*
* @return CGF
*/
CGF::CGF(const vec3& _r):
CGF::CGF(const Vec3& _r):
r(_r) {
// do nothing
}
Expand All @@ -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) {
Expand All @@ -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<double> CGF::get_grad(const vec3& r) const {
vec3 sum = vec3(0,0,0);
std::vector<double> 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);
Expand All @@ -206,7 +204,7 @@ std::vector<double> 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
*/
Expand Down Expand Up @@ -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; i<this->gtos.size(); i++) {
Expand Down
52 changes: 25 additions & 27 deletions pyqint/cgf.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,10 +22,8 @@
#pragma once

#include <iostream>
#include <Eigen/Dense>
#include <boost/math/special_functions/factorials.hpp>

typedef Eigen::Vector3d vec3;
#include "factorials.h"
#include "vec3.h"

/*
* Gaussian Type Orbital
Expand All @@ -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:
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -161,51 +159,51 @@ 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;
}

/*
* @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
* @brief Set (new) position of GTO
*
* @return void
*/
inline void set_position(const vec3& _position) {
inline void set_position(const Vec3& _position) {
this->position = _position;
}

Expand All @@ -223,7 +221,7 @@ class GTO { // Gaussian Type Orbital
class CGF { // Contracted Gaussian Function
private:
std::vector<GTO> gtos; // vector holding all gtos
vec3 r; // position of the CGF
Vec3 r; // position of the CGF

public:
/*
Expand All @@ -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{
Expand All @@ -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;
}

Expand Down Expand Up @@ -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<double> get_grad(const vec3& r) const;
std::vector<double> 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<double> 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));
}

/*
Expand Down Expand Up @@ -405,5 +403,5 @@ class CGF { // Contracted Gaussian Function
*
* @return void
*/
void set_position(const vec3 &pos);
void set_position(const Vec3 &pos);
};
28 changes: 28 additions & 0 deletions pyqint/factorials.h
View file
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit e688afb

Please sign in to comment.