Category: Quantum Mechanics and Quantum Chemistry
Blog Entry © Sunday, November 9, 2025, by James Pate Williams, Jr. Hydrogenic Wavefunctions, Radial Probability Functions, Distribution Functions, and First Moment Integrals
Blog Entry © Friday, October 31, 2025, by James Pate Williams, Jr. Quantum Mechanical Harmonic Oscillator Matrix Elements
// QMHOMatrixElements.cpp : Defines the entry point for the application.
//
#include "pch.h"
#include "framework.h"
#include "QMHOMatrixElements.h"
#include "Integration.h"
#define MAX_LOADSTRING 100
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
WCHAR buffer[32768], line[128];
int n, m;
// Forward declarations of functions included in this code module:
ATOM MyRegisterClass(HINSTANCE hInstance);
BOOL InitInstance(HINSTANCE, int);
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK About(HWND, UINT, WPARAM, LPARAM);
int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
_In_opt_ HINSTANCE hPrevInstance,
_In_ LPWSTR lpCmdLine,
_In_ int nCmdShow)
{
UNREFERENCED_PARAMETER(hPrevInstance);
UNREFERENCED_PARAMETER(lpCmdLine);
// TODO: Place code here.
// Initialize global strings
LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
LoadStringW(hInstance, IDC_QMHOMATRIXELEMENTS, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_QMHOMATRIXELEMENTS));
MSG msg;
// Main message loop:
while (GetMessage(&msg, nullptr, 0, 0))
{
if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
}
return (int) msg.wParam;
}
//
// FUNCTION: MyRegisterClass()
//
// PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
WNDCLASSEXW wcex = { };
wcex.cbSize = sizeof(WNDCLASSEX);
wcex.style = CS_HREDRAW | CS_VREDRAW;
wcex.lpfnWndProc = WndProc;
wcex.cbClsExtra = 0;
wcex.cbWndExtra = 0;
wcex.hInstance = hInstance;
wcex.hIcon = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_QMHOMATRIXELEMENTS));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_QMHOMATRIXELEMENTS);
wcex.lpszClassName = szWindowClass;
wcex.hIconSm = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));
return RegisterClassExW(&wcex);
}
//
// FUNCTION: InitInstance(HINSTANCE, int)
//
// PURPOSE: Saves instance handle and creates main window
//
// COMMENTS:
//
// In this function, we save the instance handle in a global variable and
// create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
hInst = hInstance; // Store instance handle in our global variable
HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);
if (!hWnd)
{
return FALSE;
}
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
return TRUE;
}
static double Factorial(int n)
{
// n! = n ... 2 1
double factorial = 1.0;
for (int i = 2; i <= n; i++)
{
factorial *= i;
}
return factorial;
}
static double N(int n)
{
// normilization factor
double pi = 4.0 * atan(1.0);
double factor1 = sqrt(1.0 / pi);
double factor2 = 1.0 / (pow(2.0, n) * Factorial(n));
return sqrt(factor1 * factor2);
}
static double H(double xi, int n)
{
// Hermite polynomial
if (n == 1)
{
return 2.0 * xi;
}
else if (n == 2)
{
return 4.0 * xi * xi - 2.0;
}
else if (n == 3)
{
return xi * (8.0 * xi * xi - 12.0);
}
else if (n == 4)
{
return 16.0 * pow(xi, 4.0) - 48.0 * pow(xi, 2.0) + 12.0;
}
else if (n == 5)
{
return 32.0 * pow(xi, 5.0) - 160.0 * pow(xi, 3.0) +
120.0 * xi;
}
else if (n == 6)
{
return 64.0 * pow(xi, 6.0) - 480.0 * pow(xi, 4.0) +
720.0 * xi * xi - 120.0;
}
else if (n == 7)
{
return 128.0 * pow(xi, 7.0) - 1344.0 * pow(xi, 5.0) +
3360.0 * xi * xi * xi - 1680.0 * xi;
}
else if (n == 8)
{
return 256.0 * pow(xi, 8.0) - 3584.0 * pow(xi, 6.0) +
13440.0 * pow(xi, 4.0) - 13440.0 * xi * xi + 1680.0;
}
else if (n == 9)
{
return 512.0 * pow(xi, 9.0) - 9216.0 * pow(xi, 7.0) +
48384.0 * pow(xi, 5.0) - 80640.0 * xi * xi * xi + 30240.0 * xi;
}
else if (n == 10)
{
return 1024.0 * pow(xi, 10.0) - 23040.0 * pow(xi, 8.0) +
161280.0 * pow(xi, 6.0) - 403200.0 * pow(xi, 4.0) +
302400.0 * xi * xi - 30240.0;
}
else
{
return 0.0;
}
}
static double Psi(double xi, int n)
{
return N(n) * exp(-xi * xi / 2.0) * H(xi, n);
}
static double Dx(double xi)
{
return xi * Psi(xi, n) * Psi(xi, m);
}
static double Ex(double xi)
{
return xi * xi * Psi(xi, n) * Psi(xi, m);
}
static double Fx(double xi)
{
return pow(xi, 3.0) * Psi(xi, n) * Psi(xi, m);
}
static double Gx(double xi)
{
return pow(xi, 4.0) * Psi(xi, n) * Psi(xi, m);
}
#define IDC_COMPUTE_BUTTON 1010
#define IDC_CANCEL_BUTTON 1020
#define IDC_CLEAR_BUTTON 1030
#define IDC_EDIT_MULTILINE 1040
//
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
static HFONT hFont = nullptr;
static HWND hCombo1 = nullptr;
static HWND hCombo2 = nullptr;
static HWND hCombo3 = nullptr;
static HWND hEditMultiline = nullptr;
switch (message)
{
case WM_CREATE:
CreateWindowEx(0, L"STATIC", L"n:", WS_CHILD | WS_VISIBLE,
10, 10, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
hCombo1 = CreateWindowEx(0, L"COMBOBOX", nullptr,
WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
ES_AUTOVSCROLL, 120, 10, 120, 100, hWnd, nullptr, hInst, nullptr);
CreateWindowEx(0, L"STATIC", L"m:", WS_CHILD | WS_VISIBLE,
10, 40, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
hCombo2 = CreateWindowEx(0, L"COMBOBOX", nullptr,
WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
ES_AUTOVSCROLL, 120, 40, 120, 150, hWnd, nullptr, hInst, nullptr);
hEditMultiline = CreateWindowEx(
WS_EX_CLIENTEDGE, // Extended style for sunken border
TEXT("EDIT"), // Class name
TEXT(""), // Initial text (can be blank)
WS_CHILD | WS_VISIBLE | WS_VSCROLL | ES_AUTOHSCROLL |
ES_LEFT | ES_MULTILINE | ES_AUTOVSCROLL | WS_HSCROLL | WS_VSCROLL,
250, 10, 600, 420, // Position and size
hWnd, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
CreateWindowEx(0, L"BUTTON", L"Compute", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 80, 80, 30, hWnd, (HMENU)IDC_COMPUTE_BUTTON, hInst, NULL);
CreateWindowEx(0, L"BUTTON", L"Cancel", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 120, 80, 30, hWnd, (HMENU)IDC_CANCEL_BUTTON, hInst, NULL);
CreateWindowEx(0, L"BUTTON", L"Clear", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 160, 80, 30, hWnd, (HMENU)IDC_CLEAR_BUTTON, hInst, NULL);
hFont = CreateFont(16, 0, 0, 0, FW_BOLD, FALSE, FALSE, FALSE,
UNICODE, OUT_DEFAULT_PRECIS, CLIP_DEFAULT_PRECIS,
DEFAULT_QUALITY, DEFAULT_PITCH | FF_SWISS, L"Courier New");
SendMessage(hCombo1, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hCombo2, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hEditMultiline, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"0");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"1");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"2");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"3");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"4");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"5");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"6");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"7");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"8");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"9");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"10");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"0");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"1");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"2");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"3");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"4");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"5");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"6");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"7");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"8");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"9");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"10");
SendMessage(hCombo1, CB_SETCURSEL, 0, 0);
SendMessage(hCombo2, CB_SETCURSEL, 0, 0);
//ShowWindow(hWnd, SW_SHOWMAXIMIZED);
buffer[0] = L'\0';
break;
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDC_COMPUTE_BUTTON:
{
GetWindowText(hCombo1, line, 128);
std::wstring nStr(line);
n = stoi(nStr);
GetWindowText(hCombo2, line, 128);
std::wstring mStr(line);
m = stoi(mStr);
double x0 = -10.0, x1 = 10.0;
double integ1 = Integration::SimpsonsRule(1024, x0, x1, Dx);
double integ2 = Integration::SimpsonsRule(1024, x0, x1, Ex);
double integ3 = Integration::SimpsonsRule(1024, x0, x1, Fx);
double integ4 = Integration::SimpsonsRule(1024, x0, x1, Gx);
swprintf_s(line, L"n = %d\tm = %d\t(x^1) = %+lf\r\n", n, m, integ1);
wcscat_s(buffer, line);
swprintf_s(line, L"n = %d\tm = %d\t(x^2) = %+lf\r\n", n, m, integ2);
wcscat_s(buffer, line);
swprintf_s(line, L"n = %d\tm = %d\t(x^3) = %+lf\r\n", n, m, integ3);
wcscat_s(buffer, line);
swprintf_s(line, L"n = %d\tm = %d\t(x^4) = %+lf\r\n", n, m, integ4);
wcscat_s(buffer, line);
SetWindowText(hEditMultiline, buffer);
break;
}
case IDC_CLEAR_BUTTON:
buffer[0] = '\0';
SetWindowText(hEditMultiline, buffer);
break;
case IDC_CANCEL_BUTTON:
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
break;
case WM_PAINT:
{
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
// TODO: Add any drawing code that uses hdc here...
EndPaint(hWnd, &ps);
}
break;
case WM_DESTROY:
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return 0;
}
// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
Blog Entry © Thursday, October 30, 2025, by James Pate Williams, Jr. Quantum Mechanical Harmonic Oscillator See Introduction to Quantum Mechanics with Applications to Chemistry by Linus Pauling and E. Bright Wilson, Jr. Chapter III
// QMHarmonicOscillator.cpp : Defines the entry point for the application.
//
#include "pch.h"
#include "framework.h"
#include "QMHarmonicOscillator.h"
#include "Integration.h"
#define MAX_LOADSTRING 100
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
WCHAR buffer[128];
RECT rect;
bool draw = false;
int n;
std::wstring fTitle, yTitle;
std::vector<double> xi, fx;
// Forward declarations of functions included in this code module:
ATOM MyRegisterClass(HINSTANCE hInstance);
BOOL InitInstance(HINSTANCE, int);
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK About(HWND, UINT, WPARAM, LPARAM);
int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
_In_opt_ HINSTANCE hPrevInstance,
_In_ LPWSTR lpCmdLine,
_In_ int nCmdShow)
{
UNREFERENCED_PARAMETER(hPrevInstance);
UNREFERENCED_PARAMETER(lpCmdLine);
// TODO: Place code here.
// Initialize global strings
LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
LoadStringW(hInstance, IDC_QMHARMONICOSCILLATOR, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_QMHARMONICOSCILLATOR));
MSG msg;
// Main message loop:
while (GetMessage(&msg, nullptr, 0, 0))
{
if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
}
return (int) msg.wParam;
}
//
// FUNCTION: MyRegisterClass()
//
// PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
WNDCLASSEXW wcex = { };
wcex.cbSize = sizeof(WNDCLASSEX);
wcex.style = CS_HREDRAW | CS_VREDRAW;
wcex.lpfnWndProc = WndProc;
wcex.cbClsExtra = 0;
wcex.cbWndExtra = 0;
wcex.hInstance = hInstance;
wcex.hIcon = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_QMHARMONICOSCILLATOR));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_QMHARMONICOSCILLATOR);
wcex.lpszClassName = szWindowClass;
wcex.hIconSm = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));
return RegisterClassExW(&wcex);
}
//
// FUNCTION: InitInstance(HINSTANCE, int)
//
// PURPOSE: Saves instance handle and creates main window
//
// COMMENTS:
//
// In this function, we save the instance handle in a global variable and
// create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
hInst = hInstance; // Store instance handle in our global variable
HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);
if (!hWnd)
{
return FALSE;
}
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
return TRUE;
}
#define IDC_COMPUTE_BUTTON 1010
#define IDC_CANCEL_BUTTON 1020
#define IDC_CLEAR_BUTTON 1030
static double Factorial(int n)
{
// n! = n ... 2 1
double factorial = 1.0;
for (int i = 2; i <= n; i++)
{
factorial *= i;
}
return factorial;
}
static double N(int n)
{
// normilization factor
double pi = 4.0 * atan(1.0);
double factor1 = sqrt(1.0 / pi);
double factor2 = 1.0 / (pow(2.0, n) * Factorial(n));
return sqrt(factor1 * factor2);
}
static double H(double xi, int n)
{
// Hermite polynomial
if (n == 1)
{
return 2.0 * xi;
}
else if (n == 2)
{
return 4.0 * xi * xi - 2.0;
}
else if (n == 3)
{
return xi * (8.0 * xi * xi - 12.0);
}
else if (n == 4)
{
return 16.0 * pow(xi, 4.0) - 48.0 * pow(xi, 2.0) + 12.0;
}
else if (n == 5)
{
return 32.0 * pow(xi, 5.0) - 160.0 * pow(xi, 3.0) +
120.0 * xi;
}
else if (n == 6)
{
return 64.0 * pow(xi, 6.0) - 480.0 * pow(xi, 4.0) +
720.0 * xi * xi - 120.0;
}
else if (n == 7)
{
return 128.0 * pow(xi, 7.0) - 1344.0 * pow(xi, 5.0) +
3360.0 * xi * xi * xi - 1680.0 * xi;
}
else if (n == 8)
{
return 256.0 * pow(xi, 8.0) - 3584.0 * pow(xi, 6.0) +
13440.0 * pow(xi, 4.0) - 13440.0 * xi * xi + 1680.0;
}
else if (n == 9)
{
return 512.0 * pow(xi, 9.0) - 9216.0 * pow(xi, 7.0) +
48384.0 * pow(xi, 5.0) - 80640.0 * xi * xi * xi + 30240.0 * xi;
}
else if (n == 10)
{
return 1024.0 * pow(xi, 10.0) - 23040.0 * pow(xi, 8.0) +
161280.0 * pow(xi, 6.0) - 403200.0 * pow(xi, 4.0) +
302400.0 * xi * xi - 30240.0;
}
else
{
return 0.0;
}
}
static double I(double xi, int n)
{
// first derivative of H(x) is H'(x)
if (n == 1)
{
return 2.0;
}
else if (n == 2)
{
return 8.0 * xi;
}
else if (n == 3)
{
return 24.0 * xi * xi - 12.0;
}
else if (n == 4)
{
return 64.0 * pow(xi, 3.0) - 96.0 * xi;
}
else if (n == 5)
{
return 32.0 * 5.0 * pow(xi, 4.0) - 3.0 * 160.0 * pow(xi, 2.0) +
120.0;
}
else if (n == 6)
{
return 6.0 * 64.0 * pow(xi, 5.0) - 4.0 * 480.0 * pow(xi, 3.0) +
2.0 * 720.0 * xi;
}
else if (n == 7)
{
return 7.0 * 128.0 * pow(xi, 6.0) - 5.0 * 1344.0 * pow(xi, 4.0) +
3.0 * 3360.0 * xi * xi - 1680.0;
}
else if (n == 8)
{
return 8.0 * 256.0 * pow(xi, 7.0) - 6.0 * 13584.0 * pow(xi, 5.0) +
4.0 * 13440.0 * pow(xi, 3.0) - 2.0 * 13440.0 * xi;
}
else if (n == 9)
{
return 9.0 * 512.0 * pow(xi, 8.0) - 7.0 * 9216.0 * pow(xi, 6.0) +
5.0 * 48384.0 * pow(xi, 4.0) - 3.0 * 80640.0 * xi * xi + 30240.0;
}
else if (n == 10)
{
return 10.0 * 1024.0 * pow(xi, 9.0) - 8.0 * 23040.0 * pow(xi, 7.0) +
6.0 * 161280.0 * pow(xi, 5.0) - 4.0 * 403200.0 * pow(xi, 3.0) +
2.0 * 302400.0 * xi;
}
else
{
return 0;
}
}
static double J(double xi, int n)
{
// second derivative of H(x) is H''(x)
if (n == 2)
{
return 8.0;
}
else if (n == 3)
{
return 48.0 * xi * xi;
}
else if (n == 4)
{
return 3.0 * 64.0 * pow(xi, 2.0) - 96.0;
}
else if (n == 5)
{
return 4.0 * 32.0 * 5.0 * pow(xi, 3.0) - 2.0 * 3.0 * 160.0 * xi;
}
else if (n == 6)
{
return 5.0 * 6.0 * 64.0 * pow(xi, 4.0) - 3.0 * 4.0 * 480.0 * pow(xi, 2.0) +
2.0 * 720.0;
}
else if (n == 7)
{
return 6.0 * 7.0 * 128.0 * pow(xi, 5.0) - 4.0 * 5.0 * 1344.0 * pow(xi, 3.0) +
2.0 * 3.0 * 3360.0 * xi;
}
else if (n == 8)
{
return 7.0 * 8.0 * 256.0 * pow(xi, 6.0) - 5.0 * 6.0 * 13584.0 * pow(xi, 4.0) +
3.0 * 4.0 * 13440.0 * pow(xi, 2.0) - 2.0 * 13440.0;
}
else if (n == 9)
{
return 8.0 * 9.0 * 512.0 * pow(xi, 7.0) - 6.0 * 7.0 * 9216.0 * pow(xi, 5.0) +
4.0 * 5.0 * 48384.0 * pow(xi, 3.0) - 2.0 * 3.0 * 80640.0 * xi;
}
else if (n == 10)
{
return 9.0 * 10.0 * 1024.0 * pow(xi, 8.0) - 7.0 * 8.0 * 23040.0 * pow(xi, 6.0) +
5.0 * 6.0 * 161280.0 * pow(xi, 4.0) - 3.0 * 4.0 * 403200.0 * pow(xi, 2.0) +
2.0 * 302400.0;
}
else
{
return 0;
}
}
static double Psi(double xi, int n)
{
return N(n) * exp(-xi * xi / 2.0) * H(xi, n);
}
static void FindMinMax(
std::vector<double> x,
double& min,
double& max)
{
min = x[0];
max = x[0];
for (int i = 1; i < x.size(); i++)
{
if (x[i] < min)
min = x[i];
if (x[i] > max)
max = x[i];
}
}
static void DrawTitles(
HDC hdc, RECT clientRect, const std::wstring& fTitle,
std::wstring& yTitle, int sx0, int sx1, int sy0, int sy1)
{
HFONT hCustomFont = CreateFont(
-MulDiv(10, GetDeviceCaps(hdc, LOGPIXELSY), 72), // Height in logical units
0, // Width (0 = default)
0, // Escapement
0, // Orientation
FW_BOLD, // Weight (700 = bold)
FALSE, // Italic
FALSE, // Underline
FALSE, // StrikeOut
DEFAULT_CHARSET, // Charset
OUT_DEFAULT_PRECIS, // Output precision
CLIP_DEFAULT_PRECIS, // Clipping precision
DEFAULT_QUALITY, // Quality
FIXED_PITCH | FF_MODERN,// Pitch and family
L"Courier New" // Typeface name
);
HFONT hOldFont = (HFONT)SelectObject(hdc, hCustomFont);
SIZE sz;
int width = clientRect.right - clientRect.left;
// Title: fTitle
GetTextExtentPoint32W(hdc, fTitle.c_str(), (int)fTitle.length(), &sz);
int w = sz.cx;
int h = sz.cy;
TextOutW(hdc, (width - w) / 2, h, fTitle.c_str(), (int)fTitle.length());
// x-axis title: "x"
std::wstring xTitle = L"xi";
GetTextExtentPoint32W(hdc, xTitle.c_str(), (int)xTitle.length(), &sz);
w = sz.cx;
TextOutW(hdc, sx0 + (sx1 - sx0 - w) / 2, sy1 + 2 * h, xTitle.c_str(), (int)xTitle.length());
// y-axis title: "?"
GetTextExtentPoint32W(hdc, yTitle.c_str(), (int)yTitle.length(), &sz);
w = sz.cx;
TextOutW(hdc, sx1 + w / 5, sy0 + (sy1 - sy0) / 2 - h / 2, yTitle.c_str(), (int)yTitle.length());
SelectObject(hdc, hOldFont); // Restore original font
}
static void DrawFormattedText(HDC hdc, char loctext[], RECT rect)
{
// Draw the text with formatting options
DrawTextA(hdc, loctext, (int)strlen(loctext), &rect, DT_SINGLELINE | DT_NOCLIP);
}
static double Fx1(double xi)
{
double psi = Psi(xi, n);
return xi * psi * psi;
}
static double Fx2(double xi)
{
double psi = Psi(xi, n);
return xi * xi * psi * psi;
}
static double Fx3(double xi)
{
double psi = Psi(xi, n);
return xi * xi * xi * psi * psi;
}
static double Fx4(double xi)
{
double psi = Psi(xi, n);
return pow(xi, 4.0) * psi * psi;
}
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
static HFONT hFont = nullptr;
static HWND hCombo1 = nullptr;
static HWND hCombo2 = nullptr;
static HWND hCombo3 = nullptr;
switch (message)
{
case WM_CREATE:
CreateWindowEx(0, L"STATIC", L"draw:", WS_CHILD | WS_VISIBLE,
10, 10, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
hCombo1 = CreateWindowEx(0, L"COMBOBOX", nullptr,
WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
ES_AUTOVSCROLL, 120, 10, 120, 100, hWnd, nullptr, hInst, nullptr);
CreateWindowEx(0, L"STATIC", L"average:", WS_CHILD | WS_VISIBLE,
10, 40, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
hCombo2 = CreateWindowEx(0, L"COMBOBOX", nullptr,
WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
ES_AUTOVSCROLL, 120, 40, 120, 150, hWnd, nullptr, hInst, nullptr);
CreateWindowEx(0, L"STATIC", L"n:", WS_CHILD | WS_VISIBLE,
10, 70, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
hCombo3 = CreateWindowEx(0, L"COMBOBOX", nullptr,
WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
ES_AUTOVSCROLL, 120, 70, 120, 150, hWnd, nullptr, hInst, nullptr);
CreateWindowEx(0, L"BUTTON", L"Compute", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 100, 80, 30, hWnd, (HMENU)IDC_COMPUTE_BUTTON, hInst, NULL);
CreateWindowEx(0, L"BUTTON", L"Cancel", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 140, 80, 30, hWnd, (HMENU)IDC_CANCEL_BUTTON, hInst, NULL);
CreateWindowEx(0, L"BUTTON", L"Clear", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 180, 80, 30, hWnd, (HMENU)IDC_CLEAR_BUTTON, hInst, NULL);
hFont = CreateFont(16, 0, 0, 0, FW_BOLD, FALSE, FALSE, FALSE,
UNICODE, OUT_DEFAULT_PRECIS, CLIP_DEFAULT_PRECIS,
DEFAULT_QUALITY, DEFAULT_PITCH | FF_SWISS, L"Courier New");
SendMessage(hCombo1, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hCombo2, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hCombo3, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"Hermite");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"Psi");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"|Psi^2|");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"zeros");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"x");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"x^2");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"x^3");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"x^4");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"1");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"2");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"3");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"4");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"5");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"6");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"7");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"8");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"9");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"10");
SendMessage(hCombo1, CB_SETCURSEL, 0, 0);
SendMessage(hCombo2, CB_SETCURSEL, 0, 0);
SendMessage(hCombo3, CB_SETCURSEL, 0, 0);
ShowWindow(hWnd, SW_SHOWMAXIMIZED);
buffer[0] = L'\0';
break;
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDC_COMPUTE_BUTTON:
{
if (hCombo3)
{
LRESULT result1 = SendMessage(hCombo1, CB_GETCURSEL, 0, 0);
LRESULT result2 = SendMessage(hCombo2, CB_GETCURSEL, 0, 0);
GetWindowText(hCombo3, buffer, 128);
std::wstring nStr(buffer);
n = stoi(nStr);
double x0 = -10.0, x1 = 10.0, h = (x1 - x0) / 256;
double xValue = -10.0, yValue = 0;
xi.clear();
fx.clear();
if (result1 == 0)
{
wcscpy_s(buffer, L"xi Versus Hermite Polynomial H(xi)");
fTitle = std::wstring(buffer);
wcscpy_s(buffer, L"H(xi, n)");
yTitle = std::wstring(buffer);
while (xValue <= x1)
{
xi.push_back(xValue);
yValue = H(xValue, n);
fx.push_back(yValue);
xValue += h;
}
}
else if (result1 == 1)
{
wcscpy_s(buffer, L"xi Versus Wave Function Psi(xi, n)");
fTitle = std::wstring(buffer);
wcscpy_s(buffer, L"Psi(xi, n)");
yTitle = std::wstring(buffer);
while (xValue <= x1)
{
xi.push_back(xValue);
yValue = Psi(xValue, n);
fx.push_back(yValue);
xValue += h;
}
}
else if (result1 == 2)
{
wcscpy_s(buffer, L"xi Versus |Psi^2|");
fTitle = std::wstring(buffer);
wcscpy_s(buffer, L"|Psi(xi, n)^2|");
yTitle = std::wstring(buffer);
while (xValue <= x1)
{
double psi = Psi(xValue, n);
xi.push_back(xValue);
yValue = psi * psi;
fx.push_back(yValue);
xValue += h;
}
}
GetClientRect(hWnd, &rect);
rect.left = 300;
draw = true;
InvalidateRect(hWnd, &rect, TRUE);
double integ = 0.0;
if (result2 == 0)
{
integ = Integration::SimpsonsRule(
1024, -10.0, 10.0, Fx1);
swprintf_s(buffer, L"xi * Psi^2(xi, n) = %lf", integ);
}
if (result2 == 1)
{
integ = Integration::SimpsonsRule(
1024, -10.0, 10.0, Fx2);
swprintf_s(buffer, L"xi^2 * Psi^2(xi, n) = %lf", integ);
double integ1
}
if (result2 == 2)
{
integ = Integration::SimpsonsRule(
1024, -10.0, 10.0, Fx3);
swprintf_s(buffer, L"xi^3 * Psi^2(xi, n) = %lf", integ);
}
if (result2 == 3)
{
integ = Integration::SimpsonsRule(
1024, -10.0, 10.0, Fx4);
swprintf_s(buffer, L"xi^4 * Psi^2(xi, n) = %lf", integ);
}
MessageBox(hWnd, buffer, L"Information",
MB_OK | MB_ICONINFORMATION);
}
break;
}
case IDC_CLEAR_BUTTON:
buffer[0] = '\0';
draw = false;
break;
case IDC_CANCEL_BUTTON:
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
case WM_PAINT:
{
if (draw)
{
PAINTSTRUCT ps;
POINT wPt = { };
HDC hdc = BeginPaint(hWnd, &ps);
double xMin = 0, yMin = 0;
double xMax = 0, yMax = 0;
FindMinMax(xi, xMin, xMax);
FindMinMax(fx, yMin, yMax);
double h = 0, pi = 0, plm = 0, theta = 0;
float xSpan = (float)(xMax - xMin);
float ySpan = (float)(yMax - yMin);
float width = (float)(rect.right - rect.left + 1);
float height = (float)(rect.bottom - rect.top - 32 + 1);
float sx0 = 2.0f * width / 16.0f + 300;
float sx1 = 14.0f * width / 16.0f + 300;
float sy0 = 2.0f * height / 16.0f;
float sy1 = 14.0f * height / 16.0f;
float deltaX = xSpan / 8.0f;
float deltaY = ySpan / 8.0f;
float xSlope = (sx1 - sx0) / xSpan;
float xInter = (float)(sx0 - xSlope * xMin);
float ySlope = (sy0 - sy1) / ySpan;
float yInter = (float)(sy0 - ySlope * yMax);
float px = 0, py = 0, sx = 0, sy = 0;
int i = 0;
float x = (float)xMin;
float y = (float)yMax;
DrawTitles(
hdc, rect, fTitle, yTitle,
(int)sx0, (int)sx1, (int)sy0, (int)sy1);
px = (float)xi[0];
py = (float)fx[0];
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
while (i <= 8)
{
sx = xSlope * x + xInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
LineTo(hdc, (int)sx, (int)sy1);
char numberStr[128] = "";
sprintf_s(numberStr, 128, "%5.4lf", x);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
numberStr,
(int)strlen(numberStr),
&size);
RECT textRect = { };
textRect.left = (long)(sx - size.cx / 2.0f);
textRect.right = (long)(sx + size.cx / 2.0f);
textRect.top = (long)sy1;
textRect.bottom = (long)(sy1 + size.cy / 2.0f);
DrawFormattedText(hdc, numberStr, textRect);
x += deltaX;
i++;
}
i = 0;
y = (float)yMin;
while (i <= 8)
{
sy = ySlope * y + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
LineTo(hdc, (int)sx, (int)sy);
if (i != 0)
{
char numberStr[128] = "";
sprintf_s(numberStr, 128, "%+5.3lf", y);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
numberStr,
(int)strlen(numberStr),
&size);
RECT textRect = { 0 };
textRect.left = (long)(sx0 - size.cx - size.cx / 2.0f);
textRect.right = (long)(sx0 - size.cx / 2.0f);
textRect.top = (long)(sy - size.cy / 2.0f);
textRect.bottom = (long)(sy + size.cy / 2.0f);
DrawFormattedText(hdc, numberStr, textRect);
}
y += deltaY;
i++;
}
HGDIOBJ bPenNew = NULL;
HGDIOBJ hPenOld = NULL;
bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, bPenNew);
HRGN clipRegion = CreateRectRgn(
(int)sx0, (int)sy0, // Left, Top
(int)(sx1), (int)(sy1) // Right, Bottom
);
// Apply clipping region
SelectClipRgn(hdc, clipRegion);
px = (float)xi[0];
py = (float)fx[0];
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy, &wPt);
for (size_t j = 1; j < (int)xi.size(); j++)
{
px = (float)xi[j];
py = (float)fx[j];
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
DeleteObject(bPenNew);
bPenNew = NULL;
SelectObject(hdc, hPenOld);
DeleteObject(clipRegion);
EndPaint(hWnd, &ps);
}
}
break;
case WM_DESTROY:
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return 0;
}
// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
#pragma once
class Integration
{
public:
static double TrapezoidalRule(
int n, double a, double b, double(*fx)(double));
static double SimpsonsRule(
int n, double a, double b, double(*fx)(double));
static void GenerateGaussianLegendreAbscissasAndWeights(
int n, double x[], double w[]);
static double GLIntegrate11(
int n, double x[], double w[], double(*fx)(double));
static double GLIntegrateAB(int n, double a, double b,
double x[], double w[], double(*fx)(double));
};
#include "pch.h"
#include "Integration.h"
#include "OrthogonalPolynomials.h"
#include "Zeros.h"
double Integration::TrapezoidalRule(
int n, double a, double b, double(*fx)(double))
{
double h = (b - a) / n;
double s = 0.5 * (fx(a) + fx(b));
double x = a + h;
for (int i = 1; i < n; i++)
{
s += fx(x);
x += h;
}
return h * s;
}
double Integration::SimpsonsRule(
int n, double a, double b, double(*fx)(double))
{
double h = (b - a) / n;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int i = 1; i < n; i += 2)
{
s += fx(x);
x += h2;
}
x = a + h2;
for (int i = 2; i < n; i += 2)
{
t += fx(x);
x += h2;
}
return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}
void Integration::GenerateGaussianLegendreAbscissasAndWeights(
int n, double x[], double w[])
{
char name[32] = { 0 };
OrthogonalPolynomials op;
std::vector<double> roots;
Zeros::Zero("Legendre", 0.0, 0.0, n, roots);
strcpy_s(name, sizeof(name), "Legendre");
for (int i = 0; i < n; i++)
{
double xi = roots[i];
double fd = op.g(name, xi, n);
double x2 = 1.0 - xi * xi;
x[i] = xi;
w[i] = 2.0 / (x2 * fd * fd);
}
}
double Integration::GLIntegrate11(
int n, double x[], double w[], double(*fx)(double))
{
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += w[i] * fx(x[i]);
return sum;
}
double Integration::GLIntegrateAB(int n, double a, double b,
double x[], double w[], double(*fx)(double))
{
double c = (b - a) / 2.0;
double d = (b + a) / 2.0;
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += w[i] * fx(c * x[i] + d);
return c * sum;
}
Blog Entry © Wednesday, October 29, 2025, by James Pate Williams, Jr. Hydrogenic Atomic Spectral Lines
// HydrogenSpectralLines.cpp : Defines the entry point for the application.
//
#include "pch.h"
#include "framework.h"
#include "HydrogenSpectralLines.h"
#define MAX_LOADSTRING 100
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
WCHAR buffer[32768], line[128], line1[128], line2[512];
double c = 299792458;
double h = 4.135667696e-15;
// Forward declarations of functions included in this code module:
ATOM MyRegisterClass(HINSTANCE hInstance);
BOOL InitInstance(HINSTANCE, int);
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK About(HWND, UINT, WPARAM, LPARAM);
int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
_In_opt_ HINSTANCE hPrevInstance,
_In_ LPWSTR lpCmdLine,
_In_ int nCmdShow)
{
UNREFERENCED_PARAMETER(hPrevInstance);
UNREFERENCED_PARAMETER(lpCmdLine);
// TODO: Place code here.
// Initialize global strings
LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
LoadStringW(hInstance, IDC_HYDROGENSPECTRALLINES, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_HYDROGENSPECTRALLINES));
MSG msg;
// Main message loop:
while (GetMessage(&msg, nullptr, 0, 0))
{
if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
}
return (int) msg.wParam;
}
//
// FUNCTION: MyRegisterClass()
//
// PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
WNDCLASSEXW wcex = { };
wcex.cbSize = sizeof(WNDCLASSEX);
wcex.style = CS_HREDRAW | CS_VREDRAW;
wcex.lpfnWndProc = WndProc;
wcex.cbClsExtra = 0;
wcex.cbWndExtra = 0;
wcex.hInstance = hInstance;
wcex.hIcon = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_HYDROGENSPECTRALLINES));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_HYDROGENSPECTRALLINES);
wcex.lpszClassName = szWindowClass;
wcex.hIconSm = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));
return RegisterClassExW(&wcex);
}
//
// FUNCTION: InitInstance(HINSTANCE, int)
//
// PURPOSE: Saves instance handle and creates main window
//
// COMMENTS:
//
// In this function, we save the instance handle in a global variable and
// create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
hInst = hInstance; // Store instance handle in our global variable
HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);
if (!hWnd)
{
return FALSE;
}
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
return TRUE;
}
#define IDC_COMPUTE_BUTTON 1010
#define IDC_CANCEL_BUTTON 1020
#define IDC_CLEAR_BUTTON 1030
#define IDC_EDIT_MULTILINE 1040
//
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
static HFONT hFont = nullptr;
static HWND hCombo1 = nullptr;
static HWND hCombo2 = nullptr;
static HWND hCombo3 = nullptr;
static HWND hCombo4 = nullptr;
static HWND hEditMultiline = nullptr;
switch (message)
{
case WM_CREATE:
CreateWindowEx(0, L"STATIC", L"Series:", WS_CHILD | WS_VISIBLE,
10, 10, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
hCombo1 = CreateWindowEx(0, L"COMBOBOX", nullptr,
WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
ES_AUTOVSCROLL, 120, 10, 200, 100, hWnd, nullptr, hInst, nullptr);
CreateWindowEx(0, L"STATIC", L"n2:", WS_CHILD | WS_VISIBLE,
10, 40, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
hCombo2 = CreateWindowEx(0, L"COMBOBOX", nullptr,
WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
ES_AUTOVSCROLL, 120, 40, 120, 150, hWnd, nullptr, hInst, nullptr);
CreateWindowEx(0, L"STATIC", L"atom:", WS_CHILD | WS_VISIBLE,
10, 70, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
hCombo3 = CreateWindowEx(0, L"COMBOBOX", nullptr,
WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
ES_AUTOVSCROLL, 120, 70, 120, 150, hWnd, nullptr, hInst, nullptr);
CreateWindowEx(0, L"STATIC", L"formula:", WS_CHILD | WS_VISIBLE,
10, 100, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
hCombo4 = CreateWindowEx(0, L"COMBOBOX", nullptr,
WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
ES_AUTOVSCROLL, 120, 100, 120, 150, hWnd, nullptr, hInst, nullptr);
CreateWindowEx(0, L"BUTTON", L"Compute", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 140, 80, 30, hWnd, (HMENU)IDC_COMPUTE_BUTTON, hInst, NULL);
CreateWindowEx(0, L"BUTTON", L"Cancel", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 180, 80, 30, hWnd, (HMENU)IDC_CANCEL_BUTTON, hInst, NULL);
CreateWindowEx(0, L"BUTTON", L"Clear", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 220, 80, 30, hWnd, (HMENU)IDC_CLEAR_BUTTON, hInst, NULL);
hFont = CreateFont(16, 0, 0, 0, FW_BOLD, FALSE, FALSE, FALSE,
UNICODE, OUT_DEFAULT_PRECIS, CLIP_DEFAULT_PRECIS,
DEFAULT_QUALITY, DEFAULT_PITCH | FF_SWISS, L"Courier New");
hEditMultiline = CreateWindowEx(
WS_EX_CLIENTEDGE, // Extended style for sunken border
TEXT("EDIT"), // Class name
TEXT(""), // Initial text (can be blank)
WS_CHILD | WS_VISIBLE | WS_VSCROLL | ES_AUTOHSCROLL |
ES_LEFT | ES_MULTILINE | ES_AUTOVSCROLL | WS_HSCROLL | WS_VSCROLL,
340, 10, 640, 400, // Position and size
hWnd, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
SendMessage(hCombo1, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hCombo2, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hCombo3, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hCombo4, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hEditMultiline, WM_SETFONT, (WPARAM)hFont, 0);
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"Lyman n1 = 1 n2 >= 2");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"Balmer n1 = 2 n2 >= 3");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"Paschen n1 = 3 n2 >= 4");
SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"Brackett n1 = 4 n2 >= 5");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"2");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"3");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"4");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"5");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"6");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"7");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"8");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"9");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"10");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"11");
SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"12");
SendMessage(hCombo1, CB_SETCURSEL, 0, 0);
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"hydrogen");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"helium+");
SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"lithium++");
SendMessage(hCombo4, CB_ADDSTRING, 0, (LPARAM)L"energy eV");
SendMessage(hCombo4, CB_ADDSTRING, 0, (LPARAM)L"frequency Hz");
SendMessage(hCombo4, CB_ADDSTRING, 0, (LPARAM)L"wavelength A");
SendMessage(hCombo1, CB_SETCURSEL, 0, 0);
SendMessage(hCombo2, CB_SETCURSEL, 0, 0);
SendMessage(hCombo3, CB_SETCURSEL, 0, 0);
SendMessage(hCombo4, CB_SETCURSEL, 0, 0);
buffer[0] = L'\0';
break;
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDC_COMPUTE_BUTTON:
{
int n1 = 0, n2 = 0;
WCHAR series[16] = { 0 };
if (hCombo1)
{
GetWindowText(hCombo1, line, 128);
}
if (std::wcsstr(line, L"Lyman"))
{
wcscpy_s(series, L"Lyman");
n1 = 1;
}
else if (std::wcsstr(line, L"Balmer"))
{
wcscpy_s(series, L"Balmer");
n1 = 2;
}
else if (std::wcsstr(line, L"Paschen"))
{
wcscpy_s(series, L"Paschen");
n1 = 3;
}
else if (std::wcsstr(line, L"Brackett"))
{
wcscpy_s(series, L"Brackett");
n1 = 4;
}
if (hCombo2)
{
GetWindowText(hCombo2, line, 128);
std::wstring n2Str(line);
n2 = stoi(n2Str);
}
if (n2 < n1 + 1)
{
MessageBox(hWnd, L"n2 must be greater than equal n1 + 1",
L"Warning", MB_OK | MB_ICONWARNING);
break;
}
double Z = 1.0;
GetWindowText(hCombo3, line, 128);
if (wcscmp(line, L"hydrogen") == 0)
{
Z = 1.0;
}
else if (wcscmp(line, L"helium+") == 0)
{
Z = 2.0;
}
else if (wcscmp(line, L"lithium++") == 0)
{
Z = 3.0;
}
double deltaE = -13.6 * Z * Z * (1.0 / (n1 * n1) - 1.0 / (n2 * n2));
double nu = deltaE / h;
double lambda = c / nu / 1.0e-10;
double number = deltaE;
GetWindowText(hCombo4, line1, 128);
if (wcscmp(line1, L"energy eV") == 0)
{
number = deltaE;
}
else if (wcscmp(line1, L"frequency Hz") == 0)
{
number = nu;
}
else if (wcscmp(line1, L"wavelength A") == 0)
{
number = lambda;
}
swprintf_s(line2, L"%s\t%s\t%s\t%d\t%d\t%lf\r\n",
series, line, line1, n1, n2, number);
wcscat_s(buffer, line2);
if (hEditMultiline)
{
SetWindowText(hEditMultiline, buffer);
}
break;
}
case IDC_CLEAR_BUTTON:
buffer[0] = '\0';
SetWindowText(hEditMultiline, buffer);
break;
case IDC_CANCEL_BUTTON:
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
break;
case WM_PAINT:
{
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
// TODO: Add any drawing code that uses hdc here...
EndPaint(hWnd, &ps);
}
break;
case WM_DESTROY:
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return 0;
}
// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
Blog Entry © Wednesday, September 3, 2025, by James Pate Williams, Jr. Solution of Exercise 1.11 in Modern Quantum Chemistry an Introduction to Advanced Electronic Structure Theory © 1996 by Attila Szabo and Neil S. Ostlund
// EigenVV2x2.cpp (c) Wednesday, September 3, 2025
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Solution to Exercise 1.11 in "Quantum Chemistry
// an Introduction to Advanced Electronic Structure
// Theory (c) 1996 by Attila Szabo and Neil S. Ostlund
#include "pch.h"
#include "framework.h"
#include "EigenVV2x2.h"
#define MAX_LOADSTRING 100
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
std::wstring text; // output wide string
// Forward declarations of functions included in this code module:
ATOM MyRegisterClass(HINSTANCE hInstance);
BOOL InitInstance(HINSTANCE, int);
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK About(HWND, UINT, WPARAM, LPARAM);
int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
_In_opt_ HINSTANCE hPrevInstance,
_In_ LPWSTR lpCmdLine,
_In_ int nCmdShow)
{
UNREFERENCED_PARAMETER(hPrevInstance);
UNREFERENCED_PARAMETER(lpCmdLine);
// TODO: Place code here.
// Initialize global strings
LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
LoadStringW(hInstance, IDC_EIGENVV2X2, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_EIGENVV2X2));
MSG msg;
// Main message loop:
while (GetMessage(&msg, nullptr, 0, 0))
{
if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
}
return (int) msg.wParam;
}
//
// FUNCTION: MyRegisterClass()
//
// PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
WNDCLASSEXW wcex;
wcex.cbSize = sizeof(WNDCLASSEX);
wcex.style = CS_HREDRAW | CS_VREDRAW;
wcex.lpfnWndProc = WndProc;
wcex.cbClsExtra = 0;
wcex.cbWndExtra = 0;
wcex.hInstance = hInstance;
wcex.hIcon = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_EIGENVV2X2));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_EIGENVV2X2);
wcex.lpszClassName = szWindowClass;
wcex.hIconSm = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));
return RegisterClassExW(&wcex);
}
//
// FUNCTION: InitInstance(HINSTANCE, int)
//
// PURPOSE: Saves instance handle and creates main window
//
// COMMENTS:
//
// In this function, we save the instance handle in a global variable and
// create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
hInst = hInstance; // Store instance handle in our global variable
HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);
if (!hWnd)
{
return FALSE;
}
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
return TRUE;
}
#define IDC_STATIC1 1000
#define IDC_STATIC2 1010
#define IDC_STATIC3 1020
#define IDC_STATIC4 1030
#define IDC_EDIT_A11 2000
#define IDC_EDIT_A12 2010
#define IDC_EDIT_A21 2020
#define IDC_EDIT_A22 2030
#define IDC_EDIT_MULTILINE 3000
#define IDC_BUTTON_COMPUTE 4000
#define IDC_BUTTON_CANCEL 4010
//
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
static HFONT hFont = { };
static HWND hEditMultiline = { };
static HWND hEditA11 = { };
static HWND hEditA12 = { };
static HWND hEditA21 = { };
static HWND hEditA22 = { };
switch (message)
{
case WM_CREATE:
{
hFont = CreateFont(16, 0, 0, 0, FW_NORMAL, FALSE, FALSE, FALSE,
ANSI_CHARSET, OUT_DEFAULT_PRECIS, CLIP_DEFAULT_PRECIS,
DEFAULT_QUALITY, DEFAULT_PITCH | FF_SWISS, L"Courier New Bold");
CreateWindowEx(0, L"STATIC", L"A[1][1]:", WS_CHILD | WS_VISIBLE,
10, 10, 80, 20, hWnd, (HMENU)IDC_STATIC1, hInst, NULL);
CreateWindowEx(0, L"STATIC", L"A[1][2]:", WS_CHILD | WS_VISIBLE,
10, 40, 80, 20, hWnd, (HMENU)IDC_STATIC2, hInst, NULL);
CreateWindowEx(0, L"STATIC", L"A[2][1]:", WS_CHILD | WS_VISIBLE,
10, 70, 80, 20, hWnd, (HMENU)IDC_STATIC3, hInst, NULL);
CreateWindowEx(0, L"STATIC", L"A[2][2]:", WS_CHILD | WS_VISIBLE,
10, 100, 80, 20, hWnd, (HMENU)IDC_STATIC4, hInst, NULL);
hEditA11 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
100, 10, 200, 20, hWnd, (HMENU)IDC_EDIT_A11, hInst, NULL);
hEditA12 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
100, 40, 200, 20, hWnd, (HMENU)IDC_EDIT_A12, hInst, NULL);
hEditA21 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
100, 70, 200, 20, hWnd, (HMENU)IDC_EDIT_A21, hInst, NULL);
hEditA22 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
100, 100, 200, 20, hWnd, (HMENU)IDC_EDIT_A22, hInst, NULL);
hEditMultiline = CreateWindowEx(
WS_EX_CLIENTEDGE, // Extended style for sunken border
TEXT("EDIT"), // Class name
TEXT(""), // Initial text (can be blank)
WS_CHILD | WS_VISIBLE | WS_VSCROLL | ES_AUTOHSCROLL |
ES_LEFT | ES_MULTILINE | ES_AUTOVSCROLL | WS_HSCROLL | WS_VSCROLL,
310, 10, 300, 300, // Position and size
hWnd, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
CreateWindowEx(0, L"BUTTON", L"Compute", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
10, 130, 80, 30, hWnd, (HMENU)IDC_BUTTON_COMPUTE, hInst, NULL);
CreateWindowEx(0, L"BUTTON", L"Cancel", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
220, 130, 80, 30, hWnd, (HMENU)IDC_BUTTON_CANCEL, hInst, NULL);
SendMessage(hEditMultiline, WM_SETFONT, (WPARAM)hFont, TRUE);
SetDlgItemText(hWnd, IDC_EDIT_A11, L"3");
SetDlgItemText(hWnd, IDC_EDIT_A12, L"1");
SetDlgItemText(hWnd, IDC_EDIT_A21, L"1");
SetDlgItemText(hWnd, IDC_EDIT_A22, L"3");
}
break;
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDC_BUTTON_COMPUTE:
{
WCHAR line[128] = L"";
text = L"";
// eigenvalues
std::vector<double> omega(3);
// matrix
std::vector<std::vector<double>> A(3,
std::vector<double>(3));
std::vector<std::vector<double>> B(3,
std::vector<double>(3));
// eigenvectors
std::vector<std::vector<double>> c(3,
std::vector<double>(3));
GetWindowText(hEditA11, line, 128);
std::wstring A11Str(line);
A[1][1] = std::stod(A11Str);
GetWindowText(hEditA12, line, 128);
std::wstring A12Str(line);
A[1][2] = std::stod(A12Str);
GetWindowText(hEditA21, line, 128);
std::wstring A21Str(line);
A[2][1] = std::stod(A21Str);
GetWindowText(hEditA22, line, 128);
std::wstring A22Str(line);
A[2][2] = std::stod(A22Str);
double term1 = A[1][1] + A[2][2];
double term2 = pow(A[2][2] - A[1][1], 2.0);
double term3 = 4.0 * A[1][2] * A[2][1];
// compute eigenvalues
omega[1] = 0.5 * (term1 - sqrt(term2 + term3));
omega[2] = 0.5 * (term1 + sqrt(term2 + term3));
swprintf_s(line, L"Eigenvalues:\r\n\r\n");
text += std::wstring(line);
swprintf_s(line, L"omega 1 = %13.10lf\r\n", omega[1]);
text += std::wstring(line);
swprintf_s(line, L"omega 2 = %13.10lf\r\n\r\n", omega[2]);
text += std::wstring(line);
// compute eigenvalues using a unitary transformation
// matrix A must be symmetric
double theta0 = 0.0;
if (A[1][1] == A[2][2])
{
theta0 = 0.5 * acos(0.0);
}
else
{
theta0 = 0.5 * atan(2.0 * A[1][2] /
(A[1][1] - A[2][2]));
}
// compute the eigenvalues
omega[1] = A[1][1] * pow(cos(theta0), 2.0) +
A[2][2] * pow(sin(theta0), 2.0) +
A[1][2] * sin(2.0 * theta0);
omega[2] = A[1][1] * pow(cos(theta0), 2.0) +
A[2][2] * pow(sin(theta0), 2.0) -
A[1][2] * sin(2.0 * theta0);
swprintf_s(line, L"Eigenvalues:\r\n\r\n");
text += std::wstring(line);
swprintf_s(line, L"omega 1 = %13.10lf\r\n", omega[1]);
text += std::wstring(line);
swprintf_s(line, L"omega 2 = %13.10lf\r\n\r\n", omega[2]);
text += std::wstring(line);
// compute eigenvectors
c[1][1] = cos(theta0);
c[1][2] = sin(theta0);
c[2][1] = sin(theta0);
c[2][2] = -cos(theta0);
swprintf_s(line, L"Eigenvectors:\r\n\r\n");
text += std::wstring(line);
swprintf_s(line, L"c 11 = %13.10lf\r\n", c[1][1]);
text += std::wstring(line);
swprintf_s(line, L"c 12 = %13.10lf\r\n", c[1][2]);
text += std::wstring(line);
swprintf_s(line, L"c 21 = %13.10lf\r\n", c[2][1]);
text += std::wstring(line);
swprintf_s(line, L"c 22 = %13.10lf\r\n", c[2][2]);
text += std::wstring(line);
SetWindowText(hEditMultiline, text.c_str());
break;
}
case IDC_BUTTON_CANCEL:
{
PostQuitMessage(0);
break;
}
break;
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
break;
case WM_PAINT:
{
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
// TODO: Add any drawing code that uses hdc here...
EndPaint(hWnd, &ps);
}
break;
case WM_DESTROY:
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return 0;
}
// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
Blog Entry © Tuesday, August 26, 2025, by James Pate Williams, Jr. and the Microsoft Copilot, Hydrogen-like Radial Wavefunction, Its First Derivative, and Its Probability Distribution Function
Blog Entry © Sunday, August 24, 2025, by James Pate Williams, Jr. Corrections to New Quantum Chemical Total Molecular Ground-State Energies for the Helium Hydride Cation (a Hetero Nuclear molecule) and the Hydrogen Molecule (a Homo Nuclear Molecule)
Blog Entry © Friday, August 22, 2025, by James Pate Williams, Jr. New Quantum Chemical Total Molecular Ground-State Energies for the Helium Hydride Cation (a Hetero Nuclear molecule) and the Hydrogen Molecule (a Homo Nuclear Molecule)
Blog Entry © Tuesday, August 19, 2025, by James Pate Williams, Jr., Continuation of Answers to the Exercises in Chapter 1 of Modern Quantum Chemistry an Introduction to Advanced Electronic Structure Theory by Attila Szabo and Neil S. Ostlund
Note: Later on, Tuesday, August 19, 2025, I added five C++ source code files.
#include <vector>
#include <random>
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(
int n, int row, int col,
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);
static void DblGenerateNonSingular(
double scale, double& determinant,
int n, unsigned int seed,
std::vector<std::vector<double>>& Mr);
};
#include "DblLinearAlgebra.h"
#include <iomanip>
#include <iostream>
void DblLinearAlgebra::DblPrintMatrix(
int m, int n, std::vector<std::vector<double>>& A)
{
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
std::cout << std::setprecision(6) << std::setw(9);
if (fabs(A[i][j]) > 1.0e-12)
{
std::cout << A[i][j] << ' ';
}
else
{
std::cout << 0 << ' ';
}
}
std::cout << std::endl;
}
}
void DblLinearAlgebra::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)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
C[i][j] = A[i][j] + B[i][j];
}
}
}
void DblLinearAlgebra::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)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
C[i][j] = A[i][j] - B[i][j];
}
}
}
void DblLinearAlgebra::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)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
double sum = { 0 };
for (size_t k = 0; k < p; k++)
{
sum += A[i][k] * B[k][j];
}
C[i][j] = sum;
}
}
}
void DblLinearAlgebra::DblAnticommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C)
{
std::vector<std::vector<double>> D(n,
std::vector<double>(n));
std::vector<std::vector<double>> E(n,
std::vector<double>(n));
DblMultiply(n, n, n, A, B, D);
DblMultiply(n, n, n, B, A, E);
DblAddition(n, n, D, E, C);
}
void DblLinearAlgebra::DblCommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C)
{
std::vector<std::vector<double>> D(n,
std::vector<double>(n));
std::vector<std::vector<double>> E(n,
std::vector<double>(n));
DblMultiply(n, n, n, A, B, D);
DblMultiply(n, n, n, B, A, E);
DblSubtraction(n, n, D, E, C);
}
// https://www.geeksforgeeks.org/dsa/determinant-of-a-matrix/
double getDet(std::vector<std::vector<double>>& mat, int n) {
// Base case: if the matrix is 1x1
if (n == 1) {
return mat[0][0];
}
// Base case for 2x2 matrix
if (n == 2) {
return mat[0][0] * mat[1][1] -
mat[0][1] * mat[1][0];
}
// Recursive case for larger matrices
double res = 0;
for (int col = 0; col < n; ++col) {
// Create a submatrix by removing the first
// row and the current column
std::vector<std::vector<double>> sub(n - 1,
std::vector<double>(n - 1));
for (int i = 1; i < n; ++i) {
int subcol = 0;
for (int j = 0; j < n; ++j) {
// Skip the current column
if (j == col) continue;
// Fill the submatrix
sub[i - 1LL][subcol++] = mat[i][j];
}
}
// Cofactor expansion
int sign = (col % 2 == 0) ? 1 : -1;
res += sign * mat[0][col] * getDet(sub, n - 1);
}
return res;
}
double DblLinearAlgebra::DblDeterminant(
int n, int row, int col,
std::vector<std::vector<double>>& A)
{
return getDet(A, A.size());
}
bool DblLinearAlgebra::DblGaussianElimination(
int m, int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot)
{
bool failure = false;
std::vector<double> c(m);
b.resize(n);
x.resize(n);
for (size_t i = 0; i < m; i++)
c[i] = -1;
for (size_t j = 0; j < n; j++)
{
bool found = false;
size_t i = j;
while (i < n && !found)
{
if (M[i][j] != 0)
found = true;
else
i++;
}
if (!found)
{
failure = true;
break;
}
if (i > j)
{
for (size_t l = j; l < n; l++)
{
double t = M[i][l];
M[i][l] = M[j][l];
M[j][l] = t;
t = b[i];
b[i] = b[j];
b[j] = t;
}
}
double d = 1.0 / M[j][j];
for (size_t k = j + 1; k < n; k++)
c[k] = d * M[k][j];
for (size_t k = j + 1; k < n; k++)
{
for (size_t l = j + 1; l < n; l++)
M[k][l] = M[k][l] - c[k] * M[j][l];
b[k] = b[k] - c[k] * b[j];
}
}
for (long long i = (long long)n - 1; i >= 0; i--)
{
double sum = 0;
for (size_t j = i + 1; j < n; j++)
sum += M[i][j] * x[j];
x[i] = (b[i] - sum) / M[i][i];
}
return failure;
}
bool DblLinearAlgebra::DblGaussianFactor(
int n, std::vector<std::vector<double>>& M,
std::vector<size_t>& pivot)
{
// returns false if matrix is singular
std::vector<double> d(n);
double awikod, col_max, ratio, row_max, temp;
int flag = 1;
size_t i_star, itemp;
for (size_t i = 0; i < n; i++)
{
pivot[i] = i;
row_max = 0;
for (size_t j = 0; j < n; j++)
row_max = fmax(row_max, fabs(M[i][j]));
if (row_max == 0)
{
flag = 0;
row_max = 1.0;
}
d[i] = row_max;
}
if (n <= 1) return flag != 0;
// factorization
for (size_t k = 0; k < (size_t)(n - 1LL); k++)
{
// determine pivot row the row i_star
col_max = fabs(M[k][k]) / d[k];
i_star = k;
for (size_t i = k + 1; i < n; i++)
{
awikod = fabs(M[i][k]) / d[i];
if (awikod > col_max)
{
col_max = awikod;
i_star = i;
}
}
if (col_max == 0)
flag = 0;
else
{
if (i_star > k)
{
// make k the pivot row by
// interchanging with i_star
flag *= -1;
itemp = pivot[i_star];
pivot[i_star] = pivot[k];
pivot[k] = itemp;
temp = d[i_star];
d[i_star] = d[k];
d[k] = temp;
for (size_t j = 0; j < n; j++)
{
temp = M[i_star][j];
M[i_star][j] = M[k][j];
M[k][j] = temp;
}
}
// eliminate x[k]
for (size_t i = k + 1; i < n; i++)
{
M[i][k] /= M[k][k];
ratio = M[i][k];
for (size_t j = k + 1; j < n; j++)
M[i][j] -= ratio * M[k][j];
}
}
if (M[n - 1LL][n - 1LL] == 0) flag = 0;
}
if (flag == 0)
return false;
return true;
}
bool DblLinearAlgebra::DblGaussianSolution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot)
{
if (!DblGaussianFactor(n, M, pivot))
return false;
return DblSubstitution(n, M, b, x, pivot);
}
bool DblLinearAlgebra::DblSubstitution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot)
{
double sum = 0.0;
size_t n1 = n - 1LL;
if (n == 1)
{
x[0] = b[0] / M[0][0];
return true;
}
// forward substitution
x[0] = b[pivot[0]];
for (size_t i = 1; i < n; i++)
{
double sum = 0.0;
for (size_t j = 0; j < i; j++)
sum += M[i][j] * x[j];
x[i] = b[pivot[i]] - sum;
}
// backward substitution
x[n1] /= M[n1][n1];
for (long long i = n - 2LL; i >= 0; i--)
{
double sum = 0.0;
for (size_t j = i + 1; j < n; j++)
sum += M[i][j] * x[j];
x[i] = (x[i] - sum) / M[i][i];
}
return true;
}
bool DblLinearAlgebra::DblInverse(
int n, std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& Mi)
{
std::vector<double> b(n);
std::vector<double> x(n);
std::vector<size_t> pivot(n);
std::vector<std::vector<double>> Mc(n,
std::vector<double>(n));
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
Mc[i][j] = M[i][j];
}
}
if (!DblGaussianFactor(n, Mc, pivot))
return false;
for (size_t i = 0; i < n; i++)
{
b[i] = 0;
}
for (size_t i = 0; i < n; i++)
{
b[i] = 1;
if (!DblSubstitution(n, Mc, b, x, pivot))
return false;
b[i] = 0;
for (size_t j = 0; j < n; j++)
Mi[j][i] = x[j];
}
return true;
}
void DblLinearAlgebra::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)
{
C.resize(n, std::vector<double>(n));
I.resize(n, std::vector<double>(n));
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
C[i][j] = I[i][j] = 0;
}
for (size_t i = 0; i < n; i++)
C[i][i] = I[i][i] = 1;
a[0] = 1;
for (size_t i = 1; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
for (size_t k = 0; k < n; k++)
{
double sum = 0.0;
for (size_t l = 0; l < n; l++)
sum += M[j][l] * C[l][k];
C[j][k] = sum;
}
}
double tr = 0.0;
for (size_t j = 0; j < n; j++)
tr += C[j][j];
a[i] = -tr / i;
for (size_t j = 0; j < n; j++)
{
for (size_t k = 0; k < n; k++)
C[j][k] += a[i] * I[j][k];
}
}
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
double sum = 0.0;
for (size_t k = 0; k < n; k++)
sum += M[i][k] * C[k][j];
C[i][j] = sum;
}
}
double trace = 0.0;
for (size_t i = 0; i < n; i++)
trace += C[i][i];
trace /= n;
a[n - 1LL] = -trace;
double factor = 1.0;
if ((n - 1) % 2 != 0)
factor = -1.0;
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
adjoint[i][j] = factor * C[i][j];
}
}
void DblLinearAlgebra::DblMatrixKernel(
int m, int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X,
size_t& r)
{
double D = 0.0;
std::vector <int> c(m);
std::vector <int> d(n);
r = 0;
for (size_t i = 0; i < m; i++)
c[i] = -1;
size_t j, k = 1;
Step2:
for (j = 0; j < m; j++)
{
if (M[j][k] != 0 && c[j] == -1)
break;
}
if (j == m)
{
r++;
d[k] = 0;
goto Step4;
}
D = -1.0 / M[j][k];
M[j][k] = -1;
for (size_t s = k + 1; s < n; s++)
{
M[j][s] = D * M[j][s];
for (size_t i = 0; i < m; i++)
{
if (i != j)
{
D = M[i][k];
M[i][k] = 0;
}
}
}
for (size_t s = k + 1; s < n; s++)
{
for (size_t i = 0; i < m; i++)
{
M[i][s] += D * M[j][s];
}
}
c[j] = (int)k;
d[k] = (int)j;
Step4:
if (k < n - 1)
{
k++;
goto Step2;
}
X.resize(n, std::vector<double>(n));
if (r != 0)
{
for (k = 0; k < n; k++)
{
if (d[k] == 0)
{
for (size_t i = 0; i < n; i++)
{
if (d[i] > 0)
X[k][i] = M[d[i]][k];
else if (i == k)
X[k][i] = 1;
else
X[k][i] = 0;
}
}
}
}
}
void DblLinearAlgebra::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)
{
double D = 0.0;
size_t r = 0;
std::vector<std::vector<double>> copyM(
m, std::vector<double>(n));
std::vector <int> c(m);
std::vector <int> d(n);
for (size_t i = 0; i < m; i++)
c[i] = -1;
size_t j = 0, k = 1;
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
N[i][j] = copyM[i][j] = M[i][j];
}
}
Step2:
for (size_t j = 0; j < m; j++)
{
if (copyM[j][k] != 0 && c[j] == -1)
break;
}
if (j == m)
{
r++;
d[k] = 0;
goto Step4;
}
D = -1.0 / copyM[j][k];
copyM[j][k] = -1;
for (size_t s = k + 1; s < n; s++)
{
copyM[j][s] = D * copyM[j][s];
for (size_t i = 0; i < m; i++)
{
if (i != j)
{
D = copyM[i][k];
copyM[i][k] = 0;
}
}
}
for (size_t s = k + 1; s < n; s++)
{
for (size_t i = 0; i < m; i++)
{
copyM[i][s] += D * copyM[j][s];
}
}
c[j] = (int)k;
d[k] = (int)j;
Step4:
if (k < (size_t)(n - 1LL))
{
k++;
goto Step2;
}
rank = (int)(n - r) ;
for (j = 0; j < m; j++)
{
if (c[j] != 0)
{
for (size_t i = 0; i < m; i++)
{
N[i][c[j]] = M[i][c[j]];
}
}
}
}
void DblLinearAlgebra::DblGenerateNonSingular(
double scale, double& determinant,
int n, unsigned int seed,
std::vector<std::vector<double>>& Mr)
{
bool failure = false;
std::mt19937 rng(seed);
std::uniform_real_distribution<double> dist(0.0, 1.0);
while (true)
{
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
Mr[i][j] = scale * dist(rng);
}
}
determinant = DblDeterminant(n, 0, 0, Mr);
failure = determinant == 0;
if (!failure)
return;
}
}
#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 std::complex<double> CmpDeterminant(
int n,
std::vector<std::vector<std::complex<double>>>& Ac);
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 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);
static void CmpGenerateNonSingular(
double scale, std::complex<double>& determinant,
int n, unsigned int seed,
std::vector<std::vector<std::complex<double>>>& Mc);
};
#include "CmpLinearAlgebra.h"
#include <iomanip>
#include <iostream>
#include <random>
void CmpLinearAlgebra::CmpPrintMatrix(
int m, int n,
std::vector<std::vector<std::complex<double>>>& Ac)
{
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
if (Ac[i][j]._Val[0] >= 0)
std::cout << '+';
else
std::cout << '-';
std::cout << std::setprecision(6) << std::setw(9);
if (fabs(Ac[i][j]._Val[0]) > 1.0e-12)
{
std::cout << fabs(Ac[i][j]._Val[0]) << ' ';
}
else
{
std::cout << 0 << ' ';
}
if (Ac[i][j]._Val[1] >= 0)
std::cout << '+';
else
std::cout << '-';
std::cout << std::setprecision(6) << std::setw(9);
if (fabs(Ac[i][j]._Val[1]) > 1.0e-12)
{
std::cout << fabs(Ac[i][j]._Val[1]) << "i\t";
}
else
{
std::cout << 0 << "i\t";
}
}
std::cout << std::endl;
}
}
void CmpLinearAlgebra::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)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
C[i][j] = A[i][j] + B[i][j];
}
}
}
void CmpLinearAlgebra::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)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
C[i][j] = A[i][j] - B[i][j];
}
}
}
void CmpLinearAlgebra::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)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
std::complex<double> sum = 0;
for (size_t k = 0; k < p; k++)
{
sum += A[i][k] * B[k][j];
}
C[i][j] = sum;
}
}
}
void CmpLinearAlgebra::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)
{
std::vector<std::vector<std::complex<double>>> D(n,
std::vector<std::complex<double>>(n));
std::vector<std::vector<std::complex<double>>> E(n,
std::vector<std::complex<double>>(n));
CmpMultiply(n, n, n, A, B, D);
CmpMultiply(n, n, n, B, A, E);
CmpAddition(n, n, D, E, C);
}
void CmpLinearAlgebra::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)
{
std::vector<std::vector<std::complex<double>>> D(n,
std::vector<std::complex<double>>(n));
std::vector<std::vector<std::complex<double>>> E(n,
std::vector<std::complex<double>>(n));
CmpMultiply(n, n, n, A, B, D);
CmpMultiply(n, n, n, B, A, E);
CmpSubtraction(n, n, D, E, C);
}
void CmpLinearAlgebra::CmpAdjoint(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& Ac,
std::vector<std::vector<std::complex<double>>>& Ad)
{
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
Ad[j][i] = std::conj(Ac[i][j]);
}
}
}
// https://www.geeksforgeeks.org/dsa/determinant-of-a-matrix/
std::complex<double> getDet(
std::vector<std::vector<std::complex<double>>>& mat, int n) {
// Base case: if the matrix is 1x1
if (n == 1) {
return mat[0][0];
}
// Base case for 2x2 matrix
if (n == 2) {
return mat[0][0] * mat[1][1] -
mat[0][1] * mat[1][0];
}
// Recursive case for larger matrices
std::complex<double> res = 0;
for (int col = 0; col < n; ++col) {
// Create a submatrix by removing the first
// row and the current column
std::vector<std::vector<std::complex<double>>> sub(n - 1,
std::vector<std::complex<double>>(n - 1));
for (int i = 1; i < n; ++i) {
int subcol = 0;
for (int j = 0; j < n; ++j) {
// Skip the current column
if (j == col) continue;
// Fill the submatrix
sub[i - 1LL][subcol++] = mat[i][j];
}
}
// Cofactor expansion
int sign = (col % 2 == 0) ? 1 : -1;
std::complex<double> csign(sign, 0.0);
res = res + csign * mat[0][col] * getDet(sub, n - 1);
}
return res;
}
std::complex<double> CmpLinearAlgebra::CmpDeterminant(
int n, std::vector<std::vector<std::complex<double>>>& A)
{
return getDet(A, A.size());
}
bool CmpLinearAlgebra::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)
{
bool failure = false;
std::vector<std::complex<double>> c(m);
b.resize(n);
x.resize(n);
for (size_t i = 0; i < m; i++)
c[i] = -1;
for (size_t j = 0; j < n; j++)
{
bool found = false;
size_t i = j;
while (i < n && !found)
{
if (abs(A[i][j]) != 0)
found = true;
else
i++;
}
if (!found)
{
failure = true;
break;
}
if (i > j)
{
for (size_t l = j; l < n; l++)
{
std::complex<double> t = A[i][l];
A[i][l] = A[j][l];
A[j][l] = t;
t = b[i];
b[i] = b[j];
b[j] = t;
}
}
std::complex<double> d = 1.0 / A[j][j];
for (size_t k = j + 1; k < n; k++)
c[k] = d * A[k][j];
for (size_t k = j + 1; k < n; k++)
{
for (size_t l = j + 1; l < n; l++)
A[k][l] = A[k][l] - c[k] * A[j][l];
b[k] = b[k] - c[k] * b[j];
}
}
for (long long i = (long long)n - 1; i >= 0; i--)
{
std::complex<double> sum = 0;
for (size_t j = i + 1; j < n; j++)
sum += A[i][j] * x[j];
x[i] = (b[i] - sum) / A[i][i];
}
return failure;
}
bool CmpLinearAlgebra::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)
{
std::complex<double> sum = 0;
size_t n1 = n - 1LL;
if (n == 1)
{
x[0] = b[0] / M[0][0];
return true;
}
// forward substitution
x[0] = b[pivot[0]];
for (size_t i = 1; i < n; i++)
{
std::complex<double> sum = 0;
for (size_t j = 0; j < i; j++)
sum += M[i][j] * x[j];
x[i] = b[pivot[i]] - sum;
}
// backward substitution
x[n1] /= M[n1][n1];
for (long long i = n - 2LL; i >= 0; i--)
{
std::complex<double> sum = 0;
for (size_t j = i + 1; j < n; j++)
sum += M[i][j] * x[j];
x[i] = (x[i] - sum) / M[i][i];
}
return true;
}
static std::complex<double> complex_max(
std::complex<double> a, std::complex<double> b) {
return (std::abs(a) > std::abs(b)) ? a : b;
}
bool CmpLinearAlgebra::CmpGaussianFactor(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<size_t>& pivot)
{
// returns false if matrix is singular
std::vector<std::complex<double>> d(n);
std::complex<double> awikod = 0, col_max = 0, ratio = 0, row_max = 0, temp = 0;
int flag = 1;
size_t i_star, itemp;
for (size_t i = 0; i < n; i++)
{
pivot[i] = i;
row_max = 0;
for (size_t j = 0; j < n; j++)
row_max = complex_max(row_max, abs(M[i][j]));
if (abs(row_max) == 0)
{
flag = 0;
row_max = 1;
}
d[i] = row_max;
}
if (n <= 1) return flag != 0;
// factorization
for (size_t k = 0; k < (size_t)n - 1LL; k++)
{
// determine pivot row the row i_star
col_max = abs(M[k][k]) / d[k];
i_star = k;
for (size_t i = k + 1; i < n; i++)
{
awikod = abs(M[i][k]) / d[i];
if (abs(awikod) > abs(col_max))
{
col_max = awikod;
i_star = i;
}
}
if (abs(col_max) == 0)
flag = 0;
else
{
if (i_star > k)
{
// make k the pivot row by
// interchanging with i_star
flag *= -1;
itemp = pivot[i_star];
pivot[i_star] = pivot[k];
pivot[k] = itemp;
temp = d[i_star];
d[i_star] = d[k];
d[k] = temp;
for (size_t j = 0; j < n; j++)
{
temp = M[i_star][j];
M[i_star][j] = M[k][j];
M[k][j] = temp;
}
}
// eliminate x[k]
for (size_t i = k + 1; i < n; i++)
{
M[i][k] /= M[k][k];
ratio = M[i][k];
for (size_t j = k + 1; j < n; j++)
M[i][j] -= ratio * M[k][j];
}
}
if (abs(M[n - 1LL][n - 1LL]) == 0) flag = 0;
}
if (flag == 0)
return false;
return true;
}
bool CmpLinearAlgebra::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)
{
if (!CmpGaussianFactor(n, M, pivot))
return false;
return CmpSubstitution(n, n, M, b, x, pivot);
}
bool CmpLinearAlgebra::CmpInverse(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& Mi)
{
std::vector<std::complex<double>> b(n);
std::vector<std::complex<double>> x(n);
std::vector<size_t> pivot(n);
std::vector<std::vector<std::complex<double>>> Mc(n,
std::vector<std::complex<double>>(n));
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
Mc[i][j] = M[i][j];
}
}
if (!CmpGaussianFactor(n, Mc, pivot))
return false;
for (size_t i = 0; i < n; i++)
{
b[i] = 0;
}
for (size_t i = 0; i < n; i++)
{
b[i] = 1;
if (!CmpSubstitution(n, n, Mc, b, x, pivot))
return false;
b[i] = 0;
for (size_t j = 0; j < n; j++)
Mi[j][i] = x[j];
}
return true;
}
void CmpLinearAlgebra::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)
{
C.resize(n, std::vector<std::complex<double>>(n));
I.resize(n, std::vector<std::complex<double>>(n));
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
C[i][j] = I[i][j] = 0;
}
for (size_t i = 0; i < n; i++)
C[i][i] = I[i][i] = 1;
a[0] = 1;
for (size_t i = 1; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
for (size_t k = 0; k < n; k++)
{
std::complex<double> sum = 0.0;
for (size_t l = 0; l < n; l++)
sum += M[j][l] * C[l][k];
C[j][k] = sum;
}
}
std::complex<double> tr = 0.0;
for (size_t j = 0; j < n; j++)
tr += C[j][j];
std::complex<double> ci = 0;
ci._Val[0] = (double)i;
a[i] = -tr / ci;
for (size_t j = 0; j < n; j++)
{
for (size_t k = 0; k < n; k++)
C[j][k] += a[i] * I[j][k];
}
}
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
{
std::complex<double> sum = 0.0;
for (size_t k = 0; k < n; k++)
sum += M[i][k] * C[k][j];
C[i][j] = sum;
}
}
std::complex<double> trace = 0.0;
for (size_t i = 0; i < n; i++)
trace += C[i][i];
trace /= n;
a[n - 1LL] = -trace;
std::complex<double> factor = 1.0;
if ((n - 1) % 2 != 0)
factor = -1.0;
for (size_t i = 0; i < n; i++)
{
for (size_t j = 0; j < n; j++)
adjoint[i][j] = factor * C[i][j];
}
}
void CmpLinearAlgebra::CmpMatrixKernel(
int m, int n,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& X,
size_t& r)
{
std::complex<double> D = 0;
std::vector<int> c(m);
std::vector<int> d(n);
r = 0;
for (size_t i = 0; i < m; i++)
c[i] = -1;
size_t j = 0, k = 1;
Step2:
for (j = 0; j < m; j++)
{
if (abs(M[j][k]) != 0 && c[j] == -1)
break;
}
if (j == m)
{
r++;
d[k] = 0;
goto Step4;
}
D = -1.0 / M[j][k];
M[j][k] = -1;
for (size_t s = k + 1; s < n; s++)
{
M[j][s] = D * M[j][s];
for (size_t i = 0; i < m; i++)
{
if (i != j)
{
D = M[i][k];
M[i][k] = 0;
}
}
}
for (size_t s = k + 1; s < n; s++)
{
for (size_t i = 0; i < m; i++)
{
M[i][s] += D * M[j][s];
}
}
c[j] = (int)k;
d[k] = (int)j;
Step4:
if (k < n - 1)
{
k++;
goto Step2;
}
X.resize(n, std::vector<std::complex<double>>(n));
if (r != 0)
{
for (k = 0; k < n; k++)
{
if (d[k] == 0)
{
for (size_t i = 0; i < n; i++)
{
if (d[i] > 0)
X[k][i] = M[d[i]][k];
else if (i == k)
X[k][i] = 1;
else
X[k][i] = 0;
}
}
}
}
}
void CmpLinearAlgebra::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)
{
std::complex<double> D = 0.0;
size_t r = 0;
std::vector<std::vector<std::complex<double>>> copyM(
m, std::vector<std::complex<double>>(n));
std::vector<int> c(m);
std::vector<int> d(n);
for (size_t i = 0; i < m; i++)
c[i] = -1;
size_t j = 0, k = 1;
for (size_t i = 0; i < m; i++)
{
for (size_t j = 0; j < n; j++)
{
N[i][j] = copyM[i][j] = M[i][j];
}
}
Step2:
for (size_t j = 0; j < m; j++)
{
if (abs(copyM[j][k]) != 0 && c[j] == -1)
break;
}
if (j == m)
{
r++;
d[k] = 0;
goto Step4;
}
D = -1.0 / copyM[j][k];
copyM[j][k] = -1;
for (size_t s = k + 1; s < n; s++)
{
copyM[j][s] = D * copyM[j][s];
for (size_t i = 0; i < m; i++)
{
if (i != j)
{
D = copyM[i][k];
copyM[i][k] = 0;
}
}
}
for (size_t s = k + 1; s < n; s++)
{
for (size_t i = 0; i < m; i++)
{
copyM[i][s] += D * copyM[j][s];
}
}
c[j] = (int)k;
d[k] = (int)j;
Step4:
if (k < (size_t)(n - 1LL))
{
k++;
goto Step2;
}
rank = (int)(n - r);
for (j = 0; j < m; j++)
{
if (c[j] != 0)
{
for (size_t i = 0; i < m; i++)
{
N[i][c[j]] = M[i][c[j]];
}
}
}
}
void CmpLinearAlgebra::CmpGenerateNonSingular(
double scale, std::complex<double>& cDeterminant,
int n, unsigned int seed,
std::vector<std::vector<std::complex<double>>>& Mc)
{
bool failure = false;
std::mt19937 rng(seed);
std::uniform_real_distribution<double> dist(0.0, 1.0);
while (true)
{
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
Mc[i][j]._Val[0] = scale * dist(rng);
Mc[i][j]._Val[1] = scale * dist(rng);
}
}
cDeterminant = CmpDeterminant(n, Mc);
if (cDeterminant._Val[0] != 0 || cDeterminant._Val[1] != 0)
break;
}
}
// 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.
// Program (c) Tuesday, August 19, 2025 by James Pate Williams, Jr.
#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include "DblLinearAlgebra.h"
#include "CmpLinearAlgebra.h"
int main()
{
// static data matrices
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 } };
// some array dimensions
int m = 3, n = 3, p = 3;
// a couple of 3x1 vectors
std::vector<double> br(3);
std::vector<size_t> pivot(3);
// 3x3 real matrices
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));
// a 4x4 real matrix
std::vector<std::vector<double>> Mr(4, std::vector<double>(4));
// 3x3 complex matrices
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));
std::vector<std::vector<std::complex<double>>> Mc(4,
std::vector<std::complex<double>>(4));
// copy static real matrices to dynamic matrices
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];
}
}
// copy static complex matrices to dynamic matrices
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];
}
}
// See "Elementary Numerical Analysis an
// Algorithmic Approach" (c) 1980 by S. D. Conte
// and Carl de Boor
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;
// complex matrix multiplication
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 from Szabo and Ostlund
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 flag = "
<< 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;
double rDeterminant = 0;
DblLinearAlgebra::DblGenerateNonSingular(
2.0, rDeterminant, 4, 1, Mr);
std::cout << "rDeterminant = ";
std::cout << rDeterminant << std::endl;
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 flag = " << 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);
std::complex<double> cDeterminant = 0;
CmpLinearAlgebra::CmpGenerateNonSingular(
2.0, cDeterminant, 4, 1, Mc);
std::cout << std::endl;
std::cout << "complex determinant = ";
std::cout << cDeterminant << std::endl;
double rDeterminantA = 0;
std::vector<std::vector<double>> A44r(4,
std::vector<double>(4));
DblLinearAlgebra::DblGenerateNonSingular(
2.0, rDeterminantA, 4, 2, A44r);
double rDeterminantB = 0;
std::vector<std::vector<double>> B44r(4,
std::vector<double>(4));
DblLinearAlgebra::DblGenerateNonSingular(
2.0, rDeterminantB, 4, 3, B44r);
std::cout << std::endl;
std::vector<std::vector<double>> C44r(4,
std::vector<double>(4));
DblLinearAlgebra::DblMultiply(4, 4, 4, A44r, B44r, C44r);
std::cout << "|A| = " << rDeterminantA << std::endl;
std::cout << "|B| = " << rDeterminantB << std::endl;
bool failure = false;
double rDeterminantC =
DblLinearAlgebra::DblDeterminant(4, 0, 0, C44r);
std::cout << "|AB| = " << rDeterminantC << std::endl;
std::cout << "|A||B| = " << rDeterminantA *
rDeterminantB << std::endl;
// Exercise 1.6 with 4x4 complex determinants
std::vector<std::vector<std::complex<double>>> A44c(4,
std::vector<std::complex<double>>(4));
std::cout << std::endl;
std::complex<double> cDeterminantA = 0;
CmpLinearAlgebra::CmpGenerateNonSingular(
2.0, cDeterminantA, 4, 2, A44c);
std::vector<std::vector<std::complex<double>>> B44c(4,
std::vector<std::complex<double>>(4));
std::complex<double> cDeterminantB = 0;
CmpLinearAlgebra::CmpGenerateNonSingular(
2.0, cDeterminantB, 4, 3, B44c);
std::vector<std::vector<std::complex<double>>> C44c(4,
std::vector<std::complex<double>>(4));
CmpLinearAlgebra::CmpMultiply(4, 4, 4, A44c, B44c, C44c);
std::cout << "|A| = " << cDeterminantA << std::endl;
std::cout << "|B| = " << cDeterminantB << std::endl;
failure = false;
std::complex<double> cDeterminantC =
CmpLinearAlgebra::CmpDeterminant(4, C44c);
std::cout << "|AB| = " << cDeterminantC << std::endl;
std::cout << "|A||B| = " << cDeterminantA *
cDeterminantB << std::endl;
std::cout << "\nEnter any key to halt: ";
char line[128] = "";
std::cin.getline(line, 128);
}
Numerical Explorations 