Commit 3e53ba6e authored by Jonas Schwab's avatar Jonas Schwab
Browse files

Update Tutorial/Solutions/Exercise_3/ALF/Analysis/

parent a883fd7b
! Copyright (C) 2019 The ALF project
!
! The ALF project 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, either version 3 of the License, or
! (at your option) any later version.
!
! The ALF project 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 License
! along with Foobar. If not, see http://www.gnu.org/licenses/.
!
! Under Section 7 of GPL version 3 we require you to fulfill the following additional terms:
!
! - It is our hope that this program makes a contribution to the scientific community. Being
! part of that community we feel that it is reasonable to require you to give an attribution
! back to the original authors if you have benefitted from this program.
! Guidelines for a proper citation can be found on the project's homepage
! http://alf.physik.uni-wuerzburg.de .
!
! - We require the preservation of the above copyright notice and this license in all original files.
!
! - We prohibit the misrepresentation of the origin of the original source files. To obtain
! the original source files please visit the homepage http://alf.physik.uni-wuerzburg.de .
!
! - If you make substantial changes to the program we require you to either consider contributing
! to the ALF project or to mark your material in a reasonable way as different from the original version.
Program ana
!--------------------------------------------------------------------
!> @author
!> ALF-project
!
!> @brief
!> Wrapper program for analyzing Monte Carlo bins.
!>
!> @details
!> Analyzes scalar observables, equal-time and timedisplaced correlation functions.
!>
!> Hand over any number of file names as command line argumens to this program.
!> If will analyze names ending in "_scal" as scalar observables, "_eq" asequal-time
!> correlation functions and "_tau" as timedisplaced correlation functions.
!>
!> For timedisplaced observables that are listen in "names_PH" in namelist "VAR_PH_names"
!> in file "parameters", it is assumed that they are symmetric in imaginary time around tau = beta/2
!>
!
!--------------------------------------------------------------------
use ana_mod
implicit none
Integer :: i, n, nargs
Character (len=64) :: name
Character (len=64), allocatable :: names(:)
nargs = COMMAND_ARGUMENT_COUNT()
allocate( names(nargs) )
do i = 1, nargs
CALL GET_COMMAND_ARGUMENT(i, name)
names(i) = name
enddo
do n=1, size(names)
name = names(n)
i = len(trim(name)) -4
if ( i > 1 ) then
if ( name(i:) == '_scal' ) then
print *, ''
print '(A,A)', "analyzing scalar observable ", name
call Cov_vec(name)
endif
endif
enddo
do n=1, size(names)
name = names(n)
i = len(trim(name)) -2
if ( i > 1 ) then
if ( name(i:) == '_eq' ) then
print *, ''
print '(A,A)', "analyzing equal time correlation ", name
call Cov_eq(name)
endif
endif
enddo
do n=1, size(names)
name = names(n)
i = len(trim(name)) -3
if ( i > 1 ) then
if ( name(i:) == '_tau' ) then
print *, ''
print '(A,A)', "analyzing time displaced correlation ", name
call Cov_tau(name)
endif
endif
enddo
deallocate( names )
end Program ana
! Copyright (C) 2019 The ALF project
!
! The ALF project 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, either version 3 of the License, or
! (at your option) any later version.
!
! The ALF project 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 License
! along with Foobar. If not, see http://www.gnu.org/licenses/.
!
! Under Section 7 of GPL version 3 we require you to fulfill the following additional terms:
!
! - It is our hope that this program makes a contribution to the scientific community. Being
! part of that community we feel that it is reasonable to require you to give an attribution
! back to the original authors if you have benefitted from this program.
! Guidelines for a proper citation can be found on the project's homepage
! http://alf.physik.uni-wuerzburg.de .
!
! - We require the preservation of the above copyright notice and this license in all original files.
!
! - We prohibit the misrepresentation of the origin of the original source files. To obtain
! the original source files please visit the homepage http://alf.physik.uni-wuerzburg.de .
!
! - If you make substantial changes to the program we require you to either consider contributing
! to the ALF project or to mark your material in a reasonable way as different from the original version.
module ana_mod
!--------------------------------------------------------------------
!> @author
!> ALF-project
!
!> @brief
!> Collection of routines for postprocessing the Monte Carlo bins
!
!--------------------------------------------------------------------
use iso_fortran_env, only: output_unit, error_unit
Use Errors
Use MyMats
Use Matrix
Use Lattices_v3, only: Unit_cell, Lattice, Make_lattice, Fourier_K_to_R
Use Predefined_Lattices, only: Predefined_Latt
contains
Subroutine read_vec(file, sgn, bins)
!--------------------------------------------------------------------
!> @author
!> ALF Collaboration
!>
!> @brief
!> Reads in bins of scalar observables from file
!>
!> @param [IN] file Character(len=64)
!> \verbatim
!> Name of file that gets read in
!> \endverbatim
!> @param [OUT] sgn Real(:)
!> \verbatimam
!> Sign of bins
!> \endverbatim
!> @param [OUT] bins Complex(:,:)
!> \verbatim
!> Monte Carlo bins
!> \endverbatim
!-------------------------------------------------------------------
Implicit none
Character (len=64), intent(in) :: file
Real (Kind=Kind(0.d0)), allocatable, intent(out) :: sgn(:)
Complex (Kind=Kind(0.d0)), pointer, intent(out) :: bins(:,:)
Integer :: N, N1, I, Nobs, Nbins, stat
Real (Kind=Kind(0.d0)) :: X
Complex (Kind=Kind(0.d0)), Allocatable :: tmp(:)
open(Unit=10, File=file, status="old", action='read')
read(10,*) NOBS
NOBS = NOBS-1
allocate(tmp(NOBS))
rewind(10)
Nbins = 0
do
read(10, *, iostat=stat)
if (stat /= 0) exit
Nbins = Nbins + 1
enddo
rewind(10)
allocate(bins(NOBS,Nbins))
allocate(sgn(Nbins))
do N = 1,Nbins
read(10,*) N1, (tmp(I), I=1,size(Tmp,1)),X
bins(:,N) = Tmp(:)
sgn(N) = X
enddo
close(10)
deallocate(tmp)
End Subroutine read_vec
!==============================================================================
Subroutine read_latt(file, sgn, bins, bins0, Latt, Latt_unit, dtau, Channel)
!--------------------------------------------------------------------
!> @author
!> ALF Collaboration
!>
!> @brief
!> Reads in bins of lattice-type observables (both equal time and timedisplaced) from file
!>
!> @param [IN] file Character(len=64)
!> \verbatim
!> Name of file that gets read in
!> \endverbatim
!> @param [OUT] sgn Real(:)
!> \verbatimam
!> Sign of bins
!> \endverbatim
!> @param [OUT] bins Complex(:,:,:,:,:)
!> \verbatim
!> Monte Carlo bins of correlation
!> @param [OUT] bins0 Complex(:,:)
!> \verbatim
!> Monte Carlo bins of background
!> \endverbatim
!> @param [OUT] Latt Type(Lattice)
!> \verbatim
!> Bravais lattice
!> \endverbatim
!> @param [OUT] Latt_unit Type(Unit_cell)
!> \verbatim
!> Unit cell
!> \endverbatim
!> @param [OUT] dtau Real
!> \verbatim
!> Size of imaginary time step
!> \endverbatim
!> @param [OUT] file Character(len=2)
!> \verbatim
!> MaxEnt Channel. Relevant for timedisplaced correlation.
!> \endverbatim
!-------------------------------------------------------------------
Implicit none
Character (len=64), intent(in) :: file
Real (Kind=Kind(0.d0)), allocatable, intent(out) :: sgn(:)
Complex (Kind=Kind(0.d0)), pointer , intent(out) :: bins(:,:,:,:,:)
Complex (Kind=Kind(0.d0)), pointer , intent(out) :: bins0(:,:)
Type (Lattice) , intent(out) :: Latt
Type (Unit_cell) , intent(out) :: Latt_unit
Real (Kind=Kind(0.d0)) , intent(out) :: dtau
Character (len=2) , intent(out) :: Channel
Character (len=64) :: file_aux, str_temp1
Integer, allocatable :: List(:,:), Invlist(:,:) ! For orbital structure of Unit cell
Integer :: no, no1, n, nt, nb, Ntau, Ndim, Nbins, stat, Ndim_unit
Real(Kind=Kind(0.d0)) :: X
Real(Kind=Kind(0.d0)), allocatable :: Xk_p(:,:)
Real(Kind=Kind(0.d0)) :: x_p(2), a1_p(2), a2_p(2), L1_p(2), L2_p(2)
logical :: file_exists
Integer :: L1, L2
Character (len=64) :: Model, Lattice_type
NAMELIST /VAR_Lattice/ L1, L2, Lattice_type, Model
write(file_aux, '(A,A)') trim(file), "_info"
inquire(file=file_aux, exist=file_exists)
if(file_exists) then
open(Unit=10, File=file_aux, status="old", action='read')
11 format(A22, A)
12 format(A22, I10)
13 format(A22, *(E26.17E3))
read(10, *)
read(10, 11) str_temp1, Channel
read(10, 12) str_temp1, Ntau
read(10, 13) str_temp1, dtau
read(10, *)
read(10, 12) str_temp1, Latt%N
read(10, 13) str_temp1, L1_p
read(10, 13) str_temp1, L2_p
read(10, 13) str_temp1, a1_p
read(10, 13) str_temp1, a2_p
read(10, *)
read(10, 12) str_temp1, Latt_unit%N_coord
read(10, 12) str_temp1, Latt_unit%Norb
read(10, 12) str_temp1, Ndim_unit
allocate(Latt_unit%Orb_pos_p(Latt_unit%Norb, Ndim_unit))
do no = 1, Latt_unit%Norb
read(10, 13) str_temp1, Latt_unit%Orb_pos_p(no,:)
enddo
close(10)
Call Make_Lattice(L1_p, L2_p, a1_p, a2_p, Latt)
Ndim = Latt%N*Latt_Unit%Norb
else
Channel = '--'
open(Unit=10, File='parameters', status="old", action='read')
read(10, NML=VAR_lattice)
close(10)
Call Predefined_Latt(Lattice_type, L1, L2, Ndim, List, Invlist, Latt, Latt_Unit)
open(Unit=10, File=file, status="old", action='read')
Read(10, *, iostat=stat) X, Latt_unit%Norb, Latt%N, Ntau, dtau
if (stat /= 0) then
rewind(10)
Ntau = 1
dtau = -1.d0
Read(10, *) X, Latt_unit%Norb, Latt%N
endif
close(10)
endif
! Determine the number of bins.
open(Unit=10, File=file, status="old", action='read')
nbins = 0
do
if(Ntau == 1) then
read(10, *, iostat=stat) X, Latt_unit%Norb, Latt%N
else
read(10, *, iostat=stat) X, Latt_unit%Norb, Latt%N, Ntau, dtau
endif
if (stat /= 0) exit
Do no = 1, Latt_unit%Norb
Read(10,*)
enddo
do n = 1, Latt%N
Read(10,*)
do nt = 1, Ntau
do no = 1, Latt_unit%Norb
do no1 = 1, Latt_unit%Norb
read(10,*)
enddo
enddo
enddo
enddo
nbins = nbins + 1
enddo
rewind(10)
! Allocate space
Allocate(bins(Latt%N, Ntau, Latt_unit%Norb, Latt_unit%Norb, Nbins))
Allocate(sgn(Nbins), Xk_p(2, Latt%N), bins0(Latt_unit%Norb, Nbins))
do nb = 1, nbins
if(Ntau == 1) then
read(10, *) sgn(nb), Latt_unit%Norb, Latt%N
else
read(10, *) sgn(nb), Latt_unit%Norb, Latt%N, Ntau, dtau
endif
Do no = 1, Latt_unit%Norb
Read(10,*) bins0(no,nb)
Enddo
do n = 1, Latt%N
Read(10,*) Xk_p(1,n), Xk_p(2,n)
do nt = 1, Ntau
do no = 1, Latt_unit%Norb
do no1 = 1, Latt_unit%Norb
read(10,*) bins(n,nt,no,no1,nb)
enddo
enddo
enddo
enddo
enddo
close(10)
do n = 1, Latt%N
x_p = dble(Latt%listk(n,1))*Latt%b1_p + dble(Latt%listk(n,2))*Latt%b2_p
x = (x_p(1)-Xk_p(1,n))**2 + (x_p(2)-Xk_p(2,n))**2
if ( x > 0.00001 ) then
Write(error_unit,*) "Error in read_latt: momenta do not fit", x, x_p, Xk_p(1,n)
error stop
endif
enddo
End Subroutine read_latt
!==============================================================================
Subroutine jack(func, data, N_skip, N_rebin, best, err)
!--------------------------------------------------------------------
!> @author
!> ALF Collaboration
!>
!> @brief
!> Performs jackknife error Analysis
!>
!-------------------------------------------------------------------
Implicit none
Complex (Kind=Kind(0.d0)), External :: func
Complex (Kind=Kind(0.d0)), intent(in) :: data(:,:)
Integer, intent(in) :: N_skip, N_rebin
Complex (Kind=Kind(0.d0)), intent(out) :: best, err
Complex (Kind=Kind(0.d0)), allocatable :: data_r(:,:), j(:), X(:), X2(:)
Real (Kind=Kind(0.d0)), allocatable :: Rhelp(:)
real (Kind=Kind(0.d0)) :: XM, XERR
Integer :: Nobs, Nbins, Nbins_r
Integer :: N_r, N1, N
Nobs = size(data,1)
Nbins = size(data,2)
! Skipping and Rebinning
Nbins_r = (Nbins-N_skip)/N_rebin
Allocate( data_r(Nobs,Nbins_r), X(Nobs) )
N = N_skip
Do N_r = 1, Nbins_r
X(:) = cmplx(0.D0,0.d0,kind(0.d0))
Do N1 = 1, N_rebin
N = N + 1
X(:) = X(:) + data(:,N)
enddo
data_r(:,N_r) = X(:)/dble(N_rebin)
enddo
! Creating jackknife bins
Allocate( j(Nbins_r), X2(Nobs) )
X(:) = cmplx(0.D0,0.d0,kind(0.d0))
do N = 1, Nbins_r
X(:) = X(:) + data_r(:,N)
enddo
do N = 1, Nbins_r
X2(:) = (X(:) - data_r(:,N))/dble(Nbins_r-1)
j(N) = func(X2)
enddo
! Calculate standard deviation of jackknife bins
Allocate (Rhelp(Nbins_r))
do N = 1, Nbins_r
Rhelp(N) = dble(j(N))
enddo
call errcalc(Rhelp, xm, xerr)
best = cmplx(xm , 0.d0, kind(0.D0))
err = cmplx(xerr, 0.d0, kind(0.D0))
do N = 1, Nbins_r
Rhelp(N) = aimag(j(N))
enddo
call errcalc(Rhelp, xm, xerr)
best = best + cmplx( 0.d0, xm, kind(0.D0) )
err = err + cmplx( 0.d0, xerr, kind(0.D0) )
err = err * dble(Nbins_r)
deallocate( data_r, X, j, X2, Rhelp )
End Subroutine jack
!==============================================================================
subroutine auto(func, data, N_skip, res)
!--------------------------------------------------------------------
!> @author
!> ALF Collaboration
!>
!> @brief
!> Calculates autocorrelation
!-------------------------------------------------------------------
implicit none
complex (Kind=Kind(0.d0)), External :: func
complex (Kind=Kind(0.d0)), intent(in) :: data(:,:)
Integer, intent(in) :: N_skip
Real (Kind=Kind(0.d0)), intent(inout) :: res(:)
Integer :: N_obs, N_bins, N_auto, ntau, nt
complex (Kind=Kind(0.d0)), allocatable :: data1(:), data2(:)
complex (Kind=Kind(0.d0)) :: mean, X1, X2
N_obs = size(data,1)
N_bins = size(data,2) - N_skip
N_auto = size(res)
allocate( data1(N_obs), data2(N_obs) )
do ntau = 1, N_auto
X1 = 0.0
X2 = 0.0
mean = 0.0
do nt = 1, N_bins - ntau
data1(:) = data(:,nt+N_skip)
mean = mean + func(data1)
enddo
mean = mean / dble(N_bins - ntau)
mean = func(mean)
do nt = 1, N_bins - ntau
data1(:) = data(:,nt+N_skip)
data2(:) = data(:,nt+N_skip+ntau)
X1 = X1 + (func(data1)-mean)*(func(data2)-mean)
X2 = X2 + (func(data1)-mean)*(func(data1)-mean)
enddo
! X1 = X1 / dble(N_bins - ntau)
! X2 = X2 / dble(N_bins - ntau)
Res(ntau) = dble( X1/X2 )
enddo
end subroutine auto
!==============================================================================
subroutine Cov_tau(file)
!--------------------------------------------------------------------
!> @author
!> ALF-project
!
!> @brief
!> Analysis of imaginary time displaced correlation functions.
!>
!> @param [IN] file Character(len=64)
!> \verbatim
!> Name of file that gets read in
!> \endverbatim
!> @param [IN] PartHole logical
!> \verbatimam
!> If true, it is assumed that observable is symmetric in imaginary time around tau = beta/2
!> \endverbatim
!
!--------------------------------------------------------------------
Implicit none
Character (len=64), intent(in) :: file
Character (len=64) :: name_obs
Real (Kind=Kind(0.d0)), allocatable :: sgn(:)
Complex (Kind=Kind(0.d0)), pointer :: Bins_raw(:,:,:,:,:), Bins0_raw(:,:)
Type (Lattice) :: Latt
Type (Unit_cell) :: Latt_unit
real (Kind=Kind(0.d0)) :: dtau
Character (len=2) :: Channel
Integer :: i
i = len(trim(file)) - 4
name_obs = file(:i)
call read_latt(file, sgn, bins_raw, bins0_raw, Latt, Latt_unit, dtau, Channel)
call ana_tau(name_obs, sgn, bins_raw, bins0_raw, Latt, Latt_unit, dtau, Channel)
end subroutine Cov_tau
!==============================================================================
Subroutine ana_tau(name_obs, sgn, bins_raw, bins0_raw, Latt, Latt_unit, dtau, Channel)
Implicit none
Character (len=64), intent(in) :: name_obs
Real (Kind=Kind(0.d0)), allocatable, intent(in) :: sgn(:)
Complex (Kind=Kind(0.d0)), pointer , intent(in) :: Bins_raw(:,:,:,:,:)
Complex (Kind=Kind(0.d0)), pointer , intent(in) :: Bins0_raw(:,:)
Type (Lattice) , intent(in) :: Latt
Type (Unit_cell) , intent(in) :: Latt_unit
Real (Kind=Kind(0.d0)) , intent(in) :: dtau
Character (len=2) , intent(in) :: Channel
Logical :: PartHole