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?

$\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$
    – Mam Mam
    May 15, 2023 at 13:48
  • 1
    $\begingroup$ Please, see my answer. $\endgroup$ May 15, 2023 at 16:27

1 Answer 1

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.

$\endgroup$
12
  • 1
    $\begingroup$ Yes it is. Also we can't handle large nmax*jmax with NMinimize. $\endgroup$ 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$ 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$ May 17, 2023 at 9:03
  • 1
    $\begingroup$ @MamMam You can use it but result is not the same as evaluated with StepMonitor. $\endgroup$ 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$ May 18, 2023 at 3:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.