Commit 5d0e3d2e authored by Jefferson Stafusa E. Portela's avatar Jefferson Stafusa E. Portela
Browse files

Exercise 2 of part 2 (tv model) and minor changes elsewhere.

parent 426c7c3a
......@@ -296,7 +296,7 @@ git clone git@git.physik.uni-wuerzburg.de:ALF/ALF_Tutorial.git
Here we will modify the code so as to allow for different hopping matrix elements along the $x$ and $y$ directions of a square lattice.
\exerciseitem{Modifying the hopping}
To do so we start from the module \texttt{Hamiltonian\_Hubbard\_Plain\_Vanilla\_mod.F90} found in \texttt{<ALF\_DIR>/Prog/Hamiltonians/}, which we will simply call ``\texttt{Vanilla}'', and then:
To do so we start from the module \texttt{Hamiltonian\_Hubbard\_Plain\_Vanilla\_mod.F90} found in \texttt{\$ALF\_DIR/Prog/Hamiltonians/}, which we will simply call ``\texttt{Vanilla}'', and then:
\begin{itemize}
\item Add \texttt{Ham\_Ty} to the \texttt{VAR\_Hubbard\_Plain\_Vanilla} name space in the parameter file \texttt{parameters}.
\item Declare a new variable, \texttt{Ham\_Ty}, in the module's specification (just search for the declaration of \texttt{Ham\_T} in \texttt{Vanilla}).
......@@ -356,6 +356,8 @@ The above form is readily included in the ALF since the interaction is written
\hat{H} = -t \sum_{i} \left( \hat{c}^{\dagger}_{i} \hat{c}^{\phantom\dagger}_{i+a} + \hat{c}^{\dagger}_{i+a} \hat{c}^{\phantom\dagger}_{i} \right) + V \sum_{i} \left( \hat{n}_{i} - 1/2 \right) \left( \hat{n}_{i+a} - 1/2 \right).
\end{equation}
Note that the t-V model is already implemented in ALF in the module \texttt{Hamiltonian\_tV\_mod.F90} found in \texttt{Prog/Hamiltonians/}. While it can be used for checking your own results\footnote{A short simulation of the t-V model can be conveniently run using its Jupyter notebook available in pyALF.}, you are not supposed to reproduce that implementation -- which is more general and makes use of predefined structures -- but instead to write a simpler one, based on the module \texttt{Hamiltonian\_Hubbard\_Plain\_Vanilla\_mod.F90} as detailed below.
\red{[Mention more general tV implementation in pred. strutct. (for checking correctness), as well as the tV Jupyter notebook.]}
......
......@@ -151,7 +151,10 @@
Type (Lattice), private :: Latt
Type (Unit_cell), private :: Latt_unit
Integer, private :: L1, L2
real (Kind=Kind(0.d0)), private :: Ham_T , Ham_Ty , Ham_U, Ham_chem
real (Kind=Kind(0.d0)), private :: Ham_T , Ham_U, Ham_chem
!!!!!!! Modifications for Exercise 1a
real (Kind=Kind(0.d0)), private :: Ham_Ty
!!!!!!!
real (Kind=Kind(0.d0)), private :: Dtau, Beta, Theta
Integer , private :: N_part
Character (len=64), private :: Model, Lattice_type
......@@ -187,7 +190,10 @@
NAMELIST /VAR_Lattice/ L1, L2, Lattice_type, Model
!!!!!!! Modifications for Exercise 1a
!NAMELIST /VAR_Hubbard_Plain_Vanilla/ Ham_T, ham_chem, ham_U, Dtau, Beta, Projector, Theta, Symm, N_part
NAMELIST /VAR_Hubbard_Plain_Vanilla/ Ham_T, Ham_Ty, ham_chem, ham_U, Dtau, Beta, Projector, Theta, Symm, N_part
!!!!!!!
......@@ -250,7 +256,10 @@
CALL MPI_BCAST(Projector ,1, MPI_LOGICAL , 0,Group_Comm,ierr)
CALL MPI_BCAST(Dtau ,1, MPI_REAL8 , 0,Group_Comm,ierr)
CALL MPI_BCAST(Beta ,1, MPI_REAL8 , 0,Group_Comm,ierr)
CALL MPI_BCAST(Ham_T ,1, MPI_REAL8 , 0,Group_Comm,ierr)
CALL MPI_BCAST(ham_T ,1, MPI_REAL8 , 0,Group_Comm,ierr)
!!!!!!! Modifications for Exercise 1a
CALL MPI_BCAST(Ham_Ty ,1, MPI_REAL8 , 0,Group_Comm,ierr)
!!!!!!!
CALL MPI_BCAST(ham_chem ,1, MPI_REAL8 , 0,Group_Comm,ierr)
CALL MPI_BCAST(ham_U ,1, MPI_REAL8 , 0,Group_Comm,ierr)
#endif
......@@ -350,10 +359,10 @@
If ( L2 > 1 ) then
Iy = Latt%nnlist(I,0,1)
!!!!!!! Modifications for Exercise 1a
Op_T(1,nf)%O(I, Iy) = cmplx(-Ham_T, 0.d0, kind(0.D0))
Op_T(1,nf)%O(Iy, I ) = cmplx(-Ham_T, 0.d0, kind(0.D0))
! Op_T(1,nf)%O(I, Iy) = cmplx(-Ham_Ty, 0.d0, kind(0.D0))
! Op_T(1,nf)%O(Iy, I ) = cmplx(-Ham_Ty, 0.d0, kind(0.D0))
!Op_T(1,nf)%O(I, Iy) = cmplx(-Ham_T, 0.d0, kind(0.D0))
!Op_T(1,nf)%O(Iy, I ) = cmplx(-Ham_T, 0.d0, kind(0.D0))
Op_T(1,nf)%O(I, Iy) = cmplx(-Ham_Ty, 0.d0, kind(0.D0))
Op_T(1,nf)%O(Iy, I ) = cmplx(-Ham_Ty, 0.d0, kind(0.D0))
!!!!!!!
endif
Op_T(1,nf)%O(I, I ) = cmplx(-Ham_chem, 0.d0, kind(0.D0))
......@@ -414,10 +423,10 @@
If (L2 > 1 ) Then
Iy = Latt%nnlist(I,0,1)
!!!!!!! Modifications for Exercise 1a
H0(I, Iy) = -Ham_T *(1.d0 - Delta)
H0(Iy, I ) = -Ham_T *(1.d0 - Delta)
!H0(I, Iy) = -Ham_Ty *(1.d0 - Delta)
!H0(Iy, I ) = -Ham_Ty *(1.d0 - Delta)
!H0(I, Iy) = -Ham_T *(1.d0 - Delta)
!H0(Iy, I ) = -Ham_T *(1.d0 - Delta)
H0(I, Iy) = -Ham_Ty *(1.d0 - Delta)
H0(Iy, I ) = -Ham_Ty *(1.d0 - Delta)
!!!!!!!
Endif
Enddo
......
.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
all: $(OBJS1) $(BINS)
cov_tau_ph.F90 : cov_tau.F90
cpp -w -P -DPartHole cov_tau.F90 -o cov_tau_ph.F90
%.o: %.F90
$(ALF_FC) -c $*.F90 $(ALF_FLAGS_ANA)
%.out: %.F90
$(ALF_FC) -c $(ALF_FLAGS_ANA) $<
$(ALF_FC) -o $@ $*.o $(OBJS1) $(ALF_LIB)
tidy:
rm -f $(OBJS) $(OBJS1) $(MODS) cov_tau_ph.F90
clean: tidy
rm -f $(BINS)
! Copyright (C) 2016-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 MaxEnt_Wrapper
!--------------------------------------------------------------------
!> @author
!> ALF-project
!
!> @brief
!> General wrapper for maxent. Handles particle, partile-hole, particle-particle,
!> channels as well as zero temperature. See documentation for details.
!>
!
!--------------------------------------------------------------------
Use MaxEnt_stoch_mod
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)
Implicit none
Real (Kind=Kind(0.d0)), INTENT(INOUT),allocatable :: XCOV(:,:), XQMC(:), XTAU(:)
Real (Kind=Kind(0.d0)), INTENT (IN) :: Tolerance
Integer, INTENT(IN) :: Ntau_st, Ntau_en
Integer, INTENT(INOUT) :: Ntau
end Subroutine Rescale
end Interface
Real (Kind=Kind(0.d0)), Dimension(:) , allocatable :: XQMC, XQMC_st, XTAU, Xtau_st, &
& Alpha_tot, om_bf, alp_bf, xom, A
Real (Kind=Kind(0.d0)), Dimension(:,:), allocatable :: XCOV, XCOV_st
Real (Kind=Kind(0.d0)) :: X_moments(2), Xerr_moments(2)
Real (Kind=Kind(0.d0)), External :: XKER_ph, Back_trans_ph, XKER_pp, Back_trans_pp, &
& XKER_p, Back_trans_p, XKER_T0, Back_trans_T0
Character (Len=64) :: command, File1, File2
Complex (Kind=Kind(0.d0)) :: Z
Integer :: Ngamma, Ndis, NBins, NSweeps, Nwarm, N_alpha, N_cov
Integer :: N_skip, N_rebin, Norb
Real (Kind=Kind(0.d0)) :: OM_st, OM_en, alpha_st, R, Tolerance
Logical :: Checkpoint
Character (Len=2) :: Channel
Integer :: nt, nt1, io_error, n,nw, nwp, ntau, N_alpha_1, i, nbin_qmc
Integer :: ntau_st, ntau_en, ntau_new, Ntau_old
Real (Kind=Kind(0.d0)) :: dtau, pi, xmom1, x,x1,x2, tau, omp, om, Beta,err, delta, Dom
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
NAMELIST /VAR_errors/ N_skip, N_rebin, N_cov
open(unit=30,file='parameters',status='old',action='read', iostat=io_error)
if (io_error.eq.0) then
READ(30,NML=VAR_errors)
READ(30,NML=VAR_Max_Stoch)
else
write(error_unit,*) 'No file parameters '
error stop 1
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
Allocate ( XCOV(NTAU,NTAU), XQMC(NTAU),XTAU(NTAU) )
XCOV = 0.d0
Do nt = 1,NTAU
read(10,*) xtau(nt), xqmc(nt), err
xcov(nt,nt) = err*err
Enddo
if (N_cov.eq.1) then
do nt = 1,ntau
do nt1 = 1,ntau
read(10,*) xcov(nt,nt1)
enddo
enddo
endif
close(10)
dtau = Xtau(2) - Xtau(1)
Zero = 1.D-10
pi = acos(-1.d0)
Ntau_st = 1
Ntau_en = Ntau
Select Case (Channel)
Case ("PH")
xmom1 = pi * xqmc(1)
Case ("PP")
xmom1 = 2.d0* pi * xqmc(1)
Case ("P")
xmom1 = pi * real(norb,Kind(0.d0))
! 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)
! it will be removed.
Ntau_en = Ntau - 1
Ntau_st = 1
if ( xcov(1,1) < zero ) ntau_st = 2
Case ("T0")
xmom1 = pi*xqmc(1)
Ntau_st = 1
if ( xcov(1,1) < zero ) ntau_st = 2
Case default
Write(error_unit,*) "Channel not yet implemented"
error stop 1
end Select
Ntau_old = Ntau
Call Rescale ( XCOV, XQMC,XTAU, Ntau_st, Ntau_en, Tolerance, NTAU)
Write(50,"('Data has been rescaled from Ntau ', I4,' to ', I4)") NTAU_old, Ntau
If ( Ntau <= 4 ) then
write(error_unit,*) 'Not enough data!'
error stop 1
Endif
If ( nbin_qmc > 2*Ntau .and. N_cov == 0 ) Write(50,*) 'Consider using the covariance. You seem to have enough bins'
If ( nbin_qmc < 2*Ntau .and. N_cov == 1 ) Write(50,*) 'You do not seem to have enough bins for a reliable estimate of the covariance '
! Store
Allocate ( XCOV_st(NTAU,NTAU), XQMC_st(NTAU),XTAU_st(NTAU) )
XCOV_st = XCOV
XQMC_st = XQMC
XTAU_st = XTAU
Allocate (Alpha_tot(N_alpha) )
do nt = 1,N_alpha
alpha_tot(nt) = alpha_st*(R**(nt-1))
enddo
write(50,"('First Moment, Beta ',2x,F12.6,2x,F12.6)") Xmom1, Beta
Select Case (Channel)
Case ("PH")
If (N_cov == 1 ) then
Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_ph, Back_Trans_ph, Beta, &
& Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm, N_Cov)
else
Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_ph, Back_Trans_ph, Beta, &
& Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm)
endif
! Beware: Xqmc and cov are modified in the MaxEnt_stoch call.
Case ("PP")
If (N_cov == 1 ) then
Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_pp, Back_Trans_pp, Beta, &
& Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm, N_Cov)
else
Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_pp, Back_Trans_pp, Beta, &
& Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm)
endif
! Beware: Xqmc and cov are modified in the MaxEnt_stoch call.
Case ("P")
If (N_cov == 1 ) then
Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_p, Back_Trans_p, Beta, &
& Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm, N_Cov)
else
Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_p, Back_Trans_p, Beta, &
& Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm)
endif
! Beware: Xqmc and cov are modified in the MaxEnt_stoch call.
Case ("T0")
If (N_cov == 1 ) then
Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_T0, Back_Trans_T0, Beta, &
& Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm, N_Cov)
else
Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_T0, Back_Trans_T0, Beta, &
& Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm)
endif
! Beware: Xqmc and cov are modified in the MaxEnt_stoch call.
Case default
Write(error_unit,*) "Channel not yet implemented"
error stop 1
end Select
If ( .not. Checkpoint ) then
Command = "rm dump*"
Call System (Command)
Command = "ls"
Call System (Command)
endif
Open (Unit=10,File="energies",status="unknown")
Do n = 1,N_alpha
Read(10,*) X,X1,X2
enddo
Write(50,"('Chisq at lowest temperature: 1/T, Chi^2, Error',2x,F12.6,2x,F12.6,2x,F12.6)") X,X1,X2
close(50)
Open (Unit = 10,File="Best_fit", Status ="unknown")
Allocate (om_bf(Ngamma), alp_bf(Ngamma) )
DO i = 1, Ngamma
read(10,*) om_bf(i), alp_bf(i)
Enddo
close(10)
Open (Unit = 11,File="Data_out", Status ="unknown")
do nt = 1,Ntau
X = 0.d0
tau = xtau_st(nt)
Select Case (Channel)
Case ("PH")
do i = 1,Ngamma
X = X + alp_bf(i)*Xker_ph(tau,om_bf(i), beta)
enddo
Case ("PP")
do i = 1,Ngamma
X = X + alp_bf(i)*Xker_pp(tau,om_bf(i), beta)
enddo
Case ("P")
do i = 1,Ngamma
X = X + alp_bf(i)*Xker_p(tau,om_bf(i), beta)
enddo
Case ("T0")
do i = 1,Ngamma
X = X + alp_bf(i)*Xker_T0(tau,om_bf(i), beta)
enddo
Case default
Write(error_unit,*) "Channel not yet implemented"
error stop 1
end Select
Write(11,"(F14.7,2x,F14.7,2x,F14.7,2x,F14.7)") xtau_st(nt), xqmc_st(nt), sqrt(xcov_st(nt,nt)), xmom1*X
enddo
close(11)
N_alpha_1 = N_alpha - 10
File1 ="Aom_ps"
file2 =File_i(File1,N_alpha_1)
Open(Unit=66,File=file2,status="unknown")
Allocate (xom(Ndis), A(Ndis))
do nw = 1,Ndis
read(66,*) xom(nw), A(nw), x, x1, x2
enddo
close(66)
Dom = xom(2) - xom(1)
! Compute the real frequency Green function.
delta = Dom
Open (Unit=43,File="Green", Status="unknown", action="write")
pi = acos(-1.d0)
do nw = 1,Ndis
Z = cmplx(0.d0,0.d0,Kind(0.d0))
om = xom(nw)
do nwp = 1,Ndis
omp = xom(nwp)
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
enddo
close(43)
end Program MaxEnt_Wrapper
Real (Kind=Kind(0.d0)) function XKER_ph(tau,om, beta)
Implicit None
real (Kind=Kind(0.d0)) :: tau, om, pi, beta
pi = 3.1415927
XKER_ph = (exp(-tau*om) + exp(-( beta - tau )*om ) )/( pi*(1.d0 + exp( - beta * om ) ) )
end function XKER_ph
Real (Kind=Kind(0.d0)) function XKER_pp(tau,om, beta)
Implicit None
real (Kind=Kind(0.d0)) :: tau, om, pi, beta
pi = 3.1415927
XKER_pp = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
end function XKER_pp
Real (Kind=Kind(0.d0)) function XKER_p(tau,om, beta)
Implicit None
real (Kind=Kind(0.d0)) :: tau, om, pi, beta
pi = 3.1415927
XKER_p = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
end function XKER_p
Real (Kind=Kind(0.d0)) function XKER_T0(tau,om, beta)
Implicit None
real (Kind=Kind(0.d0)) :: tau, om, pi, beta
pi = 3.1415927
XKER_T0 = exp(-tau*om) / pi
end function XKER_T0
Real (Kind=Kind(0.d0)) function Back_trans_ph(Aom, om, beta)
Implicit None
real (Kind=Kind(0.d0)) :: Aom, om, beta
Back_trans_ph = Aom/(1.d0 + exp(-beta*om) )
! This gives S(q,om) = chi(q,om)/(1 - e^(-beta om))
end function BACK_TRANS_PH
Real (Kind=Kind(0.d0)) function Back_trans_pp(Aom, om, beta)
Implicit None
real (Kind=Kind(0.d0)) :: Aom, om, beta
real (Kind=Kind(0.d0)) :: Zero
Zero = 1.D-8
if ( abs(om) < zero ) then
Back_trans_pp = beta * Aom/2.d0
else
Back_trans_pp = Aom * (1.d0 - exp(-beta*om) ) / (om *( 1.d0 + exp(-beta*om) ) )
endif
! This gives = chi(q,om)/omega
end function BACK_TRANS_PP
Real (Kind=Kind(0.d0)) function Back_trans_p(Aom, om, beta)
Implicit None
real (Kind=Kind(0.d0)) :: Aom, om, beta
Back_trans_p = Aom
end function BACK_TRANS_P
Real (Kind=Kind(0.d0)) function Back_trans_T0(Aom, om, beta)
Implicit None
real (Kind=Kind(0.d0)) :: Aom, om, beta
Back_trans_T0 = Aom
end function BACK_TRANS_T0
Subroutine Rescale ( XCOV, XQMC,XTAU, Ntau_st, Ntau_en, Tolerance, Ntau)
Implicit none
Real (Kind=Kind(0.d0)), INTENT(INOUT), allocatable :: XCOV(:,:), XQMC(:), XTAU(:)
Real (Kind=Kind(0.d0)), INTENT (IN) :: Tolerance
Integer, INTENT(IN) :: Ntau_st, Ntau_en
Integer, INTENT(INOUT) :: Ntau
!Local
Integer :: Ntau_new, nt, nt1
Integer, allocatable :: List(:)
Real (Kind=Kind(0.d0)), dimension(:,:), allocatable :: XCOV_st
Real (Kind=Kind(0.d0)), dimension(:) , allocatable :: XQMC_st, XTAU_st
! Count the number of elements
ntau_new = 0
Do nt = ntau_st,ntau_en
if ( sqrt(xcov(nt,nt))/xqmc(nt) < Tolerance .and. xqmc(nt) > 0.d0 ) then
ntau_new = ntau_new + 1
endif
Enddo
Allocate ( XCOV_st(NTAU_new,NTAU_new), XQMC_st(NTAU_new),XTAU_st(NTAU_new), List(NTAU_new) )
ntau_new = 0
Do nt = ntau_st,ntau_en
if ( sqrt(xcov(nt,nt))/xqmc(nt) < Tolerance .and. xqmc(nt) > 0.d0 ) then
ntau_new = ntau_new + 1
List(ntau_new) = nt
endif
Enddo
do nt = 1,ntau_new
XQMC_st(nt) = XQMC( List(nt) )