subroutine WriteSingHisttofile()
integer :: i, j, k, BinNo, a, b, iunit
real(dp) :: MaxFII, MinFII, BinIter, BinVal, SingExcit(NoOrbs, NoOrbs)
! <ik|jk> terms where all i, j and k are virtual
! Coulomb.
ROHistSCijkVir(:, :) = 0.0_dp
MaxFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 2, NoOcc + 1)
MinFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 2, NoOcc + 1)
do i = NoOcc + 1, NoOrbs
do k = NoOcc + 1, NoOrbs
do j = i + 1, NoOrbs
if (FourIndInts(i, k, j, k) > MaxFII) MaxFII = FourIndInts(i, k, j, k)
if (FourIndInts(i, k, j, k) < MinFII) MinFII = FourIndInts(i, k, j, k)
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
ROHistSCijkVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = NoOcc + 1, NoOrbs
do k = NoOcc + 1, NoOrbs
do j = i + 1, NoOrbs
if (.not. near_zero(FourIndInts(i, k, j, k))) then
BinNo = CEILING((FourIndInts(i, k, j, k) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSCijkVir(2, BinNo) = ROHistSCijkVir(2, BinNo) + 1.0
end if
end do
end do
end do
! Exchange.
ROHistSEijkVir(:, :) = 0.0_dp
MaxFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 1, NoOcc + 2)
MinFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 1, NoOcc + 2)
do i = NoOcc + 1, NoOrbs
do k = NoOcc + 1, NoOrbs
do j = i + 1, NoOrbs
if (FourIndInts(i, k, k, j) > MaxFII) MaxFII = FourIndInts(i, k, k, j)
if (FourIndInts(i, k, k, j) < MinFII) MinFII = FourIndInts(i, k, k, j)
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
ROHistSEijkVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = NoOcc + 1, NoOrbs
do k = NoOcc + 1, NoOrbs
do j = i + 1, NoOrbs
if (.not. near_zero(FourIndInts(i, k, k, j))) then
BinNo = CEILING((FourIndInts(i, k, k, j) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSEijkVir(2, BinNo) = ROHistSEijkVir(2, BinNo) + 1.0
end if
end do
end do
end do
! Antisymmetric.
ROHistSASijkVir(:, :) = 0.0_dp
MaxFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 2, NoOcc + 1) - FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 1, NoOcc + 2)
MinFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 2, NoOcc + 1) - FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 1, NoOcc + 2)
do i = NoOcc + 1, NoOrbs
do k = NoOcc + 1, NoOrbs
do j = i + 1, NoOrbs
if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) > MaxFII) MaxFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) < MinFII) MinFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
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
ROHistSASijkVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = NoOcc + 1, NoOrbs
do k = NoOcc + 1, NoOrbs
do j = i + 1, NoOrbs
if (.not. (FourIndInts(i, k, j, k) .isclose.FourIndInts(i, k, k, j))) then
BinNo = CEILING(((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSASijkVir(2, BinNo) = ROHistSASijkVir(2, BinNo) + 1.0
end if
end do
end do
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFSingijkVir', status='unknown')
do j = 1, 4002
if (any(.not. near_zero([ROHistSCijkVir(2, j), ROHistSEijkVir(2, j), ROHistSASijkVir(2, j)]))) then
write(iunit, '(6F20.10)') ROHistSCijkVir(1, j), ROHistSCijkVir(2, j), ROHistSEijkVir(1, j), ROHistSEijkVir(2, j),&
&ROHistSASijkVir(1, j), ROHistSASijkVir(2, j)
end if
end do
close(iunit)
end if
if ((Iteration > 1) .and. (.not. tNotConverged)) then
iunit = get_free_unit()
open(iunit, file='HistRotSingijkVir', status='unknown')
do j = 1, 4002
if (any(.not. near_zero([ROHistSCijkVir(2, j), ROHistSEijkVir(2, j), ROHistSASijkVir(2, j)]))) then
write(iunit, '(6F20.10)') ROHistSCijkVir(1, j), ROHistSCijkVir(2, j), ROHistSEijkVir(1, j), ROHistSEijkVir(2, j),&
&ROHistSASijkVir(1, j), ROHistSASijkVir(2, j)
end if
end do
close(iunit)
end if
! <ik|jk> where k is occupied, and i and j are both virtual.
! Coulomb.
ROHistSCkOcijVir(:, :) = 0.0_dp
MaxFII = FourIndInts(NoOcc + 1, 1, NoOcc + 2, 1)
MinFII = FourIndInts(NoOcc + 1, 1, NoOcc + 2, 1)
do i = NoOcc + 1, NoOrbs
do k = 1, NoOcc
do j = i + 1, NoOrbs
if (FourIndInts(i, k, j, k) > MaxFII) MaxFII = FourIndInts(i, k, j, k)
if (FourIndInts(i, k, j, k) < MinFII) MinFII = FourIndInts(i, k, j, k)
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
ROHistSCkOcijVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = NoOcc + 1, NoOrbs
do k = 1, NoOcc
do j = i + 1, NoOrbs
if (.not. near_zero(FourIndInts(i, k, j, k))) then
BinNo = CEILING((FourIndInts(i, k, j, k) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSCkOcijVir(2, BinNo) = ROHistSCkOcijVir(2, BinNo) + 1.0
end if
end do
end do
end do
! Exchange.
ROHistSEkOcijVir(:, :) = 0.0_dp
MaxFII = FourIndInts(NoOcc + 1, 1, 1, NoOcc + 2)
MinFII = FourIndInts(NoOcc + 1, 1, 1, NoOcc + 2)
do i = NoOcc + 1, NoOrbs
do k = 1, NoOcc
do j = i + 1, NoOrbs
if (FourIndInts(i, k, k, j) > MaxFII) MaxFII = FourIndInts(i, k, k, j)
if (FourIndInts(i, k, k, j) < MinFII) MinFII = FourIndInts(i, k, k, j)
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
ROHistSEkOcijVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = NoOcc + 1, NoOrbs
do k = 1, NoOcc
do j = i + 1, NoOrbs
if (.not. near_zero(FourIndInts(i, k, k, j))) then
BinNo = CEILING((FourIndInts(i, k, k, j) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSEkOcijVir(2, BinNo) = ROHistSEkOcijVir(2, BinNo) + 1.0
end if
end do
end do
end do
!antisymmetric
ROHistSASkOcijVir(:, :) = 0.0_dp
MaxFII = FourIndInts(NoOcc + 1, 1, NoOcc + 2, 1) - FourIndInts(NoOcc + 1, 1, 1, NoOcc + 2)
MinFII = FourIndInts(NoOcc + 1, 1, NoOcc + 2, 1) - FourIndInts(NoOcc + 1, 1, 1, NoOcc + 2)
do i = NoOcc + 1, NoOrbs
do k = 1, NoOcc
do j = i + 1, NoOrbs
if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) > MaxFII) MaxFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) < MinFII) MinFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
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
ROHistSASkOcijVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = NoOcc + 1, NoOrbs
do k = 1, NoOcc
do j = i + 1, NoOrbs
if (.not. (FourIndInts(i, k, j, k) .isclose.FourIndInts(i, k, k, j))) then
BinNo = CEILING(((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSASkOcijVir(2, BinNo) = ROHistSASkOcijVir(2, BinNo) + 1.0
end if
end do
end do
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFSingkOcijVir', status='unknown')
do j = 1, 4002
if (any(.not. near_zero([ROHistSCkOcijVir(2, j), ROHistSEkOcijVir(2, j), ROHistSASkOcijVir(2, j)]))) then
write(iunit, '(6F20.10)') ROHistSCkOcijVir(1, j), ROHistSCkOcijVir(2, j), ROHistSEkOcijVir(1, j), &
ROHistSEkOcijVir(2, j), ROHistSASkOcijVir(1, j), ROHistSASkOcijVir(2, j)
end if
end do
close(iunit)
end if
if ((Iteration > 1) .and. (.not. tNotConverged)) then
iunit = get_free_unit()
open(iunit, file='HistRotSingkOcijVir', status='unknown')
do j = 1, 4002
if (any(.not. near_zero([ROHistSCkOcijVir(2, j), ROHistSEkOcijVir(2, j), ROHistSASkOcijVir(2, j)]))) then
write(iunit, '(6F20.10)') ROHistSCkOcijVir(1, j), ROHistSCkOcijVir(2, j), ROHistSEkOcijVir(1, j), &
ROHistSEkOcijVir(2, j), ROHistSASkOcijVir(1, j), ROHistSASkOcijVir(2, j)
end if
end do
close(iunit)
end if
! <ik|jk> where i and k are both occupied, and j virtual.
! Coulomb.
ROHistSCikOcjVir(:, :) = 0.0_dp
MaxFII = FourIndInts(1, 1, NoOcc + 1, 1)
MinFII = FourIndInts(1, 1, NoOcc + 1, 1)
do i = 1, NoOcc
do k = 1, NoOcc
do j = NoOcc + 1, NoOrbs
if (FourIndInts(i, k, j, k) > MaxFII) MaxFII = FourIndInts(i, k, j, k)
if (FourIndInts(i, k, j, k) < MinFII) MinFII = FourIndInts(i, k, j, k)
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
ROHistSCikOcjVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = 1, NoOcc
do k = 1, NoOcc
do j = NoOcc + 1, NoOrbs
if (.not. near_zero(FourIndInts(i, k, j, k))) then
BinNo = CEILING((FourIndInts(i, k, j, k) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSCikOcjVir(2, BinNo) = ROHistSCikOcjVir(2, BinNo) + 1.0
end if
end do
end do
end do
! Exchange.
ROHistSEikOcjVir(:, :) = 0.0_dp
MaxFII = FourIndInts(1, 1, 1, NoOcc + 1)
MinFII = FourIndInts(1, 1, 1, NoOcc + 1)
do i = 1, NoOcc
do k = 1, NoOcc
do j = NoOcc + 1, NoOrbs
if (FourIndInts(i, k, k, j) > MaxFII) MaxFII = FourIndInts(i, k, k, j)
if (FourIndInts(i, k, k, j) < MinFII) MinFII = FourIndInts(i, k, k, j)
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
ROHistSEikOcjVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = 1, NoOcc
do k = 1, NoOcc
do j = NoOcc + 1, NoOrbs
if (.not. near_zero(FourIndInts(i, k, k, j))) then
BinNo = CEILING((FourIndInts(i, k, k, j) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSEikOcjVir(2, BinNo) = ROHistSEikOcjVir(2, BinNo) + 1.0
end if
end do
end do
end do
! Antisymmetrised.
ROHistSASikOcjVir(:, :) = 0.0_dp
MaxFII = FourIndInts(1, 1, NoOcc + 1, 1) - FourIndInts(1, 1, 1, NoOcc + 1)
MinFII = FourIndInts(1, 1, NoOcc + 1, 1) - FourIndInts(1, 1, 1, NoOcc + 1)
do i = 1, NoOcc
do k = 1, NoOcc
do j = NoOcc + 1, NoOrbs
if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) > MaxFII) MaxFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) < MinFII) MinFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
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
ROHistSASikOcjVir(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do i = 1, NoOcc
do k = 1, NoOcc
do j = NoOcc + 1, NoOrbs
if (.not. (FourIndInts(i, k, j, k) .isclose.FourIndInts(i, k, k, j))) then
BinNo = CEILING(((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSASikOcjVir(2, BinNo) = ROHistSASikOcjVir(2, BinNo) + 1.0
end if
end do
end do
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFSingikOcjVir', status='unknown')
do j = 1, 4002
if (any(.not. near_zero([ROHistSCikOcjVir(2, j), ROHistSEikOcjVir(2, j), ROHistSASikOcjVir(2, j)]))) then
write(iunit, '(6F20.10)') ROHistSCikOcjVir(1, j), ROHistSCikOcjVir(2, j), ROHistSEikOcjVir(1, j), &
ROHistSEikOcjVir(2, j), ROHistSASikOcjVir(1, j), ROHistSASikOcjVir(2, j)
end if
end do
close(iunit)
end if
if ((Iteration > 1) .and. (.not. tNotConverged)) then
iunit = get_free_unit()
open(iunit, file='HistRotSingikOcjVir', status='unknown')
do j = 1, 4002
if (any(.not. near_zero([ROHistSCikOcjVir(2, j), ROHistSEikOcjVir(2, j), ROHistSASikOcjVir(2, j)]))) then
write(iunit, '(6F20.10)') ROHistSCikOcjVir(1, j), ROHistSCikOcjVir(2, j), ROHistSEikOcjVir(1, j), &
ROHistSEikOcjVir(2, j), ROHistSASikOcjVir(1, j), ROHistSASikOcjVir(2, j)
end if
end do
close(iunit)
end if
! Single excitations connected to the HF determinant.
ROHistSing(:, :) = 0.0_dp
MaxFII = 0.0_dp
MinFII = 0.0_dp
do j = NoOcc + 1, NoOrbs
do i = 1, NoOcc
SingExcit(i, j) = 0.0_dp
if (i == j) cycle
a = SymLabelList2_rot(i)
b = SymLabelList2_rot(j)
do k = 1, NoOcc + 1
SingExcit(i, j) = SingExcit(i, j) + real(TMAT2D(2 * a, 2 * b), dp) + ((2 * FourIndInts(i, k, j, k)) - FourIndInts(i, k, k, j))
end do
if (SingExcit(i, j) > MaxFII) MaxFII = SingExcit(i, j)
if (SingExcit(i, j) < MinFII) MinFII = SingExcit(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
ROHistSing(1, i) = BinVal
BinVal = BinVal + BinIter
end do
do j = NoOcc + 1, NoOrbs
do i = 1, NoOcc
if (i == j) cycle
if (.not. near_zero(SingExcit(i, j))) then
BinNo = CEILING((SingExcit(i, j) - MinFII) * 4002 / (MaxFII - MinFII))
ROHistSing(2, BinNo) = ROHistSing(2, BinNo) + 1.0
end if
end do
end do
if (Iteration == 0) then
iunit = get_free_unit()
open(iunit, file='HistHFSingExcHF', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistSing(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistSing(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='HistRotSingExcHF', status='unknown')
do j = 1, 4002
if (.not. near_zero(ROHistSing(2, j))) then
do i = 1, 2
write(iunit, '(F20.10)', advance='no') ROHistSing(i, j)
end do
write(iunit, *) ''
end if
end do
close(iunit)
end if
end subroutine WriteSingHisttofile