Commit 6dab1fb7 authored by Johannes Hofmann (Laptop)'s avatar Johannes Hofmann (Laptop)
Browse files

Adding Ex. 3

parent 7491980e
......@@ -421,27 +421,63 @@ with $J_{zz} = V $ and $J_{xx} = 2t$. Hence when $V/t = 2$ we reproduce the Hei
\exercise{Adding a new observable}
\red{[Stub]}
This exercise illustrates the modifications that are required to implement a new observable, the correlation function of the bond-hopping in the 1-dimensional Hubbard chain. This observable is interesting in many different setups. For example, the model studied here exhibits an emergent $SO(4)$ symmetry that relates the anti-ferromagnetic order parameter and the bond dimerization (I. Affleck, PRL 55, 1355 (1985); I. Affleck and F. D. M. Haldane, PRB 36, 5291 (1987)). Another example, where this quantity is useful to investigate, is the 1-dimensional Su-Schrieffer-Heger model describing an electron-phonon system.
Here the task if to define a new observable, the kinetic energy correlation, given by
\exerciseitem{Applying Wick's theorem}
Here the task is to define a new equal-time observable, the kinetic energy correlation, given by
\begin{align}
&\left\langle \hat{O}_{i,\delta} \hat{O}_{j,\delta'} \right\rangle - \left\langle \hat{O}_{i,\delta} \right\rangle \left\langle \hat{O}_{j,\delta'} \right\rangle = S_O\big(i-j,\delta,\delta'\big)
&\left\langle \hat{O}_{i,\delta} \hat{O}_{j,\delta'} \right\rangle - \left\langle \hat{O}_{i,\delta} \right\rangle \left\langle \hat{O}_{j,\delta'} \right\rangle = S_O\big(i-j,\delta,\delta'\big) \label{eq:new_obs}
\end{align}
where
where $i,j$ refer to the unit cells and $\delta$ encodes the bond label. Since we a working in 1D, there is only on bond per unit cell and $\delta=ax$ such that
\begin{align}
&\hat{O}_{i,x} = \sum_{\sigma}\left( \hat{c}^\dagger_{i,\sigma}\hat{c}^{\phantom\dagger}_{i+ax,\sigma} +H.c. \right).
\end{align}
[...]\\
In the 1-D Hubbard we have emergent $SO(4)$ symmetry:
Note that the first term of Eq.~\ref{eq:new_obs} is of the generic form $\sum_{\sigma,\sigma'}\left\langle \hat{c}^\dagger_{i_1,\sigma}\hat{c}^{\phantom\dagger}_{i_2,\sigma} \hat{c}^\dagger_{j_1,\sigma'}\hat{c}^{\phantom\dagger}_{j_2,\sigma'} \right\rangle$. This expectation value can be readily decomposed into single particle Greens function by using Wick's theorem. It can be applied for a fixed field configuration $\Phi$ since the Hamiltonian $\hat{H}(\Phi)$ is then bi-linear in the fermion operators.
\begin{align}
\sum_{\sigma,\sigma'}\left\langle \hat{c}^\dagger_{i_1,\sigma}\hat{c}^{\phantom\dagger}_{i_2,\sigma} \hat{c}^\dagger_{j_1,\sigma'}\hat{c}^{\phantom\dagger}_{j_2,\sigma'} \right\rangle_\Phi & =
\sum_{\sigma,\sigma'}\left(\left\langle \hat{c}^\dagger_{i_1,\sigma}\hat{c}^{\phantom\dagger}_{i_2,\sigma}\right\rangle_\Phi
\left\langle\hat{c}^\dagger_{j_1,\sigma'}\hat{c}^{\phantom\dagger}_{j_2,\sigma'} \right\rangle_\Phi
+\left\langle \hat{c}^\dagger_{i_1,\sigma}\hat{c}^{\phantom\dagger}_{j_2,\sigma'}\right\rangle_\Phi
\left\langle\hat{c}^{\phantom\dagger}_{i_2,\sigma}\hat{c}^\dagger_{j_1,\sigma'} \right\rangle_\Phi
\right)%\\
%& =
%\mathtt{N\_SUN}^2 \left\langle \hat{c}^\dagger_{i_1,\uparrow}\hat{c}^{\phantom\dagger}_{i_2,\uparrow}\right\rangle_\Phi
%\left\langle\hat{c}^\dagger_{j_1,\uparrow}\hat{c}^{\phantom\dagger}_{j_2,\uparrow} \right\rangle_\Phi
%+\mathtt{N\_SUN}\left\langle \hat{c}^\dagger_{i_1,\uparrow}\hat{c}^{\phantom\dagger}_{j_2,\uparrow}\right\rangle_\Phi
%\left\langle\hat{c}^{\phantom\dagger}_{i_2,\uparrow}\hat{c}^\dagger_{j_1,\uparrow} \right\rangle_\Phi
\end{align}
The second term of vanishes for $\sigma\neq\sigma'$ due to flavor symmetry (Mz-decoupling used here in the vanilla version) or due to the $SU(2)$ symmetry (density decoupling available in the generic implementation of Hubbard model). The single-particle Greens function are provided in the \texttt{Obser} routine, where all equal-time observables are measured, as
\begin{align}
GRC(i,j,\sigma) & = \left\langle \hat{c}^\dagger_{i,\sigma}\hat{c}^{\phantom\dagger}_{j,\sigma}\right\rangle_\Phi\\
GR(i,j,\sigma) & = \left\langle \hat{c}^{\phantom\dagger}_{i,\sigma}\hat{c}^\dagger_{j,\sigma}\right\rangle_\Phi
\end{align}
\exerciseitem{Necessary code modifications}
In the directory \texttt{\$ALF\_DIR/Solutions/Exercise\_3} we have duplicated the ALF and commented the changes that have to be carried out to the file \texttt{Hamiltonian\_Hubbard\_Plain\_Vanilla\_mod.F90} found in \texttt{\$ALF\_DIR/Prog/Hamiltonians/}, which we here shorten to ``\texttt{Vanilla}''. The following are the essential steps to be carried out:
\begin{itemize}
\item Introduce the new observable and allocate the memory required to store the measurements. This is done in the subroutine \texttt{Alloc\_obs(Ltau)} by increasing the length of the array \texttt{Obs\_eq} appropriately and adding a new case to specify the filename in which the results are stored on disc. (You might want to revisit this section later on to add the time-displaced version of the correlation function by changing \texttt{Obs\_tau} in the same fashion.)
\item The actual measurements are taken in the subroutine \texttt{Obser(GR,Phase,Ntau)}. While \texttt{GR} is passed to the subroutine, the first lines of code already implemented construct $\mathtt{GRC}=1-\mathtt{GR}^T$. (This section does not have to be modified, but it is useful to keep this in mind for future reference when you implement a new model from scratch.)
\item The measurement of a equal-time correlation function consists of two separate parts, the connected one given by the first term in Eq.~\ref{eq:new_obs}, and the background, the second term.
\item Implement the measurement of the connected part, stored as $\mathtt{Obs\_eq}(6)$ in this example, using the Wick decomposition sketched above.
\item Keep in mind that the background $\sum_i\left\langle \hat{O}_{i,\delta} \right\rangle\neq0$ is non-vanishing and has to be measured separately (you can compare with the density correlation function), and is stored in $\mathtt{Obs\_eq}(6)\%\mathtt{Obs\_Latt0}(1)$.
\item The analysis tool will then automatically combine both contributions and evaluate Eq.~\ref{eq:new_obs} using the jackknife method to estimate the mean and error or the correlation function.
\end{itemize}
The 1-D Hubbard exhibits an emergent $SO(4)$ symmetry:
\begin{align}
\left\langle \bar{S}(r)S(0) \right\rangle &\sim \frac{(-1)^r}{r}\ln^d(r)\\
\left\langle \hat{O}_{r,x} \hat{O}_{0,x} \right\rangle - \left\langle \hat{O}_{r,x} \right\rangle \left\langle \hat{O}_{0,x} \right\rangle &\sim \frac{(-1)^r}{r}\ln^\beta(r)
\end{align}
where $d=??$ and $\beta=??$ \cite{references}.
where $d=1/2$ and $\beta=-3/2$ (T. Sato, M. Hohenadler, et.al, ArXiv:2005.08996 (2020)) and this exercise provides all the tool required to study it. Be aware the large system sizes, and therefore long run times, are required to extract the logarithmic scaling corrections. More details are discussed in the appendix of above reference (ArXiv:2005.08996).
[It should be added to the Predefined Structures.]\\
Finally, it is straight forward to implement the time-displaced version of this correlation function, following essentially the same steps as described above. The observable is now stored in \texttt{Obs\_tau}, the measurements are taken in \texttt{ObserT}, and you can find the implementation in the solutions to this exercise as well.
......
.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