Category: Hydrogen Atom

Blog Entry © Tuesday, July 29, 2025, Double and Triple Monte Carlo Integration by James Pate Williams, Jr.

Blog Entry Tuesday July 29 2025Download
#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;
}

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

Approximation of the Ground-State Total Energy of a Beryllium AtomDownload

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;
}

Blog Entry © Thursday, March 27, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Lithium (Li, Z = 3) Total Ground-State Energy Numerical Experiments

Lithium Atom ExperimentsDownload

Blog Entry © Tuesday, March 25, 2025, by James Pate Williams, Jr. Hydrogen Radial Wavefunctions and Related Functions

Hydrogen Radial Wavefunctions and Related FunctionsDownload

Revised Helium Ground-State Total Energy Computation (c) Sunday, March 23, 2025, by James Pate Williams, Jr.

Revised Computation of the Ground-State Energy of a Helium AtomDownload

Blog Entry © Saturday, November 30, 2024, by James Pate Williams, Jr. Partial Solution of the Schrödinger Equation for Hydrogen in Parabolic Coordinates

Solution of the Schrödinger Equation for Hydrogen in Parabolic CoordinatesDownload

Blog Entry © Friday, November 1, 2024, by James Pate Williams, Jr. Calculation of the Overlap Matrix for the Water Molecule (H2O) Using a Contracted Set of Gaussian Orbitals

Reference: https://content.wolfram.com/sites/19/2012/02/Ho.pdf

I reproduced most of the computations in the MATHMATICA reference. The water molecule is a planar molecule that lies in the YZ-plane.

Gaussian Overlap Integral CPP Source CodeDownload

Blog Entry (c) Monday, October 28, 2024, by James Pate Williams, Jr. Two Methods of Computing the Gaussian Type Orbital 1s Integrals (Corrected Version)

Blog Entry Monday October 24 2024Download
HFIntegrals1s CPP Source CodeDownload

Blog Entry (c) Tuesday, October 22, 2024, by James Pate Williams, Jr. Selected Exercises from “Modern Quantum Chemistry Introduction to Advanced Electronic Structure Theory” by Attila Szabo and Neil S. Ostlund (Dover Books and Kindle)

Exercises from Modern Quantum ChemistryDownload
Commutator Modern Quantum Chemistry C Source CodeDownload

Source code for the solutions to Exercise 1.18 and 1.19.

Exercise 1_18 Modern Quantum Chemistry Source CodeDownload