Eigenvalues for schroedinger eqn - extremal and project out

I’m still trying to learn both FF and numerical methods details - initially I was going to do some
EM problems but wanted to get a fun 3D geometry and had meshing issues.
So, I backed up to 2D QM issue of interest . This code has been hacked up terribly
for learning but it found the lowest eigenvalue for sch eqn by complex time
stepping ( see the reference) which is similar apparently to the “power method”
to find largest eigenvalue but this finds the ground state. It seems to work well
by projecting out the prior eigenvectors. Curious if anyone cares to comment
on any utility of this approach or how it should be done and similar with the
FF coding. Thanks.

//  run with MPI:  ff-mpirun -np 4 script.edp
// code borrowed from examples/ffddm/heat-torus-3d-surf.edp


marchywka@happy:/home/documents/cpp/proj/freefem/play$ FreeFem++ schtsmatR.edp  -dt .001  -vc .01  -nx 200 -ny 200

  author = {Michael J Marchywka},
  title = {$PROJ},
license={Knowledge sir should be free to all },
publisher={Mike Marchywka},
  filename = {$PROJ},
  url = {},
  version = {0.0.0},
  date-started = {$DATE}

Complex time steps the schrodeinger eqn to get lower eigenvalue
and project out those already found. 

[1]Solving the Schr\"{o}dinger eigenvalue problem by the imaginary time propagation technique using splitting methods with complex coefficients, Bader , Philipp[ ...] Casas , Fernando; The Journal of Chemical Physics. 2013 09/28/2013

[2]How the inter-electronic potential Ans\"{a}tze affect the bound state solutions of a planar two-electron quantum dot model, Caruso , F.[ ...] Silveira , F.; Physica E: Low-dimensional Systems and Nanostructures. 2019 2019

Did not do too bad but some out of order, 

> x=1:25;  k=1; c=pi*pi/100/100; for( i in 1:5) { for ( j in 1:5 )  { w=i*i+j*j; x[k]=c*w;k=k+1;  }} ; sort(x) 
 [1] 0.001973921 0.004934802 0.004934802 0.007895684 0.009869604 0.009869604
 [7] 0.012830486 0.012830486 0.016778327 0.016778327 0.017765288 0.019739209
[13] 0.019739209 0.024674011 0.024674011 0.025660971 0.025660971 0.028621853
[19] 0.028621853 0.031582734 0.033556655 0.033556655 0.040465378 0.040465378
[25] 0.049348022

 2854  FreeFem++ schpro2.edp  -dt .001  -vc 0  -nx 200 -ny 200
 energy 0 0.00197335
 energy 1 0.00493256
 energy 2 0.00986318
 energy 3 0.00493234
 energy 4 0.00788999
 energy 5 0.00985347
 energy 6 0.0128064
 energy 7 0.0177122
 energy 8 0.0128016
 energy 9 0.0167442
 energy 10 0.0167405
 energy 11 0.0196924
 energy 12 0.025572
 energy 13 0.0196778



load "medit"                        // HPDDM plugin
include "macro_ddm.idp"             // additional DDM functions
load "msh3"

// syntax errors due to forgotten macro are a huge problem
// and not clear from eror message !!!!!!!!!!!!!!!!!!!!
//macro init(i)i// EOM                        // scalar field initialization

macro Grad(u)[dx(u), dy(u)]// EOM    // two-dimensional gradient
func Pk = P2;                               // finite element space
// box - to + or 2x in length 
int szx = getARGV("-szx", 50);               // refinement factor
int szy = getARGV("-szy", 50);               // refinement factor
int nx =getARGV("-nx", 100);
int ny = getARGV("-ny", 100);

real dt = getARGV("-dt", 1e-4);              // time step
int itermax = getARGV("-itermax", 10000);    // number of iterations
int nstates = getARGV("-nstates", 20);       // number of energies to keep 

// error threshold to view and accept
real l2dun= getARGV("-l2dun", 1e-4);              //  
real l2view= getARGV("-l2view", 2e-4);              // 

// potential params 
int sz2 = getARGV("-sz2", 10);               // refinement factor
real vcenter= getARGV("-vc", 0.0);              //  V scale factor  
real vx= getARGV("-vx", 1.0);              // scale factor in x for V 
real vy= getARGV("-vy", 1.0);              // scale factoor in y for V

cout<<" szx "<<szx<<" szy "<<szy<<" sz2 "<<sz2<<endl;

// 2real R = 3, r=1, h=0.2;
//macro  sum(v,xx) { real s=0; for(int i=0; i<xx[].n; ++i) s=s+xx[][i]; v= s; } // EOM
//func real a(real t) { real slow=.3*1.5; real fast=1.5*1.5; real bp1=.5-1.2*sz2/sz; real bp2=1.0-bp1;real norm=fast*(bp1+1-bp2)+slow*(bp2-bp1);  if (t<bp1) return fast*t/norm; if ( t<bp2) return bp1*fast/norm+(t-bp1)*slow/norm; return bp1*fast/norm+(bp2-bp1)*slow/norm+(t-bp2)*fast/norm;  } 
func real ax(real t) { real slow=.3*1.5; real fast=1.5*1.5; real bp1=.5-1.2*sz2/szx; real bp2=1.0-bp1;real norm=fast*(bp1+1-bp2)+slow*(bp2-bp1);  if (t<bp1) return fast*t/norm; if ( t<bp2) return bp1*fast/norm+(t-bp1)*slow/norm; return bp1*fast/norm+(bp2-bp1)*slow/norm+(t-bp2)*fast/norm;  } 
func real ay(real t) { real slow=.3*1.5; real fast=1.5*1.5; real bp1=.5-1.2*sz2/szy; real bp2=1.0-bp1;real norm=fast*(bp1+1-bp2)+slow*(bp2-bp1);  if (t<bp1) return fast*t/norm; if ( t<bp2) return bp1*fast/norm+(t-bp1)*slow/norm; return bp1*fast/norm+(bp2-bp1)*slow/norm+(t-bp2)*fast/norm;  } 

//func  fx=  -szx + ax(x)*2*szx;
func  fx=  -szx + x*2*szx;
//func  fx=  -szx + 2*x*(.1*(x-.5)*(x-.5)+0*(y-.5)*(y-.5))*2*szx;
//func  fx=  -szx + 4*(x-.5)*((x-.5)*(x-.5)+0*(y-.5)*(y-.5))*2*szx*szx*szx;

//func fy=  -szy + ay(y)*2*szy;
func fy=  -szy + (y)*2*szy;
//func  fy=  -szy + 4*(y-.5)*((y-.5)*(y-.5)+0*(y-.5)*(y-.5))*2*szy*szy*szy;
cout<<" mnx mny "<<nx<<" "<<ny<<" "<<fx<<endl;

// potential function 
func real barrier( real xx, real yy)

//return vcenter*(vx*xx*xx/szx/szx+vy*yy*yy/szy/szy); 
//return -(10.9*vcenter)+vcenter/(1.0e-1+sqrt(vx*xx*xx/szx/szx+vy*yy*yy/szy/szy)); 
return 1*(vcenter)+vcenter*log(1.0e-1+sqrt(vx*xx*xx/szx/szx+vy*yy*yy/szy/szy)); 
//return 1000;
//cout<<xx<<" "<<sz2<<endl;
if ((xx<=sz2)&&(xx>= -sz2))  
if ((yy<=sz2)&&(yy>= -sz2))  return vcenter; 
return 0;

mesh Th = square(nx,ny,[fx,fy]);

fespace Vh(Th, Pk);           // local finite element space
Vh v= barrier(x,y);

real cut=.1;
real tol=.5;
Th = adaptmesh(Th, 1000.0*pow(dx(v)*dx(v)+dy(v)*dy(v),3), splitpbedge=1,  abserror=0, nbvx=1e4, cutoff=cut, err=tol, inquire=0, ratio=.01, hmin=.005,  hmax=1.5); 

// make everything with the new mesh 
v= barrier(x,y);
// not needednow
Vh xx=x;
Vh yy=y;

fespace Wh(Th, Pk);           // local finite element space
Wh psiz=1.0/(szx*szy);  
Wh psi=psiz,vpsi,wpsi=psiz; 

macro SCHM(v,w)(((Grad(w)'*Grad(v) ))) // EOM
int lowers=0;
// there has to be a better structure for this? 
real[int,int] leigens(psi[].n,nstates);
real[int] energies(nstates);

// project out lower states from current solution 
func real[int]proout(real[int] inpsi)
real[int] xxx=inpsi;
real[int] xxxf=inpsi;
for(int j=0; j<lowers; ++j)
for(int l=0; l<xxxf.n; ++l) { xxx[l]= leigens(l,j); } 
// WTF when does it tark areeays anw ehen fespace  
Wh fu; fu[]=xxxf;
Wh sht; sht[]=xxx;
real nu=int2d(Th)(sht*fu);
// this seems to be "1" now but still not sure code is right. 
real de=int2d(Th)(sht*sht); // should just do once 
real f=nu/de;
cout<<" j="<<j<<" f="<<f<<" nu="<<nu<<" de="<<de<<endl; cout.flush;
for(int l=0; l<psi[].n; ++l) { xxxf[l]= xxxf[l]-f*xxx[l];  }
} // j  
return xxxf;
} // inpsi 

for(int i=0; i<itermax; ++i)
// complex time stepping or filtering out the
// higher eigenvalues, almost opposite of the "power method"
// for eigenvalue finding 
varf schrod( psi,vpsi,solver=CG,init=1) =
int2d(Th) (   SCHM(psi,vpsi))
+int2d(Th)(  v*psi*vpsi)
+int2d(Th)( wpsi*vpsi/dt )

matrix A=schrod(Wh,Wh);
real[int] rhs=schrod(0,Wh);
//set(B,solver=LU); // nonseuqre doh 


if (false) { 
for(int jj=0; jj<lowers; ++jj)
// this may be ok for dot product not integral-
Wh test;
for(int k=0; k<psi[].n; ++k)
	{test[][k]= leigens(k,jj); } 

real nu=int2d(Th)(psi*test);
real du=int2d(Th)(test*test);
real com=nu/du;
cout<<"wtf " <<jj<<" nu "<<nu<<" du "<<du<<" eta "<<com <<endl;
for(int k=0; k<psi[].n; ++k) psi[][k]=psi[][k]-com*leigens(k,jj); 

} // jj 
} // false 


//cout<<"done solving "<<endl; cout.flush;

real  pnorm =real(sqrt(int2d(Th)(psi*conj(psi))));
real  l2 =real(sqrt(int2d(Th)((wpsi-psi)*conj(wpsi-psi))));
Wh H= -dxx(psi)-dyy(psi)+v*psi;
// going to all real loses some decay lifetime etc. 
real ei=int2d(Th)(H*psi); // do]t product? 
// this should be one after norm... 
real psi2=int2d(Th)(psi*conj(psi)); // do]t product? 
Wh Perp=H-ei*psi;
real  perp=sqrt(int2d(Th)(Perp*conj(Perp)));
cout<<" iter= "<<i<<" lowers="<<lowers<< " e ="<<ei<<" perp="<<perp<< " psi2="<<psi2<< " pnorm="<<pnorm<< " l2="<<l2<<endl;
if (l2<l2dun)
real[int] xxx=psi[];
for(int wtf=0; wtf<psi[].n; ++wtf)
	{  leigens(wtf,lowers)=xxx[wtf]; } 
// project out... 
} // l2 low enougth 

 macro myplot()cmm = "Global solution at iteration " + i, fill = 1, value = 1//
// ffglut gets sigabort after about 7 iters... 
//if ((i%3)==0){ plot(amag,wait=1,myplot); } 
if ((lowers!=0)||(l2<l2view)){ 
//NDIS(amag, " psi ");
//NDIS(psi, " psi ");
if ((l2<l2view)||((i%10)==9)){
for(int ii=0; ii<lowers; ++ii)
cout<<" energy " <<ii<<" "<<energies(ii)<<endl;
ifstream file("ioctl.txt");
int ok=0;
file>> ok;
if (ok!=0) {  medit("psi",Th,psi);  } 


} // i 

[schpro2.edp|attachment](upload://jf70tVOf1C4KgnI7Idm3ze7FTy3.edp) (9.0 KB)

it does not appear to file uploaded, fwiw,

schpro2.edp (9.0 KB)

Hi. I don’t know if it’s useful for you, but some time ago we did wrote some examples for Schrodinger eigenvalue problems using SLEPc… They are in

yes thanks I have looked at at least one or two of those and IIRC
mentioned one here to someone else concerned about
degeneracy. I can look more closely at the SLEPc althorightms.

Maybe I will run that example for a comparison and see if I can
find any reason this approach of extrmea and project may have any benefits.

SLEPc implements the power method, among many other, much more sophisticated and better, algorithms. You are kind of wasting your time reimplementing the wheel.


I would tend to agree but I was intrigues looking at the schroedinger eqn
and finding the “imaginary time” approach since it was bothering me
that I could not figure out how to use the obvious “low pass filtering”
of repeated application of H. So, I had to actually try it :slight_smile: This seems to be the opposite
of the power methid that amplifies the largest accessble eigen value.
I was also curious how many states you could project out and still
get decent eigen values. I also discovered that a dot product is not the
same as a norm.intelgral :slight_smile:

From the literature, it looks like the Dirac eqn creates a problem but I’m going
to try that in SLEPc now :slight_smile: I took the example mentioned earlier and ti seems
to convert to 3D easily ( haven’t checked the eigenvalues yet ) now
just need to put in Dirac eqn.