subroutine FindNatOrbitals() ! This routine simply takes a transformation matrix and rotates the integrals to produce a new FCIDUMP file. ! In one case the transformation matrix is read in from a file TRANSFORMMAT. ! In the other, the transformation matrix is calculated from the MP2 variational density matrix. ! MP2VDM = D2_ab = sum_ijc [ t_ij^ac ( 2 t_ij^bc - t_ji^bc ) ] ! Where : t_ij^ac = - < ab | ij > / ( E_a - E_i + E_b - Ej ) ! Ref : J. Chem. Phys. 131, 034113 (2009) - note: in Eqn 1, the cb indices are the wrong way round (should be bc). use NatOrbsMod, only: SetUpNatOrbLabels, FindNatOrbs, FillCoeffT1, DeallocateNatOrbs, PrintOccTable integer :: i, a, ierr, MinReadIn, MaxReadIn, iunit character(len=*), parameter :: this_routine = 'FindNatOrbitals' if (tUseMP2VarDenMat) write(stdout, *) '*** Transforming the HF orbitals into the MP2 approximate natural orbitals. ***' if (tFindCINatOrbs) then write(stdout, *) '*** Transforming the HF orbitals into approximate natural orbitals' write(stdout, *) 'based on the one-electron density matrix found from the wavefunction calculated above. ***' end if if (tSpinOrbs) then if (.not. tStoreSpinOrbs) then write(stdout, *) "We want to use spin orbitals - turning on tStoreSpinOrbs." tStoreSpinOrbs = .true. end if end if if (tROHF .and. tStoreSpinOrbs) call Stop_All(this_routine, "Cannot compress open shell systems into spatial " & //"orbitals when rotating, turn off ROHF.") if (tTruncRODump .and. (.not. tTruncDumpbyVal)) then NoFrozenVirt = NoTruncOrbs(1) else if (tTruncRODump) then ! If the 'number of frozen orbitals' is given as a cutoff - take NoFrozenVirt to be 0 !for all the allocation purposes - will set this later when ! we have the eigenvalues and know how many orbitals lie below it. NoFrozenVirt = 0 TruncEval = TruncEvalues(1) else NoFrozenVirt = 0 end if SpatOrbs = nBasis / 2 if (tStoreSpinOrbs) then NoOrbs = nBasis NoOcc = NEl MinReadIn = 1 MaxReadIn = nBasis if (tRotateVirtOnly) MinReadIn = NEl + 1 if (tRotateOccOnly) MaxReadIn = NEl ! If tStoreSpinOrbs ARR(:,2) is not filled, but we want to use it later, so just fill it here. do i = 1, NoOrbs ARR(BRR(i), 2) = ARR(i, 1) end do allocate(SymLabelCounts2_rot(2, 32), stat=ierr) call LogMemAlloc('SymLabelCounts2_rot', 2 * 32, 4, this_routine, SymLabelCounts2_rotTag, ierr) SymLabelCounts2_rot(:, :) = 0 ! first 8 refer to the occupied, and the second to the virtual beta spin. ! third and fourth to the occupied and virtual alpha spin. else NoOrbs = SpatOrbs NoOcc = NEl / 2 MinReadIn = 1 MaxReadIn = SpatOrbs if (tRotateVirtOnly) MinReadIn = (NEl / 2) + 1 if (tRotateOccOnly) MaxReadIn = NEl / 2 allocate(SymLabelCounts2_rot(2, 16), stat=ierr) call LogMemAlloc('SymLabelCounts2_rot', 2 * 16, 4, this_routine, SymLabelCounts2_rotTag, ierr) SymLabelCounts2_rot(:, :) = 0 ! first 8 refer to the occupied, and the second to the virtual. end if NoRotOrbs = NoOrbs call ApproxMemReq() allocate(SymLabelList2_rot(NoOrbs), stat=ierr) call LogMemAlloc('SymLabelList2_rot', NoOrbs, 4, this_routine, SymLabelList2_rotTag, ierr) SymLabelList2_rot(:) = 0 allocate(SymLabelList3_rot(NoOrbs), stat=ierr) call LogMemAlloc('SymLabelList3_rot', NoOrbs, 4, this_routine, SymLabelList3_rotTag, ierr) SymLabelList3_rot(:) = 0 allocate(SymLabelListInv_rot(NoOrbs), stat=ierr) call LogMemAlloc('SymLabelListInv_rot', NoOrbs, 4, this_routine, SymLabelListInv_rotTag, ierr) SymLabelListInv_rot(:) = 0 if (tReadInCoeff .or. tUseHFOrbs) then ! No symmetry, so no reordering of the orbitals - symlabellist just goes from 1-NoOrbs. ! When we are just reading in the coefficients and transforming, it does not matter about the ordering of the orbitals. do i = 1, NoOrbs SymLabelList2_rot(i) = i SymLabelListInv_rot(i) = i end do else if (tFindCINatOrbs .or. tUseMP2VarDenMat) then call SetupNatOrbLabels() end if ! Yet another labelling system, SymLabelList3_rot is created here. ! This indicates the label of the transformed orbital. ! In the case where we are truncating the space, the transformed orbitals are ordered according !to the size of the eigenvalues of the MP2VDM ! matrix when it is diagonalised. We wish to keep them in this order when transforming !the integrals etc, so that when we truncate the last ! NoFrozenVirt orbitals, we are removing those with the smallest MP2VDM eigenvalues (occupation numbers). ! In the case where no truncation is made however, SymLabelList3_rot is the same as SymLabelList2_rot, !so that the indexes remain the same as previously. ! This allows for the option of going straight into a spawning calc from the rotation, which is not !possible when a truncation is performed ! because of the messed up indices. if (tTruncRODump) then if (MOD(NoFrozenVirt, 2) /= 0) call Stop_All(this_routine, "Must freeze virtual spin orbitals in pairs of 2.") if (tStoreSpinOrbs) then NoRotOrbs = NoOrbs - NoFrozenVirt else NoFrozenVirt = NoFrozenVirt / 2 NoRotOrbs = NoOrbs - NoFrozenVirt end if do i = 1, NoOrbs SymLabelList3_rot(i) = i end do else do i = 1, NoOrbs SymLabelList3_rot(i) = SymLabelList2_rot(i) end do end if allocate(CoeffT1(NoOrbs, NoRotOrbs), stat=ierr) call LogMemAlloc(this_routine, NoRotOrbs * NoOrbs, 8, this_routine, CoeffT1Tag, ierr) CoeffT1(:, :) = 0.0_dp if (tSeparateOccVirt) then do i = 1, NoRotOrbs CoeffT1(i, i) = 1.0_dp end do end if if (tUEG) then call FindNatOrbs() call FillCoeffT1() else if (tReadInCoeff) then write(stdout, '(A)') " Reading in the transformation matrix from TRANSFORMMAT, and using this to rotate the HF orbitals." iunit = get_free_unit() open(iunit, file='TRANSFORMMAT', status='old') do i = 1, NoOrbs do a = 1, NoOrbs read(iunit, *) CoeffT1(a, i) end do end do close(iunit) else if (tFindCINatOrbs .or. tUseMP2VarDenMat .or. tUseHFOrbs) then if (.not. tUseHFOrbs) call FindNatOrbs() if (tUseHFOrbs) then call PrintOccTable() else call FillCoeffT1() end if end if if (tPrintRODump) then allocate(FourIndInts(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr) call LogMemAlloc('FourIndInts', (NoOrbs**4), 8, this_routine, FourIndIntsTag, ierr) ! Then, transform2ElInts write(stdout, *) 'Transforming the four index integrals' call Transform2ElIntsMemSave() write(stdout, *) 'Re-calculating the fock matrix' call CalcFOCKMatrix() write(stdout, *) 'Refilling the UMAT and TMAT2D' ! The ROFCIDUMP is also printed out in here. call RefillUMATandTMAT2D() call neci_flush(stdout) if ((tFindCINatOrbs .or. tUseMP2VarDenMat) .and. (NoDumpTruncs > 1)) call ReTruncROFciDump() if ((.not. tUseHFOrbs) .and. (.not. tReadInCoeff)) call DeallocateNatOrbs() end if if (tWriteTransMat) call WriteTransformMat() ! If a truncation is being made, the new basis will not be in the correct energetic ordering - this does not matter, as we ! never go straight into a spawning and they will be reordered when the ROFCIDUMP file is read in again. call writebasis(stdout, G1, nBasis, ARR, BRR) deallocate(CoeffT1) call LogMemDeAlloc(this_routine, CoeffT1Tag) deallocate(SymLabelList2_rot) call LogMemDeAlloc(this_routine, SymLabelList2_rotTag) deallocate(SymLabelListInv_rot) call LogMemDeAlloc(this_routine, SymLabelListInv_rotTag) if (tPrintRODump) then deallocate(FourIndInts) call LogMemDeAlloc(this_routine, FourIndIntsTag) end if end if end subroutine FindNatOrbitals