generate_random_integrals Subroutine

    subroutine generate_random_integrals(iunit, n_el, n_spat_orb, sparse, sparseT, &
            total_ms, uhf, hermitian)
        integer, intent(in) :: iunit, n_el, n_spat_orb
        real(dp), intent(in) :: sparse, sparseT
        type(SpinProj_t), intent(in) :: total_ms
        logical, optional, intent(in) :: uhf, hermitian
            !! specify if the FCIDUMP is UHF
            !! specify if the FCIDUMP is hermitian
            !!
            !! @note if uhf then assume
            !! num spin-up basis functions = num spin-down basis functions
        logical :: uhf_, hermitian_
        integer :: i, j, k, l, j_end, l_end
        real(dp) :: r, matel
        ! we get random matrix elements from the cauchy-schwartz inequalities, so
        ! only <ij|ij> are random -> random 2d matrix
        real(dp), allocatable :: umatRand(:, :)
        integer :: norb ! n_spin_orb or n_spat_orb

        def_default(hermitian_, hermitian, .true.)
        ! UHF FCIDUMPs have six sectors:
        ! two-body integrals: up-up, down-down, up-down
        ! one-body integrals: up, down
        ! nuclear repulsion
        ! delimiter: 0.0000000000000000E+00   0   0   0   0
        def_default(uhf_, uhf, .false.)

        norb = merge(2*n_spat_orb, n_spat_orb, uhf_)
        allocate(umatRand(norb, norb), source=0.0_dp)

        do i = 1, norb
            do j = 1, norb
                r = genrand_real2_dSFMT()
                if (r < sparse) then
                    umatRand(i, j) = r * r
                    if (hermitian_) then
                        umatRand(j, i) = umatRand(i, j)
                    else
                        ! have the conjugate be similar
                        umatRand(j, i) = (1 + genrand_real2_dSFMT() * 0.1) * umatRand(i, j)
                    endif
                end if
            end do
        end do

        ! write the canonical FCIDUMP header (norb here is num spatial orbitals)
        write(iunit, *) "&FCI NORB=", n_spat_orb, ",NELEC=", n_el, "MS2=", total_ms%val, ","
        write(iunit, "(A)", advance="no") "ORBSYM="

        do i = 1, n_spat_orb
            write(iunit, "(A)", advance="no") "1,"
        end do
        write(iunit, *)
        write(iunit, *) "ISYM=1,"
        write(iunit, *) "&END"
        ! generate random 4-index integrals with a given sparsity
        ! then
        ! generate random 2-index integrals with a given sparsity
        if(uhf_) then
            ! three 4-index sectors, so three calls to write_4index
            call write_4index(1, n_spat_orb, 1, n_spat_orb, hermitian_)
            write(iunit, *) 0.0000000000000000E+00, 0, 0, 0, 0 ! delimiter
            ! we can keep using the spatial orbitals as indices because of the
            ! partitioning of the dump file
            call write_4index(1, n_spat_orb, 1, n_spat_orb, hermitian_)
            write(iunit, *) 0.0000000000000000E+00, 0, 0, 0, 0
            call write_4index(1, n_spat_orb, 1, n_spat_orb, hermitian_)
            write(iunit, *) 0.0000000000000000E+00, 0, 0, 0, 0
            call write_2index(1, n_spat_orb, hermitian_)
            write(iunit, *) 0.0000000000000000E+00, 0, 0, 0, 0
            call write_2index(1, n_spat_orb, hermitian_)
            write(iunit, *) 0.0000000000000000E+00, 0, 0, 0, 0
            ! then would come the nuclear repulsion
        else ! .not. uhf_
            call write_4index(1, norb, 1, norb, hermitian_)
            call write_2index(1, n_spat_orb, hermitian_)
        endif

    contains
        subroutine write_4index(i_start, i_end, k_start, k_end, hermitian)
            ! i_end corresponds to electron 1, j_end to electron 2
            integer, intent(in) :: i_start, i_end, k_start, k_end
            logical, intent(in) :: hermitian
            integer :: k_start_
            do i = i_start, i_end
                ! j_end = i if hermitian, else i_end
                j_end = merge(i, i_end, hermitian)
                do j = 1, j_end
                    ! if the Hamiltonian *is* Hermitian, we may flip these indices
                    k_start_ = merge(i, k_start, hermitian)
                    do k = k_start_, k_end
                        l_end = merge(k, k_end, hermitian)
                        do l = 1, l_end
                            matel = sqrt(umatRand(i, j) * umatRand(k, l))
                            if (matel > eps) write(iunit, *) matel, i, j, k, l
                        end do
                    end do
                end do
            end do
        end subroutine write_4index

        subroutine write_2index(i_start, i_end, hermitian)
            integer, intent(in) :: i_start, i_end
            logical, intent(in) :: hermitian
            do i = i_start, i_end
                j_end = merge(i, i_end, hermitian)
                do j = i_start, j_end
                    r = genrand_real2_dSFMT()
                    if (r < sparseT) then
                        write(iunit, *) r, i, j, 0, 0
                    end if
                end do
            end do
        end subroutine write_2index

    end subroutine generate_random_integrals