program main ! solve 1D TDSE using variational method ! with sin-basis (fixed boundary condition). ! ! [-\frac{1}{2}\frac{d^2}{dx^2}+V(x)]\psi(x) = E\psi(x) ! implicit none integer::i,N,Nout,key,Ncpu double precision::a,b,eps double precision,allocatable::H(:,:),E(:) double precision,external::V a=-20d0 ! Left side box size b= 20d0 ! Right side box size (a < b) N=80 ! Number of sin-basis function Nout=10 ! Number of output eigen functions key=3 ! Kind of quadrature method in QUADPACK eps=1d-10 ! Epsilon of the quadrature method Ncpu=1 ! Numebr of threads in the computation allocate(H(1:N,1:N),E(1:N)); H=0d0; E=0d0 !$ call omp_set_num_threads(Ncpu) !$ call mkl_set_num_threads(Ncpu) ! Construct matrix element of Hamiltonian ! Basis function : orthnormalized sin function [a,b] ! Integral method : Adapted Gauss-legendere & Kronrod quadrature call Hamiltonian(N,H,V,a,b,key,eps) ! Diagonalize hermite matrix using lapack "dsyevd" for real hermitian matrix call Diagonalize(N,H,E) ! Normalize eigenvector call Normalize(N,H) ! Output eigenenergy to a terminal & files open(10,file="eigenenergy.d") do i=1,Nout write(6,*)i,E(i) write(10,*)i,E(i) enddo close(10) ! Output wavefunction to files call Eigenfunction_output(N,H,E,Nout,a,b) ! Output potential to a file call Potential_output(400,V,1.2*a,1.2*b) stop end program main function V(x) implicit none double precision::V,x V=0.5d0*x*x return end function V !================== subroutine Eigenfunction_output(N,H,E,Nout,a,b) implicit none integer,intent(in)::Nout,N double precision,intent(in)::a,b,H(1:N,1:N),E(1:N) ! a,b : output range, (a < b) double precision::x,dx,wfx,evec(1:N) integer::i,j,Np character(48)::fname Np=1000 dx=(b-a)/dble(Np) do j=1,Nout call intchar(j,fname) open(11,file=trim(fname)) evec(1:N)=H(1:N,j) do i=0,Np x=a+i*dx call wavefunction(x,N,evec,a,b,wfx) write(11,'(2e20.10e2,e14.7e2)')x,wfx,E(j) enddo close(11) enddo return end subroutine Eigenfunction_output subroutine Eigenenergy_output(N,E) implicit none integer,intent(in)::N double precision,intent(in)::E(1:N) integer::i,Nout Nout=10 if(N.lt.10)Nout=N do i=1,Nout write(6,*)i,E(i) enddo open(10,file="eigenenergy.d") do i=1,N write(10,*)i,E(i) enddo close(10) return end subroutine Eigenenergy_output subroutine Potential_output(N,V,a,b) implicit none integer,intent(in)::N double precision,intent(in)::a,b double precision,external::V integer::i double precision::h,x open(20,file="potential.d") h=(b-a)/dble(N) do i=0,N x=a+i*h write(20,*)x,V(x) enddo close(20) return end subroutine Potential_output subroutine wavefunction(x,N,evec,a,b,wf) implicit none integer,intent(in)::N double precision,intent(in)::a,b,x,evec(1:N) double precision,intent(out)::wf integer::i double precision::sinbasis external::sinbasis wf=0d0 do i=1,N wf=wf+evec(i)*sinbasis(i,x,a,b) enddo return end subroutine wavefunction subroutine Normalize(N,H) implicit none integer,intent(in)::N double precision,intent(inout)::H(1:N,1:N) integer::i,k double precision::s do k=1,N s=0d0 do i=1,N s=s+H(i,k)*H(i,k) enddo H(1:N,k)=H(1:N,k)/sqrt(s) enddo return end subroutine Normalize subroutine Hamiltonian(N,H,V,a,b,key,eps) implicit none integer,intent(in)::N,key double precision,intent(in)::a,b,eps double precision,intent(out)::H(1:N,1:N) double precision,external::V integer::i,j,ier double precision,parameter::pi=3.14159265358979323846264338327950288d0 double precision::c,kin,s double precision,external::qVij ! For parallel computation. H=0d0 kin=pi/(b-a) !$omp parallel do private(s,ier) do i=1,N call dqag_sk_Vij(qVij,i,i,a,b,eps,s,key,ier) H(i,i)=s enddo !$omp end parallel do do i=1,N c=i*kin H(i,i)=0.5d0*c*c+H(i,i) enddo !$omp parallel do private(s,ier) schedule(dynamic,1) do i=1,N do j=i+1,N call dqag_sk_Vij(qVij,i,j,a,b,eps,s,key,ier) H(i,j)=s enddo enddo !$omp end parallel do return end subroutine Hamiltonian function qVij(x,i,j,a,b) implicit none integer,intent(in)::i,j double precision,intent(in)::x,a,b double precision::qVij double precision::sini,sinj double precision,external::V double precision,parameter::pi=3.14159265358979323846264338327950288d0 sini=sqrt(2d0/(b-a))*sin(i*pi*(x-a)/(b-a)) sinj=sqrt(2d0/(b-a))*sin(j*pi*(x-a)/(b-a)) qVij=sini*V(x)*sinj return end function qVij function sinbasis(n,x,a,b) implicit none integer,intent(in)::n double precision,intent(in)::x,a,b ! Orthnormalized sin basis for [a,b] double precision,parameter::pi=3.14159265358979323846264338327950288d0 double precision::sinbasis sinbasis=sqrt(2d0/(b-a))*sin(n*pi*(x-a)/(b-a)) return end function sinbasis subroutine Diagonalize(N,A,ev) implicit none integer,intent(in)::N double precision,intent(inout)::A(1:N,1:N) double precision,intent(out)::ev(1:N) !for Hermite matrix !---sikinote integer::lda,lwork,liwork,info double precision,allocatable::work(:) double precision::qw(1:3) integer,allocatable::iwork(:) integer::qiw(1:3) character(1)::job,uplo job="V"; uplo="U"; lda=N !size query of dsyevd or zheevd parameter. call dsyevd(job,uplo,N,0,lda,0,qw,-1,qiw,-1,info) if(info.ne.0)then write(6,'(A)')" something wrong at dsyevd(size query) in Diagonalize" write(6,'(A)')" program stop" write(6,'(A)')" info --> ",info stop endif !just to be safe, +1 lwork=idint(qw(1))+1; liwork=qiw(1)+1 allocate(work(1:lwork),iwork(1:liwork)) work(1:lwork)=0d0; iwork(1:liwork)=0 !diagonalize. call dsyevd(job,uplo,N,A,lda,ev,work,lwork,iwork,liwork,info) if(info.ne.0)then write(6,'(A)')" something wrong at subroutine dsyevd in Diagonalize" write(6,'(A)')" program stop" write(6,'(A)')" info --> ",info stop endif deallocate(work,iwork) return end subroutine Diagonalize subroutine intchar(i,fname) implicit none integer,intent(in)::i character(*),intent(out)::fname if(i.le.9)then write(fname,'(A,i0,A)')"wf_00",i,".d" elseif(i.le.99)then write(fname,'(A,i0,A)')"wf_0",i,".d" elseif(i.le.999)then write(fname,'(A,i0,A)')"wf_",i,".d" else write(6,*)"Not supported, at subroutine filename." stop endif return end subroutine intchar