This is the continuation on my previous question Why does the minimum eigenvalue change dramatically when one basis function is added to the basis set?
Brief description of the problem: I would like to find excited states of this system: $$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. The 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 energy of the ground state $(E_{min}=-1.3803982710899914)$ was found by @Alex Trounev using variational approach and this value is in excellent agreement with data from the literature (see this answer). The energies of excited states for this system are also known from the literature: $E_{1}=-0.193746$ (the first excited state), $E_{2}=-0.07365$ (the second excited state), $E_{3}=-0.03835$ (the third excited state), etc.
To find the excited states, @Alex Trounev suggested using a penalty method described here https://arxiv.org/pdf/2009.13556.pdf .
Let's start by finding the first excited state ($|\Psi_1\rangle$):
Following the article, the first excited state can be found as follows: $|\Psi_1\rangle=\operatorname{argmin} (O[\Psi_1])$ (formula 1 in the article), where $O[\Psi_1]$ is the objective functional. The objective functional has the following form: $O[\Psi_1]=E[\Psi_1]+\lambda_0|\vec{S}_0-\vec{S}_0^*|^2$ (formula 7 in the article), where $S_0=\frac{\langle\Psi_1|\Psi_0\rangle}{\sqrt{\langle\Psi_1|\Psi_1\rangle\langle\Psi_0|\Psi_0\rangle}}$ (formula 8 in the article), and $\vec{S}_0^ * $ is a set of target overlaps (If I understand correctly, in our case it is equal to 0).
Table 1 of the article shows the scheme of actions for finding excited states. In order to get the values from this scheme, let's turn to section II B of the article.
First of all, let's write the ground state function $|\Psi_0\rangle$, which is the anchor function in this problem. It was found in the previous question.
In the code I renamed $\rho\equiv r$ for convenience.
The ground state:
In[154]:= ClearAll["Global`*"]
al = {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};
bl = {0.026284, 0.0417685, 0.057253, 0.090983, 0.124713, 0.1981855, 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};
cl = {0.006523295622643919`, -0.06054266302240082`, 0.14648250087347384`, -0.3150221502945908`, 0.416620632871218`, -0.3601131808535824`, 0.3912230593543366`,
0.012038542224830733`, 0.16818770546793724`, 0.17778604904585182`, 0.1738326674769321`, 0.12492347058872014`, 0.08226422142736733`, 0.05675211101544814`, 0.03735559092913218`, 0.026154779461694356`, 0.016669555882257724`, 0.01271934679906752`, 0.006364816098521608`,
0.00797395974846925`, -0.0001602541764108519`, 0.006595023995239419`};
nmax = Length[al];
(*the basis function*)
a[j_] := al[[j]];
b[j_] := bl[[j]];
Psi[r_, z_, j_] := Exp[-b[j]*z^2]*Exp[-a[j]*r^2];
(*the ground state function*)
c[j_] := cl[[j]];
Psi0[r_, z_] = Sum[c[j]*Psi[r, z, j], {j, 1, nmax}]
Formula 9 from the article:
(*the 1st excited state*)
(*the basis function*)
Psi1[r_, z_, j_] := Exp[-b1[j]*z^2]*Exp[-a1[j]*r^2];
(*ρ0=|Psi0(|^2)+|Psi1(|^2)(formula 9 in the article)*)
Psi02[r_, z_] :=
Sum[c[l1] *c[l2]*Psi[r, z, l1]*Psi[r, z, l2] , {l1, nmax}, {l2,
nmax}];
Psi12[r_, z_] :=
Sum[c1[l1] *c1[l2]*Psi1[r, z, l1]*Psi1[r, z, l2] , {l1, nmax}, {l2,
nmax}];
ρ0[r_, z_] := Psi02[r, z] + Psi12[r, z]
Relative normalization (formula 11 and 10 in the article):
N0 = NIntegrate[
Psi02[r, z]/ρ0[r, z]*r*2 Pi, {r, 0, Infinity}, {z, -Infinity, Infinity}]
There are 66 parameters (a1[j], b1[j] and c1[j]) in the problem (with the current choice of the excited state function structure). In order to find the derivative with respect to each parameter for relative normalization, the authors give formula 12:
Psi11[r_, z_] := Sum[c1[l] *Psi1[r, z, l] , {l, nmax}];
da1N0 = 2*
Re[NIntegrate[(
E^(-r^2 a1[1] - z^2 b1[1]) *(-r^2)*c1[1]*Psi11[r, z])/ρ0[r,
z]*r*2 Pi, {r, 0, Infinity}, {z, -Infinity, Infinity}]];
db1N0 = 2*
Re[NIntegrate[(
E^(-r^2 a1[1] - z^2 b1[1]) *(-z^2)*c1[1]*Psi11[r, z])/ρ0[r,
z]*r*2 Pi, {r, 0, Infinity}, {z, -Infinity, Infinity}]];
dc1N0 = 2*
Re[NIntegrate[(
E^(-r^2 a1[1] - z^2 b1[1]) *(-r^2)*Psi11[r, z])/ρ0[r, z]*
r*2 Pi, {r, 0, Infinity}, {z, -Infinity, Infinity}]];
A few questions:
I took a certain set of numbers as parameters of the excited wave function (
Psi11) and checked the calculations of the normalization factor and its derivatives, Mathematica gives errors. Is it due to the fact that I misunderstand the formulas from the article?The scheme shown in Table 1 (in the article) is not very clear, especially points (c)–(i). If you could write these points more simply, it would help me a lot.
I would be glad to any comments on solving the problem.
I am sorry for the not quite typical question, but in the context of my previous questions, which make up the topic series, finding a solution to this problem is of interest.
