4
$\begingroup$

There is a system that has the following Hamiltonian: $H=-\frac{1}{2}\Delta-\frac{1}{r}+\frac{B^2}{8}\rho^2-\frac{B}{2} (m + 2 m_s)$, where $r=(\rho,z,\phi), B=5, m=0, m_s=-1/2$. To find the eigenvalues, I would like to use Gaussian functions as a basis set.

The basis set: $\psi_{ij}=e^{-b_{j} z^2}e^{-a_{i} \rho^2}$ , where $a_{i}=a_1 * qa^{i-1}$, $i=1, 2 ,3,..., i_{max}$ and $b_{j}=b_1 * qb^{j-1}$, $j=1, 2 ,3,..., j_{max}$ are geometrical progressions.

I would like to minimize the sum of eigenvalues of the system (minimization in four parameters $a_1, b_1, qa, qb$).

Below is the code in which I am trying to find the eigenvalues of the system Eigenvalues[{EE[a1, b1, qa, qb], normm[a1, b1, qa, qb]}]

If I understand correctly, then the Mathematica cannot find the eigenvalues due to the complicated expressions. Is there any other way to find the sum of the eigenvalues (In order to further minimize)?

In the code I rename $\rho\equiv r$

ClearAll["Global`*"]

imax = 3; jmax = 5;
geoseq[init_, r_, n_] := init*r^(n - 1);

Psi[a1_, b1_, qa_, qb_, r_, z_, i_, j_] := 
  Exp[-geoseq[b1, qb, j]*z^2]*Exp[-geoseq[a1, qa, i]*r^2];


Kk[a1_, b1_, qa_, qb_, r_, z_, i1_, j1_, i2_, j2_] = 
  FullSimplify[
   Psi[a1, b1, qa, qb, r, z, i2, j2]*
    Laplacian[Psi[a1, b1, qa, qb, r, z, i1, j1], {r, \[Theta], z}, 
     "Cylindrical"]];
Kk1[a1_, b1_, qa_, qb_, i1_, j1_, i2_, j2_] = -1/2 2 Pi*
   Integrate[
    Kk[a1, b1, qa, qb, r, z, i1, j1, i2, j2] r, {r, 
     0, \[Infinity]}, {z, -Infinity, Infinity}, 
    Assumptions -> {a1 > 0, b1 > 0, qa > 0, qb > 0, i1 > 0, j1 > 0, 
      i2 > 0, j2 > 0}];
Kx[a1_, b1_, qa_, qb_] = 
  Table[ Kk1[a1, b1, qa, qb, i1, j1, i2, j2], {i1, 1, imax}, {i2, 1, 
    imax}, {j1, 1, jmax}, {j2, 1, jmax}];
KK[a1_, b1_, qa_, qb_] = Flatten[Kx[a1, b1, qa, qb], {{1, 3}, {2, 4}}];

B = 5; m = 0; ms = -1/2;
VC[r_, z_] := -1/Sqrt[r^2 + z^2];
Px1[a1_, b1_, qa_, qb_, i1_, j1_, i2_, j2_] = 
  2 Pi*Integrate[
    Psi[a1, b1, qa, qb, r, z, i2, 
      j2] *(VC[r, z] + 1/8 B^2 r^2 + B/2 (m + 2 ms))*
     Psi[a1, b1, qa, qb, r, z, i1, j1]*r, {r, 
     0, \[Infinity]}, {z, -Infinity, Infinity}, 
    Assumptions -> {a1 > 0, b1 > 0, qa > 0, qb > 0, i1 > 0, j1 > 0, 
      i2 > 0, j2 > 0}];
Px[a1_, b1_, qa_, qb_] = 
  Table[Px1[a1, b1, qa, qb, i1, j1, i2, j2], {i1, 1, imax}, {i2, 1, 
    imax}, {j1, 1, jmax}, {j2, 1, jmax}];
PP[a1_, b1_, qa_, qb_] = Flatten[Px[a1, b1, qa, qb], {{1, 3}, {2, 4}}];

EE[a1_, b1_, qa_, qb_] = KK[a1, b1, qa, qb] + PP[a1, b1, qa, qb];

(*normalization*)
int[a1_, b1_, qa_, qb_, i1_, j1_, i2_, j2_] = 
  2 Pi *Integrate[
    Psi[a1, b1, qa, qb, r, z, i2, j2] *
     Psi[a1, b1, qa, qb, r, z, i1, j1]  r, {r, 
     0, \[Infinity]}, {z, -Infinity, Infinity}, 
    Assumptions -> {a1 > 0, b1 > 0, qa > 0, qb > 0, i1 > 0, j1 > 0, 
      i2 > 0, j2 > 0}];
norm[a1_, b1_, qa_, qb_] = 
  Table[int[a1, b1, qa, qb, i1, j1, i2, j2], {i1, 1, imax}, {i2, 1, 
    imax}, {j1, 1, jmax}, {j2, 1, jmax}];
normm[a1_, b1_, qa_, qb_] = 
  Flatten[norm[a1, b1, qa, qb], {{1, 3}, {2, 4}}];

Eigenvalues[{EE[a1, b1, qa, qb], normm[a1, b1, qa, qb]}]
Improve this question
$\endgroup$
8
  • $\begingroup$ The sum of the eigenvalues is equal to the trace of the matrix. Is this what you need? $\endgroup$
    – Roman
    Commented May 30, 2023 at 14:09
  • $\begingroup$ @Roman, thanks! Yes, the trace of the diagonal matrix Eigenvalues[{EE[a1, b1, qa, qb], normm[a1, b1, qa, qb]}] But unfortunately Mathematica does not diagonalize this expression $\endgroup$
    – Mam Mam
    Commented May 30, 2023 at 14:23
  • 2
    $\begingroup$ What I mean is that you don't need to diagonalize. I'll write down a solution. $\endgroup$
    – Roman
    Commented May 30, 2023 at 14:41
  • 1
    $\begingroup$ @MamMam It is not clear what do you try to compute? $\endgroup$
    – Alex Trounev
    Commented May 30, 2023 at 17:12
  • $\begingroup$ @Alex Trounev, thanks for the question! I would like to get an analytic expression for the sum of the eigenvalues (which depends on 4 parameters a1, b1, qa, qb). And then minimize this expression. But the problem is that Mathematica cannot do the inversion operation for the matrices greater than imax=3, jmax=4 $\endgroup$
    – Mam Mam
    Commented May 30, 2023 at 17:25

2 Answers 2

Reset to default
6
$\begingroup$

The sum of the eigenvalues can be calculated without diagonalizing the matrix.

Example: generate random $n\times n$ Hermitian matrices:

n = 7;
A = # . ConjugateTranspose[#] &[RandomVariate[NormalDistribution[], {n, n, 2}] . {1, I}];
B = # . ConjugateTranspose[#] &[RandomVariate[NormalDistribution[], {n, n, 2}] . {1, I}];

Sum of generalized eigenvalues:

Total[Eigenvalues[{A, B}]]
(*    184.183    *)

Direct calculation without diagonalization:

Tr[A . Inverse[B]]
(*    184.183 + 4.84945*10^-13 I    *)

Chopping off the numerical-precision imaginary part:

Re[Tr[A . Inverse[B]]]
(*    184.183    *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks a lot, but in my case, the matrix normm[[a1, b1, qa, qb]] depend on four parameters and Mathematica cannot perform the inversion operation. The matrices are parameter dependent because I would like to then minimize the sum of the eigenvalues $\endgroup$
    – Mam Mam
    Commented May 30, 2023 at 15:02
  • 1
    $\begingroup$ If normm cannot be inverted, then your generalized eigenvalue problem is invalid. Invertibility of the matrix $B$ is a pre-condition. $\endgroup$
    – Roman
    Commented May 30, 2023 at 15:21
  • $\begingroup$ Thanks! I think in my case this is due to the fact that the expressions are very large. Is there a way to do the inversion of a large matrix in parts, i.e. sequentially working with its components? $\endgroup$
    – Mam Mam
    Commented May 30, 2023 at 16:11
2
$\begingroup$

First we need to test code with known parameters from my answer here. Please, pay attention that your code has been updated.

ClearAll["Global`*"]

imax = 3; jmax = 4;
geoseq[init_, r_, n_] := init*r^(n - 1);

Psi[a1_, b1_, qa_, qb_, r_, z_, i_, j_] := 
  Exp[-geoseq[b1, qb, j]*z^2]*Exp[-geoseq[a1, qa, i]*r^2];


Kk[a1_, b1_, qa_, qb_, r_, z_, i1_, j1_, i2_, j2_] = 
  1/2 FullSimplify[
    Psi[a1, b1, qa, qb, r, z, i2, j2]*
      Laplacian[Psi[a1, b1, qa, qb, r, z, i1, j1], {r, \[Theta], z}, 
       "Cylindrical"] + 
     Psi[a1, b1, qa, qb, r, z, i1, j1]*
      Laplacian[Psi[a1, b1, qa, qb, r, z, i2, j2], {r, \[Theta], z}, 
       "Cylindrical"]];
Kk1[a1_, b1_, qa_, qb_, i1_, j1_, i2_, j2_] = -1/2 2 Pi*
   Integrate[
    Kk[a1, b1, qa, qb, r, z, i1, j1, i2, j2] r, {r, 
     0, \[Infinity]}, {z, -Infinity, Infinity}, 
    Assumptions -> {a1 > 0, b1 > 0, qa > 0, qb > 0, i1 > 0, j1 > 0, 
      i2 > 0, j2 > 0}];
Kx[a1_, b1_, qa_, qb_] = 
  Table[Kk1[a1, b1, qa, qb, i1, j1, i2, j2], {i1, 1, imax}, {i2, 1, 
    imax}, {j1, 1, jmax}, {j2, 1, jmax}];
KK[a1_, b1_, qa_, qb_] = ArrayFlatten[Kx[a1, b1, qa, qb], 2];

B = 5; m = 0; ms = -1/2;
V[r_, z_] := -1/Sqrt[r^2 + z^2] + 1/8 B^2 r^2 + B/2 (m + 2 ms);
Px1[a1_, b1_, qa_, qb_, i1_, j1_, i2_, j2_] = 
  2 Pi*Integrate[
    Psi[a1, b1, qa, qb, r, z, i2, j2]*(V[r, z])*
     Psi[a1, b1, qa, qb, r, z, i1, j1]*r, {r, 
     0, \[Infinity]}, {z, -Infinity, Infinity}, 
    Assumptions -> {a1 > 0, b1 > 0, qa > 0, qb > 0, i1 > 0, j1 > 0, 
      i2 > 0, j2 > 0}];
Px[a1_, b1_, qa_, qb_] = 
  Table[Px1[a1, b1, qa, qb, i1, j1, i2, j2], {i1, 1, imax}, {i2, 1, 
    imax}, {j1, 1, jmax}, {j2, 1, jmax}];
PP[a1_, b1_, qa_, qb_] = ArrayFlatten[Px[a1, b1, qa, qb], 2];

EE[a1_, b1_, qa_, qb_] = KK[a1, b1, qa, qb] + PP[a1, b1, qa, qb];

(*normalization*)
int[a1_, b1_, qa_, qb_, i1_, j1_, i2_, j2_] = 
  2 Pi*Integrate[
    Psi[a1, b1, qa, qb, r, z, i2, j2]*
     Psi[a1, b1, qa, qb, r, z, i1, j1] r, {r, 
     0, \[Infinity]}, {z, -Infinity, Infinity}, 
    Assumptions -> {a1 > 0, b1 > 0, qa > 0, qb > 0, i1 > 0, j1 > 0, 
      i2 > 0, j2 > 0}];
norm[a1_, b1_, qa_, qb_] = 
  Table[int[a1, b1, qa, qb, i1, j1, i2, j2], {i1, 1, imax}, {i2, 1, 
    imax}, {j1, 1, jmax}, {j2, 1, jmax}];

U[a1_, b1_, qa_, qb_] = ArrayFlatten[norm[a1, b1, qa, qb], 2];


Eigenvalues[{EE[a1, b1, qa, qb], 
   U[a1, b1, qa, qb]} /. {a1 -> 1.3576595084125154`, 
   b1 -> 10.362298171005488`, qa -> 2.6994994668291183`, 
   qb -> 0.19513529311168284`}]//Chop//Sort 

(*{-1.37695, -0.146211, 1.5772, 5.09045, 6.16749, 7.94828, 13.9307, \
20.2562, 22.5456, 24.3574, 25.8672, 37.5024}*)

As we know exact value for the ground state is -1.380398866427, and in this model is about -1.37695. The first excited state (2s) =-0.193746709717, while in this model -0.146211.

Now if we increase dimension as imax = 3; jmax = 5; then we have with the same parameters out

{-1.37695, -0.187666, -0.0121764, 1.68325, 5.09161, 6.11695, 6.29631, \
8.054, 14.0402, 20.3601, 22.5507, 24.3334, 24.5055, 25.9772, 37.615}

We see, that the ground state energy is not improved, but the first exited state looks much better. Let take imax = 3; jmax = 6;, then we have

{-1.37695, -0.188485, -0.0650577, 0.0163966, 1.70671, 5.09187, \
6.11618, 6.23224, 6.32403, 8.07739, 14.0644, 20.383, 22.5519, \
24.3326, 24.4742, 24.5367, 26.0015, 37.6399}

As we can see the second excited state is about -0.0650577, while state (3 d (m = 0)) = -0.07365. For imax = 3; jmax = 7; we have

{-1.37695, -0.188638, -0.0659596, -0.0271467, 0.0227442, 1.71189, \
5.09193, 6.11603, 6.23139, 6.26886, 6.33015, 8.08256, 14.0698, \
20.3881, 22.5521, 24.3325, 24.4733, 24.5155, 24.544, 26.0069, 37.6454}

We can compare with

ground state (1 s) = -1.3793 (exact value = -1.380398866427)

the first excited state (2 s) = -0.19335 (exact value =
-0.193746709717)

the second excited state (3 d (m = 0)) = -0.07365

the third excited state (3 s) = -0.03835 Finally note, that with this model we can compute for instance imax = 15; jmax = 15; with 4 lowest states {-1.37851, -0.189091, -0.0665601, -0.0284486}.

So this is a good model even without optimization a1, b1, qa, qb. We can optimize parameters for imax=3,jmax=4as follows

var = Join[{a1, b1, qa, qb}, 
  Flatten[Table[c[i, j], {i, imax}, {j, jmax}]]]; vec = 
 Flatten[Table[
   c[i, j], {i, imax}, {j, 
    jmax}]]; 
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}]; con = 
 
 sol = 
 NMinimize[{vec . Re[EE[a1, b1, qa, qb]] . 
    vec, {vec . Re[U[a1, b1, qa, qb]] . vec == 1, a1 > 0, b1 > 0, 
    qa > 0, qb > 0}}, var]

(*Out[]= {-1.37452, {a1 -> 1.37897, b1 -> 3.48613, qa -> 2.61302, 
  qb -> 0.458814, c[1, 1] -> 0.0262767, c[1, 2] -> -0.220731, 
  c[1, 3] -> -0.0629083, c[1, 4] -> -0.399011, c[2, 1] -> -0.113762, 
  c[2, 2] -> -0.0347169, c[2, 3] -> -0.169213, c[2, 4] -> 0.128463, 
  c[3, 1] -> -0.460221, c[3, 2] -> 0.493378, c[3, 3] -> -0.201538, 
  c[3, 4] -> 0.0150695}}*)
Improve this answer
$\endgroup$
16
  • 1
    $\begingroup$ There are so many answers, but it does not matter what optimal combination we used. For instance for a1 -> 1.39282, b1 -> 19.9993, qa -> 3.4335, qb -> 0.252195 and imax = 29; jmax = 37; we have ` {-1.37924, -0.186864, -0.0663546, -0.0302142}`. Probably we should play with $n_r, n_z$ to get better results. $\endgroup$
    – Alex Trounev
    Commented Jun 1, 2023 at 7:14
  • 1
    $\begingroup$ @MamMam We can optimize parameters with respect to known values of energy {-1.3804, -0.193747, -0.07365, -0.03835}. For example, at imax=8; jmax=8 optimized parameters are {a1, b1, qa, qb}={1.2774, 21.0146, 1.82603, 0.271406}. $\endgroup$
    – Alex Trounev
    Commented Jun 1, 2023 at 11:53
  • 1
    $\begingroup$ We define these parameters once for a given combination imax=8, jmax=8. Then we can use it for any base. If you are looking for the method how to optimize parameters from the scratch, then read attentively paper Excited states in variational Monte Carlo using a penalty method :) $\endgroup$
    – Alex Trounev
    Commented Jun 1, 2023 at 15:23
  • 1
    $\begingroup$ @MamMam See last code in my answer for optimization. $\endgroup$
    – Alex Trounev
    Commented Jun 1, 2023 at 16:26
  • 1
    $\begingroup$ @MamMam Yes, you are right. It is why we need code to minimize without NMinimize. Could you explore some code in any language adopted to this problem? Then we can translate it into Mathematica .nb. $\endgroup$
    – Alex Trounev
    Commented Jun 2, 2023 at 9:12

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.