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" else 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 else 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) else 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 else 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", & "8.Sum<ij|kl>^2",& "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" else 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