From 9b261c667a8337545d1a3dfbaa416c96e6894ab0 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sat, 1 Oct 2022 19:05:36 -0400 Subject: [PATCH 1/3] Initial work. Compiles but not normalized Signed-off-by: Geoff Hutchison --- avogadro/core/gaussianset.cpp | 18 ++++ avogadro/core/gaussiansettools.cpp | 130 ++++++++++++++++++++++++++--- avogadro/core/gaussiansettools.h | 4 + 3 files changed, 140 insertions(+), 12 deletions(-) diff --git a/avogadro/core/gaussianset.cpp b/avogadro/core/gaussianset.cpp index 152b180218..ea3c2f4524 100644 --- a/avogadro/core/gaussianset.cpp +++ b/avogadro/core/gaussianset.cpp @@ -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 ; diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp index fd5a80ef44..4c7c76d99a 100644 --- a/avogadro/core/gaussiansettools.cpp +++ b/avogadro/core/gaussiansettools.cpp @@ -9,6 +9,10 @@ #include "gaussianset.h" #include "molecule.h" +#include + +using std::vector; + namespace Avogadro::Core { GaussianSetTools::GaussianSetTools(Molecule* mol) : m_molecule(mol) @@ -227,6 +231,12 @@ inline std::vector 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 ; @@ -259,7 +269,7 @@ 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]; for (unsigned int i = m_basis->gtoIndices()[moIndex]; i < m_basis->gtoIndices()[moIndex + 1]; ++i) { @@ -283,10 +293,10 @@ 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& gtoA = m_basis->gtoA(); - std::vector& gtoCN = m_basis->gtoCN(); + const vector& gtoA = m_basis->gtoA(); + const vector& 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]; for (unsigned int i = m_basis->gtoIndices()[moIndex]; i < m_basis->gtoIndices()[moIndex + 1]; ++i) { @@ -316,8 +326,8 @@ 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& gtoA = m_basis->gtoA(); - std::vector& gtoCN = m_basis->gtoCN(); + const vector& gtoA = m_basis->gtoA(); + const vector& gtoCN = m_basis->gtoCN(); // Now iterate through the D type GTOs and sum their contributions unsigned int cIndex = m_basis->cIndices()[moIndex]; @@ -356,10 +366,10 @@ 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& gtoA = m_basis->gtoA(); - std::vector& gtoCN = m_basis->gtoCN(); + const vector& gtoA = m_basis->gtoA(); + const vector& 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]; for (unsigned int i = m_basis->gtoIndices()[moIndex]; i < m_basis->gtoIndices()[moIndex + 1]; ++i) { @@ -400,10 +410,10 @@ 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& gtoA = m_basis->gtoA(); - std::vector& gtoCN = m_basis->gtoCN(); + const vector& gtoA = m_basis->gtoA(); + const vector& 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]; for (unsigned int i = m_basis->gtoIndices()[moIndex]; i < m_basis->gtoIndices()[moIndex + 1]; ++i) { @@ -456,4 +466,100 @@ final normalization values[baseIndex + i] += components[i] * componentsF[i]; } +inline void GaussianSetTools::pointG(unsigned int moIndex, const Vector3& delta, + double dr2, vector& 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& gtoA = m_basis->gtoA(); + const vector& gtoCN = m_basis->gtoCN(); + + // Now iterate through the G type GTOs and sum their contributions + unsigned int cIndex = m_basis->cIndices()[moIndex]; + 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); + 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 xyyy = delta.x() * delta.y() * delta.y() * delta.y(); + 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 xyyz = delta.x() * delta.y() * delta.y() * delta.z(); + const double xyzz = delta.x() * delta.y() * delta.z() * delta.z(); + + // 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& 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& gtoA = m_basis->gtoA(); + const vector& gtoCN = m_basis->gtoCN(); + + // Now iterate through the GTOs and sum their contributions + unsigned int cIndex = m_basis->cIndices()[moIndex]; + 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); + 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 diff --git a/avogadro/core/gaussiansettools.h b/avogadro/core/gaussiansettools.h index 32fc5cb50d..276a0d2336 100644 --- a/avogadro/core/gaussiansettools.h +++ b/avogadro/core/gaussiansettools.h @@ -124,6 +124,10 @@ class AVOGADROCORE_EXPORT GaussianSetTools std::vector& values) const; void pointF7(unsigned int index, const Vector3& delta, double dr2, std::vector& values) const; + void pointG(unsigned int index, const Vector3& delta, double dr2, + std::vector& values) const; + void pointG9(unsigned int index, const Vector3& delta, double dr2, + std::vector& values) const; // map from symmetry to angular momentum // S, SP, P, D, D5, F, F7, G, G9, etc. From 8da09287009a88ee41b85a39c886941ed08d9d99 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Thu, 28 Nov 2024 16:46:33 -0500 Subject: [PATCH 2/3] Eliminate dangling temporary warning Signed-off-by: Geoff Hutchison --- avogadro/core/gaussiansettools.cpp | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp index 4c7c76d99a..5656889b7f 100644 --- a/avogadro/core/gaussiansettools.cpp +++ b/avogadro/core/gaussiansettools.cpp @@ -271,11 +271,11 @@ inline void GaussianSetTools::pointP(unsigned int moIndex, const Vector3& delta, // 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; } } @@ -298,10 +298,11 @@ inline void GaussianSetTools::pointD(unsigned int moIndex, const Vector3& delta, // 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; } @@ -331,10 +332,11 @@ inline void GaussianSetTools::pointD5(unsigned int moIndex, // 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; } @@ -371,10 +373,11 @@ inline void GaussianSetTools::pointF(unsigned int moIndex, const Vector3& delta, // 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; } @@ -415,10 +418,11 @@ inline void GaussianSetTools::pointF7(unsigned int moIndex, // 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; } @@ -481,10 +485,11 @@ inline void GaussianSetTools::pointG(unsigned int moIndex, const Vector3& delta, // 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 - double tmpGTO = exp(-gtoA[i] * dr2); + tmpGTO = exp(-gtoA[i] * dr2); for (double& component : components) component += gtoCN[cIndex++] * tmpGTO; } @@ -496,7 +501,7 @@ inline void GaussianSetTools::pointG(unsigned int moIndex, const Vector3& delta, 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 xyyy = delta.x() * delta.y() * delta.y() * delta.y(); + 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(); @@ -504,8 +509,8 @@ inline void GaussianSetTools::pointG(unsigned int moIndex, const Vector3& delta, 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 xyyz = delta.x() * delta.y() * delta.y() * delta.z(); - const double xyzz = delta.x() * delta.y() * delta.z() * 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 @@ -535,10 +540,11 @@ inline void GaussianSetTools::pointG9(unsigned int moIndex, // 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; } From 37b95d5d3a42e664d071b9ad72b967038214aa94 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Thu, 28 Nov 2024 18:00:38 -0500 Subject: [PATCH 3/3] Added normalization -- needs testing Signed-off-by: Geoff Hutchison --- avogadro/core/gaussianset.cpp | 54 +++++++++++++++++++++++++++++++---- 1 file changed, 48 insertions(+), 6 deletions(-) diff --git a/avogadro/core/gaussianset.cpp b/avogadro/core/gaussianset.cpp index ea3c2f4524..05e66df531 100644 --- a/avogadro/core/gaussianset.cpp +++ b/avogadro/core/gaussianset.cpp @@ -452,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(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(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;