Commit 7491980e authored by Jefferson Stafusa E. Portela's avatar Jefferson Stafusa E. Portela
Browse files

Solution of exercise 2 of part 2 (tV model).

parent a4b80ffb

Too many changes to show.

To preserve performance only 248 of 248+ files are displayed.
! 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
Do nt = 1,LT_eff
Do nt1 = 1,LT_eff
Write(10,*) dble(Xcov(nt,nt1))
Enddo
Enddo
Endif
close(10)
endif
enddo
V_help = cmplx(0.d0,0.d0,kind(0.d0))
do n = 1,Nunit
do nb = 1,nbins
do nt = 1,LT_eff
V_help(nt,nb) = V_help(nt,nb) + bins(n,nt,nb)
enddo
enddo
enddo
V_help = V_help/dble(Nunit)
call COV(V_help, phase, Xcov, Xmean, N_Rebin )
write(File_out,'("g_R0")')
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 covariance
Do nt = 1,LT_eff
Do nt1 = 1,LT_eff
Write(10,*) dble(Xcov(nt,nt1))
Enddo
Enddo
Endif
close(10)
! Print susceptibilities
Open (Unit=33,File="SuscepJ" ,status="unknown")
Do n = 1,Nunit
call ERRCALCJ(Bins_chi(n,:), PhaseI, ZMean, ZERR, N_rebin )
Zmean = Zmean*dtau
Zerr = Zerr*dtau
Write(33,"(F12.6,2x,F12.6,2x,F16.8,2x,F16.8,2x,F16.8,2x,F16.8)") &
& Xk_p(1,n), Xk_p(2,n), dble(ZMean), dble(ZERR), aimag(ZMean), aimag(ZERR)
enddo
Close(33)
! Deallocate space
Deallocate ( bins, Phase,PhaseI, Xk_p, V_help, bins0)
Deallocate ( Bins_chi )
end Program Cov_tau
!!$ Interface
!!$ Integer function Rot90(n, Xk_p, Nunit)
!!$ Implicit none
!!$ Integer, INTENT(IN) :: Nunit,n
!!$ Real (Kind=Kind(0.d0)), INTENT(IN) :: Xk_p(2,Nunit)
!!$ end function Rot90
!!$ end Interface
!!$
!!$ Integer function Rot90(n, Xk_p, Nunit)
!!$
!!$ Implicit none
!!$ Integer, INTENT(IN) :: Nunit,n
!!$ Real (Kind=Kind(0.d0)), INTENT(IN) :: Xk_p(2,Nunit)
!!$
!!$ !Local
!!$ real (Kind=Kind(0.d0)) :: X1_p(2), Zero, pi, X
!!$ Integer :: m
!!$
!!$ Zero = 1.D-4
!!$ pi = acos(-1.d0)
!!$ X1_p(1) = Xk_p(2,n)
!!$ X1_p(2) = -Xk_p(1,n)
!!$ if (X1_p(1) < -pi + Zero ) X1_p(1) = X1_p(1) + 2.0*pi
!!$ if (X1_p(2) < -pi + Zero ) X1_p(2) = X1_p(2) + 2.0*pi
!!$
!!$ Rot90 = 0
!!$ Do m = 1,Nunit
!!$ X = sqrt( (X1_p(1) -Xk_p(1,m))**2 + (X1_p(2) -Xk_p(2,m))**2 )
!!$ If ( X < Zero) then
!!$ Rot90 = m
!!$ exit
!!$ endif
!!$ Enddo
!!$
!!$ end function Rot90
101 4 20.000000000000000 1
0.0000000 0.00564888 0.00037249
0.1000000 0.00564888 0.00037249
0.2000000 0.00564888 0.00037249
0.3000000 0.00564888 0.00037249
0.4000000 0.00564888 0.00037249
0.5000000 0.00564888 0.00037249
0.6000000 0.00564888 0.00037249
0.7000000 0.00564888 0.00037249
0.8000000 0.00564888 0.00037249
0.9000000 0.00564888 0.00037249
1.0000000 0.00564888 0.00037249
1.1000000 0.00564888 0.00037249
1.2000000 0.00564888 0.00037249
1.3000000 0.00564888 0.00037249
1.4000000 0.00564888 0.00037249
1.5000000 0.00564888 0.00037249
1.6000000 0.00564888 0.00037249
1.7000000 0.00564888 0.00037249
1.8000000 0.00564888 0.00037249
1.9000000 0.00564888 0.00037249
2.0000000 0.00564888 0.00037249
2.1000000 0.00564888 0.00037249
2.2000000 0.00564888 0.00037249
2.3000000 0.00564888 0.00037249
2.4000000 0.00564888 0.00037249
2.5000000 0.00564888 0.00037249
2.6000000 0.00564888 0.00037249
2.7000000 0.00564888 0.00037249
2.8000000 0.00564888 0.00037249
2.9000000 0.00564888 0.00037249
3.0000000 0.00564888 0.00037249
3.1000000 0.00564888 0.00037249
3.2000000 0.00564888 0.00037249
3.3000000 0.00564888 0.00037249
3.4000000 0.00564888 0.00037249
3.5000000 0.00564888 0.00037249
3.6000000 0.00564888 0.00037249
3.7000000 0.00564888 0.00037249
3.8000000 0.00564888 0.00037249
3.9000000 0.00564888 0.00037249
4.0000000 0.00564888 0.00037249
4.1000000 0.00564888 0.00037249
4.2000000 0.00564888 0.00037249
4.3000000 0.00564888 0.00037249
4.4000000 0.00564888 0.00037249
4.5000000 0.00564888 0.00037249
4.6000000 0.00564888 0.00037249
4.7000000 0.00564888 0.00037249
4.8000000 0.00564888 0.00037249
4.9000000 0.00564888 0.00037249
5.0000000 0.00564888 0.00037249
5.1000000 0.00564888 0.00037249
5.2000000 0.00564888 0.00037249
5.3000000 0.00564888 0.00037249
5.4000000 0.00564888 0.00037249
5.5000000 0.00564888 0.00037249
5.6000000 0.00564888 0.00037249
5.7000000 0.00564888 0.00037249
5.8000000 0.00564888 0.00037249
5.9000000 0.00564888 0.00037249
6.0000000 0.00564888 0.00037249
6.1000000 0.00564888 0.00037249
6.2000000 0.00564888 0.00037249
6.3000000 0.00564888 0.00037249
6.4000000 0.00564888 0.00037249
6.5000000 0.00564888 0.00037249
6.6000000 0.00564888 0.00037249
6.7000000 0.00564888 0.00037249
6.8000000 0.00564888 0.00037249
6.9000000 0.00564888 0.00037249
7.0000000 0.00564888 0.00037249
7.1000000 0.00564888 0.00037249
7.2000000 0.00564888 0.00037249
7.3000000 0.00564888 0.00037249
7.4000000 0.00564888 0.00037249
7.5000000 0.00564888 0.00037249
7.6000000 0.00564888 0.00037249
7.7000000 0.00564888 0.00037249
7.8000000 0.00564888 0.00037249
7.9000000 0.00564888 0.00037249
8.0000000 0.00564888 0.00037249
8.1000000 0.00564888 0.00037249
8.2000000 0.00564888 0.00037249
8.3000000 0.00564888 0.00037249
8.4000000 0.00564888 0.00037249
8.5000000 0.00564888 0.00037249
8.6000000 0.00564888 0.00037249
8.7000000 0.00564888 0.00037249
8.8000000 0.00564888 0.00037249
8.9000000 0.00564888 0.00037249
9.0000000 0.00564888 0.00037249
9.1000000 0.00564888 0.00037249
9.2000000 0.00564888 0.00037249
9.3000000 0.00564888 0.00037249
9.4000000 0.00564888 0.00037249
9.5000000 0.00564888 0.00037249
9.6000000 0.00564888 0.00037249
9.7000000 0.00564888 0.00037249
9.8000000 0.00564888 0.00037249
9.9000000 0.00564888 0.00037249
10.0000000 0.00564888 0.00037249
101 4 20.000000000000000 1
0.0000000 0.02703604 0.00056381
0.1000000 0.02558689 0.00048865
0.2000000 0.02404612 0.00046384
0.3000000 0.02266691 0.00041637
0.4000000 0.02138146 0.00038988
0.5000000 0.02011707 0.00043076
0.6000000 0.01893531 0.00043215
0.7000000 0.01789660 0.00038272
0.8000000 0.01694242 0.00030865
0.9000000 0.01596182 0.00033644
1.0000000 0.01507406 0.00036108
1.1000000 0.01427029 0.00040721
1.2000000 0.01347918 0.00039757
1.3000000 0.01277488 0.00032906
1.4000000 0.01199225 0.00031102
1.5000000 0.01136075 0.00034612
1.6000000 0.01074231 0.00037026
1.7000000 0.01005816 0.00037387
1.8000000 0.00952787 0.00036177
1.9000000 0.00899533 0.00034521
2.0000000 0.00840044 0.00037483
2.1000000 0.00787391 0.00033869
2.2000000 0.00735675 0.00027804
2.3000000 0.00688524 0.00027004
2.4000000 0.00646888 0.00024324
2.5000000 0.00610692 0.00020416
2.6000000 0.00565014 0.00026696
2.7000000 0.00529328 0.00030920
2.8000000 0.00501752 0.00041528
2.9000000 0.00474214 0.00038853
3.0000000 0.00442086 0.00045808
3.1000000 0.00415081 0.00050088
3.2000000 0.00388823 0.00053389
3.3000000 0.00372125 0.00048297
3.4000000 0.00346456 0.00047323
3.5000000 0.00331320 0.00054223
3.6000000 0.00306657 0.00059436
3.7000000 0.00282132 0.00060828
3.8000000 0.00273549 0.00056047
3.9000000 0.00260571 0.00048001
4.0000000 0.00245197 0.00042751
4.1000000 0.00229411 0.00044915
4.2000000 0.00209023 0.00046182
4.3000000 0.00198964 0.00040339
4.4000000 0.00198155 0.00041674
4.5000000 0.00188164 0.00040700
4.6000000 0.00166495 0.00041639
4.7000000 0.00164041 0.00037374
4.8000000 0.00152718 0.00038228
4.9000000 0.00151772 0.00028047
5.0000000 0.00166466 0.00020204
5.1000000 0.00188017 0.00040855
5.2000000 0.00163528 0.00031620
5.3000000 0.00120602 0.00030597
5.4000000 0.00156728 0.00036376
5.5000000 0.00158899 0.00042868
5.6000000 0.00153613 0.00054528
5.7000000 0.00145282 0.00045633
5.8000000 0.00146736 0.00052759
5.9000000 0.00139166 0.00027936
6.0000000 0.00123800 0.00024694
6.1000000 0.00116421 0.00025759
6.2000000 0.00111425 0.00025903
6.3000000 0.00070928 0.00054038
6.4000000 0.00066809 0.00043607
6.5000000 0.00049681 0.00033066
6.6000000 0.00026013 0.00047692
6.7000000 0.00005827 0.00056712
6.8000000 -0.00001694 0.00063427
6.9000000 -0.00023819 0.00079893
7.0000000 -0.00011764 0.00063917
7.1000000 -0.00042749 0.00080545
7.2000000 -0.00056104 0.00090456
7.3000000 -0.00074845 0.00100928
7.4000000 -0.00044243 0.00081248
7.5000000 -0.00040364 0.00073919
7.6000000 -0.00017409 0.00071399
7.7000000 -0.00006792 0.00071895
7.8000000 -0.00027232 0.00072780
7.9000000 -0.00023713 0.00074685
8.0000000 -0.00036099 0.00081767
8.1000000 -0.00040878 0.00079412
8.2000000 -0.00017882 0.00071644
8.3000000 -0.00014003 0.00056451
8.4000000 0.00009972 0.00020255
8.5000000 0.00011418 0.00023028
8.6000000 0.00028821 0.00007510
8.7000000 0.00023065 0.00008809
8.8000000 0.00017003 0.00012897
8.9000000 0.00006688 0.00019364
9.0000000 0.00000196 0.00020400
9.1000000 -0.00000418 0.00019048
9.2000000 -0.00007579 0.00024405
9.3000000 0.00000478 0.00020206
9.4000000 -0.00003650 0.00026808
9.5000000 0.00001283 0.00025609
9.6000000 0.00006198 0.00028659
9.7000000 0.00008960 0.00032355
9.8000000 -0.00015315 0.00043868
9.9000000 -0.00002345 0.00039626
10.0000000 0.00022370 0.00043255
101 4 20.000000000000000 1
0.0000000 0.05447215 0.00132959
0.1000000 0.04888530 0.00128299
0.2000000 0.04350947 0.00128131
0.3000000 0.03893262 0.00102061
0.4000000 0.03472979 0.00089917
0.5000000 0.03082995 0.00081492
0.6000000 0.02727093 0.00082102
0.7000000 0.02439256 0.00068971
0.8000000 0.02186980 0.00060177
0.9000000 0.01956035 0.00065817
1.0000000 0.01727521 0.00068892
1.1000000 0.01543804 0.00076754
1.2000000 0.01396098 0.00074610
1.3000000 0.01250752 0.00054616
1.4000000 0.01112890 0.00050937
1.5000000 0.00988225 0.00048212
1.6000000 0.00882137 0.00045986
1.7000000 0.00777634 0.00035345
1.8000000 0.00689779 0.00029103
1.9000000 0.00600310 0.00018762
2.0000000 0.00550746 0.00017066
2.1000000 0.00499114 0.00017591
2.2000000 0.00451968 0.00022541
2.3000000 0.00396667 0.00026069
2.4000000 0.00330035 0.00016056
2.5000000 0.00280734 0.00020435
2.6000000 0.00248310 0.00024793
2.7000000 0.00222426 0.00030891
2.8000000 0.00190277 0.00027394
2.9000000 0.00168671 0.00013605
3.0000000 0.00140177 0.00009234
3.1000000 0.00114938 0.00015680
3.2000000 0.00096309 0.00024535
3.3000000 0.00083311 0.00023592
3.4000000 0.00062827 0.00033789
3.5000000 0.00056287 0.00033206
3.6000000 0.00057207 0.00031899
3.7000000 0.00052205 0.00049389
3.8000000 0.00056458 0.00051588
3.9000000 0.00056636 0.00047025
4.0000000 0.00047716 0.00042213
4.1000000 0.00052020 0.00052663
4.2000000 0.00020283 0.00064283
4.3000000 0.00009329 0.00067813
4.4000000 0.00024959 0.00060553
4.5000000 0.00042192 0.00050570
4.6000000 0.00043264 0.00052185
4.7000000 0.00026822 0.00049353
4.8000000 0.00043508 0.00048254
4.9000000 0.00053334 0.00052446
5.0000000 0.00026729 0.00054522
5.1000000 0.00015594 0.00052036
5.2000000 0.00034095 0.00073925
5.3000000 0.00018385 0.00075678
5.4000000 0.00038371 0.00078869