! Module with some useful functions and global variables. module utils implicit none integer, parameter :: n = 400 ! Problem subdivisions real(8), parameter :: L = 20.0 ! Problem size real(8), parameter :: step = L/real(n-1,8) ! Subdivision size real(8), parameter :: pi = 4.0*atan2(1.0,1.0) real(8), parameter :: s0 = 2.0 ** 0.25 real(8), parameter :: x0 = 1.0 ! Double-well positions real(8), parameter :: v0 = 1000.0 ! Double-well height contains subroutine diago(n,d,e,vp) implicit none integer, intent(in) :: n real(8), dimension(n), intent(inout) :: d real(8), dimension(n), intent(in) :: e real(8), dimension(n,n), intent(out) :: vp integer :: info, l real(8), dimension(2*n+2) :: work vp = 0._8 ; forall(l = 1:n) vp(l,l) = 1._8 call dsteqr('V', n, d, e, vp, n, work, info) if(info == 0) then print '("Diagonalisation effectuee")' else print '("Echec de diagonalisation : Info = ",i3)', info endif end subroutine diago integer function factorial(n) implicit none integer, intent(in) :: n integer :: i factorial = 1 do i = 2, n factorial = factorial * i enddo end real(8) function hermite0(x) implicit none real(8), intent(in) :: x hermite0 = 1 end real(8) function hermite1(x) implicit none real(8), intent(in) :: x hermite1 = 2.0 * x end real(8) function hermite2(x) implicit none real(8), intent(in) :: x hermite2 = 4.0 * x**2.0 - 2.0 end real(8) function hermite3(x) implicit none real(8), intent(in) :: x hermite3 = 8.0 * x**3.0 - 12.0 * x end real(8) function hermite4(x) implicit none real(8), intent(in) :: x hermite4 = 16.0 * x**4.0 - 48.0 * x**2.0 + 12.0 end ! Returns the theoretical probability density for the harmonic potential real(8) function theoretic_harmonic_density(x,p) implicit none real(8), intent(in) :: x integer :: p theoretic_harmonic_density = & exp(-((x/s0)**2.0)) / (s0 * sqrt(pi) * (2.0 ** real(p,8)) * real(factorial(p),8)) if (p == 0) then theoretic_harmonic_density = theoretic_harmonic_density * (hermite0(x/s0) ** 2.0) else if(p == 1) then theoretic_harmonic_density = theoretic_harmonic_density * (hermite1(x/s0) ** 2.0) else if(p == 2) then theoretic_harmonic_density = theoretic_harmonic_density * (hermite2(x/s0) ** 2.0) else if(p == 3) then theoretic_harmonic_density = theoretic_harmonic_density * (hermite3(x/s0) ** 2.0) else if(p == 4) then theoretic_harmonic_density = theoretic_harmonic_density * (hermite4(x/s0) ** 2.0) endif end function theoretic_harmonic_density end module utils ! Initializes array x subroutine init_positions(x) use utils implicit none real(8), dimension(n), intent(inout) :: x integer :: i forall(i = 1:n) x(i) = real(i-1,8)*step - L/2.0 end subroutine init_positions ! Initializes array v with an infinite well potential subroutine well_potential(v,x) use utils implicit none real(8), dimension(n), intent(inout) :: v,x integer :: i forall(i = 1:n) v(i) = 0 end subroutine well_potential ! Initializes array v with a harmonic potential subroutine harmonic_potential(v,x) use utils implicit none real(8), dimension(n), intent(inout) :: v,x integer :: i forall(i = 1:n) v(i) = x(i) ** 2.0 / 2.0 end subroutine harmonic_potential ! Initializes array v with a double-well potential subroutine doublewell_potential(v,x) use utils implicit none real(8), dimension(n), intent(inout) :: v,x integer :: i forall(i = 1:n) v(i) = v0 * (((x(i) / x0) ** 2.0 - 1.0) ** 2.0) end subroutine doublewell_potential ! Initializes the diagonal d and the lower diagonal e subroutine init_diagonal(d,e,v) use utils implicit none real(8), dimension(n) , intent(inout) :: d real(8), dimension(n-1), intent(inout) :: e real(8), dimension(n) , intent(in) :: v integer :: i forall(i = 1:n) d(i) = 2.0 / (step * step) + v(i) forall(i = 1:(n-1)) e(i) = - 1.0 / (step * step) end subroutine init_diagonal ! Main program program main use utils implicit none real(8), dimension(n) :: v,x,d ! Potential, positions, H's diagonal real(8), dimension(n-1) :: e ! H's lower diagonal real(8), dimension(n,n) :: vp ! Eigen vectors integer :: i,j ! Initialization call init_positions(x) call doublewell_potential(v,x) call init_diagonal(d,e,v) ! Computation call diago(n,d,e,vp) ! Output open(1, file='results/eigenvalues.out') do i = 1, n write(1,*) (i-1), d(i), (real(i-1,8) + 0.5) * sqrt(2.0) enddo close(1) open(2, file='results/eigenvectors.out') do j = 1, n write(2,*) j, vp(j,3) enddo close(2) open(3, file='results/densities.out') do i = 1, n write(3,*) x(i),(vp(i,1) ** 2.0) * 2.0 / step + d(1),(vp(i,2) ** 2.0) * 2.0 / step + d(2) enddo close(3) end program main