CalcDerivConstr Subroutine

    subroutine CalcDerivConstr(CurrCoeff, DerivConstr)

        ! This calculates the derivative of each of the orthonormalisation
        ! constraints, l, with respect to each set of coefficients cm.

        integer :: l, i, j, a
        HElement_t(dp) :: CurrCoeff(NoOrbs, NoOrbs)
        real(dp) :: DerivConstr(NoOrbs, NoOrbs, TotNoConstraints)

        call set_timer(CalcDerivConstr_Time, 30)

        DerivConstr(:, :, :) = 0.0_dp
        do l = 1, TotNoConstraints
            i = lab(1, l)
            j = lab(2, l)
            if (i == j) then
                do a = 1, NoOrbs
                    DerivConstr(a, i, l) = CurrCoeff(a, i) * 2
                end do
            else
                do a = 1, NoOrbs
                    DerivConstr(a, j, l) = CurrCoeff(a, i)
                end do
                do a = 1, NoOrbs
                    DerivConstr(a, i, l) = CurrCoeff(a, j)
                end do
            end if
            ! DerivConstrT1 stays the same throughout the iterations
        end do

        call halt_timer(CalcDerivConstr_Time)

    end subroutine CalcDerivConstr