Skip to content

Commit b7e4a12

Browse files
ghutchisghutchis
authored
Merge pull request #1826 from ghutchis/implement-g-orbitals
Implement g orbitals
2 parents 3788209 + 37b95d5 commit b7e4a12

File tree

3 files changed

+200
-24
lines changed

3 files changed

+200
-24
lines changed

‎avogadro/core/gaussianset.cpp

+66-6
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,24 @@ unsigned int GaussianSet::addBasis(unsigned int atom, orbital type)
4444
case F7:
4545
m_numMOs += 7;
4646
break;
47+
case G:
48+
m_numMOs += 15;
49+
break;
50+
case G9:
51+
m_numMOs += 9;
52+
break;
53+
case H:
54+
m_numMOs += 21;
55+
break;
56+
case H11:
57+
m_numMOs += 11;
58+
break;
59+
case I:
60+
m_numMOs += 28;
61+
break;
62+
case I13:
63+
m_numMOs += 13;
64+
break;
4765
default:
4866
// Should never hit here
4967
;
@@ -434,12 +452,54 @@ void GaussianSet::initCalculation()
434452
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.25) * norm); //-3
435453
}
436454
} break;
437-
case G:
438-
skip = 15;
439-
break;
440-
case G9:
441-
skip = 9;
442-
break;
455+
case G: {
456+
// 16 * (2.0/pi)^0.75
457+
double norm = 11.403287525679843;
458+
double norm1 = norm / sqrt(7.0);
459+
double norm2 = norm / sqrt(35.0 / 3.0);
460+
double norm3 = norm / sqrt(35.0);
461+
m_moIndices[i] = indexMO;
462+
indexMO += 15;
463+
m_cIndices.push_back(static_cast<unsigned int>(m_gtoCN.size()));
464+
for (unsigned j = m_gtoIndices[i]; j < m_gtoIndices[i + 1]; ++j) {
465+
// molden order
466+
// xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
467+
// xxyy xxzz yyzz xxyz yyxz zzxy
468+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); // xxxx
469+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); // yyyy
470+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); // zzzz
471+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // xxxy
472+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // xxxz
473+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // yyyx
474+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // yyyz
475+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // zzzx
476+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm1); // zzzy
477+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm2); // xxyy
478+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm2); // xxzz
479+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm2); // yyzz
480+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm3); // xxyz
481+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm3); // yyxz
482+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm3); // zzxy
483+
}
484+
} break;
485+
case G9: {
486+
// 16 * (2.0/pi)^0.75
487+
double norm = 11.403287525679843;
488+
m_moIndices[i] = indexMO;
489+
indexMO += 9;
490+
m_cIndices.push_back(static_cast<unsigned int>(m_gtoCN.size()));
491+
for (unsigned j = m_gtoIndices[i]; j < m_gtoIndices[i + 1]; ++j) {
492+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); // 0
493+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //+1
494+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //-1
495+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //+2
496+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //-2
497+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //+3
498+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //-3
499+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //+4
500+
m_gtoCN.push_back(m_gtoC[j] * pow(m_gtoA[j], 2.75) * norm); //-4
501+
}
502+
} break;
443503
case H:
444504
skip = 21;
445505
break;

‎avogadro/core/gaussiansettools.cpp

+130-18
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,10 @@
99
#include "gaussianset.h"
1010
#include "molecule.h"
1111

12+
#include <iostream>
13+
14+
using std::vector;
15+
1216
namespace Avogadro::Core {
1317

1418
GaussianSetTools::GaussianSetTools(Molecule* mol) : m_molecule(mol)
@@ -227,6 +231,12 @@ inline std::vector<double> GaussianSetTools::calculateValues(
227231
case GaussianSet::F7:
228232
pointF7(i, deltas[atomIndices[i]], dr2[atomIndices[i]], values);
229233
break;
234+
case GaussianSet::G:
235+
pointG(i, deltas[atomIndices[i]], dr2[atomIndices[i]], values);
236+
break;
237+
case GaussianSet::G9:
238+
pointG9(i, deltas[atomIndices[i]], dr2[atomIndices[i]], values);
239+
break;
230240
default:
231241
// Not handled - return a zero contribution
232242
;
@@ -259,13 +269,13 @@ inline void GaussianSetTools::pointP(unsigned int moIndex, const Vector3& delta,
259269
unsigned int baseIndex = m_basis->moIndices()[moIndex];
260270
Vector3 components(Vector3::Zero());
261271

262-
// Now iterate through the P type GTOs and sum their contributions
272+
// Now iterate through the GTOs and sum their contributions
263273
unsigned int cIndex = m_basis->cIndices()[moIndex];
274+
double tmpGTO = 0.0;
264275
for (unsigned int i = m_basis->gtoIndices()[moIndex];
265276
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
266-
double tmpGTO = exp(-m_basis->gtoA()[i] * dr2);
277+
tmpGTO = exp(-m_basis->gtoA()[i] * dr2);
267278
for (unsigned int j = 0; j < 3; ++j) {
268-
// m_values[baseIndex + i] = m_basis->gtoCN()[cIndex++] * tmpGTO;
269279
components[j] += m_basis->gtoCN()[cIndex++] * tmpGTO;
270280
}
271281
}
@@ -283,15 +293,16 @@ inline void GaussianSetTools::pointD(unsigned int moIndex, const Vector3& delta,
283293

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

286-
std::vector<double>& gtoA = m_basis->gtoA();
287-
std::vector<double>& gtoCN = m_basis->gtoCN();
296+
const vector<double>& gtoA = m_basis->gtoA();
297+
const vector<double>& gtoCN = m_basis->gtoCN();
288298

289-
// Now iterate through the D type GTOs and sum their contributions
299+
// Now iterate through the GTOs and sum their contributions
290300
unsigned int cIndex = m_basis->cIndices()[moIndex];
301+
double tmpGTO = 0.0;
291302
for (unsigned int i = m_basis->gtoIndices()[moIndex];
292303
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
293304
// Calculate the common factor
294-
double tmpGTO = exp(-gtoA[i] * dr2);
305+
tmpGTO = exp(-gtoA[i] * dr2);
295306
for (double& component : components)
296307
component += gtoCN[cIndex++] * tmpGTO;
297308
}
@@ -316,15 +327,16 @@ inline void GaussianSetTools::pointD5(unsigned int moIndex,
316327
unsigned int baseIndex = m_basis->moIndices()[moIndex];
317328
double components[5] = { 0.0, 0.0, 0.0, 0.0, 0.0 };
318329

319-
std::vector<double>& gtoA = m_basis->gtoA();
320-
std::vector<double>& gtoCN = m_basis->gtoCN();
330+
const vector<double>& gtoA = m_basis->gtoA();
331+
const vector<double>& gtoCN = m_basis->gtoCN();
321332

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

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

359-
std::vector<double>& gtoA = m_basis->gtoA();
360-
std::vector<double>& gtoCN = m_basis->gtoCN();
371+
const vector<double>& gtoA = m_basis->gtoA();
372+
const vector<double>& gtoCN = m_basis->gtoCN();
361373

362-
// Now iterate through the D type GTOs and sum their contributions
374+
// Now iterate through the GTOs and sum their contributions
363375
unsigned int cIndex = m_basis->cIndices()[moIndex];
376+
double tmpGTO = 0.0;
364377
for (unsigned int i = m_basis->gtoIndices()[moIndex];
365378
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
366379
// Calculate the common factor
367-
double tmpGTO = exp(-gtoA[i] * dr2);
380+
tmpGTO = exp(-gtoA[i] * dr2);
368381
for (double& component : components)
369382
component += gtoCN[cIndex++] * tmpGTO;
370383
}
@@ -400,15 +413,16 @@ inline void GaussianSetTools::pointF7(unsigned int moIndex,
400413

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

403-
std::vector<double>& gtoA = m_basis->gtoA();
404-
std::vector<double>& gtoCN = m_basis->gtoCN();
416+
const vector<double>& gtoA = m_basis->gtoA();
417+
const vector<double>& gtoCN = m_basis->gtoCN();
405418

406-
// Now iterate through the D type GTOs and sum their contributions
419+
// Now iterate through the F type GTOs and sum their contributions
407420
unsigned int cIndex = m_basis->cIndices()[moIndex];
421+
double tmpGTO = 0.0;
408422
for (unsigned int i = m_basis->gtoIndices()[moIndex];
409423
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
410424
// Calculate the common factor
411-
double tmpGTO = exp(-gtoA[i] * dr2);
425+
tmpGTO = exp(-gtoA[i] * dr2);
412426
for (double& component : components)
413427
component += gtoCN[cIndex++] * tmpGTO;
414428
}
@@ -456,4 +470,102 @@ final normalization
456470
values[baseIndex + i] += components[i] * componentsF[i];
457471
}
458472

473+
inline void GaussianSetTools::pointG(unsigned int moIndex, const Vector3& delta,
474+
double dr2, vector<double>& values) const
475+
{
476+
// G type orbitals have 15 components and each component has a different
477+
// independent MO weighting. Many things can be cached to save time though.
478+
unsigned int baseIndex = m_basis->moIndices()[moIndex];
479+
480+
double components[15] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
481+
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
482+
483+
const vector<double>& gtoA = m_basis->gtoA();
484+
const vector<double>& gtoCN = m_basis->gtoCN();
485+
486+
// Now iterate through the G type GTOs and sum their contributions
487+
unsigned int cIndex = m_basis->cIndices()[moIndex];
488+
double tmpGTO = 0.0;
489+
for (unsigned int i = m_basis->gtoIndices()[moIndex];
490+
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
491+
// Calculate the common factor
492+
tmpGTO = exp(-gtoA[i] * dr2);
493+
for (double& component : components)
494+
component += gtoCN[cIndex++] * tmpGTO;
495+
}
496+
497+
// e.g. XXXX YYYY ZZZZ XXXY XXXZ XYYY YYYZ ZZZX ZZZY XXYY XXZZ YYZZ XXYZ XYYZ
498+
// XYZZ
499+
const double xxxx = delta.x() * delta.x() * delta.x() * delta.x();
500+
const double yyyy = delta.y() * delta.y() * delta.y() * delta.y();
501+
const double zzzz = delta.z() * delta.z() * delta.z() * delta.z();
502+
const double xxxy = delta.x() * delta.x() * delta.x() * delta.y();
503+
const double xxxz = delta.x() * delta.x() * delta.x() * delta.z();
504+
const double yyyx = delta.y() * delta.y() * delta.y() * delta.x();
505+
const double yyyz = delta.y() * delta.y() * delta.y() * delta.z();
506+
const double zzzx = delta.z() * delta.z() * delta.z() * delta.x();
507+
const double zzzy = delta.z() * delta.z() * delta.z() * delta.y();
508+
const double xxyy = delta.x() * delta.x() * delta.y() * delta.y();
509+
const double xxzz = delta.x() * delta.x() * delta.z() * delta.z();
510+
const double yyzz = delta.y() * delta.y() * delta.z() * delta.z();
511+
const double xxyz = delta.x() * delta.x() * delta.y() * delta.z();
512+
const double yyxz = delta.y() * delta.y() * delta.x() * delta.z();
513+
const double zzxy = delta.z() * delta.z() * delta.x() * delta.y();
514+
515+
// molden order
516+
// https://www.theochem.ru.nl/molden/molden_format.html
517+
// https://gau2grid.readthedocs.io/en/latest/order.html
518+
// xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
519+
// xxyy xxzz yyzz xxyz yyxz zzxy
520+
double componentsG[15] = { xxxx, yyyy, zzzz, xxxy, xxxz, yyyx, yyyz, zzzx,
521+
zzzy, xxyy, xxzz, yyzz, xxyz, yyxz, zzxy };
522+
523+
for (int i = 0; i < 15; ++i)
524+
values[baseIndex + i] += components[i] * componentsG[i];
525+
}
526+
527+
inline void GaussianSetTools::pointG9(unsigned int moIndex,
528+
const Vector3& delta, double dr2,
529+
vector<double>& values) const
530+
{
531+
// G type orbitals have 9 spherical components and each component
532+
// has a different independent MO weighting.
533+
// Many things can be cached to save time though.
534+
unsigned int baseIndex = m_basis->moIndices()[moIndex];
535+
536+
double components[9] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
537+
538+
const vector<double>& gtoA = m_basis->gtoA();
539+
const vector<double>& gtoCN = m_basis->gtoCN();
540+
541+
// Now iterate through the GTOs and sum their contributions
542+
unsigned int cIndex = m_basis->cIndices()[moIndex];
543+
double tmpGTO = 0.0;
544+
for (unsigned int i = m_basis->gtoIndices()[moIndex];
545+
i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
546+
// Calculate the common factor
547+
tmpGTO = exp(-gtoA[i] * dr2);
548+
for (double& component : components)
549+
component += gtoCN[cIndex++] * tmpGTO;
550+
}
551+
552+
double x2(delta.x() * delta.x()), y2(delta.y() * delta.y()),
553+
z2(delta.z() * delta.z());
554+
555+
double componentsG[9] = {
556+
3.0 * dr2 * dr2 - 30.0 * dr2 * z2 + 35.0 * z2 * z2 * (1.0 / 8.0),
557+
delta.x() * delta.z() * (7.0 * z2 - 3.0 * dr2) * (sqrt(5.0) / 8.0),
558+
delta.y() * delta.z() * (7.0 * z2 - 3.0 * dr2) * (sqrt(5.0) / 8.0),
559+
(x2 - y2) * (7.0 * z2 - dr2) * (sqrt(5.0) / 4.0),
560+
delta.x() * delta.y() * (7.0 * z2 - dr2) * (sqrt(5.0) / 2.0),
561+
delta.x() * delta.z() * (x2 - 3.0 * y2) * (sqrt(7.0) / 4.0),
562+
delta.y() * delta.z() * (3.0 * x2 - y2) * (sqrt(7.0) / 4.0),
563+
(x2 * x2 - 6.0 * x2 * y2 + y2 * y2) * (sqrt(35.0) / 8.0),
564+
delta.x() * delta.y() * (x2 - y2) * (sqrt(35.0) / 2.0)
565+
};
566+
567+
for (int i = 0; i < 9; ++i)
568+
values[baseIndex + i] += components[i] * componentsG[i];
569+
}
570+
459571
} // namespace Avogadro::Core

‎avogadro/core/gaussiansettools.h

+4
Original file line numberDiff line numberDiff line change
@@ -124,6 +124,10 @@ class AVOGADROCORE_EXPORT GaussianSetTools
124124
std::vector<double>& values) const;
125125
void pointF7(unsigned int index, const Vector3& delta, double dr2,
126126
std::vector<double>& values) const;
127+
void pointG(unsigned int index, const Vector3& delta, double dr2,
128+
std::vector<double>& values) const;
129+
void pointG9(unsigned int index, const Vector3& delta, double dr2,
130+
std::vector<double>& values) const;
127131

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

0 commit comments

Comments
 (0)