Skip to content

Commit

Permalink
Merge pull request #1826 from ghutchis/implement-g-orbitals
Browse files Browse the repository at this point in the history 
Implement g orbitals
  • Loading branch information
@ghutchis
ghutchis authored Nov 29, 2024
2 parents 3788209 + 37b95d5 commit b7e4a12
Show file tree
Hide file tree
Showing 3 changed files with 200 additions and 24 deletions.
72 changes: 66 additions & 6 deletions avogadro/core/gaussianset.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,24 @@ unsigned int GaussianSet::addBasis(unsigned int atom, orbital type)
case F7:
m_numMOs += 7;
break;
case G:
m_numMOs += 15;
break;
case G9:
m_numMOs += 9;
break;
case H:
m_numMOs += 21;
break;
case H11:
m_numMOs += 11;
break;
case I:
m_numMOs += 28;
break;
case I13:
m_numMOs += 13;
break;
default:
// Should never hit here
;
Expand Down Expand Up @@ -434,12 +452,54 @@ void GaussianSet::initCalculation()
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.25) * norm); //-3
}
} break;
case G:
skip = 15;
break;
case G9:
skip = 9;
break;
case G: {
// 16 * (2.0/pi)^0.75
double norm = 11.403287525679843;
double norm1 = norm / sqrt(7.0);
double norm2 = norm / sqrt(35.0 / 3.0);
double norm3 = norm / sqrt(35.0);
m_moIndices[i] = indexMO;
indexMO += 15;
m_cIndices.push_back(static_cast<unsigned int>(m_gtoCN.size()));
for (unsigned j = m_gtoIndices[i]; j < m_gtoIndices[i + 1]; ++j) {
// molden order
// xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
// xxyy xxzz yyzz xxyz yyxz zzxy
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); // xxxx
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); // yyyy
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); // zzzz
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // xxxy
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // xxxz
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // yyyx
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // yyyz
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // zzzx
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // zzzy
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm2); // xxyy
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm2); // xxzz
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm2); // yyzz
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm3); // xxyz
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm3); // yyxz
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm3); // zzxy
}
} break;
case G9: {
// 16 * (2.0/pi)^0.75
double norm = 11.403287525679843;
m_moIndices[i] = indexMO;
indexMO += 9;
m_cIndices.push_back(static_cast<unsigned int>(m_gtoCN.size()));
for (unsigned j = m_gtoIndices[i]; j < m_gtoIndices[i + 1]; ++j) {
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); // 0
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //+1
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //-1
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //+2
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //-2
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //+3
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //-3
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //+4
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //-4
}
} break;
case H:
skip = 21;
break;
Expand Down
148 changes: 130 additions & 18 deletions avogadro/core/gaussiansettools.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,10 @@
#include "gaussianset.h"
#include "molecule.h"

#include <iostream>

using std::vector;

namespace Avogadro::Core {

GaussianSetTools::GaussianSetTools(Molecule* mol) : m_molecule(mol)
Expand Down Expand Up @@ -227,6 +231,12 @@ inline std::vector<double> GaussianSetTools::calculateValues(
case GaussianSet::F7:
pointF7(i, deltas[atomIndices[i]], dr2[atomIndices[i]], values);
break;
case GaussianSet::G:
pointG(i, deltas[atomIndices[i]], dr2[atomIndices[i]], values);
break;
case GaussianSet::G9:
pointG9(i, deltas[atomIndices[i]], dr2[atomIndices[i]], values);
break;
default:
// Not handled - return a zero contribution
;
Expand Down Expand Up @@ -259,13 +269,13 @@ inline void GaussianSetTools::pointP(unsigned int moIndex, const Vector3& delta,
unsigned int baseIndex = m_basis->moIndices()[moIndex];
Vector3 components(Vector3::Zero());

// Now iterate through the P type GTOs and sum their contributions
// Now iterate through the GTOs and sum their contributions
unsigned int cIndex = m_basis->cIndices()[moIndex];
double tmpGTO = 0.0;
for (unsigned int i = m_basis->gtoIndices()[moIndex];
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
double tmpGTO = exp(-m_basis->gtoA()[i] * dr2);
tmpGTO = exp(-m_basis->gtoA()[i] * dr2);
for (unsigned int j = 0; j < 3; ++j) {
// m_values[baseIndex + i] = m_basis->gtoCN()[cIndex++] * tmpGTO;
components[j] += m_basis->gtoCN()[cIndex++] * tmpGTO;
}
}
Expand All @@ -283,15 +293,16 @@ inline void GaussianSetTools::pointD(unsigned int moIndex, const Vector3& delta,

double components[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };

std::vector<double>& gtoA = m_basis->gtoA();
std::vector<double>& gtoCN = m_basis->gtoCN();
const vector<double>& gtoA = m_basis->gtoA();
const vector<double>& gtoCN = m_basis->gtoCN();

// Now iterate through the D type GTOs and sum their contributions
// Now iterate through the GTOs and sum their contributions
unsigned int cIndex = m_basis->cIndices()[moIndex];
double tmpGTO = 0.0;
for (unsigned int i = m_basis->gtoIndices()[moIndex];
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
// Calculate the common factor
double tmpGTO = exp(-gtoA[i] * dr2);
tmpGTO = exp(-gtoA[i] * dr2);
for (double& component : components)
component += gtoCN[cIndex++] * tmpGTO;
}
Expand All @@ -316,15 +327,16 @@ inline void GaussianSetTools::pointD5(unsigned int moIndex,
unsigned int baseIndex = m_basis->moIndices()[moIndex];
double components[5] = { 0.0, 0.0, 0.0, 0.0, 0.0 };

std::vector<double>& gtoA = m_basis->gtoA();
std::vector<double>& gtoCN = m_basis->gtoCN();
const vector<double>& gtoA = m_basis->gtoA();
const vector<double>& gtoCN = m_basis->gtoCN();

// Now iterate through the D type GTOs and sum their contributions
unsigned int cIndex = m_basis->cIndices()[moIndex];
double tmpGTO = 0.0;
for (unsigned int i = m_basis->gtoIndices()[moIndex];
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
// Calculate the common factor
double tmpGTO = exp(-gtoA[i] * dr2);
tmpGTO = exp(-gtoA[i] * dr2);
for (double& component : components)
component += gtoCN[cIndex++] * tmpGTO;
}
Expand Down Expand Up @@ -356,15 +368,16 @@ inline void GaussianSetTools::pointF(unsigned int moIndex, const Vector3& delta,

double components[10] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };

std::vector<double>& gtoA = m_basis->gtoA();
std::vector<double>& gtoCN = m_basis->gtoCN();
const vector<double>& gtoA = m_basis->gtoA();
const vector<double>& gtoCN = m_basis->gtoCN();

// Now iterate through the D type GTOs and sum their contributions
// Now iterate through the GTOs and sum their contributions
unsigned int cIndex = m_basis->cIndices()[moIndex];
double tmpGTO = 0.0;
for (unsigned int i = m_basis->gtoIndices()[moIndex];
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
// Calculate the common factor
double tmpGTO = exp(-gtoA[i] * dr2);
tmpGTO = exp(-gtoA[i] * dr2);
for (double& component : components)
component += gtoCN[cIndex++] * tmpGTO;
}
Expand Down Expand Up @@ -400,15 +413,16 @@ inline void GaussianSetTools::pointF7(unsigned int moIndex,

double components[7] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };

std::vector<double>& gtoA = m_basis->gtoA();
std::vector<double>& gtoCN = m_basis->gtoCN();
const vector<double>& gtoA = m_basis->gtoA();
const vector<double>& gtoCN = m_basis->gtoCN();

// Now iterate through the D type GTOs and sum their contributions
// Now iterate through the F type GTOs and sum their contributions
unsigned int cIndex = m_basis->cIndices()[moIndex];
double tmpGTO = 0.0;
for (unsigned int i = m_basis->gtoIndices()[moIndex];
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
// Calculate the common factor
double tmpGTO = exp(-gtoA[i] * dr2);
tmpGTO = exp(-gtoA[i] * dr2);
for (double& component : components)
component += gtoCN[cIndex++] * tmpGTO;
}
Expand Down Expand Up @@ -456,4 +470,102 @@ final normalization
values[baseIndex + i] += components[i] * componentsF[i];
}

inline void GaussianSetTools::pointG(unsigned int moIndex, const Vector3& delta,
double dr2, vector<double>& values) const
{
// G type orbitals have 15 components and each component has a different
// independent MO weighting. Many things can be cached to save time though.
unsigned int baseIndex = m_basis->moIndices()[moIndex];

double components[15] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };

const vector<double>& gtoA = m_basis->gtoA();
const vector<double>& gtoCN = m_basis->gtoCN();

// Now iterate through the G type GTOs and sum their contributions
unsigned int cIndex = m_basis->cIndices()[moIndex];
double tmpGTO = 0.0;
for (unsigned int i = m_basis->gtoIndices()[moIndex];
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
// Calculate the common factor
tmpGTO = exp(-gtoA[i] * dr2);
for (double& component : components)
component += gtoCN[cIndex++] * tmpGTO;
}

// e.g. XXXX YYYY ZZZZ XXXY XXXZ XYYY YYYZ ZZZX ZZZY XXYY XXZZ YYZZ XXYZ XYYZ
// XYZZ
const double xxxx = delta.x() * delta.x() * delta.x() * delta.x();
const double yyyy = delta.y() * delta.y() * delta.y() * delta.y();
const double zzzz = delta.z() * delta.z() * delta.z() * delta.z();
const double xxxy = delta.x() * delta.x() * delta.x() * delta.y();
const double xxxz = delta.x() * delta.x() * delta.x() * delta.z();
const double yyyx = delta.y() * delta.y() * delta.y() * delta.x();
const double yyyz = delta.y() * delta.y() * delta.y() * delta.z();
const double zzzx = delta.z() * delta.z() * delta.z() * delta.x();
const double zzzy = delta.z() * delta.z() * delta.z() * delta.y();
const double xxyy = delta.x() * delta.x() * delta.y() * delta.y();
const double xxzz = delta.x() * delta.x() * delta.z() * delta.z();
const double yyzz = delta.y() * delta.y() * delta.z() * delta.z();
const double xxyz = delta.x() * delta.x() * delta.y() * delta.z();
const double yyxz = delta.y() * delta.y() * delta.x() * delta.z();
const double zzxy = delta.z() * delta.z() * delta.x() * delta.y();

// molden order
// https://www.theochem.ru.nl/molden/molden_format.html
// https://gau2grid.readthedocs.io/en/latest/order.html
// xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
// xxyy xxzz yyzz xxyz yyxz zzxy
double componentsG[15] = { xxxx, yyyy, zzzz, xxxy, xxxz, yyyx, yyyz, zzzx,
zzzy, xxyy, xxzz, yyzz, xxyz, yyxz, zzxy };

for (int i = 0; i < 15; ++i)
values[baseIndex + i] += components[i] * componentsG[i];
}

inline void GaussianSetTools::pointG9(unsigned int moIndex,
const Vector3& delta, double dr2,
vector<double>& values) const
{
// G type orbitals have 9 spherical components and each component
// has a different independent MO weighting.
// Many things can be cached to save time though.
unsigned int baseIndex = m_basis->moIndices()[moIndex];

double components[9] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };

const vector<double>& gtoA = m_basis->gtoA();
const vector<double>& gtoCN = m_basis->gtoCN();

// Now iterate through the GTOs and sum their contributions
unsigned int cIndex = m_basis->cIndices()[moIndex];
double tmpGTO = 0.0;
for (unsigned int i = m_basis->gtoIndices()[moIndex];
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
// Calculate the common factor
tmpGTO = exp(-gtoA[i] * dr2);
for (double& component : components)
component += gtoCN[cIndex++] * tmpGTO;
}

double x2(delta.x() * delta.x()), y2(delta.y() * delta.y()),
z2(delta.z() * delta.z());

double componentsG[9] = {
3.0 * dr2 * dr2 - 30.0 * dr2 * z2 + 35.0 * z2 * z2 * (1.0 / 8.0),
delta.x() * delta.z() * (7.0 * z2 - 3.0 * dr2) * (sqrt(5.0) / 8.0),
delta.y() * delta.z() * (7.0 * z2 - 3.0 * dr2) * (sqrt(5.0) / 8.0),
(x2 - y2) * (7.0 * z2 - dr2) * (sqrt(5.0) / 4.0),
delta.x() * delta.y() * (7.0 * z2 - dr2) * (sqrt(5.0) / 2.0),
delta.x() * delta.z() * (x2 - 3.0 * y2) * (sqrt(7.0) / 4.0),
delta.y() * delta.z() * (3.0 * x2 - y2) * (sqrt(7.0) / 4.0),
(x2 * x2 - 6.0 * x2 * y2 + y2 * y2) * (sqrt(35.0) / 8.0),
delta.x() * delta.y() * (x2 - y2) * (sqrt(35.0) / 2.0)
};

for (int i = 0; i < 9; ++i)
values[baseIndex + i] += components[i] * componentsG[i];
}

} // namespace Avogadro::Core
4 changes: 4 additions & 0 deletions avogadro/core/gaussiansettools.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,10 @@ class AVOGADROCORE_EXPORT GaussianSetTools
std::vector<double>& values) const;
void pointF7(unsigned int index, const Vector3& delta, double dr2,
std::vector<double>& values) const;
void pointG(unsigned int index, const Vector3& delta, double dr2,
std::vector<double>& values) const;
void pointG9(unsigned int index, const Vector3& delta, double dr2,
std::vector<double>& values) const;

// map from symmetry to angular momentum
// S, SP, P, D, D5, F, F7, G, G9, etc.
Expand Down

0 comments on commit b7e4a12

Please sign in to comment.