2
$\begingroup$

I have written a code that uses a guessed eigenvector (avecs) to construct a Hamiltonian (my hamiltonian depends on its own eigenvectors, so I can not solve it directly and I need to use iteration (self consistent field method)). then, utilizing the Eigensystem command, I can find the new values for energy and new eigenvectors.

Now I need to apply a Do loop or While loop, so that I can continue my calculations (replacing the new eigenvectors with the guessed(previous) eigenvector and constructing the new hamiltonian to find the new energies), until the difference between Abs(Etotal[i] - Etotal[i - 1]) > 0.001.

Here is a simplified part of my code, could any one help me to apply a Do/while loop?

NOS = 8;
NOE = 3;
avecs = Table[If[n == j, 1, 0], {j, 1, NOE}, {n, 1, NOS}]

interaction[n_, m_, j_, n1_, n2_] := 
  If[ n1 == n2, 
    If[n == m, 
      1/24 (1 - 3/(n^2 π^2) - 3/(n1^2 π^2))
        Subscript[a, j, n1] Subscript[a, j, n2], 
      ((1 + (-1)^(m + n)) m n)/((m - n)^2 (m + n)^2 π^2) 
        Subscript[a, j, n1] Subscript[a, j, n2]], 
     If[n == m, 
       ((1 + (-1)^(n1 + n2)) n1 n2)/((n1^2 - n2^2)^2 π^2)
         Subscript[a, j, n1] Subscript[a, j, n2] , 
       -((8 (-1 + (-1)^(m + n)) (-1 + (-1)^(n1 + n2)) m n n1 n2)/
           ((m - n)^2 (m + n)^2 (n1 - n2)^2 (n1 + n2)^2 π^4)) 
         Subscript[a,j, n1] Subscript[a, j, n2] ]] 
  /. 
    {Subscript[a, j, n1] -> avecs[[j, n1]], Subscript[a, j, n2] -> avecs[[j, n2]]}

interact = 
  Table[
    N[Sum[
       interaction[n, m, j, n1, n2], {n1, 1, NOS}, {n2, 1, NOS}, {j, 1, NOE}]], 
    {n, NOS}, {m, NOS}]   

eigensys = Eigensystem[interact] 
avecnew = Take[eigensys[[2]], -NOE]
energy = Take[eigensys[[1]], -NOE]
Etotal = Total[energy]
$\endgroup$
7
  • $\begingroup$ Do[ [whatever],{i,[startCounter],[endCounter]}], sure there is an answer out there. Look at the documentation. $\endgroup$ – user27381 Jun 23 '17 at 15:54
  • 1
    $\begingroup$ It is difficult to provide advice with so little information provided. Maybe, FixedPoint would work. Also, please include in your question the code you have tried. $\endgroup$ – bbgodfrey Jun 23 '17 at 15:55
  • $\begingroup$ Possible duplicate of Beginners problem, Do Loop, Eigenfunction iteration $\endgroup$ – user27381 Jun 23 '17 at 15:56
  • $\begingroup$ my advice for a starter, compose a fixed count Do loop first. After you validate its working and tending to converge, convert to While. ( While will run indefinitely if you make some mistake. Do and While are so similar its simple to switch ). After getting comfortable with all that look at Nest and FixedPoint $\endgroup$ – george2079 Jun 23 '17 at 16:37
  • $\begingroup$ Are you sure you don't want Abs[Etotal[i]-Etotal[i-1]] > 0.001? $\endgroup$ – m_goldberg Jun 23 '17 at 21:09
2
$\begingroup$

I suggest implementing the iteration as follows:

NOS = 8;
NOE = 3;
avecs = IdentityMatrix[NOS][[;; NOE]];

interaction[n_, m_, j_, n1_, n2_, a_] :=
  If[n1 == n2,
    If[n == m,
      1/24 (1 - 3/(n^2 π^2) - 3/(n1^2 π^2)) a[[j, n1]] a[[j, n2]], 
      ((1 + (-1)^(m + n)) m n)/((m - n)^2 (m + n)^2 π^2) a[[j, n1]] a[[j, n2]]], 
    If[n == m,
      ((1 + (-1)^(n1 + n2)) n1 n2)/((n1^2 - n2^2)^2 π^2) a[[j, n1]] a[[j, n2]],
      -((8 (-1 + (-1)^(m + n))(-1 + (-1)^(n1 + n2)) m n n1 n2)/
           ((m - n)^2 (m + n)^2 (n1 - n2)^2 (n1 + n2)^2 π^4))
        a[[j, n1]] a[[j, n2]]]]

interact[avecs_] := 
  Table[
    Sum[
      N @ interaction[n, m, j, n1, n2, avecs], 
      {n1, 1, NOS}, {n2, 1, NOS}, {j, 1, NOE}], 
    {n, NOS}, {m, NOS}]

eigensys[avecs_] := Eigensystem[interact[avecs]] // Chop

next[{prev_, avecs_}] :=
  Module[{sys},
    sys = eigensys[avecs];
    {Total[Take[First @ sys, -NOE]], Take[Last @ sys, -NOE]}]

iterResults =
  NestWhileList[next, {0, avecs}, Abs[#1[[1]] - #2[[1]]] > .001 &, 2];
eTotals = iterResults[[All, 1]]

{0, 0.168095, 0.0646216, 0.0646216}

eTotal = Last[eTotals]

0.0646216

Notes

  • I get rid of Subscript. It is really more a nuisance than a help in computation, and it is not needed for formatting output, which is its real purpose.
  • I define functions for doing the computations to take the avecs matrix as a parameter. Parameratizing with avecs makes it easier to use Mathematica's functional iterators.
  • I use NestWhileList to get the iteration history. Should you only want the final result, use NestWhile.
$\endgroup$
4
  • $\begingroup$ Thank you, I learned a lot from this $\endgroup$ – Delaram Nematollahi Jun 24 '17 at 14:55
  • $\begingroup$ m_goldberg I don't know why it gives me this error: Tag List in {{0.036032,-7.18815*10^-17,0.0284966,-5.75052*10^-18,0.00527714,-1.58433*10^-18,0.001847,-6.51987*10^-19},<<6>>,{-6.51987*10^-19,<<6>>,<<20>>}}<<1>><<6>>] is Protected. >> and also it Failed $\endgroup$ – Delaram Nematollahi Jun 24 '17 at 15:11
  • $\begingroup$ @DelaramNematollahi. I can't tell you what is going wrong without more info. I suspect you had multiple conflicting definitions of interaction and/or interact in the notebook you where you tried to evaluate my code. I suggest trying my code in a separate notebook that you open in a fresh Mathematica session. (I just now copied the code from my answer into a new notebook in a new session and it evaluated correctly) $\endgroup$ – m_goldberg Jun 24 '17 at 16:56
  • $\begingroup$ True, sorry, I am a new user. It works fine in a fresh notebook. Thank you so much :) $\endgroup$ – Delaram Nematollahi Jun 24 '17 at 17:08
2
$\begingroup$

There are several ways to iterate the function in the question to convergence. For instance, repeat the computation within a Do-loop until convergence is achieved, then terminate the loop with Throw.

NOS = 8; NOE = 3; ehist = {};
avecs = Table[If[n == j, 1, 0], {j, 1, NOE}, {n, 1, NOS}];

Catch[Do[
    interaction[n_, m_, j_, n1_, n2_] := If[ n1 == n2, 
        If[n == m, 1/24 (1 - 3/(n^2 π^2) - 3/(n1^2 π^2)) Subscript[a, j, n1]
        Subscript[a, j, n2], ((1 + (-1)^(m + n)) m n)/((m - n)^2 (m + n)^2 π^2)
        Subscript[a, j, n1] Subscript[a, j, n2]], If[n == m, ((1 + (-1)^(n1 + n2))
        n1 n2)/((n1^2 - n2^2)^2 π^2) Subscript[a, j, n1] Subscript[a, j, n2] ,
        -((8 (-1 + (-1)^(m + n)) (-1 + (-1)^(n1 + n2)) m n n1 n2)/((m - n)^2 (m + n)^2 
        (n1 - n2)^2 (n1 + n2)^2 π^4)) Subscript[a, j, n1] Subscript[a, j, n2] ]] /. 
        {Subscript[a, j, n1] -> avecs[[j, n1]], Subscript[a, j, n2] -> avecs[[j, n2]]};
    interact = Table[N[Sum[interaction[n, m, j, n1, n2], {n1, 1, NOS}, {n2, 1, NOS}, 
        {j, 1, NOE}]], {n, NOS}, {m, NOS}];
    eigensys = Eigensystem[interact] ;
    avecs = Take[eigensys[[2]], -NOE];
    energy = Take[eigensys[[1]], -NOE];
    Etotal = Total[energy];
    AppendTo[ehist, Etotal];
    If[j > 1 && Abs[ehist[[j]]/ehist[[j - 1]] - 1] < 10^-3, Throw[Etotal]], {j, 10}]]
(* 0.0646216 *)

The convergence history, if desired, is contained in ehist.

ehist
(* {0.168095, 0.0646216, 0.0646216} *)
$\endgroup$
1
  • $\begingroup$ @m_goldberg Thank you so much $\endgroup$ – Delaram Nematollahi Jun 24 '17 at 14:53

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