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,EcCCSD
  double precision              :: EcCCT
  double precision              :: get_two_e_integral,u_dot_v 

  double precision,allocatable  :: cFvv(:,:)
  double precision,allocatable  :: cFoo(:,:)
  double precision,allocatable  :: cFov(:,:)

  double precision,allocatable  :: cWoooo(:,:,:,:)
  double precision,allocatable  :: cWvvvv(:,:,:,:)
  double precision,allocatable  :: cWovvo(:,:,:,:)

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

  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


! Guess amplitudes

  allocate(tau(spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num))

! Initialization

  allocate(cFvv(spin_vir_num,spin_vir_num),cFoo(spin_occ_num,spin_occ_num),cFov(spin_occ_num,spin_vir_num), &
           cWoooo(spin_occ_num,spin_occ_num,spin_occ_num,spin_occ_num),cWvvvv(spin_vir_num,spin_vir_num,spin_vir_num,spin_vir_num),cWovvo(spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num), &
           r1(spin_occ_num,spin_vir_num),r2(spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num))

  Conv = 1d0
  nSCF = 0

  call form_tau_nc (tau)

  EcMP2 = 0.25d0*u_dot_v(OOVV,tau,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)

!   Excrement 

    nSCF = nSCF + 1

    call form_cf_nc (                                       &
        spin_fock_matrix_mo_oo,                                      &
        spin_fock_matrix_mo_ov,                                      &
        spin_fock_matrix_mo_vv,                                      &
        cFoo,cFov,cFvv)

    call form_cw_nc (tau,                                &
        cWoooo,cWovvo,cWvvvv)

!   Compute residuals

    call form_r1_nc(spin_fock_matrix_mo_ov, &
       cFoo,cFov,cFvv,r1)

    call form_r2_nc(tau,cFoo,cFov,cFvv, &
       cWoooo,cWvvvv,cWovvo,r2)

!   Check convergence 

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

!   Update 

    t1_cc(:,:)      = t1_cc(:,:)     - r1(:,:)    /delta_OV  (:,:)
    t2_cc(:,:,:,:)  = t2_cc(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
    TOUCH t1_cc t2_cc

    call form_tau_nc (tau)

!   Compute correlation energy

    call CCSD_Ec_nc(spin_fock_matrix_mo_ov,EcCCSD)

!   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(                                                        &
      cFvv,cFoo,cFov,                                                &
      cWoooo,cWvvvv,cWovvo,                                          &
      tau,                                                      &
      r1,r2)

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