# InitRotCalc Subroutine

None

## Source Code

    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