!=========================================================== ! ! File observables.f90 ! Compile: gfortran observables.f90 -o o ! Usage: ./o ! ! Read in a file with a wavefunction in the format of psi.dat: ! # E= .... ! x1 psi(x1) ! x2 psi(x2) ! ............ ! ! Outputs expectation values: ! normalization Energy

Dx Dp DxDp ! where Dx = sqrt(-^2) Dp = sqrt(-

^2) ! DxDp = Dx * Dp ! !=========================================================== program observables_expectation use, intrinsic :: iso_fortran_env implicit none integer , parameter :: d = real64 real(d), allocatable :: xstep(:),psi(:),obs(:) integer :: Nx,i,upsi real(d) :: xav, pav, x2av, p2av, Dx, Dp, DxDp real(d) :: energy,h,norm character(20) :: psifile,scratch !the first argument of the command line must be the path !to the file with the wavefunction. (GNU fortran extension...) if( iargc() /= 1) stop 'Usage: o ' call getarg(1,psifile) !If the file does not exist, we go to label 100 (stop): open(newunit=upsi,file=psifile,status='OLD',err=100) print *,"# reading wavefunction from file:", psifile !------------------------------------------------------- !We count the number of x-points in the file !we read the first comment line from the file: read(upsi,*) scratch,scratch,energy Nx=0 do read(upsi,*,end=101) Nx = Nx + 1 end do 101 continue if(mod(Nx,2) .eq. 0) Nx = Nx - 1 !we need odd # points for integrate allocate(xstep(Nx),psi(Nx),obs(Nx)) !------------------------------------------------------- !Input data: psi(x) rewind(upsi) !read again from the beginning read(upsi,*) scratch,scratch,energy !comment line in file do i = 1, Nx read(upsi,*) xstep(i),psi(i) enddo h = (xstep(Nx)-xstep(1))/(Nx-1) !------------------------------------------------------- !Calculate: !---------- norm: do i = 1,Nx obs(i) = psi(i)*psi(i) enddo norm = integrate(obs,h,Nx) !---------- : do i = 1,Nx obs(i) = xstep(i)*psi(i)*psi(i) enddo xav = integrate(obs,h,Nx)/norm !----------

/i : use quadratic polynomial at boundaries (see Derivatives.nb) !obs(1) = psi(1)*(psi(2)-psi(1))/h !naive derivative obs(1) = psi(1)*(-3.0_d*psi(1)+4.0_d*psi(2)-psi(3))/(2*h) do i = 2,Nx-1 obs(i) = psi(i)*(psi(i+1)-psi(i-1))/(2.0_d*h) enddo !obs(Nx) = psi(Nx)*(psi(Nx)-psi(Nx-1))/h !naive obs(Nx) = psi(Nx)*(psi(Nx-2)-4.0_d*psi(Nx-1)+3.0_d*psi(Nx))/(2*h) !better: O(h^2) (See Derivatives.nb) pav = -integrate(obs,h,Nx)/norm !--------- do i = 1,Nx obs(i) = xstep(i)*xstep(i)*psi(i)*psi(i) enddo x2av = integrate(obs,h,Nx)/norm !-------- !obs(1) = psi(1)*(psi(3)-2.0_d*psi(2)+psi(1))/(h*h) !naive O(h) obs(1) = psi(1 )*(2.0_d*psi(1 )-5.0_d*psi(2 )+4.0_d*psi(3 )-psi(4 ))/(h*h) !better: O(h^2) (See Derivatives.nb) do i = 2,Nx-1 obs(i) = psi(i )*(psi(i+1)-2.0_d*psi(i)+psi(i-1))/(h*h) enddo !obs(Nx) = psi(Nx)*(psi(Nx)-2.0_d*psi(Nx-1)+psi(Nx-2))/(h*h) obs(Nx) = psi(Nx)*(2.0_d*psi(Nx)-5.0_d*psi(Nx-1)+4.0_d*psi(Nx-2)-psi(Nx-3))/(h*h) !better: O(h^2) (see Derivatives.nb) p2av = -integrate(obs,h,Nx)/norm !-------- Dx Dx = sqrt(x2av - xav*xav) !-------- Dp Dp = sqrt(p2av - pav*pav) !-------- Dx . Dp DxDp = Dx*Dp !print results: print *,'# norm E

/i Dx Dp DxDp' print '(10G25.17)',norm,energy,xav,pav,x2av,p2av,Dx,Dp,DxDp stop !normal execution ends here. Error messages follow 100 stop 'Cannot open file' !----------------------------------------------------------- contains !----------------------------------------------------------- !=========================================================== ! !Simpson's rule to integrate psi(x). !For n intervals of width dx (n even) !Simpson's rule is: !int(f(x)dx) = ! (dx/3)*(f(x_0)+4 f(x_1)+2 f(x_2)+...+4 f(x_{n-1})+f(x_n)) ! !Input: Discrete values of function psi(Nx) ! Integration step dx !Returns: Integral(psi(x)psi(x) dx) !=========================================================== function integrate(psi,dx,Nx) real(d) :: integrate integer :: Nx real(d) :: psi(Nx),dx real(d) :: int integer :: i !----- zeroth order point: i = 1 int = psi(i) !----- odd order points (i=k+1 is even): do i = 2,Nx-1,2 int = int + 4.0_d*psi(i) enddo !----- even order points: do i = 3,Nx-2,2 int = int + 2.0_d*psi(i) enddo !----- last point: i = Nx int = int + psi(i) !----- measure normalization: int = int*dx/3.0_d !----- final result: integrate = int end function integrate end program observables_expectation !=========================================================== ! --------------------------------------------------------------------- ! Copyright by Konstantinos N. Anagnostopoulos (2004-2021) ! Physics Dept., National Technical University, ! konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant ! ! This program is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation, version 3 of the License. ! ! This program is distributed in the hope that it will be useful, but ! WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! General Public License for more details. ! ! You should have received a copy of the GNU General Public Liense along ! with this program. If not, see . ! -----------------------------------------------------------------------