Below is an example of a program written in fixed form layout, inherited from punched cards, where column one contains *
, !
, or C
for comments, column six may contain a continuation character, and executable statements begin at column seven. Fortran 90 introduced the Free form layout.
program sphere
c ********** Distribution of eigenvalues in a sphere of radius R
c ********** n atoms, all distances in units of 1/k
complex*16 c1
real*8 pi
integer n, nreal, lambda
integer lwork, ldvl, ldvr
parameter(c1 = (0.0,1.0))
parameter(n = 1000)
parameter(lwork = 4*n)
parameter(ldvl = 1)
parameter(ldvr = 1)
parameter(pi = 3.141592653589793238462643383279502884197)
parameter(nreal = 2000)
parameter(nlambda = 300)
real*8 akR, akR2, rho, beta
real*8 tmp1, tmp2, tmp3
real*8 x1(1:n), x2(1:n), x3(1:n)
real*8 p2d(1:nlambda,1:nlambda)
real*8 pre(1:nlambda), pim(1:nlambda)
real*8 lambdare1, lambdare2, dlambdare
real*8 lambdaim1, lambdaim2, dlambdaim
real*8 rwork(1:2*n)
complex*16 a(1:n,1:n)
complex*16 vl1(1:ldvl, 1:n), vr1(1:ldvr, 1:n)
complex*16 vl(1:n, 1:n), vr(1:n, 1:n)
complex*16 w(1:n)
complex*16 work(1:lwork)
complex*16 ctmp1, ctmp2, ctmp3, ctmp4
integer info
integer i, j, k, m, ireal, nrot
external RandomNumber
real*8 RandomNumber
c ********** Initialization of parameters
rho = 10.0D+00
akR = (2.0D+00)*pi*((3.0D+00*n/
* ((4.0D+00)*pi*rho))**((1.0D+00)/(3.0D+00)))
akR2 = akR**2
c beta = (2.8D+00)*n/((akL**2))
write(6,*) "N = ", n
write(6,*) "rho = ", rho
write(6,*) "kR = ", akR
do i = 1, nlambda
pre(i) = 0.0D+00
pim(i) = 0.0D+00
do j = 1, nlambda
p2d(i,j) = 0.0D+00
end do
end do
lambdare1 = -5.0D+00
lambdare2 = 5.0D+00
dlambdare = (lambdare2 - lambdare1)/nlambda
lambdaim1 = -1.0D+00
lambdaim2 = 7.0D+00
dlambdaim = (lambdaim2 - lambdaim1)/nlambda
evmin1 = 1000.0D+00
evmin2 = 1000.0D+00
evmax = -1000.0D+00
do ireal = 1, nreal
write(6,*) "Running", ireal, " of", nreal
c ********** Initialization of the matrix gamma
do i = 1, n
do j = 1, n
a(i,j) = 0.0D+00
end do
end do
do i = 1, n
tmp1 = (2.0D+00)*akR2
do while(tmp1.GT.akR2)
x1(i) = akR*(-1.0D+00 + (2.0D+00)*RandomNumber(idummy))
x2(i) = akR*(-1.0D+00 + (2.0D+00)*RandomNumber(idummy))
x3(i) = akR*(-1.0D+00 + (2.0D+00)*RandomNumber(idummy))
tmp1 = x1(i)**2 + x2(i)**2 + x3(i)**2
end do
end do
do i = 1, n
do j = 1, n
a(i,j) = (x1(i) - x1(j))**2 + (x2(i) - x2(j))**2 +
* (x3(i) - x3(j))**2
end do
end do
do i = 1, n
do j = 1, n
if(i.EQ.j) then
a(i,j) = 0.0D+00
else
tmp1 = cdsqrt(a(i,j))
a(i,j) = cdexp(c1*tmp1)/tmp1
c a(i,j) = c1*dsin(tmp1)/tmp1
end if
end do
end do
c ********* Using LAPACK (compile with -llapack)
c ********* Compute only the eigenvalues
c call zgeev('N', 'N', n, a, n, w, vl1, ldvl,
c * vr1, ldvr, work, lwork, rwork, info)
c ********* Compute both the eigenvalues and left and right eigenvectors
call zgeev('V', 'V', n, a, n, w, vl, n,
* vr, n, work, lwork, rwork, info)
if(info.NE.0) then
write(6,*) "Error in ZGEEV: info = ", info
pause
end if
c ********* Using IMSL
c CALL DEVLCG (n, a, n, w)
c ********* Save all eigenvalues
if(ireal.EQ.1) then
open(UNIT=12, FILE="ev.dat",
*STATUS="UNKNOWN")
else
open(UNIT=12, FILE="ev.dat",
*STATUS="UNKNOWN", ACCESS="APPEND")
end if
do i = 1, n
write(12,*) dreal(w(i)), " ", dimag(w(i))
end do
close(12)
c ********* Compute and save IPR
if(ireal.EQ.1) then
open(UNIT=12, FILE="ipr.dat",
*STATUS="UNKNOWN")
else
open(UNIT=12, FILE="ipr.dat",
*STATUS="UNKNOWN", ACCESS="APPEND")
end if
do i = 1, n
tmp1 = 0.0D+00
tmp2 = 0.0D+00
do j = 1, n
tmp1 = tmp1 + cdabs(vr(j,i))**2
tmp2 = tmp2 + cdabs(vr(j,i))**4
end do
tmp3 = tmp2/(tmp1*tmp1)
write(12,*) dreal(w(i)), " ", dimag(w(i)), " ", tmp3
end do
close(12)
c ********* Renormalize eigenfunctions to have <L|R> = 1
do i = 1, n
ctmp1 = 0.0D+00
do j = 1, n
c ctmp1 = ctmp1 + dconjg(vl(j,i))*vr(j,i)
ctmp1 = ctmp1 + vr(j,i)**2
end do
ctmp2 = (1.0D+00)/cdsqrt(ctmp1)
do j = 1, n
c vl(j,i) = vl(j,i)*dconjg(ctmp2)
vr(j,i) = vr(j,i)*ctmp2
end do
end do
c ********* Save <L|L><R|R>
if(ireal.EQ.1) then
open(UNIT=12, FILE="corr.dat",
*STATUS="UNKNOWN")
else
open(UNIT=12, FILE="corr.dat",
*STATUS="UNKNOWN", ACCESS="APPEND")
end if
do i = 1, n
tmp1 = 0.0D+00
tmp2 = 0.0D+00
do j = 1, n
tmp1 = tmp1 + cdabs(vr(j,i))**2
c tmp2 = tmp2 + cdabs(vl(j,i))**2
end do
c tmp3 = tmp1*tmp2
tmp3 = tmp1*tmp1/n
write(12,*) dreal(w(i)), " ", dimag(w(i)), " ", tmp3
end do
close(12)
c ********** Compute probability distribution
do i = 1, n
j = int((dreal(w(i)) - lambdare1)/dlambdare) + 1
if((j.GE.1).AND.(j.LE.nlambda)) then
pre(j) = pre(j) + 1.0D+00
else
write(6,*) j, " beyond 1:", nlambda
c pause
end if
end do
do i = 1, n
j = int((dimag(w(i)) - lambdaim1)/dlambdaim) + 1
if((j.GE.1).AND.(j.LE.nlambda)) then
pim(j) = pim(j) + 1.0D+00
else
write(6,*) j, " beyond 1:", nlambda
c pause
end if
end do
do i = 1, n
j = int((dreal(w(i)) - lambdare1)/dlambdare) + 1
k = int((dimag(w(i)) - lambdaim1)/dlambdaim) + 1
if((j.GE.1).AND.(j.LE.nlambda)) then
if((k.GE.1).AND.(k.LE.nlambda)) then
p2d(j,k) = p2d(j,k) + 1.0D+00
else
write(6,*) k, " beyond 1:", nlambda
c pause
end if
else
write(6,*) j, " beyond 1:", nlambda
c pause
end if
end do
end do
c ********** Normalization
tmp1 = (1.0D+00)/(n*nreal*dlambdare)
tmp2 = (1.0D+00)/(n*nreal*dlambdaim)
tmp3 = (1.0D+00)/(n*nreal*dlambdare*dlambdaim)
do i = 1, nlambda
pre(i) = tmp1*pre(i)
pim(i) = tmp2*pim(i)
do j = 1, nlambda
p2d(i,j) = tmp3*p2d(i,j)
end do
end do
c ********** Save P(lambda)
open(UNIT=7, FILE="pre.dat", STATUS="UNKNOWN")
do i = 1, nlambda
tmp1 = lambdare1 + (i - 0.5D+00)*dlambdare
write(7,*) tmp1, " ", pre(i)
end do
close(7)
open(UNIT=7, FILE="pim.dat", STATUS="UNKNOWN")
do i = 1, nlambda
tmp1 = lambdaim1 + (i - 0.5D+00)*dlambdaim
write(7,*) tmp1, " ", pim(i)
end do
close(7)
open(UNIT=7, FILE="p2d.dat", STATUS="UNKNOWN")
do i = 1, nlambda
tmp1 = lambdare1 + (i - 0.5D+00)*dlambdare
do j = 1, nlambda
tmp2 = lambdaim1 + (j - 0.5D+00)*dlambdaim
write(7,*) tmp1, " ", tmp2, " ", p2d(i,j)
end do
end do
close(7)
end
c ***************************************
real*8 function RandomNumber(idummy)
integer idummy
real*8 tmp1
c RandomNumber = rand()
call RANDOM_NUMBER(tmp1)
RandomNumber = tmp1
return
end