! ------------------------------------------------------------------------------- ! Hydrogen atom radial Schrödinger equation ! solved as a matrix equation ! ! VERSION 17.03.2006 ! ! Vesa Apaja ! ! Uses LAPACK95 ! ! Before you can compile this you need ! a) LAPACK installed ! b) BLAS installed ! c) LAPACK95 installed ! Assume these are all in $home/lib: ! f90 schro_hydrogen_lapack95.f90 -I $home/lib -L $home/lib -llapack95 -llapack -lblas ! ------------------------------------------------------------------------------------ module precision integer, parameter:: dp=kind(1.d0) end module precision program matrixschro ! LAPACK95: use la_precision, only : wp=>sp use f95_lapack, only : la_stev ,la_stevd ! use precision implicit none integer, parameter:: n=300 ! # of points integer:: i,j,jj real(dp),parameter:: h=0.2d0 ! stepsize real(dp),parameter:: r0=h ,rinf=r0+(n-1)*h ! r limits real(dp):: d(n),e(n),z(n,n),tmp,tmp2,a(n),r(n) ! ! r mesh ! forall (i=1:n) r(i) = r0+(i-1)*h ! ! Matrix elements of the tridiagonal matrix ! e = -1/h**2 ! off-diagonals d = 2/h**2 - 2/r ! diagonals ! ! LAPACK95 (real symmetric tridiagonal matrix) ! call la_stevd(d,e,z) ;print*,'uses la_stevd' !call la_stev(d,e,z) ;print*,'uses la_stev' ! print*,'first eigenvalues:' write(*,'(a5,2a15)') 'n','E','-1/n**2' write(*,'(2x,33("-"))') do i = 1, 10 write(*,'(i5,2f15.5)') i,d(i),-1/dble(i)**2 end do print*,'last eigenvalues:' do i = n-4, n print*,i,d(i) end do ! ! output eigenvector, normalised to sum(z**2)=1 ! ! z(i,j) is the i:th point of j:th eigenvector ! do i = 1, n write(55,'(1000f15.5)') r(i),(z(i,j),j=1,n) write(56,'(6f15.5)') r(i),(z(i,j),j=n-4,n) end do ! ! output eigenvalues (0 and 1 as x for plotting) ! do i = 1, n write(50,*) 0,d(i) write(50,*) 1,d(i) write(50,*) end do ! ! extras ! ------ ! orthogonality of eigenvectors ! tmp = 0.d0 do j = 1, n do jj = j+1, n tmp2 = abs(sum(z(:,j)*z(:,jj)))*h if(tmp2 > tmp) then tmp = tmp2 end if end do end do print*,'max deviation from orthogonality of eigenvectors=',tmp print*,'(Notice that this is a Riemann sum integral)' ! ! eigenvectors as a basis ! ----------------------- ! ! ! compute overlap integrals ! a = 0.d0 do j = 1, n do i = 1, n tmp = f(r(i)) a(j) = a(j) + tmp*z(i,j) end do end do ! ! demonstrate expansion ! do i = 1, n tmp = f(r(i)) tmp2 = sum(a*z(i,:)) ! sum(a(1:20)*z(i,1:20)) ! to take only part!f write(51,'(4es15.5)') r(i),tmp,tmp2,abs(tmp-tmp2) end do contains function f(r) use precision implicit none real(dp):: f real(dp),intent(in):: r ! an arbitrary function f = exp(-0.2*r)*(sin(r)+cos(5*r)) end function f end program matrixschro