InitRotCalc Subroutine

    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