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                       :: nSCF
  double precision              :: Conv
  double precision              :: EcMP2
  double precision              :: ECCSD
  double precision              :: EcCCT
  double precision              :: get_two_e_integral,u_dot_v 

! 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


! Initialization

  Conv = 1d0
  nSCF = 0

  EcMP2 = 0.25d0*u_dot_v(OOVV,tau_cc,size(OOVV))
  write(*,'(1X,A10,1X,F10.6)') 'Ec(MP2) = ',EcMP2
  write(*,'(1X,A10,1X,F10.6)') 'E (MP2) = ',EcMP2 + ERHF

!------------------------------------------------------------------------
! 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)

    nSCF = nSCF + 1

!   Check convergence 

    Conv = max(maxval(abs(r1_cc(:,:))),maxval(abs(r2_cc(:,:,:,:))))

!   Update 

    t1_cc(:,:)      = t1_cc(:,:)     - r1_cc(:,:)    /delta_OV  (:,:)
    t2_cc(:,:,:,:)  = t2_cc(:,:,:,:) - r2_cc(:,:,:,:)/delta_OOVV(:,:,:,:)
    TOUCH t1_cc t2_cc

!   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


!------------------------------------------------------------------------
! (T) correction
!------------------------------------------------------------------------
  if(doCCSDT) then
    write(*,*) "Starting (T) calculation"
    call CCSDT(EcCCT)
    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