```
subroutine InitRotCalc()
! Sets up the initial arrays to be used in the orbital rotation.
character(len=*), parameter :: this_routine = 'InitRotCalc'
integer :: i, j, Const, MinRot, MaxRot
call InitSymmArrays()
! Creates an indexing system for each of the cases with symmetry on/off, and mixing all orbitals or separating
! the occuppied from virtual.
! The arrays used in this routine are labelled with a 2 (SymLabelList2_rot and SymLabelCount2), so as to not
! mess up the spawing/FCI calcs.
do i = 1, NoOrbs
SymLabelList3_rot(i) = SymLabelList2_rot(i)
end do
! Set up constraint labels. Constraint l is the dot product of i.j.
Const = 0
do i = 1, NoOrbs
do j = i, NoOrbs
Const = Const + 1
Lab(1, Const) = i
Lab(2, Const) = j
end do
end do
! Just a check that the number of constraints labeled is the same as
! that calculated above.
write(stdout, *) 'Total number of constraints = ', TotNoConstraints
if (Const /= TotNoConstraints) then
call Stop_all(this_routine, 'ERROR in the number of constraints calculated. lmax does not equal TotNoConstraints')
end if
! Zero/initialise the arrays
! In the case where symmetry is kept, the starting transformation matrix is just the identity. Starting with a symmetric system
! means the symmetry is never broken.
! When we are breaking the symmetry, the starting transformation matrix is completely random, (and then orthonormalised).
! The ordering of the orbitals in CoeffT1 follow the ordering in SymLabelList2_rot.
CoeffT1(:, :) = 0.0_dp
if (tRotateOccOnly) then
MinRot = 1
MaxRot = NoOcc
else if (tRotateVirtOnly) then
MinRot = NoOcc + 1
MaxRot = NoOrbs
else
MinRot = 1
MaxRot = NoOrbs
end if
do i = 1, NoOrbs
CoeffT1(i, i) = 1.0_dp
end do
! If the symmetry is kept on, start with the symmetric identity matrix
! of coefficients, and it will be maintained.
! When only the occupied or virtual orbitals are rotated, the non
! rotated orbitals are still included in the transformation matrix,
! but start as the identity, and remain that way throughout.
if (lNoSymmetry) then
do i = MinRot, MaxRot
do j = MinRot, MaxRot
CoeffT1(j, i) = genrand_real2_dSFMT() * (1E-02_dp)
end do
end do
end if
! Ensures transformation matrix elements between the occupied and
! virtual orbitals are 0 (should be the case anyway though).
if (tSeparateOccVirt) call ZeroOccVirtElements(CoeffT1)
! Orthonormalise starting matrix.
call GRAMSCHMIDT_NECI(CoeffT1, NoOrbs)
! A UMATTemp is created (from the UMAT read in from the FCIDUMP) using the rotate orbs indexing system.
! i.e. in UMatTemp(1,1,1,1), the orbital involved is the first in SymLabelList2_rot.
! Doing this now, rather than using UMatInd in each transform2elint routine proved a lot faster.
DerivCoeff(:, :) = 0.0_dp
UMATTemp01(:, :, :, :) = 0.0_dp
if (((.not. tERLocalization) .and. (.not. tReadInCoeff) .and. (.not. tUseMP2VarDenMat) &
.and. (.not. tFindCINatOrbs) .and. (.not. tUseHFOrbs))&
&.or. (tERLocalization .and. tStoreSpinOrbs)) UMATTemp02(:, :, :, :) = 0.0_dp
call CopyAcrossUMAT()
call TestOrthonormality()
! With UMAT with the correct indexing and the starting coefficient,
! find the partially transformed four index integrals (and hence the
! initial potential energy), and then the initial force.
if (tERLocalization .and. (.not. tStoreSpinOrbs)) then
call Transform2ElIntsERlocal()
else
call Transform2ElInts()
end if
call FindTheForce()
if (tPrintInts) then
! This sets the initial values for the integral sums being printed.
! Values printed are then relative to these initial sums, per
! integral.
tInitIntValues = .true.
call PrintIntegrals()
tInitIntValues = .false.
end if
end subroutine InitRotCalc
```

NECI was developed by Alavi Group

© 2023

Documentation generated by FORD