Category: Quantum Mechanics and Quantum Chemistry
Blog Entry © Saturday, August 16, 2025, by James Pate Williams, Jr. Some More Elementary Quantum Chemistry
Blog Entry © Friday, August 15, 2025, by James Pate Williams, Jr. Some Elementary Quantum Chemistry
Blog Entry (c) Wednesday, August 13, 2025, by James Pate Williams, Jr. Exercises from an Online Textbook
#include <complex>
#include <vector>
class CmpLinearAlgebra
{
public:
static void CmpPrintMatrix(
int m, int n,
std::vector<std::vector<std::complex<double>>>& Ac);
static void CmpAddition(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpSubtraction(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpMultiply(
size_t m, size_t n, size_t p,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpAnticommutator(
size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpCommutator(
size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpAdjoint(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& Ac,
std::vector<std::vector<std::complex<double>>>& Ad);
static std::complex<double> CmpDeterminant(
bool& failure, int n,
std::vector<std::vector<std::complex<double>>>& A);
static bool CmpGaussianElimination(
int m, int n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot);
static bool CmpGaussianFactor(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<size_t>& pivot);
static bool CmpGaussianSolution(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot);
static bool CmpSubstitution(
int m, int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot);
static bool CmpInverse(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& Mi);
static void CmpCharPolyAndAdjoint(
int n,
std::vector<std::vector<std::complex<double>>>& C,
std::vector<std::vector<std::complex<double>>>& I,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& adjoint,
std::vector<std::complex<double>>& a);
static void CmpMatrixKernel(
int m, int n,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& X,
size_t& r);
static void CmpMatrixImage(
int m, int n,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& N,
std::vector<std::vector<std::complex<double>>>& X,
int rank);
};
#include <vector>
class DblLinearAlgebra
{
public:
static void DblPrintMatrix(
int m, int n, std::vector<std::vector<double>>& A);
static void DblAddition(
size_t m, size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblSubtraction(
size_t m, size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblMultiply(
size_t m, size_t n, size_t p,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblAnticommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblCommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static double DblDeterminant(
bool& failure, int n,
std::vector<std::vector<double>>& A);
static bool DblGaussianElimination(
int m, int n, std::vector<std::vector<double>>& A,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot);
static bool DblGaussianFactor(
int n, std::vector<std::vector<double>>& M,
std::vector<size_t>& pivot);
static bool DblGaussianSolution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot);
static bool DblSubstitution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot);
static bool DblInverse(
int n, std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& A);
static void DblCharPolyAndAdjoint(
int n,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& I,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& adjoint,
std::vector<double>& a);
static void DblMatrixKernel(
int m, int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X,
size_t& r);
static void DblMatrixImage(
int m, int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& N,
std::vector<std::vector<double>>& X,
int rank);
};
// Exercises from "Modern Quantum Chemistry An Introduction to Advanced
// Electronic Structure Theory" by Attila Szabo and Neil S. Ostlund
// https://chemistlibrary.wordpress.com/wp-content/uploads/2015/02/modern-quantum-chemistry.pdf
// Program (c) Tuesday, August 12, 2025 by James Pate Williams, Jr.
#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include "DblLinearAlgebra.h"
#include "CmpLinearAlgebra.h"
int main()
{
double AArcb[3][3] = { { 2, 3, -1 }, { 4, 4, -3 }, { -2, 3, -1 } };
double AArso[3][3] = { { 1, 1, 0 }, { 1, 2, 2 }, { 0, 2, -1 } };
double BBrso[3][3] = { { 1, -1, 1 }, { -1, 0, 0 }, { 1, 0, 1} };
double BBr[3][3] = { { 1, -1, 1 }, { -1 , 0, 0 }, { 1, 0, 1 } };
double AAcr[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
double AAci[3][3] = { { 1, 1, 2 }, { 3, 0, 1 }, { 0, 2, 4 } };
double BBcr[3][3] = { { 1, 0, 1 }, { 1 , 1, 0 }, { 0, 1, 1 } };
double BBci[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
int m = 3, n = 3, p = 3;
std::vector<double> br(3);
std::vector<size_t> pivot(3);
std::vector<std::vector<double>> Arcb(3, std::vector<double>(3));
std::vector<std::vector<double>> Arso(3, std::vector<double>(3));
std::vector<std::vector<double>> Brso(3, std::vector<double>(3));
std::vector<std::vector<double>> Br(3, std::vector<double>(3));
std::vector<std::vector<double>> Cr(3, std::vector<double>(3));
std::vector<std::vector<double>> Ai(3, std::vector<double>(3));
std::vector<std::vector<double>> Ari(3, std::vector<double>(3));
std::vector<std::vector<std::complex<double>>> Ac(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Bc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Cc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Dc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Ec(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Fc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Gc(3,
std::vector<std::complex<double>>(3));
for (int i = 0; i < m; i++)
{
for (int j = 0; j < p; j++)
{
Arcb[i][j] = AArcb[i][j];
Arso[i][j] = AArso[i][j];
Brso[i][j] = BBrso[i][j];
Ac[i][j]._Val[0] = AAcr[i][j];
Ac[i][j]._Val[1] = AAci[i][j];
}
}
for (int i = 0; i < p; i++)
{
for (int j = 0; j < n; j++)
{
Br[i][j] = BBr[i][j];
Bc[i][j]._Val[0] = BBcr[i][j];
Bc[i][j]._Val[1] = BBci[i][j];
}
}
DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Br, Cr);
std::cout << "Ar * Br = Cr Conte & de Boor" << std::endl;
DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
std::cout << std::endl;
CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
std::cout << "Ac * Bc = Cc" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
std::cout << std::endl;
// Exercise 1.2
std::cout << "Exercise 1.2 page 5 Commutator" << std::endl;
DblLinearAlgebra::DblCommutator(3, Arso, Brso, Cr);
DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
std::cout << std::endl;
std::cout << "Exercise 1.2 page 5 Anticommutator" << std::endl;
DblLinearAlgebra::DblAnticommutator(3, Arso, Brso, Cr);
DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
std::cout << std::endl;
CmpLinearAlgebra::CmpAdjoint(3, 3, Cc, Dc);
std::cout << "Exercise 1.3 page 6 Cc adjoint" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Dc);
std::cout << std::endl;
CmpLinearAlgebra::CmpAdjoint(3, 3, Ac, Ec);
CmpLinearAlgebra::CmpAdjoint(3, 3, Bc, Fc);
CmpLinearAlgebra::CmpMultiply(3, 3, 3, Fc, Ec, Gc);
std::cout << "Exercise 1.3 page 6 Bc adjoint * Ac adjoint" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Gc);
std::cout << std::endl;
std::cout << "Ar matrix" << std::endl;
DblLinearAlgebra::DblPrintMatrix(3, 3, Arcb);
bool inv = DblLinearAlgebra::DblInverse(n, Arcb, Ai);
std::cout << std::endl;
std::cout << "Ar Conte & de Boor inverse = " << inv << std::endl;
DblLinearAlgebra::DblPrintMatrix(3, 3, Ai);
std::cout << std::endl;
std::cout << "Ar * Ar inverse" << std::endl;
DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Ai, Ari);
DblLinearAlgebra::DblPrintMatrix(3, 3, Ari);
std::cout << std::endl;
std::cout << "Ac" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Ac);
std::cout << std::endl;
inv = CmpLinearAlgebra::CmpInverse(3, Ac, Bc);
std::cout << "Ac inverse = " << inv << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Bc);
CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
std::cout << std::endl;
std::cout << "Ac * AC inverse" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
}
Blog Entry © Monday, August 11, 2025, by James Pate Williams, Jr. Van der Waals Interaction Between Two Hydrogen Atoms
Blog Entry © Sunday, August 10, 2025, First-Order Perturbation Treatment of the Helium Atom by James Pate Williams, Jr.
Blog Entry © Tuesday, July 29, 2025, Double and Triple Monte Carlo Integration by James Pate Williams, Jr.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
static double randomRange(double lo, double hi)
{
return (hi - lo) * (double)rand() / RAND_MAX + lo;
}
static double integrand(double r, double w)
{
return pow(r, 4.0) * (2.0 - r) * w * w * exp(-r);
}
static double StarkEffectIntegral(double E, int N)
{
double sum = 0.0;
for (int i = 0; i <= N; i++)
{
double r = randomRange(0.0, 100.0);
double w = randomRange(-1.0, 1.0);
sum += integrand(r, w);
}
return 100.0 * 2.0 * E * sum / (16.0 * (N - 1));
}
static void firstOrderStarkEffect(double E)
{
double exact = -3.0 * E;
int N[9] = {
1000000, 2000000, 3000000, 4000000,
5000000, 6000000, 7000000, 8000000,
9000000 };
for (int n = 0; n < 9; n++)
{
int iN = N[n];
double integ = StarkEffectIntegral(E, iN);
double error = 100.0 * fabs(integ - exact) / fabs(exact);
printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
iN, integ, error);
}
printf("exact value = %13.10lf\n", exact);
}
static double ee1(int N, double R, double Z)
{
double pi = 4.0 * atan(1.0);
double sum = 0.0;
for (int i = 0; i <= N; i++)
{
double r1 = randomRange(1.0e-25, R);
double r2 = randomRange(0.0, r1);
sum += R * r1 * r1 * exp(-2.0 * Z * (r1 + r2)) * r2 * r2;
}
return 16.0 * pi * pi * sum / (N - 1);
}
static double ee2(int N, double R, double Z)
{
double pi = 4.0 * atan(1.0);
double sum = 0.0;
for (int i = 0; i <= N; i++)
{
double r1 = randomRange(1.0e-25, R);
double r2 = randomRange(r1, R);
sum += R * (R - r2) * r2 * exp(-2.0 * Z * (r1 + r2)) * r1 * r1;
}
return 16.0 * pi * pi * sum / (N - 1);
}
static void firstOrderHelium(double Z)
{
double pi = 4.0 * atan(1.0), R = 25.0;
double exact = 5.0 * pi * pi / (8.0 * pow(Z, 5.0));
int N[9] = {
1000000, 2000000, 3000000, 4000000,
5000000, 6000000, 7000000, 8000000,
9000000 };
for (int n = 0; n < 9; n++)
{
int iN = N[n];
double integ = ee1(iN, R, Z) + ee2(iN, R, Z);
double error = 100.0 * fabs(integ - exact) / fabs(exact);
printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
iN, integ, error);
}
printf("exact value = %13.10lf\n", exact);
}
int main(void)
{
firstOrderStarkEffect(2.0);
firstOrderHelium(2.0);
return 0;
}
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
static double randomRange(double lo, double hi)
{
return (hi - lo) * (double)rand() / RAND_MAX + lo;
}
static double f(double x, double y, double z)
{
return pow(sin(x), 2.0) + y * sin(z);
}
static double g(double x, double y, double z)
{
return x + y * z * z;
}
static double integral(
double x0, double x1,
double y0, double y1,
double z0, double z1,
double (*f)(double, double, double),
int N)
{
double sum = 0.0;
for (int n = 0; n <= N; n++)
{
double x = randomRange(x0, x1);
double y = randomRange(y0, y1);
double z = randomRange(z0, z1);
sum += f(x, y, z);
}
return (x1 - x0) * (y1 - y0) * (z1 - z0) *
sum / (N - 1);
}
int main(void)
{
double pi = 4.0 * atan(1.0);
double x0 = 0.0, x1 = pi;
double y0 = 0.0, y1 = 1.0;
double z0 = 0.0, z1 = pi;
double exact = 0.5 * pi * (2.0 + pi);
int N[9] = {
1000000, 2000000, 3000000, 4000000,
5000000, 6000000, 7000000, 8000000,
9000000 };
printf("integrand pow(sin(x), 2.0) + y * sin(z)\n");
printf("x = 0 to pi, y = 0 to 1, z = 0 to pi\n");
for (int n = 0; n < 9; n++)
{
int iN = N[n];
double integ = integral(
x0, x1, y0, y1, z0, z1, f, iN);
double error = 100.0 * fabs(integ - exact) / fabs(exact);
printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
iN, integ, error);
}
printf("exact value = %13.10lf\n", exact);
x0 = -1.0;
x1 = 5.0;
y0 = 2.0;
y1 = 4.0;
z0 = 0.0;
z1 = 1.0;
exact = 36.0;
printf("integrand x + y * z * z\n");
printf("x = -1 to 5, y = 2 to 4, z = 0 to 1\n");
for (int n = 0; n < 9; n++)
{
int iN = N[n];
double integ = integral(
x0, x1, y0, y1, z0, z1, g, iN);
double error = 100.0 * fabs(integ - exact) / fabs(exact);
printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
iN, integ, error);
}
printf("exact value = %13.10lf\n", exact);
return 0.0;
}
Four Methods of Numerical Double Integration: Sequential Simpson’s Rule, Multitasking Simpson’s Rule, Sequential Monte Carlo and Multitasking Monte Carlo Methods © Wednesday April 16 – 18, 2025, by James Pate Williams, Jr.
Approximation of the Ground-State Total Energy of a Beryllium Atom © Sunday, March 30 to Tuesday April 1, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Computer Science
Blog Entry © Sunday, March 29, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Slater Determinant Coefficients for Z = 2 to 4
Enter the atomic number Z (2 to 6 or 0 to quit): 2
2 1 1 + a(1)b(2)
1 0 0 - a(2)b(1)
# Even Permutations = 1
Enter the atomic number Z (2 to 6 or 0 to quit): 3
6 3 1 + a(1)b(2)c(3)
5 2 0 - a(1)b(3)c(2)
4 2 0 - a(2)b(1)c(3)
3 1 1 + a(2)b(3)c(1)
2 1 1 + a(3)b(1)c(2)
1 0 0 - a(3)b(2)c(1)
# Even Permutations = 3
Enter the atomic number Z (2 to 6 or 0 to quit): 4
24 12 0 + a(1)b(2)c(3)d(4)
23 11 1 - a(1)b(2)c(4)d(3)
22 11 1 - a(1)b(3)c(2)d(4)
21 10 0 + a(1)b(3)c(4)d(2)
20 10 0 + a(1)b(4)c(2)d(3)
19 9 1 - a(1)b(4)c(3)d(2)
18 9 1 - a(2)b(1)c(3)d(4)
17 8 0 + a(2)b(1)c(4)d(3)
16 8 0 + a(2)b(3)c(1)d(4)
15 7 1 - a(2)b(3)c(4)d(1)
14 7 1 - a(2)b(4)c(1)d(3)
13 6 0 + a(2)b(4)c(3)d(1)
12 6 0 + a(3)b(1)c(2)d(4)
11 5 1 - a(3)b(1)c(4)d(2)
10 5 1 - a(3)b(2)c(1)d(4)
9 4 0 + a(3)b(2)c(4)d(1)
8 4 0 + a(3)b(4)c(1)d(2)
7 3 1 - a(3)b(4)c(2)d(1)
6 3 1 - a(4)b(1)c(2)d(3)
5 2 0 + a(4)b(1)c(3)d(2)
4 2 0 + a(4)b(2)c(1)d(3)
3 1 1 - a(4)b(2)c(3)d(1)
2 1 1 - a(4)b(3)c(1)d(2)
1 0 0 + a(4)b(3)c(2)d(1)
# Even Permutations = 12
Enter the atomic number Z (2 to 6 or 0 to quit):
// AOPermutations.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Saturday, March 29, 2025
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Signs of the atomic orbitals in a Slater Determinant
#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
int main()
{
char alpha[] = { 'a', 'b', 'c', 'd', 'e', 'f' }, line[128] = {};
int factorial[7] = { 1, 1, 2, 6, 24, 120, 720 };
while (true)
{
int col = 0, counter = 0, row = 0, sign = 1, t = 0, Z = 0, zfact = 0;
int numberEven = 0;
std::cout << "Enter the atomic number Z (2 to 6 or 0 to quit): ";
std::cin.getline(line, 127);
std::string str(line);
Z = std::stoi(str);
if (Z == 0)
{
break;
}
if (Z < 2 || Z > 6)
{
std::cout << "Illegal Z, please try again" << std::endl;
continue;
}
zfact = factorial[Z];
std::vector<char> orb(Z);
std::vector<int> tmp(Z), vec(Z);
for (int i = 0; i < Z; i++)
{
orb[i] = alpha[i];
vec[i] = i + 1;
}
do
{
for (int i = 0; i < (int)vec.size(); i++)
{
tmp[i] = vec[i];
}
t = 0;
do
{
t++;
} while (std::next_permutation(tmp.begin(), tmp.end()));
std::cout << t << '\t' << t / 2 << '\t';
std::cout << (t / 2 & 1) << '\t';
if (Z == 2 || Z == 3)
{
if ((t / 2 & 1) == 0)
{
std::cout << "-\t";
}
else
{
std::cout << "+\t";
numberEven++;
}
}
else
{
if ((t / 2 & 1) == 1)
{
std::cout << "-\t";
}
else
{
std::cout << "+\t";
numberEven++;
}
}
for (int i = 0; i < Z; i++)
{
std::cout << orb[i] << '(' << vec[i] << ')';
}
row++;
std::cout << std::endl;
if (zfact != 2 && row == zfact)
{
std::cout << std::endl;
break;
}
row %= Z;
} while (std::next_permutation(vec.begin(), vec.end()));
std::cout << "# Even Permutations = ";
std::cout << numberEven << std::endl;
}
return 0;
}
Numerical Explorations 