Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Implement g orbitals #1826

Merged
merged 3 commits into from
Nov 29, 2024
Merged

Implement g orbitals #1826

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9b261c6
Initial work. Compiles but not normalized
ghutchis Oct 1, 2022
8da0928
Eliminate dangling temporary warning
ghutchis Nov 28, 2024
37b95d5
Added normalization -- needs testing
ghutchis Nov 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 66 additions & 6 deletions avogadro/core/gaussianset.cpp
View file
Open in desktop
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
Open in desktop
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 @@
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 @@
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;

Check notice on line 274 in avogadro/core/gaussiansettools.cpp

View check run for this annotation

Codacy Production / Codacy Static Code Analysis

avogadro/core/gaussiansettools.cpp#L274 

The scope of the variable 'tmpGTO' can be reduced.

Check notice on line 274 in avogadro/core/gaussiansettools.cpp

View check run for this annotation

Codacy Production / Codacy Static Code Analysis

avogadro/core/gaussiansettools.cpp#L274 

Variable 'tmpGTO' is assigned a value that is never used.
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 @@

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 @@
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 @@

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 @@

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 @@
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
Open in desktop
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
Loading