print_results Subroutine

  subroutine print_results(infos, bvec_mo, excitation_energy, &
                           trans, dip, spin_square, nstates)
    use precision, only: dp
    use types, only: information
    use physical_constants, only: toev => ev2htree

    implicit none

    type(information), target, intent(in) :: infos
    real(kind=dp), intent(in), dimension(:,:) :: bvec_mo
    real(kind=dp), intent(in), dimension(:) :: excitation_energy
    integer, intent(in), dimension(:,:) :: trans
    real(kind=dp), intent(in), dimension(:,:,:) :: dip
    real(kind=dp), intent(in), dimension(:) :: spin_square
    integer, intent(in) :: nstates

    integer :: istat, jstat, ij, i, j, nocca, noccb, xvec_dim, ndeex
    real(kind=dp) :: ydum, xdum, threshold, ROHF_energy, energ, f

    threshold = infos%control%conf_print_threshold
    xvec_dim = ubound(bvec_mo, 1)

    ROHF_energy = infos%mol_energy%energy
    nocca = infos%mol_prop%nelec_A
    noccb = infos%mol_prop%nelec_B

    do istat=1,nstates
      ydum = toev*excitation_energy(istat)
      write(*,'(/,1x,"State #",I4,2X,"Energy =",F12.6,1X,"eV")') istat, ydum
!     write(*,'(3x,"Symmetry of state =",4x,a)') '?a?a?'
      write(*,'(15x,"<S^2> =",1x,f9.4)') spin_square(istat)
      write(*,'(8x,"DRF",4x,"Coeff",8x,"OCC",7x,"VIR")')
      write(*,'(8x,3("-"),2x,8("-"),5x,6("-"),4x,6("-"))')
      do ij=1,xvec_dim
        i = trans(ij,1)
        j = trans(ij,2)
        xdum = bvec_mo(ij,istat)
        if (abs(xdum)>threshold) then
          write(*,'(7x,i4,1x,f9.6,6x,i4,2x,"->",2x,i4,2x)') ij,xdum,i,j
        end if
      end do
    end do

    write(*,'(/5x, "Summary table",/)')
    write(*,'(1x, "State", 6x, "Energy", 7x,"Excitation", 3x, "Excitation(eV)", &
             &2x, "<S^2>", 9x, "Transition dipole moment, a.u.",&
             &8x, "Oscillator")')
    write(*,'(11x, "Hartree", 11x, "eV", 10x, "rel. GS" &
             &18x, "X", 10x, "Y", 10x, "Z", 8x,"Abs.", 6x, "strength")')

    ndeex = 0
    do istat = 1, nstates
      if(excitation_energy(istat)<0.0_dp) ndeex = ndeex+1
    end do

    ! De-excitation
    do istat = 1, ndeex
      energ = excitation_energy(istat)-excitation_energy(1)
      f = 2.0d0 / 3.0d0 * (energ) * sum(dip(:,1,istat)**2)
      write(*,'(x, i3, 1x, f17.10, 2f13.6, 6x, &
               &f5.3, 4(1x,f10.4),2x,f10.4)') &
           istat, ROHF_energy+excitation_energy(istat), toev*excitation_energy(istat), &
           toev*energ, spin_square(istat), dip(1:3,1,istat), sqrt(sum(dip(:,1,istat)**2)), f
    end do

    ! Reference ROHF state
    write(*,'(1x, i3, 1x, f17.10, 2f13.6, 8x,&
            &"(ROHF/UHF Reference state)")') 0, ROHF_energy, 0.0_dp, -excitation_energy(1)*toev

    ! Excitation
    do istat=ndeex+1,nstates
      energ = excitation_energy(istat)-excitation_energy(1)
      f = 2.0d0 / 3.0d0 * (energ) * sum(dip(:,1,istat)**2)
      write(*,'(x, i3, 1x, f17.10, 2f13.6, 6x, &
               &f5.3, 4(1x,f10.4),2x,f10.4)') &
           istat, ROHF_energy+excitation_energy(istat), toev*excitation_energy(istat), &
           toev*energ, spin_square(istat), dip(1:3,1,istat), sqrt(sum(dip(:,1,istat)**2)), f
    end do
    write(*,*)
    write(*,"(2x,'Transition',3x,'Excitation',9x,'Transition dipole, a.u.',19x,'Oscillator',&
          &/18x,'eV',14x,'x',10x,'y',10x,'z',9x,'Abs.',7x,'strength')")
    do istat=1, nstates
       do jstat=istat+1, nstates
          energ = excitation_energy(jstat)-excitation_energy(istat)
          f = 2.0d0 / 3.0d0 * (energ) * sum(dip(:,istat,jstat)**2)
    write(*,"(3x,i0,1x,'->',1x,i0,t11,3x,f11.6,3x,3f11.4,1x,f11.4,2x,f11.4)") &
             istat,jstat,toev*energ,dip(1:3,istat,jstat), sqrt(sum(dip(:,istat,jstat)**2)), f
       enddo
    enddo
    write(*,*)

  end subroutine print_results