Working on CC

devel/cc/CC.irp.f

@@ -0,0 +1,11 @@
+program CC
+  implicit none
+  BEGIN_DOC
+  ! CC
+  ! Coupled cluster
+  ! 
+  ! This program perform coupled cluster calculation
+  ! CCSD or CCSD(T)
+  END_DOC
+  call CCSD
+end

devel/cc/CCSD.irp.f

@@ -0,0 +1,211 @@
+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

devel/cc/CCSDT.irp.f

@@ -0,0 +1,35 @@
+subroutine CCSDT(nO,nV,t1,t2,EcCCT)
+
+! Compute the (T) correction of the CCSD(T) energy
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+
+  double precision,intent(in)   :: t1(nO,nV)
+  double precision,intent(in)   :: t2(nO,nO,nV,nV)
+
+! Local variables
+
+  double precision,allocatable  :: ub(:,:,:,:,:,:)
+  double precision,allocatable  :: ubb(:,:,:,:,:,:)
+
+! Output variables
+
+  double precision,intent(out)  :: EcCCT
+
+! Memory allocation
+
+  allocate(ub(nO,nO,nO,nV,nV,nV),ubb(nO,nO,nO,nV,nV,nV))
+
+! Form CCSD(T) quantities
+
+  call form_ub(nO,nV,t1,ub)
+
+  call form_ubb(nO,nV,t2,ubb)
+
+  call form_T(nO,nV,ub,ubb,EcCCT)
+
+end subroutine CCSDT

devel/cc/EZFIO.cfg

@@ -0,0 +1,28 @@
+[energy]
+type: Threshold
+doc: Energy CC 
+interface: ezfio
+
+[cc_guess]
+type: integer
+doc: 0: MP2 amplitudes, 1: read from files
+interface: ezfio,provider,ocaml
+default: 0
+
+[cc_mode]
+type: character*(32)
+doc: mode CCD CCSD CCSD(T).
+interface: ezfio,provider,ocaml
+default:CCSD
+
+[cc_thresh]
+type: Threshold
+doc: Threshold on the convergence of the Coupled-Cluster residual. 
+interface: ezfio,provider,ocaml
+default: 1.e-10
+
+[cc_n_it_max]
+type: Strictly_positive_int
+doc: Maximum number of Coupled-Cluster iterations.
+interface: ezfio,provider,ocaml
+default: 500

devel/cc/Kronecker_delta.irp.f

@@ -0,0 +1,22 @@
+!------------------------------------------------------------------------
+function Kronecker_delta(i,j) result(delta)
+
+! Kronecker Delta
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: i,j
+
+! Output variables
+
+  double precision              :: delta
+
+  if(i == j) then
+    delta = 1d0
+  else
+    delta = 0d0
+  endif
+
+end function Kronecker_delta

devel/cc/MP2.irp.f

@@ -0,0 +1,71 @@
+subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
+
+! Perform third-order Moller-Plesset calculation
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas,nC,nO,nV,nR
+  double precision,intent(in)   :: ENuc,EHF
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas),e(nBas)
+
+! Local variables
+
+  integer                       :: i,j,a,b
+  double precision              :: eps,E2a,E2b
+
+! Output variables
+
+  double precision,intent(out)  :: EcMP2(3)
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'|    Moller-Plesset second-order calculation   |'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Compute MP2 energy
+
+  E2a = 0d0
+  E2b = 0d0
+  do i=nC+1,nO
+  do j=nC+1,nO
+    do a=nO+1,nBas-nR
+    do b=nO+1,nBas-nR
+
+      eps = e(i) + e(j) - e(a) - e(b) 
+
+!     Second-order ring diagram
+
+      E2a = E2a + ERI(i,j,a,b)*ERI(i,j,a,b)/eps
+
+!     Second-order exchange diagram
+
+      E2b = E2b + ERI(i,j,a,b)*ERI(i,j,b,a)/eps
+
+    enddo
+    enddo
+  enddo
+  enddo
+
+  EcMP2(2) = 2d0*E2a
+  EcMP2(3) = -E2b
+  EcMP2(1) = EcMP2(2) + EcMP2(3)
+
+  write(*,*)
+  write(*,'(A32)')           '-----------------------'
+  write(*,'(A32)')           ' MP2 calculation       '
+  write(*,'(A32)')           '-----------------------'
+  write(*,'(A32,1X,F16.10)') ' MP2 correlation energy',EcMP2(1)
+  write(*,'(A32,1X,F16.10)') ' Direct part           ',EcMP2(2)
+  write(*,'(A32,1X,F16.10)') ' Exchange part         ',EcMP2(3)
+  write(*,'(A32)')           '-----------------------'
+  write(*,'(A32,1X,F16.10)') ' MP2 electronic  energy',EHF + EcMP2(1)
+  write(*,'(A32,1X,F16.10)') ' MP2 total       energy',ENuc + EHF + EcMP2(1)
+  write(*,'(A32)')           '-----------------------'
+  write(*,*)
+
+end subroutine MP2

devel/cc/NEED

@@ -0,0 +1,2 @@
+hartree_fock
+mo_two_e_ints

devel/cc/README.rst

@@ -0,0 +1,4 @@
+==
+CC
+==
+

devel/cc/amplitude_guess.irp.f

@@ -0,0 +1,84 @@
+BEGIN_PROVIDER [ double precision, t1_guess, (spin_occ_num,spin_vir_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Guess amplitudes for single excitations
+ END_DOC
+
+ t1_guess(:,:)     = 0d0
+ if (cc_guess == 1) then
+    integer :: iunit
+    integer, external :: getunitandopen
+
+    character :: check
+    integer :: i, a
+    double precision :: amplitude
+
+    iunit = getunitandopen('t1','r')
+    read(iunit,*)
+    do
+      read(iunit,*,err=10) i, a, amplitude
+      i = 2*i-1
+      a = 2*a-1 - spin_occ_num
+      t1_guess(i,a) = amplitude
+    enddo
+    10 continue
+    do
+      read(iunit,*,end=20) i, a, amplitude
+      i = 2*i
+      a = 2*a - spin_occ_num
+      t1_guess(i,a) = amplitude
+    enddo
+    20 continue
+    close(iunit)
+ endif
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Guess amplitudes for double excitations
+ END_DOC
+
+ t2_guess(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
+
+ if (cc_guess == 1) then
+    integer :: iunit
+    integer, external :: getunitandopen
+    character :: check
+    integer :: i, j, a, b
+    double precision :: amplitude
+
+    iunit = getunitandopen('t2','r')
+    read(iunit,*) 
+    do
+      read(iunit,*,err=10) i, j, a, b, amplitude
+      i = 2*i-1
+      j = 2*j-1
+      a = 2*a-1 - spin_occ_num
+      b = 2*b-1 - spin_occ_num
+      t2_guess(i,j,a,b) = amplitude
+    enddo
+    10 continue
+    do
+      read(iunit,*,err=20) i, j, a, b, amplitude
+      i = 2*i
+      j = 2*j  
+      a = 2*a   - spin_occ_num
+      b = 2*b   - spin_occ_num
+      t2_guess(i,j,a,b) = amplitude
+    enddo
+    20 continue
+    do
+      read(iunit,*,end=30) i, j, a, b, amplitude
+      i = 2*i-1
+      j = 2*j
+      a = 2*a-1 - spin_occ_num
+      b = 2*b   - spin_occ_num
+      t2_guess(i,j,a,b) = amplitude
+    enddo
+    30 continue
+    close(iunit)
+ endif
+
+END_PROVIDER

devel/cc/antisymmetrize_ERI.irp.f

@@ -0,0 +1,46 @@
+subroutine antisymmetrize_ERI(ispin,nBas,ERI,db_ERI)
+
+! Antisymmetrize ERIs
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: ispin,nBas
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: i,j,k,l
+
+! Output variables
+
+  double precision,intent(out)  :: db_ERI(nBas,nBas,nBas,nBas)
+
+  if(ispin == 1) then
+
+    do l=1,nBas
+      do k=1,nBas
+        do j=1,nBas
+          do i=1,nBas
+            db_ERI(i,j,k,l) = 2d0*ERI(i,j,k,l) - ERI(i,j,l,k)
+          enddo
+        enddo
+      enddo
+    enddo
+
+  elseif(ispin == 2) then
+
+    do l=1,nBas
+      do k=1,nBas
+        do j=1,nBas
+          do i=1,nBas
+            db_ERI(i,j,k,l) = ERI(i,j,k,l) - ERI(i,j,l,k)
+          enddo
+        enddo
+      enddo
+    enddo
+
+  endif
+
+end subroutine antisymmetrize_ERI

devel/cc/denominators.irp.f

@@ -0,0 +1,70 @@
+ BEGIN_PROVIDER [ double precision, eO, (1:spin_occ_num) ]
+&BEGIN_PROVIDER [ double precision, eV, (1:spin_vir_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Eigenvalues of the spin-orbitals
+ END_DOC
+ eO(:) = spin_fock_matrix_diag_mo(1:spin_occ_num)
+ eV(:) = spin_fock_matrix_diag_mo(spin_occ_num+1:spin_mo_num)
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, delta_OV, (spin_occ_num, spin_vir_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Energy denominator for CC
+ END_DOC
+ integer :: i, a
+
+ do a=1,spin_vir_num
+   do i=1,spin_occ_num
+     delta_OV(i,a) = eV(a) - eO(i)
+   enddo
+ enddo
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, delta_OOVV, (spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Energy denominator for CC
+ END_DOC
+ integer :: i,j,a,b
+
+ do b=1,spin_vir_num
+  do a=1,spin_vir_num
+   do j=1,spin_occ_num
+    do i=1,spin_occ_num
+     delta_OOVV(i,j,a,b) = eV(a) - eO(i) + eV(b) - eO(j)
+    enddo
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, delta_OOOVVV, (spin_occ_num,spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num,spin_vir_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Energy denominator for CC
+ END_DOC
+ integer :: i,j,k,a,b,c
+
+ do c=1,spin_vir_num
+  do b=1,spin_vir_num
+   do a=1,spin_vir_num
+    do k=1,spin_occ_num
+     do j=1,spin_occ_num
+      do i=1,spin_occ_num
+       delta_OOOVVV(i,j,k,a,b,c) = eV(a) - eO(i) + eV(b) - eO(j) + eV(c) - eO(k)
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+
+

devel/cc/dimensions.irp.f

@@ -0,0 +1,20 @@
+BEGIN_PROVIDER [ integer, spin_mo_num ]
+ implicit none
+ BEGIN_DOC
+ ! Number of spin-orbitals
+ END_DOC
+ spin_mo_num = mo_num * 2
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer, spin_occ_num ]
+&BEGIN_PROVIDER [ integer, spin_vir_num ]
+ implicit none
+ BEGIN_DOC
+ ! Number of occupied/virtual spin-orbitals
+ END_DOC
+ spin_occ_num = elec_num
+ spin_vir_num = spin_mo_num - spin_occ_num
+END_PROVIDER
+
+

devel/cc/form_T.irp.f

@@ -0,0 +1,45 @@
+subroutine form_T(nO,nV,ub,ubb,EcCCT)
+
+! Compute (T) correction
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+
+  double precision,intent(in)   :: ub(nO,nO,nO,nV,nV,nV)
+  double precision,intent(in)   :: ubb(nO,nO,nO,nV,nV,nV)
+
+! Local variables
+
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+
+! Output variables
+
+  double precision,intent(out)  :: EcCCT
+
+  EcCCT = 0d0
+
+  do c=1,nV
+    do b=1,nV
+      do a=1,nV
+        do k=1,nO
+          do j=1,nO
+            do i=1,nO
+
+              EcCCT = EcCCT                                &
+                    + (ub(i,j,k,a,b,c) + ubb(i,j,k,a,b,c)) &
+                    * ubb(i,j,k,a,b,c)/delta_OOOVVV(i,j,k,a,b,c) 
+
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  EcCCT = - EcCCT/36d0
+
+end subroutine form_T

devel/cc/form_abh.irp.f

@@ -0,0 +1,97 @@
+subroutine form_abh(nO,nV,t1,tau,aoooo,bvvvv,hovvo)
+
+! Scuseria Eqs. (11),(12) and (13)
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+
+  double precision,intent(in)   :: t1(nO,nV)
+  double precision,intent(in)   :: tau(nO,nO,nV,nV)
+
+! Local variables
+
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+
+! Output variables
+
+  double precision,intent(out)  :: aoooo(nO,nO,nO,nO)
+  double precision,intent(out)  :: bvvvv(nV,nV,nV,nV)
+  double precision,intent(out)  :: hovvo(nO,nV,nV,nO)
+
+  aoooo(:,:,:,:) = OOOO(:,:,:,:)
+
+  do l=1,nO
+    do k=1,nO
+      do j=1,nO
+        do i=1,nO
+
+          do c=1,nV
+            aoooo(i,j,k,l) = aoooo(i,j,k,l) + OVOO(i,c,k,l)*t1(j,c)
+          end do
+
+          do c=1,nV
+            aoooo(i,j,k,l) = aoooo(i,j,k,l) - OVOO(j,c,k,l)*t1(i,c)
+          end do
+
+          do d=1,nV
+            do c=1,nV
+              aoooo(i,j,k,l) = aoooo(i,j,k,l) + OOVV(k,l,c,d)*tau(i,j,c,d)
+            end do
+          end do
+
+        end do
+      end do
+    end do
+  end do
+
+  bvvvv(:,:,:,:) = VVVV(:,:,:,:)
+
+  do b=1,nV 
+    do a=1,nV
+      do d=1,nV
+        do c=1,nV
+          
+          do k=1,nO
+            bvvvv(c,d,a,b) = bvvvv(c,d,a,b) - VOVV(a,k,c,d)*t1(k,b)
+          end do
+
+          do k=1,nO
+            bvvvv(c,d,a,b) = bvvvv(c,d,a,b) + VOVV(b,k,c,d)*t1(k,a)
+          end do
+
+        end do
+      end do
+    end do
+  end do
+
+  hovvo(:,:,:,:) = OVVO(:,:,:,:)
+
+  do k=1,nO 
+    do a=1,nV
+      do c=1,nV
+        do i=1,nO
+
+          do l=1,nO
+            hovvo(i,c,a,k) = hovvo(i,c,a,k) - OVOO(i,c,l,k)*t1(l,a)
+          end do
+
+          do d=1,nV
+            hovvo(i,c,a,k) = hovvo(i,c,a,k) + OVVV(k,a,c,d)*t1(i,d)
+          end do
+
+          do d=1,nV
+            do l=1,nO
+              hovvo(i,c,a,k) = hovvo(i,c,a,k) - OOVV(k,l,c,d)*tau(i,l,d,a)
+            end do
+          end do
+
+        end do
+      end do
+    end do
+  end do
+
+end subroutine form_abh

devel/cc/form_g.irp.f

@@ -0,0 +1,50 @@
+subroutine form_g(nO,nV,hvv,hoo,t1,gvv,goo)
+
+! Scuseria Eqs. (9), (10)
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+
+  double precision,intent(in)   :: hvv(nV,nV)
+  double precision,intent(in)   :: hoo(nO,nO)
+
+  double precision,intent(in)   :: t1(nO,nV)
+
+! Local variables
+
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+
+! Output variables
+
+  double precision,intent(out)  :: gvv(nV,nV)
+  double precision,intent(out)  :: goo(nO,nO)
+
+  gvv(:,:) = hvv(:,:)
+
+  do a=1,nV 
+    do c=1,nV
+      do d=1,nV
+        do k=1,nO
+          gvv(c,a) = gvv(c,a) + VOVV(a,k,c,d)*t1(k,d)
+        end do
+      end do
+    end do
+  end do
+
+  goo(:,:) = hoo(:,:)
+
+  do k=1,nO
+    do i=1,nO
+      do c=1,nV
+        do l=1,nO
+          goo(i,k) = goo(i,k) + OOOV(k,l,i,c)*t1(l,c)
+        end do
+      end do
+    end do
+  end do
+
+end subroutine form_g

devel/cc/form_h.irp.f

@@ -0,0 +1,75 @@
+subroutine form_h(nO,nV,t1,tau,hvv,hoo,hvo)
+
+! Scuseria Eqs. (5), (6) and (7)
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+
+  double precision,intent(in)   :: t1(nO,nV)
+  double precision,intent(in)   :: tau(nO,nO,nV,nV)
+
+! Local variables
+
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+
+! Output variables
+
+  double precision,intent(out)  :: hvv(nV,nV)
+  double precision,intent(out)  :: hoo(nO,nO)
+  double precision,intent(out)  :: hvo(nV,nO)
+
+  hvv(:,:) = 0d0
+
+  do b=1,nV 
+    hvv(b,b) = eV(b)
+    do a=1,nV
+      do j=1,nO
+        do k=1,nO
+          do c=1,nV
+
+            hvv(b,a) = hvv(b,a) - OOVV(j,k,b,c)*tau(j,k,a,c)
+
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  hoo(:,:) = 0d0
+
+  do i=1,nO
+    hoo(i,i) = eO(i)
+    do j=1,nO
+      do k=1,nO
+        do b=1,nV
+          do c=1,nV
+
+            hoo(i,j) = hoo(i,j) + OOVV(j,k,b,c)*tau(i,k,b,c)
+
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  hvo(:,:) = 0d0
+
+  do j=1,nO 
+    do b=1,nV
+      do c=1,nV
+        do k=1,nO
+
+          hvo(b,j) = hvo(b,j) + OOVV(j,k,b,c)*t1(k,c)
+
+        end do
+      end do
+    end do
+  end do
+
+! print*,'hvv',hvv
+
+end subroutine form_h

devel/cc/form_r1.irp.f

@@ -0,0 +1,72 @@
+subroutine form_r1(nO,nV,hvv,hoo,hvo,t1,t2,tau,r1)
+
+! Form tau in CCSD
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+
+  double precision,intent(in)   :: hvv(nV,nV)
+  double precision,intent(in)   :: hoo(nO,nO)
+  double precision,intent(in)   :: hvo(nV,nO)
+
+  double precision,intent(in)   :: t1(nO,nV)
+  double precision,intent(in)   :: t2(nO,nO,nV,nV)
+  double precision,intent(in)   :: tau(nO,nO,nV,nV)
+
+! Local variables
+
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+
+! Output variables
+
+  double precision,intent(out)  :: r1(nO,nV)
+
+  r1(:,:) = 0d0
+
+  do a=1,nV
+      do i=1,nO
+
+        do b=1,nV
+          r1(i,a) = r1(i,a) + hvv(b,a)*t1(i,b)
+        end do
+
+        do j=1,nO
+          r1(i,a) = r1(i,a) - hoo(i,j)*t1(j,a)
+        end do
+
+        do j=1,nO
+          do b=1,nV
+            r1(i,a) = r1(i,a) + hvo(b,j)*(t2(i,j,a,b) + t1(i,b)*t1(j,a))
+          end do
+        end do
+
+        do b=1,nV
+          do j=1,nO
+            r1(i,a) = r1(i,a) + OVVO(i,b,a,j)*t1(j,b)
+          end do
+        end do
+
+        do c=1,nV
+          do b=1,nV
+            do j=1,nO
+              r1(i,a) = r1(i,a) - OVVV(j,a,b,c)*tau(i,j,b,c)
+            end do
+          end do
+        end do
+
+        do b=1,nV
+          do k=1,nO
+            do j=1,nO
+              r1(i,a) = r1(i,a) - OOOV(j,k,i,b)*tau(j,k,a,b)
+            end do
+          end do
+        end do
+
+    end do
+  end do
+
+end subroutine form_r1

devel/cc/form_r2.irp.f

@@ -0,0 +1,134 @@
+subroutine form_r2(nO,nV,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
+
+! Form tau in CCSD
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+
+  double precision,intent(in)   :: gvv(nV,nV)
+  double precision,intent(in)   :: goo(nO,nO)
+  double precision,intent(in)   :: aoooo(nO,nO,nO,nO)
+  double precision,intent(in)   :: bvvvv(nV,nV,nV,nV)
+  double precision,intent(in)   :: hovvo(nO,nV,nV,nO)
+
+  double precision,intent(in)   :: t1(nO,nV)
+  double precision,intent(in)   :: t2(nO,nO,nV,nV)
+  double precision,intent(in)   :: tau(nO,nO,nV,nV)
+
+! Local variables
+
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+
+! Output variables
+
+  double precision,intent(out)  :: r2(nO,nO,nV,nV)
+
+  r2(:,:,:,:) = OOVV(:,:,:,:)
+
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
+        do i=1,nO
+
+          do k=1,nO
+            do l=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) + aoooo(i,j,k,l)*tau(k,l,a,b)
+            end do
+          end do
+
+          do d=1,nV
+            do c=1,nV
+              r2(i,j,a,b) = r2(i,j,a,b) + bvvvv(c,d,a,b)*tau(i,j,c,d)
+            end do
+          end do
+
+          do c=1,nV
+            r2(i,j,a,b) = r2(i,j,a,b) + gvv(c,a)*t2(i,j,c,b)
+          end do
+
+          do k=1,nO
+            r2(i,j,a,b) = r2(i,j,a,b) + OVOO(k,a,i,j)*t1(k,b)
+          end do
+
+          do c=1,nV
+            r2(i,j,a,b) = r2(i,j,a,b) - gvv(c,b)*t2(i,j,c,a)
+          end do
+
+          do k=1,nO
+            r2(i,j,a,b) = r2(i,j,a,b) - OVOO(k,b,i,j)*t1(k,a)
+          end do
+
+          do k=1,nO
+            r2(i,j,a,b) = r2(i,j,a,b) - goo(i,k)*t2(k,j,a,b)
+          end do
+
+          do c=1,nV
+            r2(i,j,a,b) = r2(i,j,a,b) + OVVV(j,c,b,a)*t1(i,c)
+          end do
+
+          do k=1,nO
+            r2(i,j,a,b) = r2(i,j,a,b) + goo(j,k)*t2(k,i,a,b)
+          end do
+
+          do c=1,nV
+            r2(i,j,a,b) = r2(i,j,a,b) - OVVV(i,c,b,a)*t1(j,c)
+          end do
+
+          do k=1,nO
+            do c=1,nV
+              r2(i,j,a,b) = r2(i,j,a,b) + hovvo(i,c,a,k)*t2(j,k,b,c)
+            end do
+          end do
+
+          do k=1,nO
+            do c=1,nV
+              r2(i,j,a,b) = r2(i,j,a,b) - OVVO(i,c,a,k)*t1(j,c)*t1(k,b)
+            end do
+          end do
+
+          do k=1,nO
+            do c=1,nV
+              r2(i,j,a,b) = r2(i,j,a,b) - hovvo(j,c,a,k)*t2(i,k,b,c)
+            end do
+          end do
+
+          do c=1,nV
+            do k=1,nO
+               r2(i,j,a,b) = r2(i,j,a,b) + OVVO(j,c,a,k)*t1(i,c)*t1(k,b)
+            end do
+          end do
+
+          do c=1,nV
+            do k=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) - hovvo(i,c,b,k)*t2(j,k,a,c)
+            end do
+          end do
+
+          do c=1,nV
+            do k=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) + OVVO(i,c,b,k)*t1(j,c)*t1(k,a)
+            end do
+          end do
+
+          do c=1,nV
+            do k=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) + hovvo(j,c,b,k)*t2(i,k,a,c)
+            end do
+          end do
+
+          do c=1,nV
+            do k=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) - OVVO(j,c,b,k)*t1(i,c)*t1(k,a)
+            end do
+          end do
+
+        end do
+      end do
+    end do
+  end do
+
+end subroutine form_r2

devel/cc/form_tau.irp.f

@@ -0,0 +1,34 @@
+subroutine form_tau(nO,nV,t1,t2,tau)
+
+! Form tau in CCSD
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+  double precision,intent(in)   :: t1(nO,nV)
+  double precision,intent(in)   :: t2(nO,nO,nV,nV)
+
+! Local variables
+
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+
+! Output variables
+
+  double precision,intent(out)  :: tau(nO,nO,nV,nV)
+
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
+        do i=1,nO
+
+          tau(i,j,a,b) = 0.5d0*t2(i,j,a,b) + t1(i,a)*t1(j,b)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+end subroutine form_tau

devel/cc/form_ub.irp.f

@@ -0,0 +1,46 @@
+subroutine form_ub(nO,nV,t1,ub)
+
+! Form 1st term in (T) correction
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+
+  double precision,intent(in)   :: t1(nO,nV)
+
+! Local variables
+
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+
+! Output variables
+
+  double precision,intent(out)  :: ub(nO,nO,nO,nV,nV,nV)
+
+  do c=1,nV
+    do b=1,nV
+      do a=1,nV
+        do k=1,nO
+          do j=1,nO
+            do i=1,nO
+
+              ub(i,j,k,a,b,c) = t1(i,a)*OOVV(j,k,b,c) &
+                              + t1(i,b)*OOVV(j,k,c,a) &
+                              + t1(i,c)*OOVV(j,k,a,b) &
+                              + t1(j,a)*OOVV(k,i,b,c) & 
+                              + t1(j,b)*OOVV(k,i,c,a) &
+                              + t1(j,c)*OOVV(k,i,a,b) & 
+                              + t1(k,a)*OOVV(i,j,b,c) &
+                              + t1(k,b)*OOVV(i,j,c,a) & 
+                              + t1(k,c)*OOVV(i,j,a,b)
+
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+end subroutine form_ub

devel/cc/form_ubb.irp.f

@@ -0,0 +1,64 @@
+subroutine form_ubb(nO,nV,t2,ubb)
+
+! Form 2nd term in (T) correction
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nV
+
+  double precision,intent(in)   :: t2(nO,nO,nV,nV)
+
+! Local variables
+
+  integer                       :: i,j,k,l,m
+  integer                       :: a,b,c,d,e
+
+! Output variables
+
+  double precision,intent(out)  :: ubb(nO,nO,nO,nV,nV,nV)
+
+  ubb(:,:,:,:,:,:) = 0d0
+
+  do c=1,nV
+    do b=1,nV
+      do a=1,nV
+        do k=1,nO
+          do j=1,nO
+            do i=1,nO
+
+              do e=1,nV
+                ubb(i,j,k,a,b,c) = ubb(i,j,k,a,b,c)          &
+                                 + t2(i,j,a,e)*VVVO(b,c,e,k) & 
+                                 + t2(i,j,b,e)*VVVO(c,a,e,k) & 
+                                 + t2(i,j,c,e)*VVVO(a,b,e,k) & 
+                                 + t2(k,i,a,e)*VVVO(b,c,e,j) & 
+                                 + t2(k,i,b,e)*VVVO(c,a,e,j) & 
+                                 + t2(k,i,c,e)*VVVO(a,b,e,j) & 
+                                 + t2(j,k,a,e)*VVVO(b,c,e,i) & 
+                                 + t2(j,k,b,e)*VVVO(c,a,e,i) & 
+                                 + t2(j,k,c,e)*VVVO(a,b,e,i)
+              end do
+ 
+              do m=1,nO
+                ubb(i,j,k,a,b,c) = ubb(i,j,k,a,b,c)          &
+                                 + t2(i,m,a,b)*VOOO(c,m,j,k) & 
+                                 + t2(i,m,b,c)*VOOO(a,m,j,k) & 
+                                 + t2(i,m,c,a)*VOOO(b,m,j,k) & 
+                                 + t2(j,m,a,b)*VOOO(c,m,k,i) & 
+                                 + t2(j,m,b,c)*VOOO(a,m,k,i) & 
+                                 + t2(j,m,c,a)*VOOO(b,m,k,i) & 
+                                 + t2(k,m,a,b)*VOOO(c,m,i,j) & 
+                                 + t2(k,m,b,c)*VOOO(a,m,i,j) & 
+                                 + t2(k,m,c,a)*VOOO(b,m,i,j)
+              end do
+              
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+end subroutine form_ubb

devel/cc/integral_batches.irp.f

@@ -0,0 +1,79 @@
+BEGIN_PROVIDER [ double precision, OOOO, (spin_occ_num,spin_occ_num,spin_occ_num,spin_occ_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  OOOO(:,:,:,:) = dbERI(   1:spin_occ_num   ,   1:spin_occ_num   ,   1:spin_occ_num   ,   1:spin_occ_num   )
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, OOOV, (spin_occ_num,spin_occ_num,spin_occ_num,spin_vir_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  OOOV(:,:,:,:) = dbERI(   1:spin_occ_num   ,   1:spin_occ_num   ,   1:spin_occ_num   ,spin_occ_num+1:spin_mo_num)
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, OVOO, (spin_occ_num,spin_vir_num,spin_occ_num,spin_occ_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  OVOO(:,:,:,:) = dbERI(   1:spin_occ_num   ,spin_occ_num+1:spin_mo_num,   1:spin_occ_num   ,   1:spin_occ_num   )
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, VOOO, (spin_vir_num,spin_occ_num,spin_occ_num,spin_occ_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  VOOO(:,:,:,:) = dbERI(spin_occ_num+1:spin_mo_num,   1:spin_occ_num   ,   1:spin_occ_num   ,   1:spin_occ_num   )
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, OOVV, (spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  OOVV(:,:,:,:) = dbERI(   1:spin_occ_num   ,   1:spin_occ_num   ,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num)
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, OVVO, (spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  OVVO(:,:,:,:) = dbERI(   1:spin_occ_num   ,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,   1:spin_occ_num   )
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, OVVV, (spin_occ_num,spin_vir_num,spin_vir_num,spin_vir_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  OVVV(:,:,:,:) = dbERI(   1:spin_occ_num   ,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num)
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, VOVV, (spin_vir_num,spin_occ_num,spin_vir_num,spin_vir_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  VOVV(:,:,:,:) = dbERI(spin_occ_num+1:spin_mo_num,   1:spin_occ_num   ,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num)
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, VVVO, (spin_vir_num,spin_vir_num,spin_vir_num,spin_occ_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  VVVO(:,:,:,:) = dbERI(spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,   1:spin_occ_num   )
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, VVVV, (spin_vir_num,spin_vir_num,spin_vir_num,spin_vir_num) ]
+  implicit none
+  BEGIN_DOC
+  END_DOC
+  VVVV(:,:,:,:) = dbERI(spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num)
+END_PROVIDER
+

devel/cc/save_energy.irp.f

@@ -0,0 +1,8 @@
+subroutine save_energy(E)
+  implicit none
+  BEGIN_DOC
+! Saves the energy in |EZFIO|.
+  END_DOC
+  double precision, intent(in) :: E
+  call ezfio_set_cc_energy(E)
+end

devel/cc/spatial_to_spin_ERI.irp.f

@@ -0,0 +1,52 @@
+BEGIN_PROVIDER [ double precision, dbERI, (spin_mo_num,spin_mo_num,spin_mo_num,spin_mo_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Anti-symmetrized Electron repulsion integrals in spin-orbital basis
+ END_DOC
+
+! Local variables
+
+  integer                       :: p,q,r,s
+  double precision,external     :: Kronecker_delta
+
+! Output variables
+
+  double precision, external :: get_two_e_integral
+
+  double precision :: pqrs, pqsr
+
+  PROVIDE mo_two_e_integrals_in_map
+
+    do s=1,spin_mo_num
+      do r=1,spin_mo_num
+        do q=1,spin_mo_num
+          do p=1,spin_mo_num
+
+            pqrs = Kronecker_delta(mod(p,2),mod(r,2)) &
+                 * Kronecker_delta(mod(q,2),mod(s,2)) &
+                 * get_two_e_integral(                &
+                     (p+1)/2,                         &
+                     (q+1)/2,                         &
+                     (r+1)/2,                         &
+                     (s+1)/2,                         &
+                     mo_two_e_integrals_in_map)
+
+            pqsr = Kronecker_delta(mod(p,2),mod(s,2)) &
+                 * Kronecker_delta(mod(q,2),mod(r,2)) &
+                 * get_two_e_integral(                &
+                     (p+1)/2,                         &
+                     (q+1)/2,                         &
+                     (s+1)/2,                         &
+                     (r+1)/2,                         &
+                     mo_two_e_integrals_in_map)
+
+            dbERI(p,q,r,s) = pqrs - pqsr
+          enddo
+        enddo
+      enddo
+    enddo
+
+
+END_PROVIDER
+
+

devel/cc/spatial_to_spin_MO_energy.irp.f

@@ -0,0 +1,13 @@
+BEGIN_PROVIDER [ double precision, spin_fock_matrix_diag_mo, (spin_mo_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Diagonal of the Fock matrix in the spin-orbital basis
+ END_DOC
+ integer :: p
+
+ do p=1,spin_mo_num
+    spin_fock_matrix_diag_mo(p) = fock_matrix_diag_mo((p+1)/2)
+ enddo
+
+END_PROVIDER
+

stable/utilities/extract_amplitudes.irp.f

@@ -0,0 +1,212 @@
+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
+        else if ( (s1 == 1).and.(s2 == 2) ) then
+          t2_ab(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm
+        else if ( (s1 == 2).and.(s2 == 1) ) then
+          t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm
+        else if ( (s1 == 2).and.(s2 == 2) ) then
+          t2_bb(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm
+        else
+           print *,  irp_here, ': Bug!'
+        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
+          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)
+          t2_ab(i,j,a,b) -= t1_a(i,a)*t1_b(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 + (                                                 &
+              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 + (                                                 &
+              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