program extract_amplitudes
  implicit none
  BEGIN_DOC
! Print the T1 and T2 amplitudes into a file.
  END_DOC

  read_wf = .True.
  TOUCH read_wf
  call run
end

subroutine run
  implicit none
  BEGIN_DOC
! Compute t1 and T2 amplitudes
  END_DOC
  double precision, allocatable  :: t1_a(:,:), t1_b(:,:)
  double precision, allocatable  :: t2_aa(:,:,:,:)
  double precision, allocatable  :: t2_ab(:,:,:,:)
  double precision, allocatable  :: t2_bb(:,:,:,:)
  double precision               :: phase, norm
  integer                        :: exc(0:2,2,2), degree
  integer                        :: h1,h2,p1,p2,s1,s2, istate
  integer                        :: i,j,k,l,a,b,c,d

  istate = 1

  allocate (                                                         &
      t1_a(elec_alpha_num,mo_num),                                   &
      t1_b(elec_alpha_num,mo_num),                                   &
      t2_aa(elec_alpha_num,elec_alpha_num,mo_num,mo_num),            &
      t2_ab(elec_alpha_num,elec_alpha_num,mo_num,mo_num),            &
      t2_bb(elec_alpha_num,elec_alpha_num,mo_num,mo_num) )
  
  t1_a  = 0.d0
  t1_b  = 0.d0
  t2_aa = 0.d0
  t2_ab = 0.d0
  t2_bb = 0.d0

  norm = 1.d0 / psi_coef_sorted(1,istate)

  do k=1,N_det

    call get_excitation_degree(hf_bitmask,psi_det(1,1,k),degree,N_int)

    select case (degree)

      case (1)
        call get_excitation(hf_bitmask,psi_det(1,1,k),exc,degree,phase,N_int)
        call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)

        if (h1 > elec_alpha_num) then
          print *, irp_here, "h1 > elec_alpha_num"
        endif
        if (p1 <= elec_alpha_num) then
          print *, irp_here, "p1 <= elec_alpha_num"
        endif

        if (s1 == 1) then
          t1_a(h1,p1) += phase * psi_coef(k,istate) * norm
        else if (s1 == 2) then
          t1_b(h1,p1) += phase * psi_coef(k,istate) * norm
        else
           print *,  irp_here, ': Bug!', s1, s2
        endif
        

      case (2)
        call get_excitation(hf_bitmask,psi_det(1,1,k),exc,degree,phase,N_int)
        call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)

        if (h1 > elec_alpha_num) then
          print *, irp_here, "h1 > elec_alpha_num"
        endif
        if (p1 <= elec_alpha_num) then
          print *, irp_here, "p1 <= elec_alpha_num"
        endif
        if (h2 > elec_alpha_num) then
          print *, irp_here, "h2 > elec_alpha_num"
        endif
        if (p2 <= elec_alpha_num) then
          print *, irp_here, "p2 <= elec_alpha_num"
        endif

        if ( (s1 == 1).and.(s2 == 1) ) then
          t2_aa(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm * 0.5d0
          t2_aa(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm * 0.5d0
          t2_aa(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0
          t2_aa(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm * 0.5d0
        else if ( (s1 == 1).and.(s2 == 2) ) then
          t2_ab(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm * 0.5d0
          t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0
        else if ( (s1 == 2).and.(s2 == 1) ) then
          print *,  irp_here, ': Bug!'
          stop -1
        else if ( (s1 == 2).and.(s2 == 2) ) then
          t2_bb(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm * 0.5d0
          t2_bb(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm * 0.5d0
          t2_bb(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0
          t2_bb(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm * 0.5d0
        else
           print *,  irp_here, ': Bug!'
           stop -2
        endif
        
    end select
    
  enddo

  do b=elec_alpha_num+1, mo_num
    do a=elec_alpha_num+1, mo_num
      do j=1,elec_alpha_num
        do i=1,elec_alpha_num
          if ( (i == j).or.(a == b)) then
            t2_aa(i,j,a,b) = 0.d0
            t2_bb(i,j,a,b) = 0.d0
          else
            t2_aa(i,j,a,b) -= t1_a(i,a)*t1_a(j,b)
            t2_bb(i,j,a,b) -= t1_b(i,a)*t1_b(j,b)
          endif
          t2_ab(i,j,a,b) -= 0.5d0* (t1_a(i,a)*t1_b(j,b) + t1_b(i,a)*t1_a(j,b))
        enddo
      enddo
    enddo
  enddo


  integer :: iunit, getunitandopen

  iunit = getunitandopen('t1','w')

  write(iunit,*) '# T1_a'
  do i=1,elec_alpha_num
    do a=elec_alpha_num+1, mo_num
      if (t1_a(i,a) /= 0.d0) then
        write(iunit,'(I4,X,I4,X,E20.10)') i, a, t1_a(i,a)
      endif
    enddo
  enddo

  write(iunit,*) '# T1_b'
  do i=1,elec_alpha_num
    do a=elec_alpha_num+1, mo_num
      if (t1_b(i,a) /= 0.d0) then
        write(iunit,'(I4,X,I4,X,E20.10)') i, a, t1_b(i,a)
      endif
    enddo
  enddo


  iunit = getunitandopen('t2','w')

  write(iunit,*) '# T2_aa'
  do i=1,elec_alpha_num
    do j=1,elec_alpha_num
      do a=elec_alpha_num+1, mo_num
        do b=elec_alpha_num+1, mo_num
          if (t2_aa(i,j,a,b) /= 0.d0) then
            write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_aa(i,j,a,b)
          endif
        enddo
      enddo
    enddo
  enddo


  write(iunit,*) '# T2_bb'
  do i=1,elec_alpha_num
    do j=1,elec_alpha_num
      do a=elec_alpha_num+1, mo_num
        do b=elec_alpha_num+1, mo_num
          if (t2_bb(i,j,a,b) /= 0.d0) then
            write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_bb(i,j,a,b)
          endif
        enddo
      enddo
    enddo
  enddo


  write(iunit,*) '# T2_ab'
  do i=1,elec_alpha_num
    do j=1,elec_alpha_num
      do a=elec_alpha_num+1, mo_num
        do b=elec_alpha_num+1, mo_num
          if (t2_ab(i,j,a,b) /= 0.d0) then
            write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_ab(i,j,a,b)
          endif
        enddo
      enddo
    enddo
  enddo


  double precision               :: e_cor
  double precision, external     :: get_two_e_integral
  
  e_cor = 0.d0
  do b=elec_alpha_num+1, mo_num
    do a=elec_alpha_num+1, mo_num
      do j=1,elec_alpha_num
        do i=1,elec_alpha_num
          ! The t1 terms are zero because HF
          e_cor = e_cor + 0.5d0*(                                    &
              t2_aa(i,j,a,b) + t2_bb(i,j,a,b) +                      &
              t1_a(i,a) * t1_a(j,b)           +                      &
              t1_b(i,a) * t1_b(j,b) ) * (                            &
              get_two_e_integral(i,j,a,b,mo_integrals_map) -         &
              get_two_e_integral(i,j,b,a,mo_integrals_map) )
          e_cor = e_cor + 1.0d0 * (                                  &
              t2_ab(i,j,a,b) +                                       &
              t1_a(i,a) * t1_b(j,b) ) *                              &
              get_two_e_integral(i,j,a,b,mo_integrals_map)
        enddo
      enddo
    enddo
  enddo
  e_cor = e_cor
  print '(A,F15.10)',  'HF    energy: ', hf_energy
  print '(A,F15.10)',  'corr  energy: ', e_cor
  print '(A,F15.10)',  'total energy: ', hf_energy + e_cor

  deallocate(t1_a,t1_b,t2_aa,t2_ab,t2_bb)

end