!*********************************************************** SUBROUTINE getzvecs(s,sI,zvecs,zvecsinit,nt2,nt1, & pion,ztau2x2,num) IMPLICIT INTEGER(i-n) IMPLICIT DOUBLE PRECISION(a-h,o-y) IMPLICIT COMPLEX*16(z) INCLUDE "input.f90" DIMENSION s(0:L-1,0:L-1,0:L-1,0:Lt-1) DIMENSION sI(0:L-1,0:L-1,0:L-1,0:Lt-1,1:3) DIMENSION zvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:num-1) DIMENSION zvecsinit(0:L-1,0:L-1,0:L-1,0:1,0:1,0:num-1) DIMENSION pion(0:L-1,0:L-1,0:L-1,0:Lt-1,1:3) DIMENSION zpi2x2(0:1,0:1) DIMENSION zsI2x2(0:1,0:1) DIMENSION zsSI2x2(0:1,0:1,0:1,0:1) DIMENSION ztau2x2(0:1,0:1,0:3) ! nt2 > nt1 INCLUDE "improve.f90" DO ni = 0,1; DO ns = 0,1 DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 DO npart = 0,num-1 zvecs(nx,ny,nz,nt1,ns,ni,npart) = zvecsinit(nx,ny,nz,ns,ni,npart) END DO END DO; END DO; END DO END DO; END DO DO nt = nt1+1, nt2 DO np = 0,num-1 DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 DO ni = 0,1; DO ns = 0,1 zvecs(nx,ny,nz,nt,ns,ni,np) = zvecs(nx,ny,nz,nt-1,ns,ni,np) & * (1.D0-6.D0*w0_N*h+CDSQRT(-c0*atovera*(1.D0,0.D0))*s(nx,ny,nz,nt-1)) zvecs(nx,ny,nz,nt,ns,ni,np) = zvecs(nx,ny,nz,nt,ns,ni,np) & + w1_N*h*zvecs(MOD(nx+1,L),ny,nz,nt-1,ns,ni,np) & + w1_N*h*zvecs(MOD(nx-1+L,L),ny,nz,nt-1,ns,ni,np) & + w1_N*h*zvecs(nx,MOD(ny+1,L),nz,nt-1,ns,ni,np) & + w1_N*h*zvecs(nx,MOD(ny-1+L,L),nz,nt-1,ns,ni,np) & + w1_N*h*zvecs(nx,ny,MOD(nz+1,L),nt-1,ns,ni,np) & + w1_N*h*zvecs(nx,ny,MOD(nz-1+L,L),nt-1,ns,ni,np) IF (improveN >= 1) THEN zvecs(nx,ny,nz,nt,ns,ni,np) = zvecs(nx,ny,nz,nt,ns,ni,np) & - w2_N*h*zvecs(MOD(nx+2,L),ny,nz,nt-1,ns,ni,np) & - w2_N*h*zvecs(MOD(nx-2+L,L),ny,nz,nt-1,ns,ni,np) & - w2_N*h*zvecs(nx,MOD(ny+2,L),nz,nt-1,ns,ni,np) & - w2_N*h*zvecs(nx,MOD(ny-2+L,L),nz,nt-1,ns,ni,np) & - w2_N*h*zvecs(nx,ny,MOD(nz+2,L),nt-1,ns,ni,np) & - w2_N*h*zvecs(nx,ny,MOD(nz-2+L,L),nt-1,ns,ni,np) END IF IF (improveN == 2) THEN zvecs(nx,ny,nz,nt,ns,ni,np) = zvecs(nx,ny,nz,nt,ns,ni,np) & + w3_N*h*zvecs(MOD(nx+3,L),ny,nz,nt-1,ns,ni,np) & + w3_N*h*zvecs(MOD(nx-3+L,L),ny,nz,nt-1,ns,ni,np) & + w3_N*h*zvecs(nx,MOD(ny+3,L),nz,nt-1,ns,ni,np) & + w3_N*h*zvecs(nx,MOD(ny-3+L,L),nz,nt-1,ns,ni,np) & + w3_N*h*zvecs(nx,ny,MOD(nz+3,L),nt-1,ns,ni,np) & + w3_N*h*zvecs(nx,ny,MOD(nz-3+L,L),nt-1,ns,ni,np) END IF END DO; END DO END DO; END DO; END DO END DO DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 DO nii = 0,1; DO ni = 0,1 DO nss = 0,1; DO ns = 0,1 zsSI2x2(ns,nss,ni,nii) = (0.D0,0.D0) END DO; END DO END DO; END DO zsI2x2(0,0) = sI(nx,ny,nz,nt-1,3) zsI2x2(1,1) = -sI(nx,ny,nz,nt-1,3) zsI2x2(0,1) = sI(nx,ny,nz,nt-1,1) - (0.D0,1.D0)*sI(nx,ny,nz,nt-1,2) zsI2x2(1,0) = sI(nx,ny,nz,nt-1,1) + (0.D0,1.D0)*sI(nx,ny,nz,nt-1,2) DO nii = 0,1; DO ni = 0,1 DO ns = 0,1 zsSI2x2(ns,ns,ni,nii) = zsSI2x2(ns,ns,ni,nii) & + (0.D0,1.D0)*CDSQRT(cI*atovera*(1.D0,0.D0))*zsI2x2(ni,nii) END DO END DO; END DO DO iso = 1,3 pi1 = 0.D0 & + o1/2.D0*pion(MOD(nx+1,L),ny,nz,nt-1,iso) & - o1/2.D0*pion(MOD(nx-1+L,L),ny,nz,nt-1,iso) & - o2/2.D0*pion(MOD(nx+2,L),ny,nz,nt-1,iso) & + o2/2.D0*pion(MOD(nx-2+L,L),ny,nz,nt-1,iso) & + o3/2.D0*pion(MOD(nx+3,L),ny,nz,nt-1,iso) & - o3/2.D0*pion(MOD(nx-3+L,L),ny,nz,nt-1,iso) pi2 = 0.D0 & + o1/2.D0*pion(nx,MOD(ny+1,L),nz,nt-1,iso) & - o1/2.D0*pion(nx,MOD(ny-1+L,L),nz,nt-1,iso) & - o2/2.D0*pion(nx,MOD(ny+2,L),nz,nt-1,iso) & + o2/2.D0*pion(nx,MOD(ny-2+L,L),nz,nt-1,iso) & + o3/2.D0*pion(nx,MOD(ny+3,L),nz,nt-1,iso) & - o3/2.D0*pion(nx,MOD(ny-3+L,L),nz,nt-1,iso) pi3 = 0.D0 & + o1/2.D0*pion(nx,ny,MOD(nz+1,L),nt-1,iso) & - o1/2.D0*pion(nx,ny,MOD(nz-1+L,L),nt-1,iso) & - o2/2.D0*pion(nx,ny,MOD(nz+2,L),nt-1,iso) & + o2/2.D0*pion(nx,ny,MOD(nz-2+L,L),nt-1,iso) & + o3/2.D0*pion(nx,ny,MOD(nz+3,L),nt-1,iso) & - o3/2.D0*pion(nx,ny,MOD(nz-3+L,L),nt-1,iso) zpi2x2(0,0) = pi3 zpi2x2(1,1) = -pi3 zpi2x2(0,1) = pi1 - (0.D0,1.D0)*pi2 zpi2x2(1,0) = pi1 + (0.D0,1.D0)*pi2 DO nii = 0,1; DO ni = 0,1 DO nss = 0,1; DO ns = 0,1 zsSI2x2(ns,nss,ni,nii) = zsSI2x2(ns,nss,ni,nii) & - gA*atovera/(2.D0*fpi*SQRT(qpi3))*zpi2x2(ns,nss) & * ztau2x2(ni,nii,iso) END DO; END DO END DO; END DO END DO DO np = 0,num-1 DO nii = 0,1; DO ni = 0,1 DO nss = 0,1; DO ns = 0,1 zvecs(nx,ny,nz,nt,ns,ni,np) = zvecs(nx,ny,nz,nt,ns,ni,np) & + zsSI2x2(ns,nss,ni,nii)*zvecs(nx,ny,nz,nt-1,nss,nii,np) END DO; END DO END DO; END DO END DO END DO; END DO; END DO END DO END SUBROUTINE getzvecs !*********************************************************** SUBROUTINE getzdualvecs(s,sI,zdualvecs,zdualvecsinit, & nt2,nt1,pion,ztau2x2,num) IMPLICIT INTEGER(i-n) IMPLICIT DOUBLE PRECISION(a-h,o-y) IMPLICIT COMPLEX*16(z) INCLUDE "input.f90" DIMENSION s(0:L-1,0:L-1,0:L-1,0:Lt-1) DIMENSION sI(0:L-1,0:L-1,0:L-1,0:Lt-1,1:3) DIMENSION zdualvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:num-1) DIMENSION zdualvecsinit(0:L-1,0:L-1,0:L-1,0:1,0:1,0:num-1) DIMENSION pion(0:L-1,0:L-1,0:L-1,0:Lt-1,1:3) DIMENSION zpi2x2(0:1,0:1) DIMENSION zsI2x2(0:1,0:1) DIMENSION zsSI2x2(0:1,0:1,0:1,0:1) DIMENSION ztau2x2(0:1,0:1,0:9) ! nt2 > nt1 INCLUDE "improve.f90" DO ni = 0,1; DO ns = 0,1 DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 DO npart = 0,num-1 zdualvecs(nx,ny,nz,nt2,ns,ni,npart) = zdualvecsinit(nx,ny,nz,ns,ni,npart) END DO END DO; END DO; END DO END DO; END DO DO nt = nt2,nt1+1,-1 DO np = 0,num-1 DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 DO ni = 0,1; DO ns = 0,1 zdualvecs(nx,ny,nz,nt-1,ns,ni,np) = zdualvecs(nx,ny,nz,nt,ns,ni,np) & * (1.D0-6.D0*w0_N*h+CDSQRT(-c0*atovera*(1.D0,0.D0))*s(nx,ny,nz,nt-1)) zdualvecs(nx,ny,nz,nt-1,ns,ni,np) = zdualvecs(nx,ny,nz,nt-1,ns,ni,np) & + w1_N*h*zdualvecs(MOD(nx+1,L),ny,nz,nt,ns,ni,np) & + w1_N*h*zdualvecs(MOD(nx-1+L,L),ny,nz,nt,ns,ni,np) & + w1_N*h*zdualvecs(nx,MOD(ny+1,L),nz,nt,ns,ni,np) & + w1_N*h*zdualvecs(nx,MOD(ny-1+L,L),nz,nt,ns,ni,np) & + w1_N*h*zdualvecs(nx,ny,MOD(nz+1,L),nt,ns,ni,np) & + w1_N*h*zdualvecs(nx,ny,MOD(nz-1+L,L),nt,ns,ni,np) IF (improveN >= 1) THEN zdualvecs(nx,ny,nz,nt-1,ns,ni,np) = zdualvecs(nx,ny,nz,nt-1,ns,ni,np) & - w2_N*h*zdualvecs(MOD(nx+2,L),ny,nz,nt,ns,ni,np) & - w2_N*h*zdualvecs(MOD(nx-2+L,L),ny,nz,nt,ns,ni,np) & - w2_N*h*zdualvecs(nx,MOD(ny+2,L),nz,nt,ns,ni,np) & - w2_N*h*zdualvecs(nx,MOD(ny-2+L,L),nz,nt,ns,ni,np) & - w2_N*h*zdualvecs(nx,ny,MOD(nz+2,L),nt,ns,ni,np) & - w2_N*h*zdualvecs(nx,ny,MOD(nz-2+L,L),nt,ns,ni,np) END IF IF (improveN == 2) THEN zdualvecs(nx,ny,nz,nt-1,ns,ni,np) = zdualvecs(nx,ny,nz,nt-1,ns,ni,np) & + w3_N*h*zdualvecs(MOD(nx+3,L),ny,nz,nt,ns,ni,np) & + w3_N*h*zdualvecs(MOD(nx-3+L,L),ny,nz,nt,ns,ni,np) & + w3_N*h*zdualvecs(nx,MOD(ny+3,L),nz,nt,ns,ni,np) & + w3_N*h*zdualvecs(nx,MOD(ny-3+L,L),nz,nt,ns,ni,np) & + w3_N*h*zdualvecs(nx,ny,MOD(nz+3,L),nt,ns,ni,np) & + w3_N*h*zdualvecs(nx,ny,MOD(nz-3+L,L),nt,ns,ni,np) END IF END DO; END DO END DO; END DO; END DO END DO DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 DO nii = 0,1; DO ni = 0,1 DO nss = 0,1; DO ns = 0,1 zsSI2x2(ns,nss,ni,nii) = (0.D0,0.D0) END DO; END DO END DO; END DO zsI2x2(0,0) = sI(nx,ny,nz,nt-1,3) zsI2x2(1,1) = -sI(nx,ny,nz,nt-1,3) zsI2x2(0,1) = sI(nx,ny,nz,nt-1,1) - (0.D0,1.D0)*sI(nx,ny,nz,nt-1,2) zsI2x2(1,0) = sI(nx,ny,nz,nt-1,1) + (0.D0,1.D0)*sI(nx,ny,nz,nt-1,2) DO nii = 0,1; DO ni = 0,1 DO ns = 0,1 zsSI2x2(ns,ns,ni,nii) = zsSI2x2(ns,ns,ni,nii) & + (0.D0,1.D0)*CDSQRT(cI*atovera*(1.D0,0.D0))*zsI2x2(ni,nii) END DO END DO; END DO DO iso = 1,3 pi1 = 0.D0 & + o1/2.D0*pion(MOD(nx+1,L),ny,nz,nt-1,iso) & - o1/2.D0*pion(MOD(nx-1+L,L),ny,nz,nt-1,iso) & - o2/2.D0*pion(MOD(nx+2,L),ny,nz,nt-1,iso) & + o2/2.D0*pion(MOD(nx-2+L,L),ny,nz,nt-1,iso) & + o3/2.D0*pion(MOD(nx+3,L),ny,nz,nt-1,iso) & - o3/2.D0*pion(MOD(nx-3+L,L),ny,nz,nt-1,iso) pi2 = 0.D0 & + o1/2.D0*pion(nx,MOD(ny+1,L),nz,nt-1,iso) & - o1/2.D0*pion(nx,MOD(ny-1+L,L),nz,nt-1,iso) & - o2/2.D0*pion(nx,MOD(ny+2,L),nz,nt-1,iso) & + o2/2.D0*pion(nx,MOD(ny-2+L,L),nz,nt-1,iso) & + o3/2.D0*pion(nx,MOD(ny+3,L),nz,nt-1,iso) & - o3/2.D0*pion(nx,MOD(ny-3+L,L),nz,nt-1,iso) pi3 = 0.D0 & + o1/2.D0*pion(nx,ny,MOD(nz+1,L),nt-1,iso) & - o1/2.D0*pion(nx,ny,MOD(nz-1+L,L),nt-1,iso) & - o2/2.D0*pion(nx,ny,MOD(nz+2,L),nt-1,iso) & + o2/2.D0*pion(nx,ny,MOD(nz-2+L,L),nt-1,iso) & + o3/2.D0*pion(nx,ny,MOD(nz+3,L),nt-1,iso) & - o3/2.D0*pion(nx,ny,MOD(nz-3+L,L),nt-1,iso) zpi2x2(0,0) = pi3 zpi2x2(1,1) = -pi3 zpi2x2(0,1) = pi1 - (0.D0,1.D0)*pi2 zpi2x2(1,0) = pi1 + (0.D0,1.D0)*pi2 DO nii = 0,1; DO ni = 0,1 DO nss = 0,1; DO ns = 0,1 zsSI2x2(ns,nss,ni,nii) = zsSI2x2(ns,nss,ni,nii) & - gA*atovera/(2.D0*fpi*SQRT(qpi3))*zpi2x2(ns,nss) & * ztau2x2(ni,nii,iso) END DO; END DO END DO; END DO END DO DO np = 0,num-1 DO nii = 0,1; DO ni = 0,1 DO nss = 0,1; DO ns = 0,1 zdualvecs(nx,ny,nz,nt-1,nss,nii,np) = zdualvecs(nx,ny,nz,nt-1,nss,nii,np) & + zsSI2x2(ns,nss,ni,nii)*zdualvecs(nx,ny,nz,nt,ns,ni,np) END DO; END DO END DO; END DO END DO END DO; END DO; END DO END DO END SUBROUTINE getzdualvecs !********************************************************** SUBROUTINE getmidoverlap(zvecs,zdualvecs,zdeter,aldeterabs,nt1,nt2) IMPLICIT integer(i-n) IMPLICIT double precision(a-h,o-y) IMPLICIT complex*16(z) INCLUDE "input.f90" DIMENSION zvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:n_f-1) DIMENSION zdualvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:n_f-1) DIMENSION zwave(0:L-1,0:L-1,0:L-1,0:1,0:1,0:n_f-1) DIMENSION zcorrmatrix(0:n_f-1,0:n_f-1) DIMENSION indx(0:n_f-1) DO np1 = 0,n_f-1; DO np2 = 0,n_f-1 zcorrmatrix(np2,np1) = 0.D0 DO ni = 0,1; DO ns = 0,1 DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 zcorrmatrix(np2,np1) = & zcorrmatrix(np2,np1) + & zdualvecs(nx,ny,nz,nt2,ns,ni,np2) & *zvecs(nx,ny,nz,nt1,ns,ni,np1)/L**3 END DO; END DO; END DO END DO; END DO END DO; END DO CALL zludcmp_ratio(zcorrmatrix,zdeter,aldeterabs,n_f) END SUBROUTINE getmidoverlap !********************************************************** SUBROUTINE getinvcorr(zvecs,zdualvecs,zldeter,zcorrmatrix,zcorrinv,nt) IMPLICIT integer(i-n) IMPLICIT double precision(a-h,o-y) IMPLICIT complex*16(z) INCLUDE "input.f90" DIMENSION zvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:n_f-1) DIMENSION zdualvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:n_f-1) DIMENSION zcorrmatrix(0:n_f-1,0:n_f-1) DIMENSION zcorrinv(0:n_f-1,0:n_f-1) DO np1 = 0,n_f-1; DO np2 = 0,n_f-1 zcorrmatrix(np2,np1) = 0.D0 DO ni = 0,1; DO ns = 0,1 DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 zcorrmatrix(np2,np1) = & zcorrmatrix(np2,np1) + & zdualvecs(nx,ny,nz,nt,ns,ni,np2)* & zvecs(nx,ny,nz,nt,ns,ni,np1)/L**3 END DO; END DO; END DO END DO; END DO END DO; END DO CALL zdoinv_log(zcorrmatrix,zcorrinv,zldeter,n_f) END SUBROUTINE getinvcorr ! ********************************************************** SUBROUTINE dV(zvecs,zdualvecs,ztau2x2,zcorrinv,zdVall) IMPLICIT integer(i-n) IMPLICIT double precision(a-h,o-y) IMPLICIT complex*16(z) INCLUDE "input.f90" DIMENSION zvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:n_f-1) DIMENSION zdualvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:n_f-1) DIMENSION ztau2x2(0:1,0:1,0:3) DIMENSION zcorrinv(0:n_f-1,0:n_f-1) DIMENSION zdVall(0:L-1,0:L-1,0:L-1,0:Lt-1,0:3,0:3) DIMENSION zmat(0:1,0:1,0:1,0:1) DO nt = 0,Lt-1 DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 DO nii = 0,1; DO ni = 0,1; DO nss = 0,1; DO ns = 0,1 zmat(ns,nss,ni,nii) = 0.D0 DO np2 = 0,n_f-1; DO np1 = 0,n_f-1 zmat(ns,nss,ni,nii) = & zmat(ns,nss,ni,nii) & + zdualvecs(nx,ny,nz,nt+1, & ns,ni,np1) & *zvecs(nx,ny,nz,nt,nss,nii,np2) & *zcorrinv(np2,np1) END DO; END DO END DO; END DO; END DO; END DO DO iso = 0,3; DO nsvec = 0,3 zdVall(nx,ny,nz,nt,nsvec,iso) = 0.D0 DO nii = 0,1; DO ni = 0,1; DO nss = 0,1; DO ns = 0,1 zdVall(nx,ny,nz,nt,nsvec,iso) = & zdVall(nx,ny,nz,nt,nsvec,iso) & + zmat(ns,nss,ni,nii) & *ztau2x2(ns,nss,nsvec) & *ztau2x2(ni,nii,iso) END DO; END DO; END DO; END DO END DO; END DO END DO; END DO; END DO END DO END SUBROUTINE dV ! ********************************************************** SUBROUTINE generalmidoverlap(zvecs,zdualvecs,zdeter,aldeterabs,nwhich1,nt1,nwhich2,nt2) IMPLICIT INTEGER(i-n) IMPLICIT DOUBLE PRECISION(a-h,o-y) IMPLICIT COMPLEX*16(z) INCLUDE "input.f90" DIMENSION zvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:n_f-1) DIMENSION zdualvecs(0:L-1,0:L-1,0:L-1,0:Lt,0:1,0:1,0:n_f-1) DIMENSION zcorrmatrix(0:n_f-1,0:n_f-1) DIMENSION nwhich1(0:n_f-1) DIMENSION nwhich2(0:n_f-1) DIMENSION indx(0:n_f-1) DO np1 = 0,n_f-1; DO np2 = 0,n_f-1 zcorrmatrix(np2,np1) = 0.D0 DO ni = 0,1; DO ns = 0,1 DO nz = 0,L-1; DO ny = 0,L-1; DO nx = 0,L-1 zcorrmatrix(np2,np1) = zcorrmatrix(np2,np1) & + zdualvecs(nx,ny,nz,nt2,ns,ni,nwhich2(np2)) & * zvecs(nx,ny,nz,nt1,ns,ni,nwhich1(np1))/L**3 END DO; END DO; END DO END DO; END DO END DO; END DO CALL zludcmp_ratio(zcorrmatrix,zdeter,aldeterabs,n_f) END SUBROUTINE generalmidoverlap ! **********************************************************