subroutine WriteDoubHisttofile()
integer :: i, j, k, l, BinNo, iunit
real(dp) :: MaxFII, MinFII, BinIter, OnePartOrbEnValue, BinVal
! Histogramming all coulomb terms <ij|ij> where i<j, and i and j are
! both virtual. In reality we are looking at i = <j, but the ERhistograms
! will show the i = j terms.
if (tROHistVirtCoulomb) then
ROHistDCijOcklVir(:, :) = 0.0_dp
MinFII = FourIndInts(1, 2, NoOcc + 1, NoOcc + 2)
MaxFII = FourIndInts(1, 2, NoOcc + 1, NoOcc + 2)
do i = 1, NoOcc
do j = 1, NoOcc
do k = NoOcc + 1, NoOrbs
do l = NoOcc + 1, NoOrbs
if (FourIndInts(i, j, k, l) < MinFII) MinFII = FourIndInts(i, j, k, l)
if (FourIndInts(i, j, k, l) > MaxFII) MaxFII = FourIndInts(i, j, k, l)
end do
end do
end do
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistDCijOcklVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = 1, NoOcc
do j = 1, NoOcc
do k = NoOcc + 1, NoOrbs
do l = NoOcc + 1, NoOrbs
if (.not. near_zero(FourIndInts(i, j, k, l))) then
BinNo = CEILING((FourIndInts(i, j, k, l) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistDCijOcklVir(2, BinNo) = ROHistDCijOcklVir(2, BinNo) + 1.0
end if
end do
end do
end do
end do
! Antisymmetric.
ROHistASijOcklVir(:, :) = 0.0_dp
MinFII = FourIndInts(1, 2, NoOcc + 1, NoOcc + 2) - FourIndInts(1, 2, NoOcc + 2, NoOcc + 1)
MaxFII = FourIndInts(1, 2, NoOcc + 1, NoOcc + 2) - FourIndInts(1, 2, NoOcc + 2, NoOcc + 1)
do i = 1, NoOcc
do j = 1, NoOcc
do k = NoOcc + 1, NoOrbs
do l = NoOcc + 1, NoOrbs
if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) < MinFII) MinFII = (FourIndInts(i, j, k, l) - &
FourIndInts(i, j, l, k))
if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) > MaxFII) MaxFII = (FourIndInts(i, j, k, l) - &
FourIndInts(i, j, l, k))
end do
end do
end do
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistASijOcklVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = 1, NoOcc
do j = 1, NoOcc
do k = NoOcc + 1, NoOrbs
do l = NoOcc + 1, NoOrbs
if (.not. (FourIndInts(i, j, k, l) .isclose.FourIndInts(i, j, l, k))) then
BinNo = CEILING(((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistASijOcklVir(2, BinNo) = ROHistASijOcklVir(2, BinNo) + 1.0
end if
end do
end do
end do
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFDoubijOcklVir', status='unknown')
do j = 1, 4002
if (.not. (near_zero(ROHistDCijOcklVir(2, j)) .and. near_zero(ROHistASijOcklVir(2, j)))) then
write(iunit, '(4F20.10)') ROHistDCijOcklVir(1, j), ROHistDCijOcklVir(2, j), &
ROHistASijOcklVir(1, j), ROHistASijOcklVir(2, j)
end if
end do
close(iunit)
end if
if ((.not. tNotConverged) .and. (Iteration > 1)) then
iunit = get_free_unit()
open(iunit, file='HistRotDoubijOcklVir', status='unknown')
do j = 1, 4002
if (.not. (near_zero(ROHistDCijOcklVir(2, j)) .and. near_zero(ROHistASijOcklVir(2, j)))) then
write(iunit, '(4F20.10)') ROHistDCijOcklVir(1, j), ROHistDCijOcklVir(2, j), &
ROHistASijOcklVir(1, j), ROHistASijOcklVir(2, j)
end if
end do
close(iunit)
end if
ROHistDCijklVir(:, :) = 0.0_dp
MinFII = FourIndInts(NoOrbs - 1, NoOrbs, NoOrbs - 1, NoOrbs)
MaxFII = FourIndInts(NoOrbs - 1, NoOrbs, NoOrbs - 1, NoOrbs)
do i = NoOcc + 1, NoOrbs
do j = NoOcc + 1, NoOrbs
do k = i + 1, NoOrbs
do l = j + 1, NoOrbs
if (FourIndInts(i, j, k, l) < MinFII) MinFII = FourIndInts(i, j, k, l)
if (FourIndInts(i, j, k, l) > MaxFII) MaxFII = FourIndInts(i, j, k, l)
end do
end do
end do
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistDCijklVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = NoOcc + 1, NoOrbs
do j = NoOcc + 1, NoOrbs
do k = i + 1, NoOrbs
do l = j + 1, NoOrbs
if (.not. near_zero(FourIndInts(i, j, k, l))) then
BinNo = CEILING((FourIndInts(i, j, k, l) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistDCijklVir(2, BinNo) = ROHistDCijklVir(2, BinNo) + 1.0
end if
end do
end do
end do
end do
! Antisymmetric.
ROHistASijklVir(:, :) = 0.0_dp
MinFII = FourIndInts(NoOrbs - 3, NoOrbs - 2, NoOrbs - 1, NoOrbs) - FourIndInts(NoOrbs - 3, NoOrbs - 2, NoOrbs, NoOrbs - 1)
MaxFII = FourIndInts(NoOrbs - 3, NoOrbs - 2, NoOrbs - 1, NoOrbs) - FourIndInts(NoOrbs - 3, NoOrbs - 2, NoOrbs, NoOrbs - 1)
do i = NoOcc + 1, NoOrbs
do j = NoOcc + 1, NoOrbs
do k = i + 1, NoOrbs
do l = j + 1, NoOrbs
if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) < MinFII) MinFII = (FourIndInts(i, j, k, l) - &
FourIndInts(i, j, l, k))
if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) > MaxFII) MaxFII = (FourIndInts(i, j, k, l) - &
FourIndInts(i, j, l, k))
end do
end do
end do
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistASijklVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = NoOcc + 1, NoOrbs
do j = NoOcc + 1, NoOrbs
do k = i + 1, NoOrbs
do l = j + 1, NoOrbs
if (.not. (FourIndInts(i, j, k, l) .isclose.FourIndInts(i, j, l, k))) then
BinNo = CEILING(((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistASijklVir(2, BinNo) = ROHistASijklVir(2, BinNo) + 1.0
end if
end do
end do
end do
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFDoubijklVirt', status='unknown')
do j = 1, 4002
if (.not. (near_zero(ROHistDCijklVir(2, j)) .and. near_zero(ROHistASijklVir(2, j)))) then
write(iunit, '(4F20.10)') ROHistDCijklVir(1, j), ROHistDCijklVir(2, j), &
ROHistASijklVir(1, j), ROHistASijklVir(2, j)
end if
end do
close(iunit)
end if
if ((.not. tNotConverged) .and. (Iteration > 1)) then
iunit = get_free_unit()
open(iunit, file='HistRotDoubijklVirt', status='unknown')
do j = 1, 4002
if (.not. (near_zero(ROHistDCijklVir(2, j)) .and. near_zero(ROHistASijklVir(2, j)))) then
write(iunit, '(4F20.10)') ROHistDCijklVir(1, j), ROHistDCijklVir(2, j), ROHistASijklVir(1, j), &
ROHistASijklVir(2, j)
end if
end do
close(iunit)
end if
end if
! Histogramming all one particle orbital energies (occupied and
! virtual) even though we are not changing occupied. Would like to
! see HOMO-LUMO gap etc.
if (tROHistOneElInts) then
ROHistHijVirt(:, :) = 0.0_dp
MinFII = TMAT2DRot(NoOcc + 1, NoOcc + 2)
MaxFII = TMAT2DRot(NoOcc + 1, NoOcc + 2)
do i = NoOcc + 1, NoOrbs
do j = i + 1, NoOrbs
if (TMAT2DRot(i, j) < MinFII) MinFII = TMAT2DRot(i, j)
if (TMAT2DRot(i, j) > MaxFII) MaxFII = TMAT2DRot(i, j)
end do
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistHijVirt(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = NoOcc + 1, NoOrbs
do j = i + 1, NoOrbs
if (.not. near_zero(TMAT2DRot(i, j))) then
BinNo = CEILING((TMAT2DRot(i, j) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistHijVirt(2, BinNo) = ROHistHijVirt(2, BinNo) + 1.0
end if
end do
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFHijVirt', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistHijVirt(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistHijVirt(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
if ((.not. tNotConverged) .and. (Iteration > 1)) then
iunit = get_free_unit()
open(iunit, file='HistRotHijVirt', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistHijVirt(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistHijVirt(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
ROHistHijOccVirt(:, :) = 0.0_dp
MinFII = TMAT2DRot(1, NoOcc + 1)
MaxFII = TMAT2DRot(1, NoOcc + 1)
do i = 1, NoOcc
do j = NoOcc + 1, NoOrbs
if (TMAT2DRot(i, j) < MinFII) MinFII = TMAT2DRot(i, j)
if (TMAT2DRot(i, j) > MaxFII) MaxFII = TMAT2DRot(i, j)
end do
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistHijOccVirt(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = 1, NoOcc
do j = NoOcc + 1, NoOrbs
if (.not. near_zero(TMAT2DRot(i, j))) then
BinNo = CEILING((TMAT2DRot(i, j) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistHijOccVirt(2, BinNo) = ROHistHijOccVirt(2, BinNo) + 1.0
end if
end do
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFHijOccVirt', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistHijOccVirt(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistHijOccVirt(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
if ((.not. tNotConverged) .and. (Iteration > 1)) then
iunit = get_free_unit()
open(iunit, file='HistRotHijOccVirt', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistHijOccVirt(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistHijOccVirt(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
ROHistHii(:, :) = 0.0_dp
MinFII = TMAT2DRot(1, 1)
MaxFII = TMAT2DRot(1, 1)
do i = 1, NoOrbs
if (TMAT2DRot(i, i) < MinFII) MinFII = TMAT2DRot(i, i)
if (TMAT2DRot(i, i) > MaxFII) MaxFII = TMAT2DRot(i, i)
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistHii(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = 1, NoOrbs
BinNo = CEILING((TMAT2DRot(i, i) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistHii(2, BinNo) = ROHistHii(2, BinNo) + 1.0
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFHii', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistHii(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistHii(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
if ((.not. tNotConverged) .and. (Iteration > 1)) then
iunit = get_free_unit()
open(iunit, file='HistRotHii', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistHii(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistHii(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
end if
if (tROHistOnePartOrbEn) then
ROHistOnePartOrbEn(:, :) = 0.0_dp
MaxFII = 0.0_dp
MinFII = 0.0_dp
do i = 1, NoOrbs
OnePartOrbEnValue = 0.0_dp
OnePartOrbEnValue = OnePartOrbEnValue + TMAT2DRot(i, i)
do j = 1, NoOcc
OnePartOrbEnValue = OnePartOrbEnValue + (2 * FourIndInts(i, j, i, j)) - FourIndInts(i, j, j, i)
end do
if (OnePartOrbEnValue > MaxFII) MaxFII = OnePartOrbEnValue
if (OnePartOrbEnValue < MinFII) MinFII = OnePartOrbEnValue
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistOnePartOrbEn(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = 1, NoOrbs
OnePartOrbEnValue = 0.0_dp
OnePartOrbEnValue = OnePartOrbEnValue + TMAT2DRot(i, i)
do j = 1, NoOcc
OnePartOrbEnValue = OnePartOrbEnValue + (2 * FourIndInts(i, j, i, j)) - FourIndInts(i, j, j, i)
end do
BinNo = CEILING((OnePartOrbEnValue - MinFII) * 4002 / (MaxFII - MinFII))
ROHistOnePartOrbEn(2, BinNo) = ROHistOnePartOrbEn(2, BinNo) + 1.0
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFOnePartOrbEn', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistOnePartOrbEn(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistOnePartOrbEn(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
if ((Iteration > 1) .and. (.not. tNotConverged)) then
iunit = get_free_unit()
open(iunit, file='HistRotOnePartOrbEn', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistOnePartOrbEn(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistOnePartOrbEn(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
end if
if (tROHistDoubExc) then
ROHistDoubExc(:, :) = 0.0_dp
MaxFII = 0.0_dp
MinFII = 0.0_dp
do l = NoOcc + 1, NoOrbs
do j = 1, NoOcc
do k = NoOcc + 1, NoOrbs
do i = 1, NoOcc
if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) > MaxFII) then
MaxFII = (FourIndInts(i, j, k, l)) - FourIndInts(i, j, l, k)
end if
if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) < MinFII) then
MinFII = (FourIndInts(i, j, k, l)) - FourIndInts(i, j, l, k)
end if
end do
end do
end do
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistDoubExc(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do l = NoOcc + 1, NoOrbs
do j = 1, NoOcc
do k = NoOcc + 1, NoOrbs
do i = 1, NoOcc
if (.not. (FourIndInts(i, j, k, l) .isclose.FourIndInts(i, j, l, k))) then
BinNo = CEILING(((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistDoubExc(2, BinNo) = ROHistDoubExc(2, BinNo) + 1.0
end if
end do
end do
end do
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFDoubExc', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistDoubExc(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistDoubExc(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
if ((Iteration > 1) .and. (.not. tNotConverged)) then
iunit = get_free_unit()
open(iunit, file='HistRotDoubExc', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistDoubExc(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistDoubExc(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
end if
if (tROHistER) then
ROHistER(:, :) = 0.0_dp
MaxFII = 0.0_dp
MinFII = 0.0_dp
do i = 1, NoOrbs
if (FourIndInts(i, i, i, i) > MaxFII) MaxFII = FourIndInts(i, i, i, i)
if (FourIndInts(i, i, i, i) < MinFII) MinFII = FourIndInts(i, i, i, i)
end do
BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
MaxFII = MaxFII + BinIter
MinFII = MinFII - BinIter
BinVal = MinFII
do i = 1, 4002
ROHistER(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = 1, NoOrbs
BinNo = CEILING((FourIndInts(i, i, i, i) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistER(2, BinNo) = ROHistER(2, BinNo) + 1.0
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHF-ER', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistER(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistER(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
if ((Iteration > 1) .and. (.not. tNotConverged)) then
iunit = get_free_unit()
open(iunit, file='HistRot-ER', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistER(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistER(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
end if
end subroutine WriteDoubHisttofile