C.. At this point, we transform the UMAT and TMAT into a new UMAT and
C.. TMAT and Ecore with the frozen orbitals factored in
C..
C.. a,b are frozen spinorbitals
C.. E’core = Ecore+sum_a t_aa + sum_(a<b) (
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer | :: | NHG | ||||
| integer | :: | NBASIS | ||||
| real(kind=dp) | :: | UMAT(*) |
C.. was (NHG/ISS,NHG/ISS,NHG/ISS,NHG/ISS) |
|||
| real(kind=dp) | :: | UMAT2(*) | ||||
| real(kind=dp) | :: | ECORE |
C.. was (NBASIS/ISS,NBASIS/ISS,NBASIS/ISS,NBASIS/ISS) |
|||
| type(BasisFN) | :: | G1(NHG) | ||||
| integer | :: | nBasisMax(5,*) | ||||
| integer | :: | ISS | ||||
| integer | :: | BRR(NHG) | ||||
| integer | :: | NFROZEN | ||||
| integer | :: | NTFROZEN | ||||
| integer | :: | NFROZENIN | ||||
| integer | :: | NTFROZENIN | ||||
| integer | :: | NEL |
SUBROUTINE IntFREEZEBASIS(NHG, NBASIS, UMAT, UMAT2, ECORE, &
& G1, NBASISMAX, ISS, BRR, NFROZEN, NTFROZEN, NFROZENIN, NTFROZENIN, NEL)
use SystemData, only: Symmetry, BasisFN, arr, tagarr
use OneEInts, only: GetPropIntEl, GetTMATEl, TMATSYM2, TMAT2D2, PropCore, &
OneEPropInts2, OneEPropInts, tOneElecDiag, NewTMatInd, &
GetNEWTMATEl, tCPMDSymTMat, SetupTMAT2, SWAPTMAT, &
SwapOneEPropInts, SetupPropInts2, FieldCore, OneEFieldInts, &
OneEFieldInts2, SetupFieldInts2, SwapOneEFieldInts
USE UMatCache, only: FreezeTransfer, UMatCacheData
Use UMatCache, only: FreezeUMatCache, CreateInvBrr2, FreezeUMat2D, SetupUMatTransTable
use LoggingData, only: tCalcPropEst, iNumPropToEst
use UMatCache, only: GTID
use global_utilities
use CalcData, only: nFields_it, tCalcWithField
use sym_mod, only: getsym, SetupFREEZEALLSYM, FREEZESYMLABELS
use util_mod, only: NECI_ICOPY
INTEGER NHG, NBASIS, nBasisMax(5, *), ISS
TYPE(BASISFN) G1(NHG), KSYM
HElement_t(dp) UMAT(*)
!!C.. was (NHG/ISS,NHG/ISS,NHG/ISS,NHG/ISS)
HElement_t(dp) UMAT2(*)
real(dp) ECORE
!!C.. was (NBASIS/ISS,NBASIS/ISS,NBASIS/ISS,NBASIS/ISS)
real(dp) ARR2(NBASIS, 2)
INTEGER NFROZEN, BRR(NHG), BRR2(NBASIS), GG(NHG)
TYPE(BASISFN) G2(NHG)
INTEGER NTFROZEN, NFROZENIN, NTFROZENIN, frozen_count, list_sz
INTEGER BLOCKMINW, BLOCKMAXW, FROZENBELOWW, BLOCKMINY, BLOCKMAXY, FROZENBELOWY
INTEGER BLOCKMINX, BLOCKMAXX, FROZENBELOWX, BLOCKMINZ, BLOCKMAXZ, FROZENBELOWZ
INTEGER I, J, K, L, A, B, IP, JP, IDI, IDJ, IDK, IDL, W, X, Y, Z
INTEGER IB, JB, KB, LB, AB, BB, IPB, JPB, KPB, LPB
INTEGER IDIP, IDJP, IDKP, IDLP
INTEGER IDA, IDB
INTEGER iSize, ierr
INTEGER NEL
logical :: frozen_occ, frozen_virt
! TYPE(Symmetry) KSYM
character(*), parameter :: this_routine = 'IntFreezeBasis'
! IF(tHub.or.tUEG) THEN
! CALL Stop_All("IntFreezeBasis","Freezing does not currently work with the hubbard model/UEG.")
! end if
!!C.. Just check to see if we're not in the middle of a degenerate set with the same sym
IF(NFROZEN > 0) THEN
IF(ABS(ARR(NFROZEN, 1) - ARR(NFROZEN + 1, 1)) < 1.0e-6_dp .AND. &
& G1(BRR(NFROZEN))%SYM%s == G1(BRR(NFROZEN + 1))%SYM%s) THEN
call stop_all(this_routine, "Cannot freeze in the middle of a degenerate set")
ELSE IF(ABS(ARR(NFROZEN, 1) - ARR(NFROZEN + 1, 1)) < 1.0e-6_dp) THEN
write(stdout, '(a)') 'WARNING: Freezing in the middle of a degenerate set.'
write(stdout, '(a)') 'This should only be done for debugging purposes.'
end if
end if
IF(NTFROZEN > 0) THEN
IF(ABS(ARR(NHG - NTFROZEN, 1) - ARR(NHG - NTFROZEN + 1, 1)) < 1.0e-6_dp &
& .AND. G1(BRR(NHG - NTFROZEN))%SYM%s &
& == G1(BRR(NHG - NTFROZEN + 1))%SYM%s) THEN
call stop_all(this_routine, "Cannot freeze in the middle of a degenerate virtual set")
ELSE IF(ABS(ARR(NHG - NTFROZEN, 1) - ARR(NHG - NTFROZEN + 1, 1)) < 1.0e-6_dp) THEN
write(stdout, '(a)') 'WARNING: Freezing in the middle of a degenerate set.'
write(stdout, '(a)') 'This should only be done for debugging purposes.'
end if
end if
IF(NFROZENIN > 0) THEN
IF(ABS(ARR(NEL - NFROZENIN, 1) - ARR(NEL - NFROZENIN + 1, 1)) < 1.0e-6_dp .AND. &
& G1(BRR(NEL - NFROZENIN))%SYM%s == G1(BRR(NEL - NFROZENIN + 1))%SYM%s) THEN
call stop_all(this_routine, "Cannot freeze in the middle of a degenerate set")
ELSE IF(ABS(ARR(NEL - NFROZENIN, 1) - ARR(NEL - NFROZENIN + 1, 1)) < 1.0e-6_dp) THEN
write(stdout, '(a)') 'WARNING: Freezing in the middle of a degenerate set.'
write(stdout, '(a)') 'This should only be done for debugging purposes.'
end if
end if
IF(NTFROZENIN > 0) THEN
IF(ABS(ARR(NEL + NTFROZENIN, 1) - ARR(NEL + NTFROZENIN + 1, 1)) < 1.0e-6_dp &
& .AND. G1(BRR(NEL + NTFROZENIN))%SYM%s &
& == G1(BRR(NEL + NTFROZENIN + 1))%SYM%s) THEN
call stop_all(this_routine, "Cannot freeze in the middle of a degenerate virtual set")
ELSE IF(ABS(ARR(NEL + NTFROZENIN, 1) - ARR(NEL + NTFROZENIN + 1, 1)) < 1.0e-6_dp) THEN
write(stdout, '(a)') 'WARNING: Freezing in the middle of a degenerate set.'
write(stdout, '(a)') 'This should only be done for debugging purposes.'
end if
end if
!!C.. At this point, we transform the UMAT and TMAT into a new UMAT and
!!C.. TMAT and Ecore with the frozen orbitals factored in
!!C..
!!C.. a,b are frozen spinorbitals
!!C.. E'core = Ecore+sum_a t_aa + sum_(a<b) (<ab|ab>-<ab|ba>)
!!C.. t'_ii = t_ii+ sum_a ( <ai|ai> - <ai|ia> )
!!C.. NHG contains the old number of orbitals
!!C.. NBASIS contains the new
!!C.. We first need to work out where each of the current orbitals will
!!C.. end up in the new set
! Allocate the frozen orbital list for later use
! n.b. These are for frozen occupied orbs. Frozen virts are chucked.
list_sz = nfrozen + nfrozenin !ntfrozen + ntfrozenin
if(list_sz /= 0 .or. ntfrozen + ntfrozenin /= 0) then
allocate(frozen_orb_list(list_sz), frozen_orb_reverse_map(nbasis), &
stat=ierr)
call LogMemAlloc("frozen_orb_list", list_sz, sizeof_int, &
this_routine, tagFrozen)
call LogMemAlloc("frozen_orb_reverse_map", nbasis, sizeof_int, &
this_routine, tagFrozenMap)
end if
K = 0
frozen_count = 0
DO I = 1, NHG
frozen_occ = .false.
frozen_virt = .false.
DO J = 1, NFROZEN
! if (any(brr(1:nfrozen) == i))
!C.. if orb I is to be frozen, L will become 0
IF(BRR(J) == I) frozen_occ = .true.
end do
DO J = NEL - NFROZENIN + 1, NEL
IF(BRR(J) == I) frozen_occ = .true.
end do
DO J = NEL + 1, NEL + NTFROZENIN
IF(BRR(J) == I) frozen_virt = .true.
end do
DO J = NBASIS + NFROZEN + NFROZENIN + NTFROZENIN + 1, NHG
!C.. if orb I is to be frozen, L will become 0
IF(BRR(J) == I) frozen_virt = .true.
end do
if(frozen_occ) then
GG(I) = 0
frozen_count = frozen_count + 1
frozen_orb_list(frozen_count) = i
else if(frozen_virt) then
GG(I) = 0
ELSE
!C.. we've got an orb which is not to be frozen
K = k + 1
!C.. GG(I) is the new position in G of the (old) orb I
GG(I) = K
!C.. copy the eigenvalue table to the new location
G2(K) = G1(I)
! Store the reverse mapping as well
frozen_orb_reverse_map(k) = i
end if
end do
if(frozen_count /= list_sz .or. k /= nbasis) &
call stop_all(this_routine, 'Incorrect number of orbitals frozen')
!C.. Now construct the new BRR and ARR
! DO I=1,NBASIS
!Need to run through the remaining orbitals in 2 lots, the occupied and virtual, because
!each are being shifted by different amounts. The occupied are only affected by the low energy
!frozen orbitals, but the virtuals need to also account for the inner frozen orbitals.
DO W = 1, 2
IF(W == 1) THEN
BLOCKMINW = 1
BLOCKMAXW = NEL - NFROZEN - NFROZENIN
FROZENBELOWW = NFROZEN
else if(W == 2) THEN
BLOCKMINW = NEL - NFROZEN - NFROZENIN + 1
BLOCKMAXW = NBASIS
FROZENBELOWW = NFROZEN + NFROZENIN + NTFROZENIN
end if
DO I = BLOCKMINW, BLOCKMAXW
BRR2(I) = GG(BRR(I + FROZENBELOWW))
ARR2(I, 1) = ARR(I + FROZENBELOWW, 1)
end do
end do
DO I = 1, NHG
IF(GG(I) /= 0) ARR2(GG(I), 2) = ARR(I, 2)
end do
CALL SetupTMAT2(NBASIS, 2, iSize)
if(tCalcPropEst) call SetupPropInts2(NBASIS)
if (tCalcWithField) call SetupFieldInts2(NBASIS, nFields_it)
!C.. First deal with Ecore
!Adding the energy of the occupied orbitals to the core energy.
!Need to do this for both the low energy frozen and inner frozen orbitals.
DO A = 1, NFROZEN
AB = BRR(A)
! Ecore' = Ecore + sum_a <a|h|a> where a is a frozen spin orbital
! TMATEL is the one electron integrals <a|h|a>.
ECORE = ECORE + GetTMATEl(AB, AB)
! Ecore' = Ecore + sum a<b (<ab|ab> - <ab|ba>)
DO B = A + 1, NFROZEN
BB = BRR(B)
IDA = GTID(AB)
IDB = GTID(BB)
!C.. No sign problems from permuations here as all perms even
ECORE = ECORE + get_umat_el(IDA, IDB, IDA, IDB)
!C.. If we have spin-independent integrals, or
!C.. if the spins are the same
IF(G1(AB)%MS == G1(BB)%MS) &
& ECORE = ECORE - get_umat_el(IDA, IDB, IDB, IDA)
end do
!The sum over b runs over all frozen orbitals > a, so the inner frozen orbitals too.
DO B = NEL - NFROZENIN + 1, NEL
BB = BRR(B)
IDA = GTID(AB)
IDB = GTID(BB)
!C.. No sign problems from permuations here as all perms even
ECORE = ECORE + get_umat_el(IDA, IDB, IDA, IDB)
!C.. If we have spin-independent integrals, or
!C.. if the spins are the same
IF(G1(AB)%MS == G1(BB)%MS) &
& ECORE = ECORE - get_umat_el(IDA, IDB, IDB, IDA)
end do
end do
!Need to also account for when a is the frozen inner orbitals, but b > a, so b only runs over the frozen
!inner.
DO A = NEL - NFROZENIN + 1, NEL
AB = BRR(A)
ECORE = ECORE + GetTMATEl(AB, AB)
DO B = A + 1, NEL
BB = BRR(B)
IDA = GTID(AB)
IDB = GTID(BB)
!C.. No sign problems from permuations here as all perms even
ECORE = ECORE + get_umat_el(IDA, IDB, IDA, IDB)
!C.. If we have spin-independent integrals, or
!C.. if the spins are the same
IF(G1(AB)%MS == G1(BB)%MS) &
& ECORE = ECORE - get_umat_el(IDA, IDB, IDB, IDA)
end do
end do
! Now dealing with the zero body part of the property integrals if needed
IF(tCalcPropEst) then
write(stdout, *) 'PropCore before freezing:', PropCore
DO A = 1, NFROZEN
AB = BRR(A)
! Ecore' = Ecore + sum_a <a|h|a> where a is a frozen spin orbital
! TMATEL is the one electron integrals <a|h|a>.
DO B = 1, iNumPropToEst
PropCore(B) = PropCore(B) + GetPropIntEl(AB, AB, B)
write(stdout, *) '1', PropCore(B), AB, B, GetPropIntEl(AB, AB, B)
end do
end do
! Now dealing with the zero body part of the field integrals if needed
write(stdout,*) 'FieldCore before freezing:', FieldCore
IF(tCalcWithField) then
DO A=1,NFROZEN
AB=BRR(A)
! Ecore' = Ecore + sum_a <a|h|a> where a is a frozen spin orbital
! TMATEL is the one electron integrals <a|h|a>.
DO B=1, nFields_it
FieldCore(B)=FieldCore(B)+OneEFieldInts(AB,AB,B)
write(stdout,*) '1', FieldCore(B), AB, B, OneEFieldInts(AB,AB,B)
ENDDO
ENDDO
!Need to also account for when a is the frozen inner orbitals
DO A=NEL-NFROZENIN+1,NEL
AB=BRR(A)
DO B=1, nFields_it
FieldCore(B)=FieldCore(B)+OneEFieldInts(AB,AB,B)
write(stdout,*) '2', FieldCore(B), AB, B, OneEFieldInts(AB,AB,B)
ENDDO
ENDDO
ENDIF
write(stdout,*) 'FieldCore after freezing:', FieldCore
!Need to also account for when a is the frozen inner orbitals
DO A = NEL - NFROZENIN + 1, NEL
AB = BRR(A)
DO B = 1, iNumPropToEst
PropCore(B) = PropCore(B) + GetPropIntEl(AB, AB, B)
write(stdout, *) '2', PropCore(B), AB, B, GetPropIntEl(AB, AB, B)
end do
end do
write(stdout, *) 'PropCore after freezing:', PropCore
end if
!C.. now deal with the new TMAT
FREEZETRANSFER = .true.
!First the low energy frozen orbitals.
!t'_ii = t_ii+ sum_a ( <ai|ai> - <ai|ia> )
!Again need to do this for the remaining occupied, and then the remaining virtual separately.
!The above i runs over all orbitals, whereas a is only over the occupied virtuals.
DO W = 1, 2
IF(W == 1) THEN
BLOCKMINW = 1
BLOCKMAXW = NEL - NFROZEN - NFROZENIN
FROZENBELOWW = NFROZEN
else if(W == 2) THEN
BLOCKMINW = NEL - NFROZEN - NFROZENIN + 1
BLOCKMAXW = NBASIS
FROZENBELOWW = NFROZEN + NFROZENIN + NTFROZENIN
end if
DO I = BLOCKMINW, BLOCKMAXW
IP = I + FROZENBELOWW
IB = BRR(IP)
IPB = GG(IB)
IDI = GTID(IB)
!I and J give the indexes of the TMAT. This bit accounts for the off-diagonal terms which must be copied accross.
DO Y = 1, 2
IF(Y == 1) THEN
BLOCKMINY = 1
BLOCKMAXY = NEL - NFROZEN - NFROZENIN
FROZENBELOWY = NFROZEN
else if(Y == 2) THEN
BLOCKMINY = NEL - NFROZEN - NFROZENIN + 1
BLOCKMAXY = NBASIS
FROZENBELOWY = NFROZEN + NFROZENIN + NTFROZENIN
end if
DO J = BLOCKMINY, BLOCKMAXY
JP = J + FROZENBELOWY
JB = BRR(JP)
JPB = GG(JB)
IDJ = GTID(JB)
IF(tCPMDSymTMat) THEN
TMATSYM2(NEWTMATInd(IPB, JPB)) = GetTMATEl(IB, JB)
ELSE
IF(IPB == 0 .or. JPB == 0) THEN
! write(stdout,*) 'W',W,'I',I,'J',J,'IPB',IPB,'JPB',JPB
! CALL neci_flush(stdout)
! CALL Stop_All("","here 01")
end if
if(tOneElecDiag) then
if((IB == JB) .and. (IPB == JPB)) then
TMAT2D2(IPB, 1) = GetTMATEl(IB, JB)
end if
else
TMAT2D2(IPB, JPB) = GetTMATEl(IB, JB)
end if
end if
DO A = 1, NFROZEN
AB = BRR(A)
IDA = GTID(AB)
!C.. SGN takes into account permutationnness.
!C SGN=1
!C IF(IB.GT.AB) SGN=-SGN
!C IF(JB.GT.AB) SGN=-SGN
IF(G1(IB)%MS == G1(JB)%MS) THEN
IF(tCPMDSymTMat) THEN
TMATSYM2(NEWTMATInd(IPB, JPB)) = &
& GetNEWTMATEl(IPB, JPB) + get_umat_el(IDA, IDI, IDA, IDJ)
ELSE
! IF(IPB.eq.0.or.JPB.eq.0) CALL Stop_All("","here 02")
if(tOneElecDiag) then
if(IPB == JPB) then
TMAT2D2(IPB, 1) = TMAT2D2(IPB, 1) + get_umat_el(IDA, IDI, IDA, IDJ)
else
if(abs(get_umat_el(IDA, IDI, IDA, IDJ)) > 1.0e-8_dp) then
call stop_all("", "Error here in freezing for UEG")
end if
end if
else
TMAT2D2(IPB, JPB) = TMAT2D2(IPB, JPB) + get_umat_el(IDA, IDI, IDA, IDJ)
end if
end if
end if
!C.. If we have spin-independent integrals, ISS.EQ.2.OR
!C.. if the spins are the same
IF(G1(IB)%MS == G1(AB)%MS .AND. G1(AB)%MS == G1(JB)%MS) THEN
IF(tCPMDSymTMat) THEN
TMATSYM2(NEWTMATInd(IPB, JPB)) = GetNEWTMATEl(IPB, JPB) &
& - get_umat_el(IDA, IDI, IDJ, IDA)
ELSE
if(tOneElecDiag) then
if(IPB == JPB) then
TMAT2D2(IPB, 1) = GetNEWTMATEl(IPB, JPB) &
& - get_umat_el(IDA, IDI, IDJ, IDA)
else
if(abs(get_umat_el(IDA, IDI, IDJ, IDA)) > 1.0e-8_dp) then
call stop_all("", "Error here in freezing for UEG")
end if
end if
else
! IF(IPB.eq.0.or.JPB.eq.0) CALL Stop_All("","here 03")
TMAT2D2(IPB, JPB) = GetNEWTMATEl(IPB, JPB) &
& - get_umat_el(IDA, IDI, IDJ, IDA)
end if
end if
end if
end do
DO A = NEL - NFROZENIN + 1, NEL
AB = BRR(A)
IDA = GTID(AB)
!C.. SGN takes into account permutationnness.
!C SGN=1
!C IF(IB.GT.AB) SGN=-SGN
!C IF(JB.GT.AB) SGN=-SGN
IF(G1(IB)%MS == G1(JB)%MS) THEN
IF(tCPMDSymTMat) THEN
TMATSYM2(NEWTMATInd(IPB, JPB)) = &
& GetNEWTMATEl(IPB, JPB) + get_umat_el(IDA, IDI, IDA, IDJ)
ELSE
if(tOneElecDiag) then
if(IPB == JPB) then
TMAT2D2(IPB, 1) = TMAT2D2(IPB, 1) + get_umat_el(IDA, IDI, IDA, IDJ)
else
if(abs(get_umat_el(IDA, IDI, IDA, IDJ)) > 1.0e-8_dp) then
call stop_all("", "Error here in freezing for UEG")
end if
end if
else
! IF(IPB.eq.0.or.JPB.eq.0) CALL Stop_All("","here 04")
TMAT2D2(IPB, JPB) = TMAT2D2(IPB, JPB) + get_umat_el(IDA, IDI, IDA, IDJ)
end if
end if
end if
!C.. If we have spin-independent integrals, ISS.EQ.2.OR
!C.. if the spins are the same
IF(G1(IB)%MS == G1(AB)%MS .AND. G1(AB)%MS == G1(JB)%MS) THEN
IF(tCPMDSymTMat) THEN
TMATSYM2(NEWTMATInd(IPB, JPB)) = GetNEWTMATEl(IPB, JPB) &
& - get_umat_el(IDA, IDI, IDJ, IDA)
ELSE
if(tOneElecDiag) then
if(IPB == JPB) then
TMAT2D2(IPB, 1) = GetNEWTMATEl(IPB, JPB) &
& - get_umat_el(IDA, IDI, IDJ, IDA)
else
if(abs(get_umat_el(IDA, IDI, IDJ, IDA)) > 1.0e-8_dp) then
call stop_all("", "Error here in freezing for UEG")
end if
end if
else
! IF(IPB.eq.0.or.JPB.eq.0) CALL Stop_All("","here 05")
TMAT2D2(IPB, JPB) = GetNEWTMATEl(IPB, JPB) &
& - get_umat_el(IDA, IDI, IDJ, IDA)
end if
end if
end if
end do
! write(stdout,*) "T",TMAT(IB,JB),I,J,TMAT2(IPB,JPB)
! IF(abs(TMAT(IPB,JPB)).gt.1.0e-9_dp) write(16,*) I,J,TMAT2(IPB,JPB)
end do
end do
end do
end do
! Reorganize the one-body integrals, no corrections are needed for the one-body integrals of the property integrals as long as corresponding pertubation operator does not have any two-body components.
IF(tCalcPropEst) then
DO W = 1, 2
IF(W == 1) THEN
BLOCKMINW = 1
BLOCKMAXW = NEL - NFROZEN - NFROZENIN
FROZENBELOWW = NFROZEN
else if(W == 2) THEN
BLOCKMINW = NEL - NFROZEN - NFROZENIN + 1
BLOCKMAXW = NBASIS
FROZENBELOWW = NFROZEN + NFROZENIN + NTFROZENIN
end if
DO I = BLOCKMINW, BLOCKMAXW
IP = I + FROZENBELOWW
IB = BRR(IP)
IPB = GG(IB)
DO Y = 1, 2
IF(Y == 1) THEN
BLOCKMINY = 1
BLOCKMAXY = NEL - NFROZEN - NFROZENIN
FROZENBELOWY = NFROZEN
else if(Y == 2) THEN
BLOCKMINY = NEL - NFROZEN - NFROZENIN + 1
BLOCKMAXY = NBASIS
FROZENBELOWY = NFROZEN + NFROZENIN + NTFROZENIN
end if
DO J = BLOCKMINY, BLOCKMAXY
JP = J + FROZENBELOWY
JB = BRR(JP)
JPB = GG(JB)
IF(tCPMDSymTMat) THEN
call stop_all("IntFreezeBasis", "Not implemented for tCPMD")
ELSE
! IF(IPB.eq.0.or.JPB.eq.0) THEN
! write(stdout,*) 'W',W,'I',I,'J',J,'IPB',IPB,'JPB',JPB
! CALL neci_flush(stdout)
! CALL Stop_All("","here 01")
! end if
if(tOneElecDiag) then
call stop_all("IntFreezeBasis", "Not implemented for tOneElecDiag")
else
OneEPropInts2(IPB, JPB, :) = OneEPropInts(IB, JB, :)
end if
end if
! write(stdout,*) "T",TMAT(IB,JB),I,J,TMAT2(IPB,JPB)
! IF(abs(TMAT(IPB,JPB)).gt.1.0e-9_dp) write(16,*) I,J,TMAT2(IPB,JPB)
end do
end do
end do
end do
end if
if (tCalcWithField) then
DO W=1,2
IF(W.eq.1) THEN
BLOCKMINW=1
BLOCKMAXW=NEL-NFROZEN-NFROZENIN
FROZENBELOWW=NFROZEN
ELSEIF(W.eq.2) THEN
BLOCKMINW=NEL-NFROZEN-NFROZENIN+1
BLOCKMAXW=NBASIS
FROZENBELOWW=NFROZEN+NFROZENIN+NTFROZENIN
ENDIF
DO I=BLOCKMINW,BLOCKMAXW
IP=I+FROZENBELOWW
IB=BRR(IP)
IPB=GG(IB)
DO Y=1,2
IF(Y.eq.1) THEN
BLOCKMINY=1
BLOCKMAXY=NEL-NFROZEN-NFROZENIN
FROZENBELOWY=NFROZEN
ELSEIF(Y.eq.2) THEN
BLOCKMINY=NEL-NFROZEN-NFROZENIN+1
BLOCKMAXY=NBASIS
FROZENBELOWY=NFROZEN+NFROZENIN+NTFROZENIN
ENDIF
DO J=BLOCKMINY,BLOCKMAXY
JP=J+FROZENBELOWY
JB=BRR(JP)
JPB=GG(JB)
IF(tCPMDSymTMat) THEN
call stop_all("IntFreezeBasis","Not implemented for tCPMD")
ELSE
! IF(IPB.eq.0.or.JPB.eq.0) THEN
! WRITE(stdout,*) 'W',W,'I',I,'J',J,'IPB',IPB,'JPB',JPB
! CALL neci_flush(stdout)
! CALL Stop_All("","here 01")
! ENDIF
if(tOneElecDiag) then
call stop_all("IntFreezeBasis","Not implemented for tOneElecDiag")
else
OneEFieldInts2(IPB,JPB,:)=OneEFieldInts(IB,JB,:)
endif
ENDIF
! WRITE(stdout,*) "T",TMAT(IB,JB),I,J,TMAT2(IPB,JPB)
! IF(abs(TMAT(IPB,JPB)).gt.1.0e-9_dp) WRITE(16,*) I,J,TMAT2(IPB,JPB)
ENDDO
ENDDO
ENDDO
ENDDO
endif
IF(NBASISMAX(1, 3) >= 0 .AND. ISS /= 0) THEN
!CC Only do the below if we've a stored UMAT
!C.. Now copy the relevant matrix elements of UMAT across
!C.. the primed (...P) are the new versions
DO W = 1, 2
IF(W == 1) THEN
BLOCKMINW = 1
BLOCKMAXW = NEL - NFROZEN - NFROZENIN
FROZENBELOWW = NFROZEN
ELSEIF(W == 2) THEN
BLOCKMINW = NEL - NFROZEN - NFROZENIN + 1
BLOCKMAXW = NBASIS
FROZENBELOWW = NFROZEN + NFROZENIN + NTFROZENIN
ENDIF
DO I = BLOCKMINW, BLOCKMAXW
IB = BRR(I + FROZENBELOWW)
IPB = GG(IB)
IF(ISS /= 0 .OR. G1(I)%MS == 1) THEN
IDI = GTID(IB)
IDIP = GTID(IPB)
DO X = 1, 2
IF(X == 1) THEN
BLOCKMINX = 1
BLOCKMAXX = NEL - NFROZEN - NFROZENIN
FROZENBELOWX = NFROZEN
ELSEIF(X == 2) THEN
BLOCKMINX = NEL - NFROZEN - NFROZENIN + 1
BLOCKMAXX = NBASIS
FROZENBELOWX = NFROZEN + NFROZENIN + NTFROZENIN
ENDIF
DO J = BLOCKMINX, BLOCKMAXX
JB = BRR(J + FROZENBELOWX)
JPB = GG(JB)
IF(ISS /= 0 .OR. G1(I)%MS == 1) THEN
IDJ = GTID(JB)
IDJP = GTID(JPB)
DO Y = 1, 2
IF(Y == 1) THEN
BLOCKMINY = 1
BLOCKMAXY = NEL - NFROZEN - NFROZENIN
FROZENBELOWY = NFROZEN
ELSEIF(Y == 2) THEN
BLOCKMINY = NEL - NFROZEN - NFROZENIN + 1
BLOCKMAXY = NBASIS
FROZENBELOWY = NFROZEN + NFROZENIN + NTFROZENIN
ENDIF
DO K = BLOCKMINY, BLOCKMAXY
KB = BRR(K + FROZENBELOWY)
KPB = GG(KB)
IF(ISS /= 0 .OR. G1(I)%MS == 1) THEN
IDK = GTID(KB)
IDKP = GTID(KPB)
DO Z = 1, 2
IF(Z == 1) THEN
BLOCKMINZ = 1
BLOCKMAXZ = NEL - NFROZEN - NFROZENIN
FROZENBELOWZ = NFROZEN
ELSEIF(Z == 2) THEN
BLOCKMINZ = NEL - NFROZEN - NFROZENIN + 1
BLOCKMAXZ = NBASIS
FROZENBELOWZ = NFROZEN + NFROZENIN + NTFROZENIN
ENDIF
DO L = BLOCKMINZ, BLOCKMAXZ
IF((K * (K - 1)) / 2 + I >= (L * (L - 1)) / 2 + J) THEN
LB = BRR(L + FROZENBELOWZ)
LPB = GG(LB)
IF(ISS /= 0 .OR. G1(I)%MS == 1) THEN
IDL = GTID(LB)
IDLP = GTID(LPB)
UMAT2(UMat2Ind(IDIP, IDJP, IDKP, IDLP)) = &
& UMAT(UMatInd(IDI, IDJ, IDK, IDL))
ENDIF
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
CALL neci_flush(11)
CALL neci_flush(12)
ELSEIF(Associated(UMatCacheData)) THEN
!.. We've a UMAT2D and a UMATCACHE. Go and Freeze them
!C.. NHG contains the old number of orbitals
!C.. NBASIS contains the new
!C.. GG(I) is the new position in G of the (old) orb I
CALL FreezeUMatCache(GG, NHG, NBASIS)
end if
IF(ISS == 0) CALL SetupUMatTransTable(GG, nHG, nBasis)
! do i=1,NBASIS
! write(stdout,*) i,G2(i)%Sym%S
! end do
IF(NFROZENIN > 0) THEN
! CALL GETSYM(BRR,NFROZEN,BRR((NEL-NFROZENIN+1):NEL),G1,NBASISMAX,KSym)
!This is slightly dodgey... I have commented out the above routine that finds the symmetry of the frozen orbitals.
!There is already a test to check we are not freezing in the middle of degenracies, so the symmetry should always come
!out as 0 anyway, unless we are breaking up a set of symmetry irreps somehow ... I think.
write(stdout, *) "WARNING: Setting the symmetry of the frozen orbitals to 0. This will be incorrect &
&if orbitals are frozen in the middle of a degenerate set of the same symmetery irrep."
KSym%Sym%S = 0
CALL SetupFREEZEALLSYM(KSym)
else if(NFROZEN > 0) THEN
CALL GETSYM(BRR, NFROZEN, G1, NBASISMAX, KSym)
! write(stdout,*) '************'
! write(stdout,*) 'KSym',KSym
CALL SetupFREEZEALLSYM(KSym)
END IF
CALL FREEZESYMLABELS(NHG, NBASIS, GG, .false.)
do i = 1, NBASIS
G1(i) = G2(i)
end do
FREEZETRANSFER = .false.
!C.. Copy the new BRR and ARR over the old ones
CALL SWAPTMAT()
IF(tCalcPropEst) call SwapOneEPropInts(nBasis, iNumPropToEst)
IF(tCalcWithField) call SwapOneEFieldInts(nBasis,nFields_it)
deallocate(arr)
LogDealloc(tagarr)
allocate(arr(nBasis, 2), stat=ierr)
LogAlloc(ierr, 'Arr', 2 * nBasis, 8, tagArr)
CALL NECI_ICOPY(NBASIS, BRR2, 1, BRR, 1)
CALL DCOPY(NBASIS * 2, ARR2, 1, ARR, 1)
!C.. Now reset the total number of orbitals
NHG = NBASIS
Call DetFreezeBasis(GG)
RETURN
end subroutine intfreezebasis