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