crunge.f
c
Index
c file runge.f
c
c Testing Runge Kutta
c
c For Time-Dependent Schroedinger Equation
c
c Including the Jacobi diagonalization
c
c Written for NORDPLUS SCHOOL March 1993
c Modified for Stockholm Course March 1995
c
c Based on standard routines used by
c Atomic Physics Group in Bergen
c
c Collected and tested by L. Kocbach
c
c (based in part on standard functions and some routines
c by e.g. J.P. Hansen, IBM Scientific Package (1972)
c etc..... )
c
c commons:
complex vc(32,32)
real akof(32,32),vv(32,32)
c
c akof(*,*) and vv(*,*) define the matrix elements
c see vcalc()
c
complex y(32)
dimension vcc(32,32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /blocka/ vcc
common /block4/ vc, akof, vv,velo ,vdiag, rr
common /block5/ itest,iout,nout
c
open(11,file='values',status='old')
open(13,file='eigenvals')
open(14,file='diagonals')
open(15,file='probabs')
open(16,file='adiaprob')
c
iout=0
read(11,*) nout
read(11,*) itest
read(11,*) velo
read(11,*) tstep
read(11,*) nstat
write(*,*) nout , ' nout '
write(*,*) itest, ' itest '
write(*,*) velo, ' velo '
write(*,*) tstep, ' tstep dt '
write(*,*) nstat, ' no of states, '
c
c Example case: exp(-a(i,j)*t*t) Read from 'values'
c
do 13 i = 1,nstat
13 read(11,*) vdiag(i)
c
do 3 i = 1,nstat
do 4 j = i,nstat
read(11,*) vv(i,j),akof(i,j)
write(*,*) vv(i,j),akof(i,j) ,' V a ',i,j
vv(j,i) = vv(i,j)
akof(j,i) = akof(i,j)
4 continue
3 continue
c
c Starting Time-loop
c
do 6 i = 1,nstat
y(i) = cmplx(0.0,0.0)
6 continue
c
y(1) = cmplx(1.0,0.0)
c
tmin=-10.0/velo
t=tmin
tmax=-tmin
c
c Adjusting the values to avoid BIG FILES
c
c allpoi all points
c prinpo printed points
c fmaxpo maximum printed
c
fmaxpo=300.0
c
allpoi=(tmax - tmin)/tstep
prinpo=allpoi/float(nout)
if(prinpo.gt.fmaxpo) nout=ifix(allpoi/fmaxpo)
c
c
c Time-loop
c
1000 continue
ttt=t
call runge(y,nstat,ttt,tstep)
call output(nstat,y,t)
t=t+tstep
if (t.lt.tmax) goto 1000
END
c
subroutine runge(y,mm,t,dt)
c -----------------------------
c
c performs Runge Kutta 4-order for LINEAR systems
c for TDSE : complex amplitudes necesssary
c
c input y(t) , output y(t+dt)
c matrix elements calculated in vcalc
c
integer mm
real t,dt
complex y(32),k1(32),k2(32),k3(32),k4(32)
complex res(32)
c commons:
complex vc(32,32)
real akof(32,32),vv(32,32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /block4/ vc, akof, vv,velo ,vdiag,rr
c
call prdkt(k1,y,dt,mm)
do 10 i = 1,mm
res(i) = y(i)+ k1(i)*0.5
10 continue
call vcalc(mm,t+0.5*dt,dt/2.0)
call prdkt(k2,res,dt,mm)
do 11 i = 1,mm
res(i) = y(i)+ k2(i)*0.5
11 continue
call prdkt(k3,res,dt,mm)
do 12 i = 1,mm
res(i) = y(i)+ k3(i)
12 continue
t = t + dt
call vcalc(mm,t,dt/2.0)
call prdkt(k4,res,dt,mm)
c
c iteration step:
c
do 20 i = 1,mm
y(i) = y(i) + k1(i)/6 + k2(i)/3 + k3(i)/3 + k4(i)/6
20 continue
end
c
subroutine prdkt(ki,y,dt,mm)
c -----------------------------
c
c returns ki = coup.matrix * y - vector * dt
c
real dt
complex ki(*),y(*)
integer mm
c commons:
complex vc(32,32)
real akof(32,32),vv(32,32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /block4/ vc, akof, vv,velo ,vdiag, rr
do 10 i = 1,mm
ki(i) = cmplx(0.,0.)
do 15 j = 1,mm
ki(i) = ki(i) + vc(i,j)*y(j)
15 continue
ki(i) = ki(i)*dt
10 continue
end
c
subroutine vcalc(nstat,t,dtx)
c -----------------------------
c
c evaluates coupling matrix and multiplies by i
c
real dtx,drx,t,tt,velo
c commons:
complex vc(32,32)
real akof(32,32),vv(32,32)
dimension vcc(32,32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /blocka/ vcc
common /block4/ vc, akof, vv,velo ,vdiag,rr
c
c Build the matrix
c
c Example case: exp(-a(i,j)*t*t)
c
tt=t+dtx
rr=velo*tt
c
do 3 i = 1,nstat
do 4 j = i,nstat
vcc(i,j) = vv(i,j)*exp(-akof(i,j)*rr*rr)
if (abs(vcc(i,j)) .lt. 0.000001)
. vcc(i,j) = 0.000001
vc(i,j) = cmplx(1.0,0.0) * vcc(i,j)
vc(j,i) = vc(i,j)
vcc(j,i) = vcc(i,j)
4 continue
vcc(i,i) = vcc(i,i) + vdiag(i)
vc(i,i) = vc(i,i) + cmplx(vdiag(i),0.0)
3 continue
c
c Multiply matrix by i
c
do 8 i = 1,nstat
do 9 j = 1,nstat
vc(i,j) = cmplx(0.0,-1.0) * vc(i,j)
9 continue
8 continue
return
end
c
subroutine output(nn,yy,tt)
c -----------------------------
c
complex yy(32)
c commons:
complex vc(32,32)
complex yadiab(32),zsum
real akof(32,32),vv(32,32), ach(32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /block4/ vc, akof, vv,velo ,vdiag, rr
common /block5/ itest,iout,nout
c
if(iout.eq.0) then
sum=0.0
do 333 k=1,nn
bach=cabs(yy(k))
ach(k) = bach*bach
sum=sum+ach(k)
333 continue
call exJaco(nn)
c
write(15,802) tt,(ach(l),l=1,nn), rr
c
1111 format('time: ',f12.5,' ........................', f12.8)
1112 format(f12.5,'... follow amps in adiabatic basis .', f12.8)
2222 format(' ',2f12.8,' ',f12.8)
4444 format(' ....... ', 2f12.8)
c
c Transforming to the adiabatic basis
c Using the trafo(*,*)
c
zsum=cmplx(0.0,0.0)
c
do 71 k=1,nn
zsum=cmplx(0.0,0.0)
do 73 ii=1,nn
zsum=zsum+trafo(k,ii)*yy(ii)
73 continue
yadiab(k)=zsum
71 continue
c
sum=0.0
do 334 k=1,nn
bach=cabs(yadiab(k))
ach(k) = bach*bach
sum=sum+ach(k)
334 continue
write(16,902) tt ,(ach(l),l=1,nn) ,rr
802 format(f10.3,' ',5f15.6)
902 format(' ',f10.3,' ',5f15.6)
endif
iout=iout+1
if (iout.eq.nout) iout=0
c
return
end
c
c
c jacob.f
c
c
subroutine exJaco(n)
c --------------------
c
c This is modified main program of the
c test of the diagonalisation routine
c
c See also jacob.f
c
c the eigenvectors from jacob are arranged
c in columns
c v1 v2 v3
c v1 v2 v3
c v1 v2 v3
c eigenvalues are returned on diagonal of a( )
c
c the matrix a is destroyed by jacob
c
c i-th vector(k) = s (k,i)
c i-th eigenvalue is a(i,i)
c...............................................
dimension vcc(32,32)
dimension a(32,32),s(32,32)
dimension akof(32,32),vv(32,32)
dimension vinter1(32,32),vinter2(32,32)
dimension v(32)
dimension vdiag(32)
complex vc(32,32)
dimension trafo(32,32)
common /blockb/ trafo
common /blocka/ vcc
common /block4/ vc, akof, vv,velo ,vdiag, rr
common /block5/ itest,iout,nout
c
c instead of entering, the matrix a(*,*)
c is obtained from vcc(*,*)
c constructed in subroutine vcalc( )
c
do 3 k=1,n
do 33 l=1,n
33 a(k,l)=vcc(k,l)
3 continue
501 format()
601 format(' matrix to be diagonalized')
602 format(' ',10f10.3)
632 format(f9.5,' ',10f9.5)
call jacob(n,a,s)
c
c trafo(*,*) is inverted s(*,*)
c
do 61 k=1,n
do 62 l=1,n
trafo(k,l)=s(l,k)
62 continue
61 continue
c
609 format(' vectors in columns ')
610 format(' eigenvalues for R= ',f9.5)
619 format('____________________________________________________')
do 300 k=1,n
300 v(k)=a(k,k)
c
c writing to file 13 ( see the filenames at open()
c ( at the beginning of main )
c
write(13,632) rr, (v(i),i=1,n)
write(14,632) rr, (vcc(i,i),i=1,n)
c
c
if(itest.eq.1) then
c
c test of diagonals - (itest==1)
c
do 71 k=1,n
do 72 l=1,n
suma=0.0
do 73 ii=1,n
suma=suma+vcc(k,ii)*s(ii,l)
73 continue
vinter1(k,l)=suma
72 continue
71 continue
do 91 k=1,n
do 92 l=1,n
suma=0.0
do 93 ii=1,n
suma=suma+s(ii,k)*vinter1(ii,l)
93 continue
vinter2(k,l)=suma
92 continue
91 continue
endif
c
return
end
c
c This is the diagonalisation routine
c
c See also jacob.f
c
subroutine jacob(nn,a,s)
c -------------------------
c
c
integer p,q,flag
real nu,nuf
dimension a(32,32),s(32,32)
n=iabs(nn)
if(nn) 4,18,1
1 do 3 i=1,n
do 2 k=1,n
2 s(i,k)=0.0
3 s(i,i)=1.0
4 flag=0
if(n-1) 41,18,41
41 z=0.
do 5 q=2,n
m=q-1
do 5 p=1,m
5 z=z+a(p,q)*a(p,q)
nu=1.414*sqrt(z)
z=nu/float(n)
if(z-1.e-27) 18,18,6
6 nuf=z*1.0e-10
7 nu=nu/float(n)
8 do 15 q=2,n
m=q-1
do 15 p=1,m
if(abs(a(p,q))-nu) 15,9,9
9 flag=1
z=0.5*(a(p,p)-a(q,q))
zz=-a(p,q)
sn=sign(1.,z*zz)
zz=zz*zz/(zz*zz+z*z)
z=sqrt(zz)
sn=sn*z/sqrt(2.*(1.+sqrt(1.-zz)))
snsq=sn*sn
cssq=1.0-snsq
cs=sqrt(cssq)
sncs=sn*cs
do 14 i=1,n
if((i-p)*(i-q)) 11,12,11
11 ipr=min0(i,p)
ipc=max0(i,p)
iqr=min0(i,q)
iqc=max0(i,q)
z=a(ipr,ipc)
zz=a(iqr,iqc)
a(ipr,ipc)=z*cs-zz*sn
a(iqr,iqc)=z*sn+zz*cs
12 if(nn) 14,18,13
13 z=s(i,p)
zz=s(i,q)
s(i,p)=z*cs-zz*sn
s(i,q)=z*sn+zz*cs
14 continue
z=a(p,p)
zz=a(q,q)
a(p,p)=z*cssq+zz*snsq-2.0*a(p,q)*sncs
a(q,q)=z+zz-a(p,p)
a(p,q)=(z-zz)*sncs+a(p,q)*(cssq-snsq)
15 continue
if(flag) 18,17,16
16 flag=0
go to 8
17 if(nuf-nu) 7,18,18
18 return
end