WriteTransformMat Subroutine

    subroutine WriteTransformMat()

        integer :: w, x, i, a, b, iunit

! This file is printed to be used to produce cube files from QChem.
! Line 1 is the coefficients of HF spatial orbitals 1 2 3 ... which form transformed orbital 1 etc.

        iunit = get_free_unit()
        open(iunit, file='MOTRANSFORM', FORM='UNFORMATTED', access='direct', recl=8)
        ! Need to put this back into the original order.

        x = 0
        if (tStoreSpinOrbs) then
            do i = 1, NoOrbs - 1, 2
                ! SymLabelList2_rot(i) gives the orbital label (from Dalton or QChem) corresponding to our
                ! label i.
                ! SymLabelListInv_rot(j) therefore gives the label used in CoeffT1 corresponding to the
                ! Qchem/Dalton label j.

                do a = 1, NoOrbs - 1, 2
                    b = SymLabelListInv_rot(a)
                    write(iunit, rec=x) CoeffT1(b, i)
                    ! a/b are the original (HF) orbitals, and i/j the transformed
                end do
            end do
            do i = 2, NoOrbs, 2
                do a = 2, NoOrbs, 2
                    b = SymLabelListInv_rot(a)
                    write(iunit, rec=x) CoeffT1(b, i)
                    ! a/b are the original (HF) orbitals, and i/j the transformed
                end do
            end do
        else
            w = 1
            x = 1   !keep a counter of record number
            do while (w <= 2)
                do i = 1, SpatOrbs
                    ! SymLabelList2_rot(i) gives the orbital label (from Dalton or QChem) corresponding to our
                    ! label i.
                    ! SymLabelListInv_rot(j) therefore gives the label used in CoeffT1 corresponding to the
                    ! Qchem/Dalton label j.
                    do a = 1, SpatOrbs
                        b = SymLabelListInv_rot(a)
                        write(iunit, rec=x) CoeffT1(b, i)
                        x = x + 1
                        ! a/b are the original (HF) orbitals, and i/j the transformed
                    end do
                end do
                w = w + 1
                ! print the whole matrix twice, once for alpha spin, once for beta.
            end do
        end if
        close(iunit)

        open(iunit, file='MOTRANSFORM02')
        ! Need to put this back into the original order.
        w = 1
        x = 1   !keep a counter of record number
        do while (w <= 2)
            do i = 1, SpatOrbs
                ! SymLabelList2_rot(i) gives the orbital label (from Dalton or QChem) corresponding to our
                ! label i.
                ! SymLabelListInv_rot(j) therefore gives the label used in CoeffT1 corresponding to the
                ! Qchem/Dalton label j.
                do a = 1, SpatOrbs
                    b = SymLabelListInv_rot(a)
                    write(iunit, '(F20.10)', advance='no') CoeffT1(b, i)
                    x = x + 1
                    ! a/b are the original (HF) orbitals, and i/j the transformed
                end do
                write(iunit, *) ''
            end do
            w = w + 1
            ! print the whole matrix twice, once for alpha spin, once for beta.
        end do
        close(iunit)

        open(iunit, file='TRANSFORMMAT', status='unknown')
        do i = 1, NoOrbs
            do a = 1, NoOrbs
                b = SymLabelListInv_rot(a)
                write(iunit, *) CoeffT1(b, i)
            end do
        end do
        call neci_flush(iunit)
        close(iunit)

    end subroutine WriteTransformMat