InitLocalOrbs Subroutine

public subroutine InitLocalOrbs()




Source Code

Source Code

    subroutine InitLocalOrbs()

        character(len=*), parameter :: this_routine = 'InitLocalOrbs'
        integer :: ierr

        ! Writing to output which PE is being maximised/minimised.
        write(stdout, *) '*****'
        if (tERLocalization) then
            write(stdout, *) "Calculating new molecular orbitals based on Edmiston-Reudenberg localisation,"
            write(stdout, *) "i.e. maximisation of the <ii|ii> integrals..."
            write(stdout, *) "*****"
        end if
        if (tVirtCoulombMax) then
            write(stdout, *) "Calculating new molecular orbitals based on maximisation of the sum of the"
            write(stdout, *) "<ij|ij> integrals, where i and j are both virtuals..."
            write(stdout, *) "*****"
        end if
        if (tOffDiagSqrdMin) then
            write(stdout, *) "Calculating new molecular orbitals based on mimimisation "
            write(stdout, *) "of <ij|kl>^2 integrals..."
            write(stdout, *) "*****"
        end if
        if (tOffDiagMin) then
            write(stdout, *) "Calculating new molecular orbitals based on mimimisation "
            write(stdout, *) "of <ij|kl> integrals..."
            write(stdout, *) "*****"
        end if
        if (tDoubExcMin) then
            write(stdout, *) "Calculating new molecular orbitals based on mimimisation "
            write(stdout, *) "of the double excitation hamiltonian elements."
            write(stdout, *) "*****"
        end if
        if (tOnePartOrbEnMax) then
            write(stdout, *) "Calculating new molecular orbitals based on maximisation "
            write(stdout, *) "of the virtual one particle orbital energies."
            write(stdout, *) "*****"
        else if (tMaxHLGap) then
            ! This will transform all the orbitals within a particlar group to
            ! have the same diagonal fock matrix element.
            write(stdout, *) "Transforming orbitals based on equating their diagonal fock matrix elements."
            write(stdout, *) "*****"
        end if

        ! Writing out which orthonormalisation method is being used...
        if (tLagrange) then
            if (tShake) then
                call neci_flush(stdout)
                call Stop_All(this_routine, "ERROR. Both LAGRANGE and SHAKE keywords present in the input. &
                & These two orthonormalisation methods clash.")
            end if
            write(stdout, *) "Using a Lagrange multiplier to attempt to rotate orbitals in a way to maintain orthonormality"
        else if (tShake) then
            write(stdout, *) "Using the shake algorithm to iteratively find lambdas which maintain "
            write(stdout, *) "orthonormalisation with rotation"
            write(stdout, *) "Explicity reorthonormalizing orbitals after each rotation."
        end if

        ! Check for a few possible errors.
        if (.not. TwoCycleSymGens) then
            call neci_flush(stdout)
            call Stop_All(this_routine, "ERROR. TwoCycleSymGens is false.  Symmetry is not abelian.")
        end if
        if ((tRotateOccOnly .or. tRotateVirtOnly) .and. (.not. tSeparateOccVirt)) then
            tSeparateOccVirt = .true.
            write(stdout, *) "NOTE. Cannot rotate only occupied or virtual without first separating them."
            write(stdout, *) "SEPARATEOCCVIRT keyword is being turned on."
        end if
        if ((tOffDiagSqrdMax .and. tOffDiagSqrdMin) .or. (tOffDiagMax .and. tOffDiagMin)) then
            call neci_flush(stdout)
            call Stop_All(this_routine, "ERROR. Cannot both maximise and minimise off diagonal elements simultaneously")
        end if
        if (tOnePartOrbEnMax .and. (.not. tSeparateOccVirt)) then
            call neci_flush(stdout)
            call Stop_All(this_routine, &
                          "ERROR. Cannot currently maximise the one particle orbital energies without separating occupied and virtual.")
        end if
        write(stdout, *) "*****"

        ! Zero values.
        OrthoNorm = 0.0_dp
        ERPotEnergy = 0.0_dp
        PotEnergy = 0.0_dp
        Force = 0.0_dp
        TwoEInts = 0.0_dp
        PEInts = 0.0_dp
        PEOrtho = 0.0_dp
        ForceInts = 0.0_dp
        DistCs = 0.0_dp
        DistLs = 0.0_dp
        LambdaMag = 0.0_dp
        SpatOrbs = nBasis / 2
        if (tStoreSpinOrbs) then
            NoOrbs = nBasis
            NoOcc = NEl
            NoOrbs = SpatOrbs
            NoOcc = NEl / 2
        end if
        NoRotOrbs = NoOrbs
        Iteration = 0
        OrthoForce = 0.0_dp
        ShakeIterInput = ShakeIterMax
        TotNoConstraints = (NoOrbs * (NoOrbs + 1)) / 2

        ! When maximising the one particle orbital energies, choose the zero
        ! value (Epsilon min).
        if (tRotateVirtOnly .and. tOnePartOrbEnMax) then
            EpsilonMin = ARR(NEl + 1, 1)
            write(stdout, *) 'Taking EpsilonMin to be the LUMO of the HF orbitals...'
            write(stdout, *) 'EpsilonMin  =  ', EpsilonMin
        else if (tOnePartOrbEnMax) then
            EpsilonMin = ChemPot
            write(stdout, *) 'Taking EpsilonMin to be the chemical potential (midway between HF HOMO and LUMO)...'
            write(stdout, *) 'therefore EpsilonMin  =  ', EpsilonMin
        end if

        ! Set timed routine names.
        Rotation_Time%timer_name = 'RotateTime'
        Shake_Time%timer_name = 'ShakeTime'
        FullShake_Time%timer_name = 'FullShakeTime'
        Findtheforce_Time%timer_name = 'FindtheForceTime'
        Transform2ElInts_Time%timer_name = 'Transform2ElIntsTime'
        findandusetheforce_time%timer_name = 'Findandusetheforce'
        CalcDerivConstr_Time%timer_name = 'CalcDerivConstr'
        TestOrthoConver_Time%timer_name = 'TestOrthoConver'

        ! Allocate memory.

        allocate(CoeffT1(NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('CoeffT1', NoOrbs**2, 8, this_routine, CoeffT1Tag, ierr)
        allocate(CoeffCorT2(NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('CoeffCorT2', NoOrbs**2, 8, this_routine, CoeffCorT2Tag, ierr)
        allocate(CoeffUncorT2(NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('CoeffUncT2', NoOrbs**2, 8, this_routine, CoeffUncorT2Tag, ierr)
        CoeffUncorT2(:, :) = 0.0_dp

        allocate(DerivCoeff(NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('DerivCoeff', NoOrbs**2, 8, this_routine, DerivCoeffTag, ierr)

        allocate(DiagTMAT2Dfull(NoOrbs - (NoOcc)), stat=ierr)
        call LogMemAlloc('DiagTMAT2Dfull', (NoOrbs - (NoOcc)), 8, this_routine, DiagTMAT2DfullTag, ierr)
        allocate(UMATTemp01(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('UMATTemp01', NoOrbs**4, 8, this_routine, UMATTemp01Tag, ierr)
        allocate(TwoIndInts01(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('TwoIndInts01', NoOrbs**4, 8, this_routine, TwoIndInts01Tag, ierr)
        allocate(ThreeIndInts02(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('ThreeIndInts02', NoOrbs**4, 8, this_routine, ThreeIndInts02Tag, ierr)
        allocate(FourIndInts(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('FourIndInts', NoOrbs**4, 8, this_routine, FourIndIntsTag, ierr)
        allocate(FourIndInts02(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('FourIndInts02', NoOrbs**4, 8, this_routine, FourIndInts02Tag, ierr)

        ! Partially transformed temporary arrays.
        if (tERLocalization .and. (.not. tStoreSpinOrbs)) then
            allocate(TwoIndIntsER(NoOrbs, NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('TwoIndIntsER', NoOrbs**3, 8, this_routine, TwoIndIntsERTag, ierr)
            allocate(ThreeIndInts01ER(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('ThreeIndInts01ER', NoOrbs**2, 8, this_routine, ThreeIndInts01ERTag, ierr)
            allocate(ThreeIndInts02ER(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('ThreeIndInts02ER', NoOrbs**2, 8, this_routine, ThreeIndInts02ERTag, ierr)
            allocate(FourIndIntsER(NoOrbs), stat=ierr)
            call LogMemAlloc('FourIndIntsER', NoOrbs, 8, this_routine, FourIndIntsERTag, ierr)
            allocate(TMAT2DTemp(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('TMAT2DTemp', NoOrbs**2, 8, this_routine, TMAT2DTempTag, ierr)
            allocate(TMAT2DPartRot01(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('TMAT2DPartRot01', NoOrbs**2, 8, this_routine, TMAT2DPartRot01Tag, ierr)
            allocate(TMAT2DPartRot02(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('TMAT2DPartRot02', NoOrbs**2, 8, this_routine, TMAT2DPartRot02Tag, ierr)
            allocate(TMAT2DRot(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('TMAT2DRot', NoOrbs**2, 8, this_routine, TMAT2DRotTag, ierr)
            allocate(UMATTemp02(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('UMATTemp02', NoOrbs**4, 8, this_routine, UMATTemp02Tag, ierr)

            ! Partially transformed combined arrays.
            allocate(TwoIndInts02(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('TwoIndInts02', NoOrbs**4, 8, this_routine, TwoIndInts02Tag, ierr)
            allocate(ThreeIndInts01(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('ThreeIndInts01', NoOrbs**4, 8, this_routine, ThreeIndInts01Tag, ierr)
            allocate(ThreeIndInts03(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('ThreeIndInts03', NoOrbs**4, 8, this_routine, ThreeIndInts03Tag, ierr)
            allocate(ThreeIndInts04(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('ThreeIndInts04', NoOrbs**4, 8, this_routine, ThreeIndInts04Tag, ierr)
        end if

        ! Allocate according to orthonormalisation method being used.
        if (tLagrange) then
            allocate(Lambdas(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('Lambdas', NoOrbs**2, 8, this_routine, LambdasTag, ierr)
            allocate(DerivLambda(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('DerivLambda', NoOrbs**2, 8, this_routine, DerivLambdaTag, ierr)
            Lambdas(:, :) = 0.0_dp
            DerivLambda(:, :) = 0.0_dp
        end if

        if (tShake) then
            allocate(ShakeLambda(TotNoConstraints), stat=ierr)
            call LogMemAlloc('ShakeLambda', TotNoConstraints, 8, this_routine, ShakeLambdaTag, ierr)
            ShakeLambda(:) = 0.0_dp
            allocate(ShakeLambdaNew(TotNoConstraints), stat=ierr)
            call LogMemAlloc('ShakeLambdaNew', TotNoConstraints, 8, this_routine, ShakeLambdaNewTag, ierr)
            ShakeLambdaNew(:) = 0.0_dp
            allocate(Constraint(TotNoConstraints), stat=ierr)
            call LogMemAlloc('Constraint', TotNoConstraints, 8, this_routine, ConstraintTag, ierr)
            allocate(ConstraintCor(TotNoConstraints), stat=ierr)
            call LogMemAlloc('ConstraintCor', TotNoConstraints, 8, this_routine, ConstraintCorTag, ierr)
            allocate(DerivConstrT1(NoOrbs, NoOrbs, TotNoConstraints), stat=ierr)
            call LogMemAlloc('DerivConstrT1', NoOrbs * TotNoConstraints * NoOrbs, 8, this_routine, DerivConstrT1Tag, ierr)
            DerivConstrT1(:, :, :) = 0.0_dp
            allocate(DerivConstrT2(NoOrbs, NoOrbs, TotNoConstraints), stat=ierr)
            call LogMemAlloc('DerivConstrT2', NoOrbs * TotNoConstraints * NoOrbs, 8, this_routine, DerivConstrT2Tag, ierr)
            allocate(ForceCorrect(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('ForceCorrect', NoOrbs**2, 8, this_routine, ForceCorrectTag, ierr)
            allocate(Correction(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('Correction', NoOrbs**2, 8, this_routine, CorrectionTag, ierr)
            if (tShakeApprox) then
                allocate(DerivConstrT1T2Diag(TotNoConstraints), stat=ierr)
                call LogMemAlloc('DerivConstrT1T2Diag', TotNoConstraints, 8, this_routine, DerivConstrT1T2DiagTag, ierr)
                DerivConstrT1T2Diag(:) = 0.0_dp
                allocate(DerivConstrT1T2(TotNoConstraints, TotNoConstraints), stat=ierr)
                call LogMemAlloc('DerivConstrT1T2', TotNoConstraints**2, 8, this_routine, DerivConstrT1T2Tag, ierr)
            end if
        end if

        ! Indexing arrays.
        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

        allocate(Lab(2, TotNoConstraints), stat=ierr)
        call LogMemAlloc('Lab', 2 * TotNoConstraints, 4, this_routine, LabTag, ierr)
        Lab(:, :) = 0

        ! Do any initial calculations, and set up starting values for arrays
        ! used in rotation.
        call InitRotCalc()

        ! Write out the headings for the results file.
        transform_unit = get_free_unit()
        open(transform_unit, file='Transform', status='unknown')
        if (tLagrange) then
            write(transform_unit, "(A12, 11A18)") "# Iteration", "2.PotEnergy", "3.PEInts", "4.PEOrtho", "5.Force", "6.ForceInts", &
            "7.OrthoForce", "8.Sum<ij|kl>^2",&
                        &"9.OrthoNormCondition", "10.DistMovedbyCs", "11.DistMovedByLs", "12.LambdaMag"
            write(stdout, "(A12, 11A19)") "Iteration", "2.PotEnergy", "3.PEInts", "4.PEOrtho", "5.Force", "6.ForceInts", "7.OrthoForce", &
                    &"9.OrthoNormCondition", "10.DistMovedbyCs", "11.DistMovedbyLs", "12.LambdaMag"
        else if (tERLocalization .and. tHijSqrdMin) then
            write(transform_unit, "(A12, 7A24)") "# Iteration", "2.ERPotEnergy", "3.HijSqrdPotEnergy", "4.PotEnergy", "5.Force", &
                "6.Totalcorrforce", "7.OrthoNormCondition", "8.DistMovedbyCs"
            write(stdout, "(A12, 7A24)") "# Iteration", "2.ERPotEnergy", "3.HijSqrdPotEnergy", "4.PotEnergy", "5.Force", &
                "6.Totalcorrforce", "7.OrthoNormCondition", "8.DistMovedbyCs"
        else if (tERLocalization) then
            write(transform_unit, "(A12, 5A24)") "# Iteration", "2.Sum_i<ii|ii>", "3.Force", "4.TotCorrForce", &
                "5.OrthoNormCondition", "6.DistMovedbyCs"
            write(stdout, "(A12, 5A24)") "Iteration", "2.Sum_i<ii|ii>", "3.Force", "4.TotCorrForce", "5.OrthoNormCondition", "6.DistMovedbyCs"
            write(transform_unit, "(A12, 5A24)") "# Iteration", "2.PotEnergy", "3.Force", "4.Totalcorrforce", &
                "5.OrthoNormCondition", "6.DistMovedbyCs"
            write(stdout, "(A12, 5A24)") "Iteration", "2.PotEnergy", "3.Force", "4.TotCorrForce", "5.OrthoNormCondition", "6.DistMovedbyCs"
        end if

    end subroutine InitLocalOrbs