subroutine CCSD

! CCSD module

  implicit none

! Input variables (provided by IRP)

  integer                       :: maxscf
  double precision              :: thresh

  logical                       :: doCCSDT
  integer                       :: nBas,nEl
  double precision              :: ERHF

! Local variables

  integer                       :: p,q,r,s
  double precision              :: start_CCSDT,end_CCSDT,t_CCSDT
  integer                       :: nBas2
  integer                       :: nO
  integer                       :: nV
  integer                       :: nSCF
  double precision              :: Conv
  double precision              :: EcMP2
  double precision              :: ECCSD,EcCCSD
  double precision              :: EcCCT
  double precision              :: get_two_e_integral,u_dot_v 

  double precision,allocatable  :: hvv(:,:)
  double precision,allocatable  :: hoo(:,:)
  double precision,allocatable  :: hvo(:,:)
  double precision,allocatable  :: gvv(:,:)
  double precision,allocatable  :: goo(:,:)
  double precision,allocatable  :: aoooo(:,:,:,:)
  double precision,allocatable  :: bvvvv(:,:,:,:)
  double precision,allocatable  :: hovvo(:,:,:,:)

  double precision,allocatable  :: r1(:,:)
  double precision,allocatable  :: r2(:,:,:,:)

  double precision,allocatable  :: t1(:,:)
  double precision,allocatable  :: t2(:,:,:,:)
  double precision,allocatable  :: tau(:,:,:,:)

! Hello world

  write(*,*)
  write(*,*)'**************************************'
  write(*,*)'|         CCSD calculation           |'
  write(*,*)'**************************************'
  write(*,*)

! IRP init  
  provide cc_mode
  maxSCF=cc_n_it_max
  thresh=cc_thresh

  doCCSDT = trim(cc_mode)=='CCSD(T)'

  nBas=mo_num
  nEl=elec_num
  ERHF=hf_energy

! Spatial to spin orbitals

  nBas2 = spin_mo_num


! Define occupied and virtual spaces

  nO = spin_occ_num
  nV = spin_vir_num

! Guess amplitudes

  allocate(t1(nO,nV),t2(nO,nO,nV,nV),tau(nO,nO,nV,nV))

  t1(:,:)     = t1_guess(:,:)
  t2(:,:,:,:) = t2_guess(:,:,:,:)

  call form_tau(nO,nV,t1,t2,tau)

  EcMP2 = 0.5d0*u_dot_v(pack(OOVV,.true.),pack(tau,.true.),size(OOVV))
  write(*,'(1X,A10,1X,F10.6)') 'Ec(MP2) = ',EcMP2

! Initialization

  allocate(hvv(nV,nV),hoo(nO,nO),hvo(nV,nO), &
           gvv(nV,nV),goo(nO,nO), & 
           aoooo(nO,nO,nO,nO),bvvvv(nV,nV,nV,nV),hovvo(nO,nV,nV,nO), &
           r1(nO,nV),r2(nO,nO,nV,nV))

  Conv = 1d0
  nSCF = 0

!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
  write(*,*)
  write(*,*)'----------------------------------------------------'
  write(*,*)'| CCSD calculation                                  |'
  write(*,*)'----------------------------------------------------'
  write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
            '|','#','|','E(CCSD)','|','Ec(CCSD)','|','Conv','|'
  write(*,*)'----------------------------------------------------'

  do while(Conv > thresh .and. nSCF < maxSCF)

!   Increment 

    nSCF = nSCF + 1

!   Scuseria Eqs. (5), (6) and (7)

    call form_h(nO,nV,t1,tau,hvv,hoo,hvo)

!   Scuseria Eqs. (9), (10), (11), (12) and (13)

    call form_g(nO,nV,hvv,hoo,t1,gvv,goo)

    call form_abh(nO,nV,t1,tau,aoooo,bvvvv,hovvo)

!   Compute residuals

    call form_r1(nO,nV,hvv,hoo,hvo,t1,t2,tau,r1)

    call form_r2(nO,nV,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)

!   Check convergence 

    Conv = max(maxval(abs(r1(:,:))),maxval(abs(r2(:,:,:,:))))

!   Update 

    t1(:,:)      = t1(:,:)     - r1(:,:)    /delta_OV  (:,:)
    t2(:,:,:,:)  = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)

    call form_tau(nO,nV,t1,t2,tau)
 
!   Compute correlation energy

    EcCCSD = 0.5d0*u_dot_v(pack(OOVV,.true.),pack(tau,.true.),size(OOVV))

!   Dump results

    ECCSD = ERHF + EcCCSD

    write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
      '|',nSCF,'|',ECCSD,'|',EcCCSD,'|',Conv,'|'

  end do
  write(*,*)'----------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------

! Did it actually converge?

  if(nSCF == maxSCF) then

    write(*,*)
    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
    write(*,*)'                 Convergence failed                 '
    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
    write(*,*)

    stop

  end if

! Deallocate memory

  deallocate(hvv,hoo,hvo,         &
             gvv,goo,             &
             aoooo,bvvvv,hovvo,   &
             tau,                 &
             r1,r2)

!------------------------------------------------------------------------
! (T) correction
!------------------------------------------------------------------------
  if(doCCSDT) then
    write(*,*) "Starting (T) calculation"
!   call cpu_time(start_CCSDT)V
    call CCSDT(nO,nV,t1,t2,EcCCT)
!   call cpu_time(end_CCSDT)
    call write_time(6)

!    t_CCSDT = end_CCSDT - start_CCSDT
!    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for (T) = ',t_CCSDT,' seconds'
     write(*,*)

    write(*,*)
    write(*,*)'----------------------------------------------------'
    write(*,*)'                 CCSDT(T) energy                    '
    write(*,*)'----------------------------------------------------'
    write(*,'(1X,A20,1X,F15.10)')' E(CCSD(T))  = ',ECCSD  + EcCCT
    write(*,'(1X,A20,1X,F10.6)') ' Ec(CCSD(T)) = ',EcCCSD + EcCCT
    write(*,*)'----------------------------------------------------'
    write(*,*)

    call save_energy(ECCSD  + EcCCT)

  else

    call save_energy(ECCSD)

  end if

end subroutine CCSD