Subroutine DetCalcInit
Use global_utilities
Use IntegralsData, only: NFROZEN
use SystemData, only: lms, lms2, nBasis, nBasisMax, nEl, SymRestrict
use SystemData, only: Alat, arr, brr, boa, box, coa, ecore, g1, Beta
use SystemData, only: tParity, tSpn, Symmetry, STot, NullBasisFn, tUHF, tMolpro
use sym_mod
use LoggingData, only: tLogDets
use HElem
use util_mod, only: get_free_unit, NECI_ICOPY
Type(BasisFn) ISym
integer i, j, ii, iunit
integer ierr, norb
integer nDetTot
character(25), parameter :: this_routine = 'DetCalcInit'
IF (.NOT. TCALCHMAT) THEN
write(stdout, *) "Not storing the H matrix."
IF (TENERGY .AND. .NOT. TBLOCK) THEN
write(stdout, *) "Cannot calculate energies without blocking the Hamiltonian."
TENERGY = .FALSE.
end if
IF (TENERGY .AND. NBLK /= 0) THEN
!C.. We're doing a Lanczos without calculating the H mat
write(stdout, *) "Cannot perform Lanczos without Hamiltonian"
TENERGY = .FALSE.
end if
end if
! If we want to have UMat2D, and it isn't yet filled in, generate it
! here. All of the integrals setup/freezing etc is done...
if (tDeferred_Umat2d .and. .not. tUMat2D) then
ASSERT(.not. btest(nbasis, 0))
! And fill in the array
call SetupUMat2d_dense(nBasis)
end if
!Copied Specdet information from Calc.F, so if inspect is present, but no determinant/csf specified, it will still run.
if (TSPECDET) then
if (.not. associated(specdet)) then
!specdet not allocated. Allocate it and copy fdet
allocate(specdet(nel))
write(stdout, *) "TSPECDET set, but not allocated. using FDET"
CALL NECI_ICOPY(NEL, FDET, 1, SPECDET, 1)
else if (.not. ISVALIDDET(SPECDET, NEL)) then
write(stdout, *) "TSPECDET set, but invalid. using FDET"
CALL NECI_ICOPY(NEL, FDET, 1, SPECDET, 1)
end if
end if
!C IF(TCALCHMAT.OR.NPATHS.NE.0.OR.DETINV.GT.0.OR.TBLOCK) THEN
iExcitLevel = ICILEVEL
IF (tFindDets) THEN
!C..Need to determine the determinants
IF (iExcitLevel /= 0) THEN
write(stdout, *) "Performing truncated CI at level ", iExcitLevel
IF (TSPECDET) THEN
write(stdout, *) "Using SPECDET:"
call write_det(stdout, SPECDET, .true.)!
CALL NECI_ICOPY(NEL, SPECDET, 1, FDET, 1)
ELSE
write(stdout, *) "Using Fermi DET:"
call write_det(stdout, FDET, .true.)
end if
!C.. if we're doing a truncated CI expansion
CALL GENEXCIT(FDET, iExcitLevel, NBASIS, NEL, reshape([0], [1, 1]), &
reshape([0], [1, 1]), NDET, 1, G1, .TRUE., NBASISMAX, .TRUE.)
write(stdout, *) "NDET out of GENEXCIT ", NDET
!C.. We need to add in the FDET
NDET = NDET + 1
II = NDET
NBLOCKS = 1
else if (TBLOCK) THEN
write(stdout, *) "Determining determinants and blocks."
IF (TPARITY) THEN
write(stdout, *) "Using symmetry restriction:"
CALL writeallsym(stdout, SymRestrict, .TRUE.)
end if
IF (TSPN) THEN
write(stdout, *) "Using spin restriction:", LMS
end if
if (tUHF .and. tMolpro) then
!When breaking spin symmetry in molpro, it is important to occupy alpha orbs preferentially
CALL GNDTS_BLK(NEL, nBasis, BRR, NBASISMAX, NMRKS, .TRUE., &
& NDET, G1, II, NBLOCKSTARTS, NBLOCKS, TSPN, -LMS2, TPARITY, &
& SymRestrict, IFDET,.NOT. TREAD, NDETTOT, BLOCKSYM)
else
CALL GNDTS_BLK(NEL, nBasis, BRR, NBASISMAX, NMRKS, .TRUE., &
& NDET, G1, II, NBLOCKSTARTS, NBLOCKS, TSPN, LMS2, TPARITY, &
& SymRestrict, IFDET,.NOT. TREAD, NDETTOT, BLOCKSYM)
end if
write(stdout, *) "NBLOCKS:", NBLOCKS
write(stdout, *) "Determining determinants."
IF (TPARITY) THEN
write(stdout, *) "Using symmetry restriction:"
CALL writeallsym(stdout, SymRestrict, .TRUE.)
end if
IF (TSPN) THEN
write(stdout, *) "Using spin restriction:", LMS
end if
if (tUHF .and. tMolpro) then
!When breaking spin symmetry in molpro, it is important to occupy alpha orbs preferentially
CALL GNDTS(NEL, nBasis, BRR, NBASISMAX, NMRKS, .TRUE., G1, TSPN, -LMS, TPARITY, SymRestrict, II, IFDET)
else
CALL GNDTS(NEL, nBasis, BRR, NBASISMAX, NMRKS, .TRUE., G1, TSPN, LMS, TPARITY, SymRestrict, II, IFDET)
end if
NBLOCKS = 1
NDET = II
end if
!C..
IF (II == 0) THEN
write(stdout, *) "No determinants found. Cannot continue"
call stop_all(this_routine, "No determinants found. Cannot continue")
end if
!C.. NEL now only includes active electrons
write(stdout, *) "Number of determinants found to be: ", II
write(stdout, *) "Allocating initial memory for calculation of energy..."
CALL neci_flush(stdout)
allocate(NMrks(nEl, II), stat=ierr)
LogAlloc(ierr, 'NMRKS', NEL * II, 4, tagNMRKS)
NMRKS(1:NEL, 1:II) = 0
allocate(NBLOCKSTARTS(NBLOCKS + 1), stat=ierr)
call LogMemAlloc('NBLOCKSTARTS', NBLOCKS + 1, 4, this_routine, tagNBLOCKSTARTS, ierr)
NBLOCKSTARTS(1:NBLOCKS + 1) = 0
allocate(BlockSym(NBLOCKS + 1), stat=ierr)
LogAlloc(ierr, 'BLOCKSYM', NBLOCKS + 1, BasisFNSizeB, tagBlockSym)
BLOCKSYM(1:NBLOCKS) = NullBasisFn
!C..
NDET = II
IF (iExcitLevel /= 0) THEN
!C.. Use HAMIL to temporarily hold a list of excitation levels
CALL NECI_ICOPY(NEL, FDET, 1, NMRKS, 1)
allocate(Hamil(II), stat=ierr)
LogAlloc(ierr, 'HAMIL', II, HElement_t_sizeB, tagHamil)
NDET = 0
CALL GENEXCIT(&
FDET, iExcitLevel, NBASIS, NEL, &
NMRKS(:1, :2), int(HAMIL), NDET, 1, G1, .TRUE., NBASISMAX, .FALSE.)
Deallocate(Hamil)
LogDealloc(tagHamil)
NDET = NDET + 1
NBLOCKSTARTS(1) = 1
NBLOCKSTARTS(2) = II + 1
IFDET = 1
else if (TBLOCK) THEN
if (tUHF .and. tMolpro) then
!When breaking spin symmetry in molpro, it is important to occupy alpha orbs preferentially
CALL GNDTS_BLK(NEL, nBasis, BRR, NBASISMAX, NMRKS, .FALSE., NDET, G1, II, NBLOCKSTARTS, NBLOCKS, TSPN, -LMS2, TPARITY, &
& SymRestrict, IFDET,.NOT. TREAD, NDETTOT, BLOCKSYM)
else
CALL GNDTS_BLK(NEL, nBasis, BRR, NBASISMAX, NMRKS, .FALSE., NDET, G1, II, NBLOCKSTARTS, NBLOCKS, TSPN, LMS2, TPARITY, &
& SymRestrict, IFDET,.NOT. TREAD, NDETTOT, BLOCKSYM)
end if
ELSE
if (tUHF .and. tMolpro) then
!When breaking spin symmetry in molpro, it is important to occupy alpha orbs preferentially
CALL GNDTS(NEL, nBasis, BRR, NBASISMAX, NMRKS, .FALSE., G1, TSPN, -LMS, TPARITY, SymRestrict, II, IFDET)
else
CALL GNDTS(NEL, nBasis, BRR, NBASISMAX, NMRKS, .FALSE., G1, TSPN, LMS, TPARITY, SymRestrict, II, IFDET)
end if
NBLOCKSTARTS(1) = 1
NBLOCKSTARTS(2) = II + 1
end if
if (tLogDets) THEN
iunit = get_free_unit()
open(iunit, FILE='DETS', STATUS='UNKNOWN')
DO I = 1, NDET
call write_det(iunit, NMRKS(:, I), .false.)
CALL GETSYM(NMRKS(:, I), NEL, G1, NBASISMAX, ISYM)
CALL WRITESYM(iunit, ISym%Sym, .TRUE.)
end do
close(iunit)
end if
!Update: 14.03.2018, K.Ghanem
!FDET has already been assigned in DetPreFreezeInit.
!No idea why it is overwirtten here.
!Therefore, I comment out the following code and hope for the best.
!Instead, I look for FDET in the list of determinants NMRKS and assign the index to IFDET.
!C.. Now generate the fermi determiant
!C.. Work out the fermi det
IFDET = 0
DO I = 1, NDET
IF (ALL(NMRKS(:, I) == FDET)) THEN
IFDET = I
Exit
END IF
END DO
IF (IFDET == 0) call stop_all("DetCalcInit", "Fermi determinant is not found in NMRKS!")
write(stdout, *) ' NUMBER OF SYMMETRY UNIQUE DETS ', NDET
!C write(stdout,*) ' TOTAL NUMBER OF DETS.' , NDETTOT
IF (NEVAL == 0) THEN
write(stdout, *) 'NEVAL=0. Setting NEVAL=NDET'
NEVAL = NDET
end if
IF (NEVAL > NDET) THEN
write(stdout, *) 'NEVAL>NDET. Setting NEVAL=NDET'
NEVAL = NDET
end if
IF (ABS(DETINV) > NDET) THEN
write(stdout, *) 'DETINV=', DETINV, '>NDET=', NDET
write(stdout, *) 'Setting DETINV to 0'
DETINV = 0
end if
CALL neci_flush(stdout)
!C ==----------------------------------------------------------------==
!C..Set up memory for c's, nrow and the label
IF (TCALCHMAT) THEN
write(stdout, *) "CK Size", NDET * NEVAL * HElement_t_size
allocate(CkN(nDet, nEval), stat=ierr)
LogAlloc(ierr, 'CKN', nDet * nEval, HElement_t_sizeB, tagCKN)
CKN = (0.0_dp)
allocate(Ck(nDet, nEval), stat=ierr)
LogAlloc(ierr, 'CK', nDet * nEval, HElement_t_sizeB, tagCK)
CK = (0.0_dp)
allocate(W(nEval), stat=ierr)
LogAlloc(ierr, 'W', nEval, 8, tagW)
W = 0.0_dp
end if
IF (TREAD) THEN
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex with tREAD")
#else
CALL READ_PSI(BOX, BOA, COA, NDET, NEVAL, NBASISMAX, NEL, CK, W)
#endif
end if
end if
End Subroutine DetCalcInit