# 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
// NBPROC 4
//
// 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

@software{mjm_cpp_\$PROJ,
author = {Michael J Marchywka},
title = {\$PROJ},
abstract=(),
institution={},
license={Knowledge sir should be free to all },
publisher={Mike Marchywka},
email={marchywka@hotmail.com},
authorid={orcid.org/0000-0001-9237-455X},
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
https://arxiv.org/pdf/1304.6845.pdf

[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
https://arxiv.org/pdf/1809.08326.pdf

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

marchywka@happy:/home/documents/cpp/proj/freefem/play\$

*/

include "macro_ddm.idp"             // additional DDM functions
//https://doc.freefem.org/documentation/BEM.html

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

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
cout<<vcenter<<endl;
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;
x=1;
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);

medit("grid",Th);
// make everything with the new mesh
v= barrier(x,y);
medit("potential",Th,v);
// 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;

//////////////////////////////////////////////////
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 )
+on(1,2,3,4,psi=0)
;

matrix A=schrod(Wh,Wh);
real[int] rhs=schrod(0,Wh);
//psi[]=A^-1*rhs;
//set(B,solver=GMRES);
//set(A,solver=GMRES);
//set(A,solver=LU);
set(A,solver=CG);
//set(B,solver=CG,eps=.1);
//set(B,solver=CG);
//set(B,solver=CG);
//set(B,solver=LU); // nonseuqre doh
//psi[]=B^-1*rhsa;

psi[]=A^-1*rhs;

psi[]=proout(psi[]);
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.flush;

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

real  pnorm =real(sqrt(int2d(Th)(psi*conj(psi))));
psi=psi/pnorm;
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;
wpsi=psi;
if (l2<l2dun)
{
real[int] xxx=psi[];
for(int wtf=0; wtf<psi[].n; ++wtf)
{  leigens(wtf,lowers)=xxx[wtf]; }
energies(lowers)=ei;
++lowers;
wpsi=psiz;
/////////////////////////////
// project out...
///////////////////////////////
wpsi[]=proout(wpsi[]);
psi=wpsi;
} // 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;
cout.flush;
ifstream file("ioctl.txt");
int ok=0;
file>> ok;
if (ok!=0) {  medit("psi",Th,psi);  }
}

}

} // i

``````

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
`your_install_path/share/FreeFEM/4.12/examples/hpddm/*schrodinger*`

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.

Thanks.

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 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

From the literature, it looks like the Dirac eqn creates a problem but I’m going
to try that in SLEPc now 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.