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

This question arose during the discussion of the previous questions: Why does the minimum eigenvalue change dramatically when one basis function is added to the basis set? and Finding excited states using the condition of wave functions orthogonality

Brief description of the problem:

Consider the following system with the 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).
In order to solve this problem by the matrix method, it is necessary to choose a set of basis functions. To describe ground and excited states with the z component of the angular momentum m = 0 (1s, 2s, 3d (m = 0), 3s ...), a Gaussian basis set is well suited. Such basis functions have the following form: $$\psi_j=e^{-b_{j} z^2}e^{-a_{j} \rho^2}$$, where $$a_{j}$$ and $$b_{j}$$ are parameters. The task is reduced to the correct finding of parameters.

The parameters describing the ground state (1s) with high accuracy were found in this work (see here https://arxiv.org/abs/1709.05553 or here https://pubs.aip.org/aip/jcp/article/147/24/244108/195534/Accurate-and-balanced-anisotropic-Gaussian-type). But the basis presented in this article does not describe excited states (2s, 3d (m = 0), 3s ...) well, since there are not enough basis functions to describe them. The basis set parameters presented in this article has a complex structure, so increasing the basis set functions in this case is a difficult task.

It seems that a more general way to find the basis set of parameters is presented in this article (see here https://drive.google.com/file/d/1q74hAn0UAdNd8DtPkoCsdYPkr61-xhr2/view?usp=sharing or here https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.2220930140). The sets of parameters that the authors use are geometric progressions:
$$\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. About finding these parameters, the authors write the following: "the parameters $$\alpha_1$$, $$\alpha_N$$, $$\beta_1$$ and $$\beta_M$$ were varied in a four dimensional optimization search which minimized the sum of the eigen-values of interest" and "We have performed the calculations presented in this paper with a basis with ten Gaussians in the $$\rho$$ direction and twelve Gaiissians in the z direction, for a total of 120 different basis functions.". I really don't understand what that means?

Energy values from the articles:
ground state (1s) = - 1.3793 (exact value = -1.380398866427)
the first excited state (2s) = -0.19335 (exact value = -0.193746709717)
the second excited state (3d (m=0)) = -0.07365
the third excited state (3s) = -0.03835

All these states should have the same structure of the basis set function ($$\psi_j=e^{-b_{j} z^2}e^{-a_{j} \rho^2}$$) due to the m=0 (the z component of the angular momentum).

I have a few questions:
I really can't understand what the authors mean in varied optimization search of the parameters $$\alpha_1$$, $$\alpha_N$$, $$\beta_1$$ and $$\beta_M$$, could you please explain to me what that means?

Why do the parameters ($$\alpha_1$$, $$\alpha_N$$, $$\beta_1$$ and $$\beta_M$$) have different indexes, shouldn't their number be the same? Since they are paired in the basis function.

I am well aware that this question may not be very typical, but it should be considered in the context of the previous two questions, since this question is a continuation of the reasoning that is demonstrated there. On the other hand, a possible answer will be demonstrated in Mathematica, which has optimization and minimization tools.

I would like to draw attention to the fact that to get the energy from the Aldrich, C., & Greene, R. L article it's needed the value from the table 1 to subtract the field value and divide by 2 (this is due to the fact that the authors of this article do not take into account the spin). For example for 2s state: value from the table 1 for magnetic field = 5 is equal 4.6133 so the energy of the 2s state = (4.6133 - 5)/2 = -0,19335

• @Alex Trounev, your feedback is needed ) May 11 at 21:47
• If $r=(\rho,z,\phi), what do you means by 1/r May 13 at 14:24 • @Roland F, thanks for the comment!$r$is just a coordinate,$\rho$is a polar coordinate ($1/r = 1/\sqrt{\rho^2+z^2}$).$1/r$is Coulomb's low May 14 at 9:31 ## 2 Answers To compute exited states $$\Psi_i$$ with some basis functions $$\Phi_j$$ we need to solve equation $$H\Psi=\lambda U\Psi$$, where $$H_{ij}=\langle\Psi_i|\hat{H}|\Psi_j\rangle$$, and $$U_{ij}=\langle\Phi_i|\Phi_j\rangle$$. Last matrix is not inversible in general case, so we have some numerical problem here. For example, exited states of the ground state computed here can be calculated with Eigensystem[] as follows a0j = {1.25, 1.25003, 1.25006, 1.251393, 1.252726, 1.2587955, 1.264865, 1.260499, 1.303186, 1.546232, 2.337185, 4.507307, 9.426298, 20.270036, 43.977048, 95.675215, 208.32642, 453.736619, 988.32236, 2152.803005, 4689.357171, 10214.647096}; b0j = {0.026284, 0.0417685, 0.057253, 0.09098300000000001, 0.124713, 0.19818550000000001, 0.271658, 0.417292, 0.61588, 1.008925, 1.949878, 4.24735, 9.251849, 20.15297, 43.898488, 95.622497, 208.291042, 453.712878, 988.306428, 2152.792314, 4689.349997, 10214.642282}; nmax = Length[a0j]; Psi[r_, z_, j_] := Exp[-b[j]*z^2]*Exp[-a[j]*r^2]; (*kinetic energy*) Kk[r_, z_, n1_, n2_] := FullSimplify[ Psi[r, z, n2]* Laplacian[Psi[r, z, n1], {r, \[Theta], z}, "Cylindrical"]]; ss = Integrate[ Kk[r, z, l1, l2] r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, Assumptions -> {a[l1] > 0, b[l1] > 0, a[l2] > 0, b[l2] > 0}]; Kx = -1/2 2 Pi Sum[c[l1] c[l2] ss, {l1, nmax}, {l2, nmax}]; (*potential energy*) VP1[r_, z_] := -1/Sqrt[r^2 + z^2] + 25/8 r^2 - 5/2; Px = Integrate[ Psi[r, z, l2]*VP1[r, z]*Psi[r, z, l1]*r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, Assumptions -> {a[l1] > 0, b[l1] > 0, a[l2] > 0, b[l2] > 0}]; Px = 2 Pi Sum[c[l1] c[l2] Px, {l1, nmax}, {l2, nmax}]; int = Integrate[ Psi[r, z, l1] Psi[r, z, l2] r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, Assumptions -> {a[l1] > 0, b[l1] > 0, a[l2] > 0, b[l2] > 0}]; norm = {2 Pi Sum[c[l1] c[l2] int, {l1, nmax}, {l2, nmax}] == 1}; cons = Join[Table[a[n] == a0j[[n]], {n, nmax}], Table[b[n] == b0j[[n]], {n, nmax}], norm]; mat = (Kx + Px) /. Join[Table[a[n] -> a0j[[n]], {n, nmax}], Table[b[n] -> b0j[[n]], {n, nmax}]] // Chop; matrix = Last@CoefficientArrays[mat, Table[c[n], {n, nmax}], "Symmetric" -> True] U = Table[ 2 Pi int /. Join[Table[a[n] -> a0j[[n]], {n, nmax}], Table[b[n] -> b0j[[n]], {n, nmax}]], {l1, nmax}, {l2, nmax}]; mat = Inverse[U] . matrix; {val, vec} = Eigensystem[mat];  Therefore, the ground state and exited states are given by val // Sort (*Out[]= {-1.3804, -0.193692, -0.0640854, 0.0576317, 0.278105, \ 0.641224, 1.27678, 2.35197, 4.03324, 5.46049, 9.06797, 16.7551, \ 33.9461, 71.7555, 154.926, 337.651, 738.55, 1617.97, 3553.12, \ 7850.89, 17666., 42031.}*)  First and second state are in a good agreement with ground state (1s) = -1.380398866427, and the first excited state (2s) = -0.193746709717. But all other eigenvalues are differ from exact values. Apparently we need to extend our basis with constrain on U so that Inverse[U] exists. Update 1. To implement main idea from the paper "Hydrogen-Like Systems in Arbitrary Magnetic Fields. A Variational Approach" we use code nmax = 5; Psi[r_, z_, j_] := Exp[-b[j]*z^2]*Exp[-a[j]*r^2]; (*kinetic energy*) Kk[r_, z_, n1_, n2_] := FullSimplify[ Psi[r, z, n2]* Laplacian[Psi[r, z, n1], {r, \[Theta], z}, "Cylindrical"]]; ss = Integrate[ Kk[r, z, l1, l2] r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, Assumptions -> {a[l1] > 0, b[l1] > 0, a[l2] > 0, b[l2] > 0}]; Kx = -1/2 2 Pi Sum[c[l1] c[l2] ss, {l1, nmax}, {l2, nmax}]; (*potential energy*) VP1[r_, z_] := -1/Sqrt[r^2 + z^2] + 25/8 r^2 - 5/2; Px = Integrate[ Psi[r, z, l2]*VP1[r, z]*Psi[r, z, l1]*r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, Assumptions -> {a[l1] > 0, b[l1] > 0, a[l2] > 0, b[l2] > 0}]; Px = 2 Pi Sum[c[l1] c[l2] Px, {l1, nmax}, {l2, nmax}]; int = Integrate[ Psi[r, z, l1] Psi[r, z, l2] r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, Assumptions -> {a[l1] > 0, b[l1] > 0, a[l2] > 0, b[l2] > 0}]; norm = 2 Pi Sum[c[l1] c[l2] int, {l1, nmax}, {l2, nmax}]; U = Table[2 Pi int, {l1, nmax}, {l2, nmax}]; var = Join[{a1, b1, qa, qb}, Table[c[i], {i, nmax}]]; a = a1; b = b1; Do[{a[i] = a[i - 1] qa, b[i] = b[i - 1] qb}, {i, 2, nmax}]; sol = NMinimize[{Kx + Px, {norm == 1., a1 > 0, b1 > 0, qa > 0, qb > 0}}, var] (*Out[]= {-1.37514 + 0. I, {a1 -> 1.24159, b1 -> 0.455892, qa -> 2.12986, qb -> 3.77828, c -> 0.49116, c -> 0.318804, c -> 0.127836, c -> 0.0669821, c -> 0.0506232}}*)  Note, that we compute the ground state as sol[]=-1.37514, and it is not so bad compared to exact value -1.380398866427. Other eigenvalues computed as follows abq = Take[sol[], 4]; matrix = Last@CoefficientArrays[(Kx + Px) /. abq, Table[c[n], {n, nmax}], "Symmetric" -> True] mat = Inverse[U /. abq] . matrix; {val, vec} = Eigensystem[mat]; val // Sort (*Out[]= {-1.37514, 2.20185, 11.2994, 43.117, 159.438}*)  The problem of this method is that even for nmax=8 we have badly conditioning matrix U. • Well, I have tested method from the paper linked above and got message about U that Result for Inverse of badly conditioned matrix may contain significant numerical errors. Therefore we need computer codes written for a CDC Cyber 74 computer or something like that. May 12 at 16:34 • @MamMam Please, see Update 1 to my answer. May 13 at 8:59 • @MamMam Yes it is, matrix is$H_{ij}\$. May 13 at 12:01
• @MamMam Actually they used very advanced code with basis functions 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]];`. If you like, I can show you my code with implementation this idea. May 13 at 13:55
• Well, then please open new topic. :) May 13 at 14:32

Maxwell: There are three forms of fields: |E|=|B| light, |E|>|B| electric fields with a bit of movement of charges. |E|<|B| magetic fields of strong neutral current or worse spin fields on surfaces, with a bit of extra charges. They are easly discriminated by their 4-vector potential and these sectors are Lorentz invariant.

The Hamiltonian of the H-atom in a magnetic field has the following forms with B M = B_z L_z the constant magnetic energy of the angular momentum in state M

Cylindrical coordinates

$$H = -r^{-1} \ \partial_r \left(\ r \ \partial_ r \psi \right) - \partial_{zz}\psi + B M \ \psi + C^2 \psi \ r^2 - \frac{1}{\sqrt{r^2+z^2}} \psi$$

Spherical coordinates

$$H = -\frac{1}{r^2} \partial_r \left( r^2 \partial_r \psi\right) - \frac{1}{r^2 \sin \theta} \partial_{ \theta} (\sin \theta \partial_\theta \psi) - \frac{1}{r} \psi \ + \ B M \psi + C^2 r^2 \psi$$

This is a problematic Hamiltonian for a simple reason: Its good for small magnetic fields everywhere smaller than the electric field. It's ok if the central charge is not a point but is closed at $$r=0$$ by the oscillator potential of a constant charge density in the nucleus.

In a strong magnetic field, if there is a region where B > E eg on the surface of a neutron star, the electric field can locally transformed away yielding a magnetic field of helix form. The effect is known from northern aura, the wave functions concentrates along any magnetic line between energetic boundaries eg by a surface.

In this case, with $$L_z = \pm M \hbar$$ and no limit by $$L^2$$ , there is no lower bound. By emission of a photon and aquiring one more wavelength on the circle the energy can go down to -oo.

Its essential therefore to exactly fixing the numerical constants first. The rest can be seen in

Liberman Hydrogen Atom

• Thank you very much! It's really useful information May 13 at 20:34
• This is nice explanation (+1). May 14 at 2:07