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