4
$\begingroup$

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.

$\endgroup$
9
  • 1
    $\begingroup$ @Susi Lehtola, could you please help with this issue $\endgroup$
    – Mam Mam
    May 30, 2023 at 22:35
  • $\begingroup$ Please provide the error message in a code block: mattermodeling.meta.stackexchange.com/q/417/5 $\endgroup$ Jun 1, 2023 at 18:25
  • $\begingroup$ @Nike Dattani, thanks! Wolfram Mathematica does not give errors, I have used 12 version $\endgroup$
    – Mam Mam
    Jun 4, 2023 at 11:42
  • 1
    $\begingroup$ When you said that it doesn't work for more than 20 basis functions, you mean it just takes too long? How long? $\endgroup$ Jun 4, 2023 at 12:57
  • $\begingroup$ @Nike Dattani, yes, I meant that it takes a very long time due to the fact that the analytic expressions become very large. I can't say the exact time as I didn't expect the final result. I have been waiting for more than a several hours (this is for 20 basis functions). I need to use about 120 functions, so the time will increase dramatically each time with an increase in the basis set $\endgroup$
    – Mam Mam
    Jun 4, 2023 at 14:11

0

You must log in to answer this question.