# Help with Mathematica for optimizing a Gaussian basis set for a hydrogen atom in a magnetic field

This is a follow-up to my previous question: Optimal Gaussian basis set for hydrogen atom in magnetic field

## Brief description of the problem

I would like to find the ground and excited states of the system:

$$H=-\frac{1}{2}\Delta-\frac{1}{r}+\frac{B^2}{8}\rho^2-\frac{B}{2} (m + 2 m_s),\tag{1}$$ in which $$r=(\rho,z,\phi), B=5, m=0, m_s=-1/2$$.

## Method

To find the eigenvalues I would like to use Gaussian functions: $$\psi_{ij}=e^{-b_{j} z^2}e^{-a_{i} \rho^2},\tag{2}$$

in which,

\begin{align} a_{i}&=a_1 qa^{i-1}, i=1, 2 ,3,..., i_{\textrm{max}},\tag{3}\\b_{j}&=b_1 qb^{j-1}, j=1, 2 ,3,..., j_{\textrm{max}},\tag{4} \end{align}

are geometrical progressions.

## Code

I used Mathematica (in the code I rename $$\alpha_i \equiv a_i$$, $$\beta_j \equiv b_j$$, $$\rho\equiv r$$):

ClearAll["Global*"]

imax = 4; jmax = 5;

B = 5; lz = m = 0; ms = -1/2;

VP1[r_, z_] := -1/Sqrt[r^2 + z^2];

Psi[r_, z_, i_, j_] := Exp[-b[j]*z^2]*Exp[-a[i]*r^2];

a[1] = a1;
b[1] = b1;
Do[a[i] = a[i - 1] qa, {i, 2, imax}];
Do[b[j] = b[j - 1] qb, {j, 2, jmax}];

(*kinetic energy*)
Kk = FullSimplify[
Psi[r, z, i2, j2] *
Laplacian[Psi[r, z, i1, j1] , {r, \[Theta], z}, "Cylindrical"]];
Kk1 = -1/2 2 Pi*
Integrate[Kk r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity},
Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0,
a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];
Kx = Table[
Kk1, {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}];
KK[a1_, b1_, qa_, qb_] = Flatten[Kx, {{1, 3}, {2, 4}}];

(*potential energy + A kinetic energy*)
Px1 = 2 Pi*
Integrate[
Psi[r, z, i2, j2] *(VP1[r, z] + 1/8 B^2 r^2 + B/2 (m + 2 ms))*
Psi[r, z, i1, j1]*r, {r, 0, \[Infinity]}, {z, -Infinity,
Infinity},
Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0,
a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];
Px = Table[
Px1, {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}];
PP = Flatten[Px, {{1, 3}, {2, 4}}];

EE = KK + PP;

(*normalization*)
int = 2 Pi*
Integrate[
Psi[r, z, i2, j2] *Psi[r, z, i1, j1]  r, {r,
0, \[Infinity]}, {z, -Infinity, Infinity},
Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0,
a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];
norm = Table[
int, {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}];
BB = Flatten[norm, {{1, 3}, {2, 4}}];

IBB = Inverse[BB](**)

(*the sum of eigenvalues to minimize*)
Tr[EE.IBB]


## Problem

The code cannot handle more than 20 basis functions, but the paper from the 1970s was able to do 100+ basis functions. I will be glad to see any comments. Please use any programming language convenient for you, I will try to figure it out.

• Please provide the error message in a code` block: mattermodeling.meta.stackexchange.com/q/417/5 Commented Jun 1, 2023 at 18:25