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

Exercise 2 of part 2 (tv model) done, but untested.

parent 5d0e3d2e
...@@ -296,7 +296,7 @@ git clone git@git.physik.uni-wuerzburg.de:ALF/ALF_Tutorial.git ...@@ -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. 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} \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 here shorten to ``\texttt{Vanilla}'', and then:
\begin{itemize} \begin{itemize}
\item Add \texttt{Ham\_Ty} to the \texttt{VAR\_Hubbard\_Plain\_Vanilla} name space in the parameter file \texttt{parameters}. \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}). \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}).
...@@ -358,29 +358,24 @@ The above form is readily included in the ALF since the interaction is written ...@@ -358,29 +358,24 @@ The above form is readily included in the ALF since the interaction is written
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. 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.]}
\exerciseitem{Define new model} \exerciseitem{Define new model}
\red{[To be updated:]}\\ In the directory \texttt{\$ALF\_DIR/Solutions/Exercise\_2} 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:
In the directory \texttt{Solutions/Exercise\_3} we have duplicated the ALF and commented the changes that have to be carried out to the file
\texttt{Hamiltonian\_Examples.f90 } in the \texttt{Prog} directory so as to include the \texttt{t\_V} model.
Here are the steps to be carried out.
\begin{itemize} \begin{itemize}
\item Add the t-V name space in the parameter file so as to read in the appropriate variables. \item Add the \texttt{VAR\_t\_V} name space in the file \texttt{parameters} and set the necessary variables -- or simply rename the \texttt{VAR\_Hubbard\_Plain\_Vanilla} name space to \texttt{VAR\_t\_V} and, within it, \texttt{Ham\_U} to \texttt{Ham\_Vint}. (Ignore the name space \texttt{VAR\_tV}, which is used by the general implementation mentioned above.)
\item Declare new variables in the \texttt{Hamiltonian\_Examples.f90 } file. \item Declare a new variable, \texttt{Ham\_Vint}, in \texttt{Vanilla}'s specification.
\item In the \texttt{Ham\_set} subroutine of the file \texttt{Hamiltonian\_Examples.f90} set and read in the parameters for the new model, \item Add the \texttt{VAR\_t\_V} name space (\texttt{namelist}) declaration at the \texttt{Ham\_Set} subroutine of \texttt{Vanilla}, containing the same variables the name space contains in \texttt{parameters}, and read it in.
\texttt{t\_V}. For this model \texttt{NF=1} and \texttt{N\_SUN=1} since we are working with spinless fermions \item Still in the \texttt{Ham\_set} subroutine of \texttt{Vanilla}: set \texttt{NF=1}, since we are working with spinless fermions; change the MPI broadcast call for \texttt{Ham\_U} to broadcast \texttt{Ham\_Vint} instead; and change similarly the output to the info file.
\item In the \texttt{Ham\_V} subroutine you will have to add the new interaction. \item In the \texttt{Ham\_V} subroutine you have to add the new interaction. For a given bond at a given time slice, we need to decouple the interaction:
For a given bond at a given time-slice, we need to decouple the interaction:
\begin{equation} \begin{equation}
e^{\Delta \tau \frac{V}{2} \left( \hat{c}^{\dagger}_{i} \hat{c}^{\phantom\dagger}_{i+a} + \hat{c}^{\dagger}_{i+a} \hat{c}^{\phantom\dagger}_{i} \right)^2} = e^{\Delta \tau \frac{V}{2} \left( \hat{c}^{\dagger}_{i} \hat{c}^{\phantom\dagger}_{i+a} + \hat{c}^{\dagger}_{i+a} \hat{c}^{\phantom\dagger}_{i} \right)^2} =
\sum_{l= \pm1, \pm 2} \gamma_l e^{ \sqrt{\Delta \tau \frac{V}{2}} \eta_l \left( \hat{c}^{\dagger}_{i} \hat{c}^{\phantom\dagger}_{i+a} + \hat{c}^{\dagger}_{i+a} \hat{c}^{\phantom\dagger}_{i} \right) } = \sum_{l= \pm1, \pm 2} \gamma_l e^{ g \eta_l \left( \hat{c}^{\dagger}_{i}, \hat{c}^{\dagger}_{i+a} \right) O \sum_{l= \pm1, \pm 2} \gamma_l e^{ \sqrt{\Delta \tau \frac{V}{2}} \eta_l \left( \hat{c}^{\dagger}_{i} \hat{c}^{\phantom\dagger}_{i+a} + \hat{c}^{\dagger}_{i+a} \hat{c}^{\phantom\dagger}_{i} \right) } = \sum_{l= \pm1, \pm 2} \gamma_l e^{ g \eta_l \left( \hat{c}^{\dagger}_{i}, \hat{c}^{\dagger}_{i+a} \right) O
\left(\hat{c}^{\phantom\dagger}_{i}, \hat{c}^{\phantom\dagger}_{i+a} \right)^{T} } \left(\hat{c}^{\phantom\dagger}_{i}, \hat{c}^{\phantom\dagger}_{i+a} \right)^{T} }.
\end{equation} \end{equation}
Here is how this translates in the code. Here is how this translates in the code.
\begin{lstlisting}[style=fortran] \begin{lstlisting}[style=fortran]
Allocate(Op_V(Latt%N,N_FL)) Allocate(Op_V(Latt%N,N_FL))
......
...@@ -40,6 +40,7 @@ ...@@ -40,6 +40,7 @@
Use Errors Use Errors
Use MyMats Use MyMats
Use Matrix Use Matrix
use iso_fortran_env, only: output_unit, error_unit
Implicit none Implicit none
Integer :: Nunit, Norb, N_auto Integer :: Nunit, Norb, N_auto
Integer :: no, no1, n, nbins, n_skip, nb, NT, NT1, Lt, N_rebin, N_cov, ierr, N_Back Integer :: no, no1, n, nbins, n_skip, nb, NT, NT1, Lt, N_rebin, N_cov, ierr, N_Back
...@@ -61,8 +62,8 @@ ...@@ -61,8 +62,8 @@
N_auto = 0 N_auto = 0
OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr) OPEN(UNIT=5,FILE='parameters',STATUS='old',ACTION='read',IOSTAT=ierr)
IF (ierr /= 0) THEN IF (ierr /= 0) THEN
WRITE(*,*) 'unable to open <parameters>',ierr WRITE(error_unit,*) 'unable to open <parameters>',ierr
STOP error stop 1
END IF END IF
READ(5,NML=VAR_errors) READ(5,NML=VAR_errors)
CLOSE(5) CLOSE(5)
...@@ -93,8 +94,8 @@ ...@@ -93,8 +94,8 @@
nbins = Nbins - n_skip nbins = Nbins - n_skip
Write(6,*) "Effective # of bins: ", Nbins Write(6,*) "Effective # of bins: ", Nbins
if(Nbins <= 1) then if(Nbins <= 1) then
write (*,*) "Effective # of bins smaller than 2. Analysis impossible!" write (error_unit,*) "Effective # of bins smaller than 2. Analysis impossible!"
stop 1 error stop 1
endif endif
if (mod(Lt-1,2) == 0 ) then if (mod(Lt-1,2) == 0 ) then
Lt_eff = (Lt -1 ) /2 + 1 Lt_eff = (Lt -1 ) /2 + 1
......
...@@ -152,6 +152,9 @@ ...@@ -152,6 +152,9 @@
Type (Unit_cell), private :: Latt_unit Type (Unit_cell), private :: Latt_unit
Integer, private :: L1, L2 Integer, private :: L1, L2
real (Kind=Kind(0.d0)), private :: Ham_T , ham_U, Ham_chem real (Kind=Kind(0.d0)), private :: Ham_T , ham_U, Ham_chem
!!!!! Modifications for Exercise 2
real (Kind=Kind(0.d0)), private :: Ham_Vint
!!!!!
real (Kind=Kind(0.d0)), private :: Dtau, Beta, Theta real (Kind=Kind(0.d0)), private :: Dtau, Beta, Theta
Integer , private :: N_part Integer , private :: N_part
Character (len=64), private :: Model, Lattice_type Character (len=64), private :: Model, Lattice_type
...@@ -189,6 +192,9 @@ ...@@ -189,6 +192,9 @@
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_chem, ham_U, Dtau, Beta, Projector, Theta, Symm, N_part
!!!!! Modifications for Exercise 2
NAMELIST /VAR_t_V/ Ham_T, Ham_chem, Ham_Vint, Dtau, Beta, Projector, Theta, Symm
!!!!!
#ifdef MPI #ifdef MPI
...@@ -226,14 +232,24 @@ ...@@ -226,14 +232,24 @@
Write(6,*) 'For one-dimensional lattices set L2=1' Write(6,*) 'For one-dimensional lattices set L2=1'
stop stop
endif endif
READ(5,NML=VAR_Hubbard_Plain_Vanilla) !!!!! Modifications for Exercise 2
!READ(5,NML=VAR_Hubbard_Plain_Vanilla)
READ(5,NML=VAR_t_V)
!!!!!
CLOSE(5) CLOSE(5)
Ltrot = nint(beta/dtau) Ltrot = nint(beta/dtau)
if (Projector) Thtrot = nint(theta/dtau) if (Projector) Thtrot = nint(theta/dtau)
Ltrot = Ltrot+2*Thtrot Ltrot = Ltrot+2*Thtrot
N_SUN = 1 N_SUN = 1
N_FL = 2 !!!!! Modifications for Exercise 2
!N_FL = 2
N_FL = 1
If (L2 /= 1) then
Write(6,*) "The t_V model is implemented only for L2 = 1"
Stop
Endif
!!!!!
#ifdef MPI #ifdef MPI
Endif Endif
...@@ -252,7 +268,10 @@ ...@@ -252,7 +268,10 @@
CALL MPI_BCAST(Beta ,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)
CALL MPI_BCAST(ham_chem ,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) !!!!! Modifications for Exercise 2
!CALL MPI_BCAST(ham_U ,1, MPI_REAL8 , 0,Group_Comm,ierr)
CALL MPI_BCAST(Ham_Vint ,1, MPI_REAL8 , 0,Group_Comm,ierr)
!!!!!
#endif #endif
! Setup the Bravais lattice ! Setup the Bravais lattice
...@@ -287,7 +306,10 @@ ...@@ -287,7 +306,10 @@
endif endif
Write(50,*) 'dtau,Ltrot_eff: ', dtau,Ltrot Write(50,*) 'dtau,Ltrot_eff: ', dtau,Ltrot
Write(50,*) 't : ', Ham_T Write(50,*) 't : ', Ham_T
Write(50,*) 'Ham_U : ', Ham_U !!!!! Modifications for Exercise 2
!Write(50,*) 'Ham_U : ', Ham_U
Write(50,*) 'Ham_Vint : ', Ham_Vint
!!!!!
Write(50,*) 'Ham_chem : ', Ham_chem Write(50,*) 'Ham_chem : ', Ham_chem
Close(50) Close(50)
#ifdef MPI #ifdef MPI
...@@ -457,13 +479,21 @@ ...@@ -457,13 +479,21 @@
Integer :: nf, I Integer :: nf, I
Real (Kind=Kind(0.d0)) :: X Real (Kind=Kind(0.d0)) :: X
!!!!! Modifications for Exercise 2
Integer :: i2
!!!!!
Allocate(Op_V(Ndim,N_FL)) Allocate(Op_V(Ndim,N_FL))
ham_U = Ham_Vint ! for testing
do nf = 1,N_FL do nf = 1,N_FL
do i = 1, Ndim do i = 1, Ndim
!!!!! Modifications for Exercise 2
Call Op_make(Op_V(i,nf), 1) Call Op_make(Op_V(i,nf), 1)
!Call Op_make(Op_V(i,nf), 2)
!!!!!
enddo enddo
enddo enddo
...@@ -471,11 +501,21 @@ ...@@ -471,11 +501,21 @@
X = 1.d0 X = 1.d0
if (nf == 2) X = -1.d0 if (nf == 2) X = -1.d0
Do i = 1,Ndim Do i = 1,Ndim
Op_V(i,nf)%P(1) = I !!!!! Modifications for Exercise 2
Op_V(i,nf)%O(1,1) = cmplx(1.d0, 0.d0, kind(0.D0)) !Op_V(i,nf)%P(1) = i
Op_V(i,nf)%g = X*SQRT(CMPLX(DTAU*ham_U/2.d0, 0.D0, kind(0.D0))) !Op_V(i,nf)%O(1,1) = cmplx(1.d0, 0.d0, kind(0.D0))
Op_V(i,nf)%alpha = cmplx(0.d0, 0.d0, kind(0.D0)) !Op_V(i,nf)%g = X*SQRT(CMPLX(DTAU*ham_U/2.d0, 0.D0, kind(0.D0)))
!Op_V(i,nf)%alpha = cmplx(0.d0, 0.d0, kind(0.D0))
!Op_V(i,nf)%type = 2
i2 = Latt%nnlist(i,1,0)
Op_V(i,nf)%P(1) = i
Op_V(i,nf)%P(2) = i2
Op_V(i,nf)%O(1,2) = cmplx(1.d0 ,0.d0, kind(0.d0))
Op_V(i,nf)%O(2,1) = cmplx(1.d0 ,0.d0, kind(0.d0))
Op_V(i,nf)%g = sqrt(cmplx(Dtau*Ham_Vint/2.d0, 0.d0, kind(0.d0)))
Op_V(i,nf)%alpha = cmplx(0d0 ,0.d0, kind(0.d0))
Op_V(i,nf)%type = 2 Op_V(i,nf)%type = 2
!!!!!
Call Op_set( Op_V(i,nf) ) Call Op_set( Op_V(i,nf) )
Enddo Enddo
Enddo Enddo
...@@ -600,6 +640,9 @@ ...@@ -600,6 +640,9 @@
Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, Z, ZP,ZS, ZZ, ZXY, ZDen Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, Z, ZP,ZS, ZZ, ZXY, ZDen
Integer :: I,J, imj, nf, Ix, Iy Integer :: I,J, imj, nf, Ix, Iy
Real (Kind=Kind(0.d0)) :: X Real (Kind=Kind(0.d0)) :: X
!!!!! Modifications for Exercise 2
Integer ::I1, J1, no_I, no_J
!!!!!
ZP = PHASE/Real(Phase, kind(0.D0)) ZP = PHASE/Real(Phase, kind(0.D0))
ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0)))
...@@ -626,8 +669,11 @@ ...@@ -626,8 +669,11 @@
Zkin = Zkin* dble(N_SUN) Zkin = Zkin* dble(N_SUN)
Do I = 1,Latt%N Do I = 1,Latt%N
Ix = Latt%nnlist(I,1,0) Ix = Latt%nnlist(I,1,0)
Zkin = Zkin + GRC(I,Ix,1) + GRC(Ix,I,1) & !!!!! Modifications for Exercise 2
& + GRC(I,Ix,2) + GRC(Ix,I,2) !Zkin = Zkin + GRC(I,Ix,1) + GRC(Ix,I,1) &
! & + GRC(I,Ix,2) + GRC(Ix,I,2)
Zkin = Zkin + GRC(I,Ix,1) + GRC(Ix,I,1)
!!!!!
Enddo Enddo
If (L2 > 1) then If (L2 > 1) then
Do I = 1,Latt%N Do I = 1,Latt%N
...@@ -642,7 +688,11 @@ ...@@ -642,7 +688,11 @@
ZPot = cmplx(0.d0, 0.d0, kind(0.D0)) ZPot = cmplx(0.d0, 0.d0, kind(0.D0))
Do I = 1,Ndim Do I = 1,Ndim
ZPot = ZPot + Grc(i,i,1) * Grc(i,i, 2) !!!!! Modifications for Exercise 2
!ZPot = ZPot + Grc(i,i,1) * Grc(i,i, 2)
i1 = Latt%nnlist(i,1,0)
ZPot = ZPot + Grc(i,i,1) * Grc(i1,i1, 1) + Grc(i,i1,1)*Gr(i,i1,1)
!!!!!
Enddo Enddo
Zpot = Zpot*ham_U Zpot = Zpot*ham_U
Obs_scal(2)%Obs_vec(1) = Obs_scal(2)%Obs_vec(1) + Zpot * ZP*ZS Obs_scal(2)%Obs_vec(1) = Obs_scal(2)%Obs_vec(1) + Zpot * ZP*ZS
...@@ -650,7 +700,10 @@ ...@@ -650,7 +700,10 @@
Zrho = cmplx(0.d0,0.d0, kind(0.D0)) Zrho = cmplx(0.d0,0.d0, kind(0.D0))
Do I = 1,Ndim Do I = 1,Ndim
Zrho = Zrho + Grc(i,i,1) + Grc(i,i,2) !!!!! Modifications for Exercise 2
!Zrho = Zrho + Grc(i,i,1) + Grc(i,i,2)
Zrho = Zrho + Grc(i,i,1)
!!!!!
enddo enddo
Obs_scal(3)%Obs_vec(1) = Obs_scal(3)%Obs_vec(1) + Zrho * ZP*ZS Obs_scal(3)%Obs_vec(1) = Obs_scal(3)%Obs_vec(1) + Zrho * ZP*ZS
Obs_scal(4)%Obs_vec(1) = Obs_scal(4)%Obs_vec(1) + (Zkin + Zpot)*ZP*ZS Obs_scal(4)%Obs_vec(1) = Obs_scal(4)%Obs_vec(1) + (Zkin + Zpot)*ZP*ZS
...@@ -660,25 +713,47 @@ ...@@ -660,25 +713,47 @@
Obs_eq(I)%Ave_sign = Obs_eq(I)%Ave_sign + Real(ZS,kind(0.d0)) Obs_eq(I)%Ave_sign = Obs_eq(I)%Ave_sign + Real(ZS,kind(0.d0))
Enddo Enddo
Do I = 1,Latt%N !!!!! Modifications for Exercise 2
Do J = 1,Latt%N !Do I = 1,Latt%N
imj = latt%imj(I,J) ! Do J = 1,Latt%N
ZXY = GRC(I,J,1) * GR(I,J,2) + GRC(I,J,2) * GR(I,J,1) ! imj = latt%imj(I,J)
ZZ = GRC(I,J,1) * GR(I,J,1) + GRC(I,J,2) * GR(I,J,2) + & ! ZXY = GRC(I,J,1) * GR(I,J,2) + GRC(I,J,2) * GR(I,J,1)
(GRC(I,I,2) - GRC(I,I,1))*(GRC(J,J,2) - GRC(J,J,1)) ! ZZ = GRC(I,J,1) * GR(I,J,1) + GRC(I,J,2) * GR(I,J,2) + &
! (GRC(I,I,2) - GRC(I,I,1))*(GRC(J,J,2) - GRC(J,J,1))
ZDen = (GRC(I,I,1) + GRC(I,I,2)) * (GRC(I,I,1) + GRC(I,I,2)) + &
& GRC(I,J,1) * GR(I,J,1) + GRC(I,J,2) * GR(I,J,2) ! ZDen = (GRC(I,I,1) + GRC(I,I,2)) * (GRC(I,I,1) + GRC(I,I,2)) + &
Obs_eq(1)%Obs_Latt(imj,1,1,1) = Obs_eq(1)%Obs_Latt(imj,1,1,1) + (GRC(I,J,1) + GRC(I,J,2))*ZP*ZS ! & GRC(I,J,1) * GR(I,J,1) + GRC(I,J,2) * GR(I,J,2)
Obs_eq(2)%Obs_Latt(imj,1,1,1) = Obs_eq(2)%Obs_Latt(imj,1,1,1) + ZZ *ZP*ZS ! Obs_eq(1)%Obs_Latt(imj,1,1,1) = Obs_eq(1)%Obs_Latt(imj,1,1,1) + (GRC(I,J,1) + GRC(I,J,2))*ZP*ZS
Obs_eq(3)%Obs_Latt(imj,1,1,1) = Obs_eq(3)%Obs_Latt(imj,1,1,1) + ZXY *ZP*ZS ! Obs_eq(2)%Obs_Latt(imj,1,1,1) = Obs_eq(2)%Obs_Latt(imj,1,1,1) + ZZ *ZP*ZS
Obs_eq(4)%Obs_Latt(imj,1,1,1) = Obs_eq(4)%Obs_Latt(imj,1,1,1) + (2.d0*ZXY + ZZ)*ZP*ZS/3.d0 ! Obs_eq(3)%Obs_Latt(imj,1,1,1) = Obs_eq(3)%Obs_Latt(imj,1,1,1) + ZXY *ZP*ZS
Obs_eq(5)%Obs_Latt(imj,1,1,1) = Obs_eq(5)%Obs_Latt(imj,1,1,1) + ZDen * ZP * ZS ! Obs_eq(4)%Obs_Latt(imj,1,1,1) = Obs_eq(4)%Obs_Latt(imj,1,1,1) + (2.d0*ZXY + ZZ)*ZP*ZS/3.d0
! Obs_eq(5)%Obs_Latt(imj,1,1,1) = Obs_eq(5)%Obs_Latt(imj,1,1,1) + ZDen * ZP * ZS
! enddo
! Obs_eq(5)%Obs_Latt0(1) = Obs_eq(5)%Obs_Latt0(1) + (GRC(I,I,1) + GRC(I,I,2)) * ZP*ZS
!enddo
Z = cmplx(dble(N_SUN), 0.d0, kind(0.D0))
Do I1 = 1,Ndim
I = I1 !List(I1,1)
no_I = 1 !List(I1,2)
Do J1 = 1,Ndim
J = J1 !List(J1,1)
no_J = 1 !List(J1,2)
imj = latt%imj(I,J)
! Green
Obs_eq(1)%Obs_Latt(imj,1,no_I,no_J) = Obs_eq(1)%Obs_Latt(imj,1,no_I,no_J) + &
& Z * GRC(I1,J1,1) * ZP*ZS ! Green
Obs_eq(2)%Obs_Latt(imj,1,no_I,no_J) = Obs_eq(2)%Obs_Latt(imj,1,no_I,no_J) + &
& Z * GRC(I1,J1,1) * GR(I1,J1,1) * ZP*ZS! SpinZ
Obs_eq(3)%Obs_Latt(imj,1,no_I,no_J) = Obs_eq(3)%Obs_Latt(imj,1,no_I,no_J) + &
& Z * GRC(I1,J1,1) * GR(I1,J1,1) * ZP*ZS ! SpinXY
Obs_eq(4)%Obs_Latt(imj,1,no_I,no_J) = Obs_eq(4)%Obs_Latt(imj,1,no_I,no_J) + &
& ( GRC(I1,I1,1) * GRC(J1,J1,1) * Z + &
& GRC(I1,J1,1) * GR(I1,J1,1 ) ) * Z * ZP*ZS ! Den
enddo enddo
Obs_eq(5)%Obs_Latt0(1) = Obs_eq(5)%Obs_Latt0(1) + (GRC(I,I,1) + GRC(I,I,2)) * ZP*ZS Obs_eq(4)%Obs_Latt0(no_I) = Obs_eq(4)%Obs_Latt0(no_I) + Z * GRC(I1,I1,1) * ZP * ZS
enddo enddo
!!!!!
...@@ -732,30 +807,53 @@ ...@@ -732,30 +807,53 @@
Obs_tau(I)%Ave_sign = Obs_tau(I)%Ave_sign + Real(ZS,kind(0.d0)) Obs_tau(I)%Ave_sign = Obs_tau(I)%Ave_sign + Real(ZS,kind(0.d0))
Enddo Enddo
Endif Endif
Do I = 1,Latt%N !!!!! Modifications for Exercise 2
Do J = 1,Latt%N ! Do I = 1,Latt%N
imj = latt%imj(I,J) ! Do J = 1,Latt%N
! imj = latt%imj(I,J)
ZZ = (GTT(I,I,1) - GTT(I,I,2) ) * ( G00(J,J,1) - G00(J,J,2) ) & !
& - G0T(J,I,1) * GT0(I,J,1) - G0T(J,I,2) * GT0(I,J,2) ! ZZ = (GTT(I,I,1) - GTT(I,I,2) ) * ( G00(J,J,1) - G00(J,J,2) ) &
ZXY = - G0T(J,I,1) * GT0(I,J,2) - G0T(J,I,2) * GT0(I,J,1) ! & - G0T(J,I,1) * GT0(I,J,1) - G0T(J,I,2) * GT0(I,J,2)
! ZXY = - G0T(J,I,1) * GT0(I,J,2) - G0T(J,I,2) * GT0(I,J,1)
!
ZDen = (cmplx(2.d0,0.d0,kind(0.d0)) - GTT(I,I,1) - GTT(I,I,2) ) * & !
& (cmplx(2.d0,0.d0,kind(0.d0)) - G00(J,J,1) - G00(J,J,2) ) & ! ZDen = (cmplx(2.d0,0.d0,kind(0.d0)) - GTT(I,I,1) - GTT(I,I,2) ) * &
& -G0T(J,I,1) * GT0(I,J,1) - G0T(J,I,2) * GT0(I,J,2) ! & (cmplx(2.d0,0.d0,kind(0.d0)) - G00(J,J,1) - G00(J,J,2) ) &
! & -G0T(J,I,1) * GT0(I,J,1) - G0T(J,I,2) * GT0(I,J,2)
Obs_tau(1)%Obs_Latt(imj,NT+1,1,1) = Obs_tau(1)%Obs_Latt(imj,NT+1,1,1) + (GT0(I,J,1) + GT0(I,J,2))*ZP*ZS !
Obs_tau(2)%Obs_Latt(imj,NT+1,1,1) = Obs_tau(2)%Obs_Latt(imj,NT+1,1,1) + ZZ *ZP*ZS ! Obs_tau(1)%Obs_Latt(imj,NT+1,1,1) = Obs_tau(1)%Obs_Latt(imj,NT+1,1,1) + (GT0(I,J,1) + GT0(I,J,2))*ZP*ZS
Obs_tau(3)%Obs_Latt(imj,NT+1,1,1) = Obs_tau(3)%Obs_Latt(imj,NT+1,1,1) + ZXY *ZP*ZS ! Obs_tau(2)%Obs_Latt(imj,NT+1,1,1) = Obs_tau(2)%Obs_Latt(imj,NT+1,1,1) + ZZ *ZP*ZS
Obs_tau(4)%Obs_Latt(imj,NT+1,1,1) = Obs_tau(4)%Obs_Latt(imj,NT+1,1,1) + (2.d0*ZXY + ZZ)*ZP*ZS/3.d0