SUBROUTINE WriteHistogram()
INTEGER :: i,IterRead, io1, io2, io3, Tot_No_Unique_Dets,ierr,val,j
integer :: posn, FinalPop
real(dp) :: ShiftRead,AllERead,NumParts,DDOT
real(dp),dimension(lenof_sign) :: norm,norm1,norm2,norm3
real(dp) :: norm_c,norm1_c,norm2_c,norm3_c
real(dp) :: AvVarEnergy, VarEnergy, GroundE_Ever, ProjGroundE
real(dp) , allocatable :: CKN(:),HOrderedHist(:),HOrderedInstHist(:)
integer , allocatable :: ExpandedWalkerDets(:,:)
CHARACTER(len=22) :: abstr,abstr2
LOGICAL :: exists
character(len=*), parameter :: t_r='WriteHistogram'
!This will open a file called SpawnHist-"Iter" on unit number 17.
abstr = 'SpawnHist-'//str(Iter)
IF(iProcIndex.eq.0) THEN
write(stdout,*) "Writing out the average wavevector up to iteration number: ", Iter
CALL neci_flush(stdout)
end if
IF(iProcIndex.eq.0) THEN
AllHistogram(:,:)=0.0_dp
AllInstHist(:,:)=0.0_dp
AllAvAnnihil(:,:)=0.0_dp
AllInstAnnihil(:,:)=0.0_dp
end if
CALL MPIReduce(Histogram,MPI_SUM,AllHistogram)
CALL MPIReduce(InstHist,MPI_SUM,AllInstHist)
CALL MPIReduce(InstAnnihil,MPI_SUM,AllInstAnnihil)
CALL MPIReduce(AvAnnihil,MPI_SUM,AllAvAnnihil)
BeforeNormHist(:) = AllHistogram(1,:)
IF(iProcIndex.eq.0) THEN
norm=0.0_dp
norm1=0.0_dp
norm2=0.0_dp
norm3=0.0_dp
do i=1,Det
do j=1,lenof_sign
norm(j)=norm(j)+AllHistogram(j,i)**2
norm1(j)=norm1(j)+AllInstHist(j,i)**2
norm2(j)=norm2(j)+AllInstAnnihil(j,i)**2
norm3(j)=norm3(j)+AllAvAnnihil(j,i)**2
end do
end do
#ifdef CMPLX_
norm_c=SQRT(sum(norm))
norm1_c=SQRT(sum(norm1))
norm2_c=SQRT(sum(norm2))
norm3_c=SQRT(sum(norm3))
#else
norm=SQRT(norm)
norm1=SQRT(norm1)
norm2=SQRT(norm2)
norm3=SQRT(norm3)
#endif
do i=1,Det
do j=1,lenof_sign
#ifdef CMPLX_
AllHistogram(j,i)=AllHistogram(j,i)/norm_c
AllInstHist(j,i)=AllInstHist(j,i)/norm1_c
IF(.not. near_zero(norm2_c)) THEN
AllInstAnnihil(j,i)=AllInstAnnihil(j,i)/norm2_c
end if
IF(.not. near_zero(norm3_c)) THEN
AllAvAnnihil(j,i)=AllAvAnnihil(j,i)/norm3_c
end if
#else
AllHistogram(j,i)=AllHistogram(j,i)/norm(j)
AllInstHist(j,i)=AllInstHist(j,i)/norm1(j)
IF(.not. near_zero(norm2(j))) THEN
AllInstAnnihil(j,i)=AllInstAnnihil(1,i)/norm2(j)
end if
IF(.not. near_zero(norm3(j))) THEN
AllAvAnnihil(j,i)=AllAvAnnihil(j,i)/norm3(j)
end if
#endif
end do
end do
io1 = get_free_unit()
open(io1,FILE=abstr,STATUS='UNKNOWN')
abstr = 'Energies-'//str(Iter - iWriteHistEvery)
abstr2 = 'Energies-'//str(Iter)
io2 = get_free_unit()
open(io2,FILE=abstr2,STATUS='UNKNOWN')
INQUIRE(FILE=abstr,EXIST=exists)
IF(exists) THEN
io3 = get_free_unit()
open(io3,FILE=abstr,STATUS='OLD',POSITION='REWIND',ACTION='READ')
do while(.true.)
read(io3,"(I13,3G25.16)",END=99) IterRead,ShiftRead,AllERead,NumParts
write(io2,"(I13,3G25.16)") IterRead,ShiftRead,AllERead,NumParts
end do
99 CONTINUE
#ifdef CMPLX_
IF(near_zero(AllHFOut(1))) then
write(io2,"(I13,3G25.16)") Iter,DiagSft,AllERead,SUM(AllTotPartsLastOutput)
ELSE
write(io2,"(I13,3G25.16)") Iter,DiagSft,AllENumOut/AllHFOut,SUM(AllTotPartsLastOutput)
end if
#else
IF(near_zero(AllHFOut(1))) THEN
write(io2,"(I13,3G25.16)") Iter,DiagSft(1),AllERead,AllTotPartsLastOutput(1)
ELSE
write(io2,"(I13,3G25.16)") Iter,DiagSft(1),AllENumOut(1)/AllHFOut(1),AllTotPartsLastOutput(1)
end if
#endif
close(io2)
close(io3)
ELSE
open(io2,FILE=abstr2,STATUS='UNKNOWN')
#ifdef CMPLX_
write(io2,"(I13,3G25.16)") Iter,DiagSft,AllENumOut/AllHFOut,SUM(AllTotPartsLastOutput)
#else
write(io2,"(I13,3G25.16)") Iter,DiagSft(1),AllENumOut(1)/AllHFOut(1),AllTotPartsLastOutput(1)
#endif
close(io2)
end if
norm=0.0_dp
norm1=0.0_dp
Tot_No_Unique_Dets = 0
do i=1,Det
posn=binary_search_ilut(CurrentDets(:,1:TotWalkers), FCIDETS(:,i), NifD+1)
if (posn.lt.0) then
FinalPop = 0
else
FinalPop = int(CurrentDets(NifD+1,posn))
end if
do j=1,lenof_sign
norm(j)=norm(j)+AllHistogram(j,i)**2
norm1(j)=norm1(j)+AllAvAnnihil(j,i)**2
end do
IF(lenof_sign.eq.1) THEN
write(io1,"(I13,6G25.16,I13,G25.16)") i, &
AllHistogram(1,i), norm, AllInstHist(1,i), &
AllInstAnnihil(1,i), AllAvAnnihil(1,i), norm1, &
FinalPop, BeforeNormHist(i)
ELSE
#ifdef CMPLX_
write(io1,"(I13,6G25.16)") i, AllHistogram(1,i), sum(norm), &
AllInstHist(1,i), AllInstAnnihil(1,i), &
AllAvAnnihil(1,i), sum(norm1)
#else
write(io1,"(I13,6G25.16)") i, AllHistogram(1,i), norm(1), &
AllInstHist(1,i), AllInstAnnihil(1,i), &
AllAvAnnihil(1,i), norm1(1)
#endif
end if
IF(.not. near_zero(AllHistogram(1,i))) Tot_No_Unique_Dets = Tot_No_Unique_Dets + 1
end do
if(tCalcVariationalEnergy) then
!Calculate the variational FCIMC energy
!Check that Det eq NDet
if(Det.ne.NDet) call stop_all(t_r,'Error')
allocate(CKN(1:Det))
CKN = 0.0_dp
!We need to reorder the vectors so that they are in the same order as the hamiltonian
allocate(HOrderedHist(1:Det))
HOrderedHist(:) = AllHistogram(1,:)
allocate(HOrderedInstHist(1:Det))
HOrderedInstHist(:) = AllInstHist(1,:)
call sort(ReIndex(1:Det),HOrderedHist(1:Det),HOrderedInstHist(1:Det))
#ifdef CMPLX_
call stop_all(t_r, "not implemented for complex")
#else
call my_hpsi(Det,1,NROW,LAB,HAMIL,HOrderedHist,CKN,.true.)
#endif
AvVarEnergy = DDOT(Det,HOrderedHist,1,CKN,1)
CKN = 0.0_dp
#ifdef CMPLX_
call stop_all(t_r, "not implemented for complex")
#else
call my_hpsi(Det,1,NROW,LAB,HAMIL,HOrderedInstHist,CKN,.true.)
#endif
VarEnergy = DDOT(Det,HOrderedInstHist,1,CKN,1)
deallocate(CKN)
deallocate(HOrderedHist,HOrderedInstHist)
if(tDiagAllSpaceEver) then
allocate(ExpandedWalkerDets(NEl,Tot_No_Unique_Dets),stat=ierr)
val=1
do i=1,Det
if(.not. near_zero(AllHistogram(1,i))) then
call decode_bit_det(ExpandedWalkerDets(:,val),FCIDets(0:NIfTot,i))
val=val+1
end if
end do
if((val-1).ne.Tot_No_Unique_Dets) call stop_all(t_r,'Wrong counting')
call LanczosFindGroundE(ExpandedWalkerDets,Tot_No_Unique_Dets,GroundE_Ever,ProjGroundE,.true.)
if(abs(GroundE_Ever-ProjGroundE).gt.1.0e-7_dp) call stop_all(t_r,'Why do these not agree?!')
write(Tot_Unique_Dets_Unit,"(2I14,3G25.15)") Iter,Tot_No_Unique_Dets,AvVarEnergy,VarEnergy,GroundE_Ever
else
write(Tot_Unique_Dets_Unit,"(2I14,2G25.15)") Iter,Tot_No_Unique_Dets,AvVarEnergy,VarEnergy
end if
else
if(tDiagAllSpaceEver) then
allocate(ExpandedWalkerDets(NEl,Tot_No_Unique_Dets),stat=ierr)
val=1
do i=1,Det
if(.not. near_zero(AllHistogram(1,i))) then
call decode_bit_det(ExpandedWalkerDets(:,val),FCIDets(0:NIfTot,i))
val=val+1
end if
end do
if((val-1).ne.Tot_No_Unique_Dets) call stop_all(t_r,'Wrong counting')
call LanczosFindGroundE(ExpandedWalkerDets,Tot_No_Unique_Dets,GroundE_Ever,ProjGroundE,.true.)
if(abs(GroundE_Ever-ProjGroundE).gt.1.0e-7_dp) call stop_all(t_r,'Why do these not agree?!')
write(Tot_Unique_Dets_Unit,"(2I14,3G25.15)") Iter,Tot_No_Unique_Dets,GroundE_Ever
else
write(Tot_Unique_Dets_Unit,"(2I14)") Iter, Tot_No_Unique_Dets
end if
end if
close(io1)
end if
InstHist(:,:)=0.0_dp
InstAnnihil(:,:)=0.0_dp
END SUBROUTINE WriteHistogram