Commit e23f723a authored by Jonas Schwab's avatar Jonas Schwab

Update Tutorial/Solutions/Exercise_1/ALF/Analysis/

parent 5cc4b256
.PHONY : all tidy clean
BINS= cov_scal.out cov_tau.out cov_tau_ph.out Max_SAC.out cov_eq.out
OBJS= cov_scal.o cov_tau.o cov_tau_ph.o Max_SAC.o cov_eq.o
OBJS1 = Predefined_Latt_mod.o
MODS=predefined_lattices.mod
BINS=cov_eq.out cov_scal.out cov_tau.out cov_tau_ph.out Max_SAC.out ana.out #repack_latt.out
OBJS=cov_eq.o cov_scal.o cov_tau.o cov_tau_ph.o Max_SAC.o ana.o
OBJS1=Predefined_Latt_mod.o ana_mod.o
MODS=predefined_lattices.mod ana_mod.mod
all: $(OBJS1) $(BINS)
cov_tau_ph.F90 : cov_tau.F90
cpp -w -P -DPartHole cov_tau.F90 -o cov_tau_ph.F90
all: $(BINS)
Predefined_Latt_mod.o: ../Prog/Predefined_Latt_mod.F90
$(ALF_FC) -c -o $@ $(ALF_FLAGS_ANA) $<
cov_tau_ph.o: cov_tau.F90
$(ALF_FC) -c -o $@ $(ALF_FLAGS_ANA) -DPartHole $<
%.o: %.F90
$(ALF_FC) -c $*.F90 $(ALF_FLAGS_ANA)
$(ALF_FC) -c -o $@ $(ALF_FLAGS_ANA) $<
%.out: %.F90
$(ALF_FC) -c $(ALF_FLAGS_ANA) $<
$(ALF_FC) -o $@ $*.o $(OBJS1) $(ALF_LIB)
%.out: $(OBJS1) %.o
$(ALF_FC) -o $@ $^ $(ALF_LIB)
tidy:
rm -f $(OBJS) $(OBJS1) $(MODS) cov_tau_ph.F90
......
! Copyright (C) 2016-2019 The ALF project
! Copyright (C) 2016-2020 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
......@@ -42,9 +42,8 @@
!--------------------------------------------------------------------
Use MaxEnt_stoch_mod
use iso_fortran_env, only: output_unit, error_unit
Implicit None
!Implicit Real (Kind=Kind(0.d0)) (A-G,O-Z)
!Implicit Integer (H-N)
Interface
Subroutine Rescale ( XCOV, XQMC,XTAU, Ntau_st, Ntau_en, Tolerance, Ntau)
......@@ -67,7 +66,7 @@
Integer :: Ngamma, Ndis, NBins, NSweeps, Nwarm, N_alpha, N_cov
Integer :: N_skip, N_rebin, Norb
Integer :: N_skip, N_rebin, N_Back, N_auto, Norb
Real (Kind=Kind(0.d0)) :: OM_st, OM_en, alpha_st, R, Tolerance
Logical :: Checkpoint
Character (Len=2) :: Channel
......@@ -78,9 +77,9 @@
Real (Kind=Kind(0.d0)) :: Zero
NAMELIST /VAR_Max_Stoch/ Ngamma, Ndis, NBins, NSweeps, Nwarm, N_alpha, &
& OM_st, OM_en, alpha_st, R, Checkpoint, Channel, Tolerance
& OM_st, OM_en, alpha_st, R, Checkpoint, Tolerance
NAMELIST /VAR_errors/ N_skip, N_rebin, N_cov
NAMELIST /VAR_errors/ N_skip, N_rebin, N_cov, N_Back, N_auto
open(unit=30,file='parameters',status='old',action='read', iostat=io_error)
if (io_error.eq.0) then
......@@ -92,26 +91,9 @@
endif
close(30)
Open(unit=50,File='Info_MaxEnt',Status="unknown")
write(50,*) 'Channel :: ', Channel
If (Channel == "ph" ) then
Write(50,*) 'Om_start is set to zero. PH channel corresponds to symmetric data '
Om_st = 0.d0
endif
Write(50, "('Covariance :: ',I2)") N_cov
Write(50, "('Checkpoint :: ',L )") Checkpoint
Write(50, "('Om_st, Om_en :: ',2x,F12.6,2x,F12.6)") Om_st, Om_en
Write(50, "('Delta Om :: ',2x,F12.6)") (Om_en - Om_st)/real(Ndis,kind(0.d0))
Write(50, "('Bins, Sweeps, Warm :: ',2x,I4,2x,I4,2x,I4)") NBins, NSweeps, Nwarm
If (N_alpha <= 10 ) then
Write(error_unit,*) 'Not enough temperatures: N_alpha has to be bigger than 10'
error stop 1
Endif
Write(50, "('N_Alpha, Alpha_st,R:: ',2x,I4,F12.6,2x,F12.6)") N_alpha, alpha_st, R
open (unit=10,File="g_dat", status="unknown")
read(10,*) ntau, nbin_qmc, Beta, Norb
read(10,*) ntau, nbin_qmc, Beta, Norb, Channel
Allocate ( XCOV(NTAU,NTAU), XQMC(NTAU),XTAU(NTAU) )
XCOV = 0.d0
Do nt = 1,NTAU
......@@ -128,6 +110,23 @@
close(10)
dtau = Xtau(2) - Xtau(1)
Open(unit=50,File='Info_MaxEnt',Status="unknown")
write(50,*) 'Channel :: ', Channel
If (Channel == "ph" ) then
Write(50,*) 'Om_start is set to zero. PH channel corresponds to symmetric data '
Om_st = 0.d0
endif
Write(50, "('Covariance :: ',I2)") N_cov
Write(50, "('Checkpoint :: ',L1)") Checkpoint
Write(50, "('Om_st, Om_en :: ',2x,F12.6,2x,F12.6)") Om_st, Om_en
Write(50, "('Delta Om :: ',2x,F12.6)") (Om_en - Om_st)/real(Ndis,kind(0.d0))
Write(50, "('Bins, Sweeps, Warm :: ',2x,I4,2x,I4,2x,I4)") NBins, NSweeps, Nwarm
If (N_alpha <= 10 ) then
Write(error_unit,*) 'Not enough temperatures: N_alpha has to be bigger than 10'
error stop 1
Endif
Write(50, "('N_Alpha, Alpha_st,R:: ',2x,I4,F12.6,2x,F12.6)") N_alpha, alpha_st, R
Zero = 1.D-10
pi = acos(-1.d0)
Ntau_st = 1
......@@ -138,7 +137,7 @@
Case ("PP")
xmom1 = 2.d0* pi * xqmc(1)
Case ("P")
xmom1 = pi * real(norb,Kind(0.d0))
xmom1 = pi * ( xqmc(1) + xqmc(ntau) )
! Remove the tau = beta point from the data since it is correlated
! due to the sum rule with the tau=0 data point. Also if the tau = 0
! data point has no fluctations (due to particle-hole symmetry for instance)
......@@ -222,9 +221,9 @@
If ( .not. Checkpoint ) then
Command = "rm dump*"
Call System (Command)
Call EXECUTE_COMMAND_LINE(Command)
Command = "ls"
Call System (Command)
Call EXECUTE_COMMAND_LINE(Command)
endif
Open (Unit=10,File="energies",status="unknown")
......@@ -295,7 +294,7 @@
Z = Z + A(nwp)/cmplx( om - omp, delta, kind(0.d0))
enddo
Z = Z * dom
write(43,"('X'2x,F14.7,2x,F16.8,2x,F16.8)" ) xom(nw), dble(Z), -Aimag(Z)/pi
write(43,"('X',2x,F14.7,2x,F16.8,2x,F16.8)") xom(nw), dble(Z), -Aimag(Z)/pi
enddo
close(43)
......
! 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
This diff is collapsed.
......@@ -44,7 +44,7 @@
Use MyMats
Use Matrix
Use Lattices_v3
Use Predefined_Lattices
Use Predefined_Lattices
use iso_fortran_env, only: output_unit, error_unit
Implicit none
......@@ -76,6 +76,8 @@
NAMELIST /VAR_errors/ n_skip, N_rebin, N_Cov, N_Back, N_auto
N_skip = 1
N_rebin = 1
N_Back = 1
N_auto = 0
OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr)
......@@ -114,7 +116,7 @@
Close(10)
Write(6,*) "# of bins: ", Nbins
nbins = Nbins - n_skip
Write(6,*) "Effective # of bins: ", Nbins
Write(6,*) "Effective # of bins: ", Nbins/N_rebin
N_auto=min(N_auto,Nbins/3)
if(Nbins <= 1) then
write (error_unit,*) "Effective # of bins smaller than 2. Analysis impossible!"
......
......@@ -61,6 +61,8 @@
NAMELIST /VAR_errors/ n_skip, N_rebin, N_Cov, N_Back, N_auto
N_skip = 1
N_rebin = 1
N_auto=0
OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr)
IF (ierr /= 0) THEN
......@@ -109,7 +111,7 @@
endif
OPEN (UNIT=21, FILE='Var_scalJ', STATUS='unknown')
WRITE(21,*) 'Effective number of bins, and bins: ', Nbins_eff, Nbins
WRITE(21,*) 'Effective number of bins, and bins: ', Nbins_eff/N_rebin, Nbins
ALLOCATE (EN(Nbins_eff), SIGN(Nbins_eff))
DO IOBS = 1,NOBS
WRITE(21,*)
......
......@@ -63,13 +63,14 @@
Real (Kind=Kind(0.d0)), allocatable :: Xk_p(:,:)
Character (len=64) :: File_out
NAMELIST /VAR_errors/ n_skip, N_rebin, N_Cov, N_Back, N_auto
NAMELIST /VAR_errors/ N_skip, N_rebin, N_Cov, N_Back, N_auto
N_skip = 1
N_rebin = 1
N_Back = 1
N_auto = 0
N_Cov = 0
OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr)
IF (ierr /= 0) THEN
WRITE(error_unit,*) 'unable to open <parameters>',ierr
......
! Copyright (C) 2016 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:
!
! 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:
!
! - 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:
!
! - 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 Cov_tau
!--------------------------------------------------------------------
!> @author
!> ALF-project
!
!> @brief
!> Analysis of imaginary time displaced correlation functions.
!
!--------------------------------------------------------------------
Use Errors
Use MyMats
Use Matrix
use iso_fortran_env, only: output_unit, error_unit
Implicit none
Integer :: Nunit, Norb, N_auto
Integer :: no, no1, n, nbins, n_skip, nb, NT, NT1, Lt, N_rebin, N_cov, ierr, N_Back
Integer :: Lt_eff
real (Kind=Kind(0.d0)):: X, Y, dtau, X_diag
Complex (Kind=Kind(0.d0)), allocatable :: Xmean(:), Xcov(:,:)
Complex (Kind=Kind(0.d0)) :: Zmean, Zerr
Complex (Kind=Kind(0.d0)) :: Z, Z_diag
Real (Kind=Kind(0.d0)) :: Zero=1.D-8
Real (Kind=Kind(0.d0)), allocatable :: Phase(:)
Complex (Kind=Kind(0.d0)), allocatable :: PhaseI(:)
Complex (Kind=Kind(0.d0)), allocatable :: Bins(:,:,:), Bins_chi(:,:), OneBin(:,:,:)
Complex (Kind=Kind(0.d0)), allocatable :: Bins0(:,:)
Complex (Kind=Kind(0.d0)), allocatable :: V_help(:,:)
Real (Kind=Kind(0.d0)), allocatable :: Xk_p(:,:)
Character (len=64) :: File_out
NAMELIST /VAR_errors/ n_skip, N_rebin, N_Cov, N_Back, N_auto
N_Back = 1
N_auto = 0
OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr)
IF (ierr /= 0) THEN
WRITE(error_unit,*) 'unable to open <parameters>',ierr
error stop 1
END IF
READ(5,NML=VAR_errors)
CLOSE(5)
! Determine the number of bins.
Open ( Unit=10, File="intau", status="unknown" )
nbins = 0
do
Read(10,*,End=10) X,Norb,Nunit, LT, dtau
Do no = 1,Norb
Read(10,*) Z
enddo
do n = 1,Nunit
Read(10,*) X,Y
do nt = 1,LT
do no = 1,Norb
do no1 = 1,Norb
read(10,*) Z
enddo
enddo
enddo
enddo
!Write(6,*) nbins
nbins = nbins + 1
enddo
10 continue
Close(10)
Write(6,*) "# of bins: ", Nbins
nbins = Nbins - n_skip
Write(6,*) "Effective # of bins: ", Nbins
if(Nbins <= 1) then
write (error_unit,*) "Effective # of bins smaller than 2. Analysis impossible!"
error stop 1
endif
if (mod(Lt-1,2) == 0 ) then
Lt_eff = (Lt -1 ) /2 + 1
else
Lt_eff = Lt/2
endif
! Allocate space
Allocate ( bins(Nunit,Lt_eff,Nbins), Phase(Nbins),PhaseI(Nbins), Xk_p(2,Nunit), &
& V_help(Lt_eff,Nbins), bins0(Nbins,Norb))
Allocate ( Bins_chi(Nunit,Nbins) )
Allocate ( OneBin(Lt,Norb,Norb) )
Allocate (Xmean(Lt_eff), Xcov(Lt_eff,Lt_eff))
bins = 0.d0
bins0 = cmplx(0.d0,0.d0,Kind(0.d0))
Open ( Unit=10, File="intau", status="unknown" )
do nb = 1, nbins + n_skip
if (nb > n_skip ) then
Read(10,*,End=10) Phase(nb-n_skip),no,no1,n, X
PhaseI(nb-n_skip) = cmplx(Phase(nb-n_skip),0.d0,Kind(0.d0))
Z_diag = cmplx(0.d0,0.d0,kind(0.d0))
Do no = 1,Norb
Read(10,*) Z
If ( N_Back == 1 ) bins0(nb-n_skip,no) = Z
Z_diag = Z_diag + bins0(nb-n_skip,no)*bins0(nb-n_skip,no)/Phase(nb-n_skip)
Enddo
do n = 1,Nunit
Read(10,*) Xk_p(1,n), Xk_p(2,n)
! Read
do nt = 1,Lt
do no = 1,norb
do no1 = 1,Norb
read(10,*) OneBin(nt,no,no1)
enddo
enddo
enddo
if ( sqrt(Xk_p(1,n)**2 + Xk_p(2,n)**2) < 1.D-6 ) then
Do nt = 1,Lt
do no = 1,norb
do no1 = 1,Norb
OneBin(nt,no,no1) = OneBin(nt,no,no1) - &
& bins0(nb-n_skip,no)*bins0(nb-n_skip,no1)*cmplx(dble(Nunit),0.d0,kind(0.d0))&
& /Phase(nb-n_skip)
enddo
enddo
enddo
endif
! tau = nt*dtau nt = 0,Ltrot
! beta - tau = (lt-1)*dtau - nt*dtau = (lt -1 - nt)*dtau
! --> g(nt + 1 ) = g(lt -1 -nt +1) --> g(nt + 1 ) = g(lt -nt ) --> g(nt ) = g(lt -nt +1 )
do nt = 1,Lt_eff
do no = 1,Norb
bins(n,nt,nb-n_skip) = bins(n,nt,nb-n_skip) + &
& ( OneBin(nt,no,no) + OneBin(Lt - nt + 1,no,no) ) / cmplx(2.d0,0.d0,Kind(0.d0))
enddo
enddo
Z = cmplx(0.d0,0.d0,kind(0.d0))
Do nt = 1,Lt_eff -1
do no = 1,Norb
do no1 = 1,Norb
Z = Z + cmplx(0.5d0,0.d0,Kind(0.d0)) * ( OneBin(nt,no,no1) + Onebin(nt+1,no,no1) )
enddo
enddo
enddo
Z = Z*cmplx(2.d0,0.d0,Kind(0.d0))
Bins_chi(N,Nb-n_skip) = Z
enddo
else
Read(10,*,End=10) X,no,no1,n,Y
do no = 1,Norb
Read(10,*) Z
enddo
do n = 1,Nunit
Read(10,*) X,Y
do nt = 1,LT
do no = 1,Norb
do no1 = 1,Norb
read(10,*) Z
enddo
enddo
enddo
enddo
endif
enddo
close(10)
!do n = 1,Nbins
! Write(6,*) Phase(n)
!Enddo
do n = 1,Nunit
if ( Xk_p(1,n) >= -zero .and. XK_p(2,n) >= -zero ) then
call COV(bins(n,:,:), phase, Xcov, Xmean, N_rebin )
write(File_out,'("g_",F4.2,"_",F4.2)') Xk_p(1,n), Xk_p(2,n)
Open (Unit=10,File=File_out,status="unknown")
Write(10,*) Lt_eff, nbins/N_rebin, real(lt-1,kind(0.d0))*dtau, Norb
do nt = 1, Lt_eff
Write(10,"(F14.7,2x,F16.8,2x,F16.8)") &
& dble(nt-1)*dtau, dble(Xmean(nt)), sqrt(abs(dble(Xcov(nt,nt))))
enddo
If (N_cov == 1) Then ! print covarariance