1
$\begingroup$

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?

Improve this question
$\endgroup$
2
  • $\begingroup$ @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. $\endgroup$ Commented May 15, 2023 at 13:48
  • 1
    $\begingroup$ Please, see my answer. $\endgroup$ Commented May 15, 2023 at 16:27

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
12
  • 1
    $\begingroup$ Yes it is. Also we can't handle large nmax*jmax with NMinimize. $\endgroup$ Commented May 16, 2023 at 2:31
  • 1
    $\begingroup$ 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] $\endgroup$ Commented May 17, 2023 at 3:49
  • 1
    $\begingroup$ 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). $\endgroup$ Commented May 17, 2023 at 9:03
  • 1
    $\begingroup$ @MamMam You can use it but result is not the same as evaluated with StepMonitor. $\endgroup$ Commented May 17, 2023 at 11:36
  • 1
    $\begingroup$ @MamMam Yes it is. This is well know mathematical problem. You can ask about this problem on math.stackexchange.com $\endgroup$ Commented May 18, 2023 at 3:16

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.