# Optimal basis set of Gaussian functions for describing a quantum system (part 2)

It's a part 2 of this question Optimal basis set of Gaussian functions for describing a quantum system (part 1)

There, the answer was given to a question related to finding geometric progressions of the coefficients ($$\alpha_i$$ and $$\beta_i$$) of a Gaussian basis set in the general case using minimization. But there was still the second question regarding why the authors of the article use different indexes for the parameters $$N$$ and $$M$$ for $$\alpha$$ and $$\beta$$ respectively.

Brief description of the problem: The system has the following Hamiltonian $$H=-\frac{1}{2}\Delta-\frac{1}{r}+\frac{25}{8}\rho^2-5/2$$, where $$r=(\rho,z,\phi)$$ is a coordinate in the cylindrical system. Physically, this Hamiltonian describes a hydrogen atom in a magnetic field equal to 5 (in dimensionless units).

To find the ground and excited states of the system by the matrix method, one can use the Gaussian set of basis functions $$\psi_j=e^{-b_{j} z^2}e^{-a_{j} \rho^2}$$, where $$a_{j}$$ and $$b_{j}$$ are parameters. This approach is described in the article C. Aldrich & R.L. Greene (https://drive.google.com/file/d/1q74hAn0UAdNd8DtPkoCsdYPkr61-xhr2/view).

The authors use as sets of parameters geometric progressions in the following form:

$$\alpha_j=\alpha_{j-1}(\frac{\alpha_N}{\alpha_1})^{1/(N-1)}$$

$$\beta_j=\alpha_{j-1}(\frac{\beta_M}{\beta_1})^{1/(M-1)}$$,
where $$\alpha_1$$, $$\alpha_N$$, $$\beta_1$$ and $$\beta_M$$ are the first and the last parameters.

Why do authors use different indices for parameters (N and M for $$\alpha$$ and $$\beta$$ respectively)? Aren't they paired in the basis function?

• @Alex Trounev, could you please explain, why did the authors use different indexes for the sequences α and β in the article, N and M respectively? Could you show the code with implementation this idea. Commented May 15, 2023 at 13:48

In the paper "Hydrogen-Like Systems in Arbitrary Magnetic Fields. A Variational Approach" the method to optimize basis functions has been proposed. They described basis functions dependent on 4 quantum numbers in a form

$$\psi_{ijq}=z^q e^{-a_i r^2} e^{-b_j z^2}, \psi_m=r^{|m|}e^{i m\theta}$$

where $$a_i,b_j$$ are parameters to be optimized. The variational wave functions in terms of a set of basis functions can be written as follows $$\Psi_{mq}=\sum_{ij}\psi_{ijq}\psi_m$$

The code to optimize energy is not so differ from what we discussed here:

nmax = 3; jmax = 4;
B = 5; Lz = m = 0; q = 0; ms = -1/2;
VP1[r_, z_] := -1/Sqrt[r^2 + z^2] + 1/8 B^2 r^2 + B/2 (m + 2 ms);

Psi[r_, z_, i_, j_, q_] := z^q Exp[-b[j]*z^2]*Exp[-a[i]*r^2];
Psim[r_, \[Theta]_, m_] := r^Abs[m] Exp[I m \[Theta]];
(*kinetic energy*)
Kk = 1/2 FullSimplify[
Psi[r, z, i2, j2, q] Psim[r, \[Theta], -m]*
Laplacian[
Psi[r, z, i1, j1, q] Psim[r, \[Theta], m], {r, \[Theta], z},
"Cylindrical"] +
Psi[r, z, i1, j1, q] Psim[r, \[Theta], m]*
Laplacian[
Psi[r, z, i2, j2, q] Psim[r, \[Theta], -m], {r, \[Theta], z},
"Cylindrical"]];

ss = 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 = - 1/2 2 Pi Sum[
c[i1, j1] c[i2, j2] ss, {i1, nmax}, {i2, nmax}, {j1, jmax}, {j2,
jmax}];

(*potential energy*)
Px = Integrate[
Psi[r, z, i2, j2, q] Psim[r, \[Theta], -m]*VP1[r, z]*
Psi[r, z, i1, j1, q] Psim[r, \[Theta], m]*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 = 2 Pi Sum[
c[i1, j1] c[i2, j2] Px, {i1, nmax}, {i2, nmax}, {j1, jmax}, {j2,
jmax}];

int = Integrate[
Psi[r, z, i2, j2, q] Psim[r, \[Theta], -m]*
Psi[r, z, i1, j1, q] Psim[r, \[Theta], m] 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 = {2 Pi Sum[
c[i1, j1] c[i2, j2] int, {i1, nmax}, {i2, nmax}, {j1,
jmax}, {j2, jmax}] == 1};

U = 2 Pi ArrayFlatten[
Table[ int, {i1, nmax}, {i2, nmax}, {j1, jmax}, {j2, jmax}], 2];

var = Join[{a1, b1, qa, qb},
Flatten[Table[c[i, j], {i, nmax}, {j, jmax}]]]; a[1] = a1;
b[1] = b1; Do[a[i] = a[i - 1] qa, {i, 2, nmax}]; Do[
b[j] = b[j - 1] qb, {j, 2, jmax}]; con =
Join[norm, {a1 > 0, b1 > 0, qa > 0, qb > 0}];

sol = NMinimize[{Kx + Px, con}, var]

(*Out[]= {-1.37834 + 0. I, {a1 -> 1.39282, b1 -> 19.9993,
qa -> 3.4335, qb -> 0.252195, c[1, 1] -> -0.00902664,
c[1, 2] -> 0.0164657, c[1, 3] -> -0.317789, c[1, 4] -> -0.375104,
c[2, 1] -> 0.0377631, c[2, 2] -> -0.209612, c[2, 3] -> -0.11075,
c[2, 4] -> 0.0505215, c[3, 1] -> -0.158495, c[3, 2] -> 0.0175225,
c[3, 3] -> 0.0210703, c[3, 4] -> -0.0143571}}*)


Please, pay attention that we use 3 a and 4 b to prepare system of $$3\times 4=12$$ basis functions. Eigenvalues can be computed as follows

abq = Take[sol[[2]], 4]

Hij =
Last@CoefficientArrays[(Kx + Px) /. abq,
Flatten[Table[c[i, j], {i, nmax}, {j, jmax}]], "Symmetric" -> True]

mat = Inverse[U /. abq] . Hij; {val, vec} = Eigensystem[mat];

val // Sort // Chop

(*Out[]= {-1.37834, 0.542836, 5.79876, 6.72021, 7.97674, 13.9054, \
32.4115, 34.3706, 35.4432, 40.6693, 41.9055, 67.6762}*)


The ground state is -1.37834, and it is not bad compared to exact value -1.380398866427.

• Yes it is. Also we can't handle large nmax*jmax with NMinimize. Commented May 16, 2023 at 2:31
• Use pts = Reap[ NMinimize[{Kx + Px, con}, var, StepMonitor :> Sow[var]]][[2, 1]]; and ListPlot[Table[pts[[All, i]], {i, 4}], PlotLegends -> Take[var, 4], PlotRange -> All] Commented May 17, 2023 at 3:49
• You can try, but unfortunately, EvaluationMonitor is not working with NMinimize as expected it is why we use StepMonitor (see Options for NMinimize for details). Commented May 17, 2023 at 9:03
• @MamMam You can use it but result is not the same as evaluated with StepMonitor. Commented May 17, 2023 at 11:36
• @MamMam Yes it is. This is well know mathematical problem. You can ask about this problem on math.stackexchange.com Commented May 18, 2023 at 3:16