Merge pull request #263 from Ydrnan/dev-stable-cc

utils_cc, ccsd, ccsd(t) + tests

src/ccsd/80.ccsd_spin.bats

@@ -0,0 +1,225 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
+
+
+function run() {
+  thresh1=1e-6
+  thresh2=1e-6
+  test_exe scf || skip
+  qp set_file $1
+  qp edit --check
+  #qp run scf
+  qp set_frozen_core
+  qp set utils_cc cc_par_t true
+  qp set utils_cc cc_thresh_conv 1e-12
+  file="$(echo $1 | sed 's/.ezfio//g')"
+  qp run ccsd_spin_orb | tee $file.ccsd.out
+  energy1="$(grep 'E(CCSD)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
+  energy2="$(grep 'E(T)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
+  #rm $file.ccsd.out
+  eq $energy1 $2 $thresh1
+  eq $energy2 $3 $thresh2
+}
+
+@test "b2_stretched" {
+run  b2_stretched.ezfio -49.136487344382 -0.003497589175
+}
+
+@test "be" {
+run  be.ezfio -14.623559003577 -0.000230982022
+}
+
+@test "c2h2" {
+run  c2h2.ezfio -12.394008897618 -0.010790491561
+}
+
+@test "ch4" {
+run  ch4.ezfio -40.390721785799 -0.004476100282
+}
+
+@test "clf" {
+run  clf.ezfio -559.186562904081 -0.006577143392
+}
+
+@test "clo" {
+run  clo.ezfio -534.564874409332 -0.007584571424
+}
+
+@test "co2" {
+run  co2.ezfio -188.129602527766 -0.018040668885
+}
+
+@test "dhno" {
+run  dhno.ezfio -130.816650109473 -0.012197331453
+}
+
+@test "f2" {
+run  f2.ezfio -199.287826338097 -0.017592872692
+}
+
+@test "f" {
+run  f.ezfio -99.616644511121 -0.003624525307
+}
+
+@test "h2o2" {
+run  h2o2.ezfio -151.182552729963 -0.009511682086
+}
+
+@test "h2o" {
+run  h2o.ezfio -76.237710276526 -0.003001800577
+}
+
+@test "h2s" {
+run  h2s.ezfio -398.861214015390 -0.003300559757
+}
+
+@test "h3coh" {
+run  h3coh.ezfio -115.221296424969 -0.003566171432
+}
+
+@test "hbo" {
+run  hbo.ezfio -100.213539770415 -0.006851489212
+}
+
+@test "hcn" {
+run  hcn.ezfio -93.190247992657 -0.013418135043
+}
+
+@test "hco" {
+run  hco.ezfio -113.405413962350 -0.007973455337
+}
+
+@test "lif" {
+run  lif.ezfio -107.270402903250 -0.007742969005
+}
+
+@test "n2" {
+run  n2.ezfio -109.355358930472 -0.018477744342
+}
+
+@test "n2h4" {
+run  n2h4.ezfio -111.556885923139 -0.009048077008
+}
+
+@test "nh3" {
+run  nh3.ezfio -56.465503060954 -0.007638273755
+}
+
+@test "oh" {
+run  oh.ezfio -75.614606132774 -0.004126661739
+}
+
+@test "sih2_3b1" {
+run  sih2_3b1.ezfio -290.016780973072 -0.000497825874
+}
+
+@test "sih3" {
+run  sih3.ezfio -5.575343504534 -0.002094123268
+}
+
+@test "so" {
+run  so.ezfio -26.035945178665 -0.010594351274
+}
+
+#@test "b2_stretched" {
+#run  b2_stretched.ezfio -49.136487344382 -49.139984933557
+#}
+#
+#@test "be" {
+#run  be.ezfio -14.623559003577 -14.623789985599
+#}
+#
+#@test "c2h2" {
+#run  c2h2.ezfio -12.394008897618 -12.404799389179
+#}
+#
+#@test "ch4" {
+#run  ch4.ezfio -40.390721784961 -40.395197884406
+#}
+#
+#@test "clf" {
+#run  clf.ezfio -559.186562906072 -559.193140046904
+#}
+#
+#@test "clo" {
+#run  clo.ezfio -534.564874409333 -534.572458980757
+#}
+#
+#@test "co2" {
+#run  co2.ezfio -188.129602511724 -188.147643198675
+#}
+#
+#@test "dhno" {
+#run  dhno.ezfio -130.816650109473 -130.828847440925
+#}
+#
+#@test "f2" {
+#run  f2.ezfio -199.287826338097 -199.305419210789
+#}
+#
+#@test "f" {
+#run  f.ezfio -99.616644511120 -99.620269036428
+#}
+#
+#@test "h2o2" {
+#run  h2o2.ezfio -151.182552729963 -151.192064412049
+#}
+#
+#@test "h2o" {
+#run  h2o.ezfio -76.237710276526 -76.240712077103
+#}
+#
+#@test "h2s" {
+#run  h2s.ezfio -398.861214015416 -398.864514575146
+#}
+#
+#@test "h3coh" {
+#run  h3coh.ezfio -115.221296424969 -115.224862596401
+#}
+#
+#@test "hbo" {
+#run  hbo.ezfio -100.213539770415 -100.220391259627
+#}
+#
+#@test "hcn" {
+#run  hcn.ezfio -93.190247983000 -93.203666131216
+#}
+#
+#@test "hco" {
+#run  hco.ezfio -113.405413962350 -113.413387417687
+#}
+#
+#@test "lif" {
+#run  lif.ezfio -107.270402903211 -107.278145872216
+#}
+#
+#@test "n2" {
+#run  n2.ezfio -109.355358930472 -109.373836674814
+#}
+#
+#@test "n2h4" {
+#run  n2h4.ezfio -111.556885922642 -111.565934000556
+#}
+#
+#@test "nh3" {
+#run  nh3.ezfio -56.465503060954 -56.473141334709
+#}
+#
+#@test "oh" {
+#run  oh.ezfio -75.614606131897 -75.618732794235
+#}
+#
+#@test "sih2_3b1" {
+#run  sih2_3b1.ezfio -290.016780973071 -290.017278798946
+#}
+#
+#@test "sih3" {
+#run  sih3.ezfio -5.575343504534 -5.577437627802
+#}
+#
+#@test "so" {
+#run  so.ezfio -26.035945181998 -26.046539528491
+#}
+

src/ccsd/81.ccsd_space.bats

@@ -0,0 +1,225 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
+
+
+function run() {
+  thresh1=1e-6
+  thresh2=1e-6
+  test_exe scf || skip
+  qp set_file $1
+  qp edit --check
+  #qp run scf
+  qp set_frozen_core
+  qp set utils_cc cc_par_t true
+  qp set utils_cc cc_thresh_conv 1e-12
+  file="$(echo $1 | sed 's/.ezfio//g')"
+  qp run ccsd_space_orb | tee $file.ccsd.out
+  energy1="$(grep 'E(CCSD)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
+  energy2="$(grep 'E(T)' $file.ccsd.out | tail -n 1 | awk '{printf $3}')"
+  #rm $file.ccsd.out
+  eq $energy1 $2 $thresh1
+  eq $energy2 $3 $thresh2
+}
+
+@test "b2_stretched" {
+run  b2_stretched.ezfio -49.136487344382 -0.003497589175
+}
+
+@test "be" {
+run  be.ezfio -14.623559003577 -0.000230982022
+}
+
+@test "c2h2" {
+run  c2h2.ezfio -12.394008897618 -0.010790491561
+}
+
+@test "ch4" {
+run  ch4.ezfio -40.390721785799 -0.004476100282
+}
+
+@test "clf" {
+run  clf.ezfio -559.186562904081 -0.006577143392
+}
+
+#@test "clo" {
+#run  clo.ezfio -534.564874409332 -0.007584571424
+#}
+
+@test "co2" {
+run  co2.ezfio -188.129602527766 -0.018040668885
+}
+
+#@test "dhno" {
+#run  dhno.ezfio -130.816650109473 -0.012197331453
+#}
+
+@test "f2" {
+run  f2.ezfio -199.287826338097 -0.017592872692
+}
+
+#@test "f" {
+#run  f.ezfio -99.616644511121 -0.003624525307
+#}
+
+@test "h2o2" {
+run  h2o2.ezfio -151.182552729963 -0.009511682086
+}
+
+@test "h2o" {
+run  h2o.ezfio -76.237710276526 -0.003001800577
+}
+
+@test "h2s" {
+run  h2s.ezfio -398.861214015390 -0.003300559757
+}
+
+@test "h3coh" {
+run  h3coh.ezfio -115.221296424969 -0.003566171432
+}
+
+@test "hbo" {
+run  hbo.ezfio -100.213539770415 -0.006851489212
+}
+
+@test "hcn" {
+run  hcn.ezfio -93.190247992657 -0.013418135043
+}
+
+#@test "hco" {
+#run  hco.ezfio -113.405413962350 -0.007973455337
+#}
+
+@test "lif" {
+run  lif.ezfio -107.270402903250 -0.007742969005
+}
+
+@test "n2" {
+run  n2.ezfio -109.355358930472 -0.018477744342
+}
+
+@test "n2h4" {
+run  n2h4.ezfio -111.556885923139 -0.009048077008
+}
+
+@test "nh3" {
+run  nh3.ezfio -56.465503060954 -0.007638273755
+}
+
+#@test "oh" {
+#run  oh.ezfio -75.614606132774 -0.004126661739
+#}
+
+#@test "sih2_3b1" {
+#run  sih2_3b1.ezfio -290.016780973072 -0.000497825874
+#}
+
+#@test "sih3" {
+#run  sih3.ezfio -5.575343504534 -0.002094123268
+#}
+
+#@test "so" {
+#run  so.ezfio -26.035945178665 -0.010594351274
+#}
+
+#@test "b2_stretched" {
+#run  b2_stretched.ezfio -49.136487344382 -49.139984933557
+#}
+#
+#@test "be" {
+#run  be.ezfio -14.623559003577 -14.623789985599
+#}
+#
+#@test "c2h2" {
+#run  c2h2.ezfio -12.394008897618 -12.404799389179
+#}
+#
+#@test "ch4" {
+#run  ch4.ezfio -40.390721784961 -40.395197884406
+#}
+#
+#@test "clf" {
+#run  clf.ezfio -559.186562906072 -559.193140046904
+#}
+#
+##@test "clo" {
+##run  clo.ezfio -534.564874409333 -534.572458980757
+##}
+#
+#@test "co2" {
+#run  co2.ezfio -188.129602511724 -188.147643198675
+#}
+#
+##@test "dhno" {
+##run  dhno.ezfio -130.816650109473 -130.828847440925
+##}
+#
+#@test "f2" {
+#run  f2.ezfio -199.287826338097 -199.305419210789
+#}
+#
+##@test "f" {
+##run  f.ezfio -99.616644511120 -99.620269036428
+##}
+#
+#@test "h2o2" {
+#run  h2o2.ezfio -151.182552729963 -151.192064412049
+#}
+#
+#@test "h2o" {
+#run  h2o.ezfio -76.237710276526 -76.240712077103
+#}
+#
+#@test "h2s" {
+#run  h2s.ezfio -398.861214015416 -398.864514575146
+#}
+#
+#@test "h3coh" {
+#run  h3coh.ezfio -115.221296424969 -115.224862596401
+#}
+#
+#@test "hbo" {
+#run  hbo.ezfio -100.213539770415 -100.220391259627
+#}
+#
+#@test "hcn" {
+#run  hcn.ezfio -93.190247983000 -93.203666131216
+#}
+#
+##@test "hco" {
+##run  hco.ezfio -113.405413962350 -113.413387417687
+##}
+#
+#@test "lif" {
+#run  lif.ezfio -107.270402903211 -107.278145872216
+#}
+#
+#@test "n2" {
+#run  n2.ezfio -109.355358930472 -109.373836674814
+#}
+#
+#@test "n2h4" {
+#run  n2h4.ezfio -111.556885922642 -111.565934000556
+#}
+#
+#@test "nh3" {
+#run  nh3.ezfio -56.465503060954 -56.473141334709
+#}
+#
+##@test "oh" {
+##run  oh.ezfio -75.614606131897 -75.618732794235
+##}
+#
+##@test "sih2_3b1" {
+##run  sih2_3b1.ezfio -290.016780973071 -290.017278798946
+##}
+#
+##@test "sih3" {
+##run  sih3.ezfio -5.575343504534 -5.577437627802
+##}
+#
+##@test "so" {
+##run  so.ezfio -26.035945181998 -26.046539528491
+##}
+

src/ccsd/NEED

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

src/ccsd/README.md

@@ -0,0 +1,31 @@
+# CCSD in spin orbitals and spatial orbitals
+
+CCSD and CCSD(T) in spin orbitals for open and closed shell systems.
+CCSD and CCSD(T) in spatial orbitals for closed shell systems.
+
+## Calculations
+The program will automatically choose the version in spin or spatial orbitals
+To run the general program:  
+```
+qp run ccsd
+```
+Nevertheless, you can enforce the run in spin orbitals with
+```
+qp run ccsd_spin_orb
+```
+
+## Settings
+The settings can be changed with:
+```
+qp set utils_cc cc_#param #val
+```
+For more informations on the settings, look at the module utils_cc and its documentation. 
+
+## Org files
+The org files are stored in the directory org in order to avoid overwriting on user changes.
+The org files can be modified, to export the change to the source code, run
+```
+./TANGLE_org_mode.sh and
+mv *.irp.f ../.
+```
+

src/ccsd/ccsd.irp.f

@@ -0,0 +1,18 @@
+program ccsd
+
+  implicit none
+
+  BEGIN_DOC
+  ! CCSD program
+  END_DOC
+
+  read_wf = .True.
+  touch read_wf
+
+  if (.not. cc_ref_is_open_shell) then
+    call run_ccsd_space_orb
+  else 
+    call run_ccsd_spin_orb
+  endif
+
+end

src/ccsd/ccsd_space_orb.irp.f

@@ -0,0 +1,12 @@
+! Code
+
+program ccsd
+
+  implicit none
+
+  read_wf = .True.
+  touch read_wf
+
+  call run_ccsd_space_orb
+  
+end

src/ccsd/ccsd_spin_orb.irp.f

@@ -0,0 +1,16 @@
+! Prog
+
+program ccsd
+
+  implicit none
+
+  BEGIN_DOC
+  ! CCSD in spin orbitals
+  END_DOC
+
+  read_wf = .True.
+  touch read_wf
+
+  call run_ccsd_spin_orb
+  
+end

src/ccsd/ccsd_t_space_orb.irp.f

@@ -0,0 +1,412 @@
+! Dumb way
+
+subroutine ccsd_par_t_space(nO,nV,t1,t2,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO, nV)
+  double precision, intent(in)  :: t2(nO, nO, nV, nV)
+  double precision, intent(out) :: energy
+  
+  double precision, allocatable :: W(:,:,:,:,:,:)
+  double precision, allocatable :: V(:,:,:,:,:,:)
+  integer :: i,j,k,a,b,c
+  
+  allocate(W(nO,nO,nO,nV,nV,nV))
+  allocate(V(nO,nO,nO,nV,nV,nV))
+
+  call form_w(nO,nV,t2,W)   
+  call form_v(nO,nV,t1,W,V)
+
+  energy = 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
+              energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c))  !delta_ooovvv(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  
+  energy = energy / 3d0
+  
+  deallocate(V,W)
+end
+
+subroutine form_w(nO,nV,t2,W)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t2(nO, nO, nV, nV)
+  double precision, intent(out) :: W(nO, nO, nO, nV, nV, nV)
+  
+  integer :: i,j,k,l,a,b,c,d
+
+  W = 0d0
+  do c = 1, nV
+    print*,'W:',c,'/',nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+
+              do d = 1, nV
+                W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+                ! chem (bd|ai)
+                ! phys <ba|di>
+                + cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+                + cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+                + cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+                + cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+                + cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+                + cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+              enddo
+
+              do l = 1, nO
+                W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+                ! chem (ck|jl)
+                ! phys <cj|kl>
+                - cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
+                - cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
+                - cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
+                - cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
+                - cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
+                - cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
+              enddo
+
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end
+
+subroutine form_v(nO,nV,t1,w,v)
+
+implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO, nV)
+  double precision, intent(in)  :: W(nO, nO, nO, nV, nV, nV)
+  double precision, intent(out) :: V(nO, nO, nO, nV, nV, nV)
+
+  integer :: i,j,k,a,b,c
+
+  V = 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
+              V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+              + cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+              + cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+              + cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end
+
+! Main
+
+subroutine ccsd_par_t_space_v2(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO,nV), f_o(nO), f_v(nV)
+  double precision, intent(in)  :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO)
+  double precision, intent(out) :: energy
+  
+  double precision, allocatable :: W(:,:,:,:,:,:)
+  double precision, allocatable :: V(:,:,:,:,:,:)
+  double precision, allocatable :: W_ijk(:,:,:), V_ijk(:,:,:)
+  double precision, allocatable :: X_vvvo(:,:,:,:), X_ovoo(:,:,:,:), X_vvoo(:,:,:,:)
+  double precision, allocatable :: T_vvoo(:,:,:,:), T_ovvo(:,:,:,:), T_vo(:,:)
+  integer                       :: i,j,k,l,a,b,c,d
+  double precision              :: e,ta,tb, delta, delta_ijk
+ 
+  !allocate(W(nV,nV,nV,nO,nO,nO))
+  !allocate(V(nV,nV,nV,nO,nO,nO))
+  allocate(W_ijk(nV,nV,nV), V_ijk(nV,nV,nV))
+  allocate(X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO), X_vvoo(nV,nV,nO,nO))
+  allocate(T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO), T_vo(nV,nO))
+
+  ! Temporary arrays
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,T_vvoo,T_ovvo,T_vo,X_vvvo,X_ovoo,X_vvoo, &
+  !$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
+  !$OMP PRIVATE(a,b,c,d,i,j,k,l) &
+  !$OMP DEFAULT(NONE)
+  
+  !v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+  !X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j)
+  
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do a = 1, nV
+      do b = 1, nV
+        do d = 1, nV
+          X_vvvo(d,b,a,i) = v_vvvo(b,a,d,i)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do d = 1, nV
+          T_vvoo(d,c,k,j) = t2(k,j,c,d)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+ 
+  !v_vooo(c,j,k,l) * t2(i,l,a,b) &
+  !X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
+
+  !$OMP DO collapse(3)
+  do k = 1, nO
+    do j = 1, nO
+      do c = 1, nV
+        do l = 1, nO
+           X_ovoo(l,c,j,k) = v_vooo(c,j,k,l)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do b = 1, nV
+      do a = 1, nV
+        do l = 1, nO
+          T_ovvo(l,a,b,i) = t2(i,l,a,b)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+                     
+  !v_vvoo(b,c,j,k) * t1(i,a) &
+  !X_vvoo(b,c,k,j) * T1_vo(a,i) &
+  
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do b = 1, nV
+          X_vvoo(b,c,k,j) = v_vvoo(b,c,j,k)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(1)
+  do i = 1, nO
+    do a = 1, nV
+      T_vo(a,i) = t1(i,a)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call wall_time(ta)
+  energy = 0d0
+  do i = 1, nO
+    do j = 1, nO
+      do k = 1, nO
+        delta_ijk = f_o(i) + f_o(j) + f_o(k)
+        call form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_ijk)
+        call form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,W_ijk,V_ijk)
+        !$OMP PARALLEL &
+        !$OMP SHARED(energy,nV,i,j,k,W_ijk,V_ijk,f_o,f_v,delta_ijk) &
+        !$OMP PRIVATE(a,b,c,e,delta) &
+        !$OMP DEFAULT(NONE)
+        e = 0d0
+        !$OMP DO
+        do c = 1, nV
+          do b = 1, nV
+            do a = 1, nV
+              delta = 1d0 / (delta_ijk - f_v(a) - f_v(b) - f_v(c))
+              !energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c))  !delta_ooovvv(i,j,k,a,b,c)
+              e = e + (4d0 * W_ijk(a,b,c) + W_ijk(b,c,a) + W_ijk(c,a,b)) &
+                       * (V_ijk(a,b,c) - V_ijk(c,b,a)) * delta
+            enddo
+          enddo
+        enddo
+        !$OMP END DO
+        !$OMP CRITICAL
+        energy = energy + e
+        !$OMP END CRITICAL
+        !$OMP END PARALLEL
+      enddo
+    enddo
+    call wall_time(tb)
+    write(*,'(F12.2,A5,F12.2,A2)') dble(i)/dble(nO)*100d0, '% in ', tb - ta, ' s'
+  enddo
+  
+  energy = energy / 3d0
+
+  deallocate(W_ijk,V_ijk,X_vvvo,X_ovoo,T_vvoo,T_ovvo,T_vo)
+  !deallocate(V,W)
+end
+
+! W_ijk
+
+subroutine form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV,i,j,k
+  !double precision, intent(in) :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO)
+  double precision, intent(in)  :: X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO)
+  double precision, intent(out) :: W(nV,nV,nV)!,nO,nO,nO)
+  
+  integer :: l,a,b,c,d
+
+  !W = 0d0
+  !do i = 1, nO
+  !  do j = 1, nO
+  !    do k = 1, nO
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W) &
+  !$OMP PRIVATE(a,b,c,d,l) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        W(a,b,c) = 0d0
+
+        do d = 1, nV
+          !W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+          W(a,b,c) = W(a,b,c) &
+          ! chem (bd|ai)
+          ! phys <ba|di>
+          !+ cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+          !+ cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+          !+ cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+          !+ cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+          !+ cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+          !+ cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+          + X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) &
+          + X_vvvo(d,c,a,i) * T_vvoo(d,b,j,k) & ! bc kj
+          + X_vvvo(d,a,c,k) * T_vvoo(d,b,j,i) & ! prev ac ik
+          + X_vvvo(d,b,c,k) * T_vvoo(d,a,i,j) & ! prev ab ij
+          + X_vvvo(d,c,b,j) * T_vvoo(d,a,i,k) & ! prev bc kj
+          + X_vvvo(d,a,b,j) * T_vvoo(d,c,k,i) ! prev ac ik
+        enddo
+        
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+         
+        do l = 1, nO
+          !W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+          W(a,b,c) = W(a,b,c) &
+          ! chem (ck|jl)
+          ! phys <cj|kl>
+          !- cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
+          !- cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
+          !- cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
+          !- cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
+          !- cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
+          !- cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
+          - X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
+          - X_ovoo(l,b,k,j) * T_ovvo(l,a,c,i) & ! bc kj
+          - X_ovoo(l,b,i,j) * T_ovvo(l,c,a,k) & ! prev ac ik
+          - X_ovoo(l,a,j,i) * T_ovvo(l,c,b,k) & ! prev ab ij
+          - X_ovoo(l,a,k,i) * T_ovvo(l,b,c,j) & ! prev bc kj
+          - X_ovoo(l,c,i,k) * T_ovvo(l,b,a,j) ! prev ac ik
+        enddo
+
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+  !    enddo
+  !  enddo
+  !enddo
+
+end
+
+! V_ijk
+
+subroutine form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,w,v)
+
+implicit none
+
+  integer, intent(in)           :: nO,nV,i,j,k
+  !double precision, intent(in)  :: t1(nO,nV)
+  double precision, intent(in)  :: T_vo(nV,nO)
+  double precision, intent(in)  :: X_vvoo(nV,nV,nO,nO)
+  double precision, intent(in)  :: W(nV,nV,nV)!,nO,nO,nO)
+  double precision, intent(out) :: V(nV,nV,nV)!,nO,nO,nO)
+
+  integer :: a,b,c
+
+  !V = 0d0
+  !do i = 1, nO
+  !  do j = 1, nO
+  !    do k = 1, nO
+  
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,i,j,k,T_vo,X_vvoo,W,V) &
+  !$OMP PRIVATE(a,b,c) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        !V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+        V(a,b,c) = W(a,b,c) &
+        !+ cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+        !+ cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+        !+ cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+        + X_vvoo(b,c,k,j) * T_vo(a,i) &
+        + X_vvoo(a,c,k,i) * T_vo(b,j) &
+        + X_vvoo(a,b,j,i) * T_vo(c,k)
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+  !    enddo
+  !  enddo
+  !enddo
+
+end

src/ccsd/ccsd_t_spin_orb.irp.f

@@ -0,0 +1,376 @@
+! v1
+
+subroutine ccsd_par_t_spin(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,v_vvvo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+  double precision, intent(in)  :: v_ooov(nO,nO,nO,nV)
+  double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO), v_vvvo(nV,nV,nV,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: t3(:,:,:,:,:,:), s(:,:)
+  double precision              :: e_t, e_st, e_dt, delta_abc, delta
+  integer                       :: i,j,k,l,m,a,b,c,d,e
+
+  allocate(t3(nO,nO,nO,nV,nV,nV), s(nO,nV))
+
+  t3 = 0d0
+  
+  ! T3
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              do e = 1, nV
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+                  + t2(j,k,a,e) * v_vvvo(b,c,e,i) &
+                  - t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
+                  - t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
+                  - t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
+                  - t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+                  + t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+                  + t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+                  + t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+                  + t2(j,i,c,e) * v_vvvo(b,a,e,k)   ! + P(ik) P(ac)
+              enddo
+              do m = 1, nO
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+                  + t2(m,i,b,c) * v_ooov(j,k,m,a) &
+                  - t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
+                  - t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
+                  - t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
+                  - t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+                  + t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+                  + t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+                  + t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+                  + t2(m,k,b,a) * v_ooov(j,i,m,c)   ! + P(ik) P(ac)
+              enddo
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) * (1d0 / delta)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  
+
+  ! E_T
+  e_t = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              e_t = e_t + t3(i,j,k,a,b,c) * delta * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  e_t = e_t / 36d0
+
+  ! E_ST
+  s = 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
+              s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  e_st = 0d0
+  do a = 1, nV
+    do i = 1, nO
+      e_st = e_st + s(i,a) * t1(i,a)
+    enddo
+  enddo
+  e_st = e_st * 0.25d0
+
+  ! E_DT
+  e_dt = 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
+              e_dt = e_dt + t2(i,j,a,b) * f_ov(k,c) * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  e_dt = e_dt * 0.25d0
+
+  ! (T)
+  !print*,e_t,e_st,e_dt
+  energy = e_t + e_st + e_dt
+
+  deallocate(t3,s)
+  
+end
+
+! v2
+
+subroutine ccsd_par_t_spin_v2(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+  double precision, intent(in)  :: v_ooov(nO,nO,nO,nV)
+  double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: t3_bc(:,:,:,:), s(:,:), e_t(:), e_dt(:)
+  double precision, allocatable :: A_vovv(:,:,:,:), v_vvvo(:,:,:,:)
+  double precision, allocatable :: T_voov(:,:,:,:), B_ooov(:,:,:,:)
+  double precision              :: e_st, delta_abc, delta, ta, tb
+  integer                       :: i,j,k,l,m,a,b,c,d,e
+
+  allocate(t3_bc(nO,nO,nO,nV), s(nO,nV), e_t(nV), e_dt(nV))
+  allocate(A_vovv(nV,nO,nV,nV),v_vvvo(nV,nV,nV,nO),T_voov(nV,nO,nO,nV),B_ooov(nO,nO,nO,nV))
+
+  call gen_v_spin(cc_nV_m,cc_nV_m,cc_nV_m,cc_nO_m, &
+       cc_nV_S,cc_nV_S,cc_nV_S,cc_nO_S, &
+       cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, &
+       nV,nV,nV,nO, v_vvvo)
+  
+  ! Init 
+  s = 0d0
+  e_t = 0d0
+  e_st = 0d0
+  e_dt = 0d0
+
+  call wall_time(ta)
+  !$OMP PARALLEL &
+  !$OMP PRIVATE(i,j,k,m,a,b,c,e) &
+  !$OMP SHARED(A_vovv,ta,tb,t3_bc,s,e_t,e_st,e_dt,t2,v_vvvo,v_ooov, &
+  !$OMP v_vvoo,f_o,f_v,f_ov,delta,delta_abc,nO,nV,T_voov,B_ooov) &
+  !$OMP DEFAULT(NONE)
+
+  !$OMP DO collapse(3)
+  do c = 1, nV
+    do b = 1, nV
+      do i = 1, nO
+        do e = 1, nV
+          A_vovv(e,i,b,c) = v_vvvo(b,c,e,i)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$omp do collapse(3)
+  do a = 1, nV
+    do k = 1, nO
+      do j = 1, nO
+        do e = 1, nV
+          T_voov(e,j,k,a) = t2(j,k,a,e)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$omp end do nowait
+
+  !$omp do collapse(3)
+  do a = 1, nV
+    do k = 1, nO
+      do j = 1, nO
+        do m = 1, nO
+          B_ooov(m,j,k,a) = v_ooov(j,k,m,a)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$omp end do
+
+  do c = 1, nV
+    do b = 1, nV
+
+      ! T3(:,:,:,:,b,c)
+      ! Init
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              t3_bc(i,j,k,a) = 0d0
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              do e = 1, nV
+                t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
+                   !+ t2(j,k,a,e) * v_vvvo(b,c,e,i) &
+                   !- t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
+                   !- t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
+                   !- t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
+                   !- t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+                   !+ t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+                   !+ t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+                   !+ t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+                   !+ t2(j,i,c,e) * v_vvvo(b,a,e,k)   ! + P(ik) P(ac)
+                   + T_voov(e,j,k,a) * A_vovv(e,i,b,c) &
+                   - T_voov(e,i,k,a) * A_vovv(e,j,b,c) & ! - P(ij)
+                   - T_voov(e,j,i,a) * A_vovv(e,k,b,c) & ! - P(ik)
+                   - T_voov(e,j,k,b) * A_vovv(e,i,a,c) & ! - P(ab)
+                   - T_voov(e,j,k,c) * A_vovv(e,i,b,a) & ! - P(ac)
+                   + T_voov(e,i,k,b) * A_vovv(e,j,a,c) & ! + P(ij) P(ab)
+                   + T_voov(e,i,k,c) * A_vovv(e,j,b,a) & ! + P(ij) P(ac)
+                   + T_voov(e,j,i,b) * A_vovv(e,k,a,c) & ! + P(ik) P(ab)
+                   + T_voov(e,j,i,c) * A_vovv(e,k,b,a)   ! + P(ik) P(ac)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              do m = 1, nO
+                t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
+                   !+ t2(m,i,b,c) * v_ooov(j,k,m,a) &
+                   !- t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
+                   !- t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
+                   !- t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
+                   !- t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+                   !+ t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+                   !+ t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+                   !+ t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+                   !+ t2(m,k,b,a) * v_ooov(j,i,m,c)   ! + P(ik) P(ac)
+                   + t2(m,i,b,c) * B_ooov(m,j,k,a) &
+                   - t2(m,j,b,c) * B_ooov(m,i,k,a) & ! - P(ij)
+                   - t2(m,k,b,c) * B_ooov(m,j,i,a) & ! - P(ik)
+                   - t2(m,i,a,c) * B_ooov(m,j,k,b) & ! - P(ab)
+                   - t2(m,i,b,a) * B_ooov(m,j,k,c) & ! - P(ac)
+                   + t2(m,j,a,c) * B_ooov(m,i,k,b) & ! + P(ij) P(ab)
+                   + t2(m,j,b,a) * B_ooov(m,i,k,c) & ! + P(ij) P(ac)
+                   + t2(m,k,a,c) * B_ooov(m,j,i,b) & ! + P(ik) P(ab)
+                   + t2(m,k,b,a) * B_ooov(m,j,i,c)   ! + P(ik) P(ac)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+
+      !$OMP DO
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+               delta = f_o(i) + f_o(j) + f_o(k) - delta_abc 
+               t3_bc(i,j,k,a) = t3_bc(i,j,k,a) * (1d0 / delta)
+             enddo
+           enddo
+         enddo
+       enddo
+      !$OMP END DO
+
+      ! E_T
+      !$OMP DO
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              e_t(a) = e_t(a) + t3_bc(i,j,k,a) * delta * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO nowait
+
+      ! E_ST
+      !$OMP DO
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO nowait
+
+      ! E_DT
+      !$OMP DO
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              e_dt(a) = e_dt(a) + t2(i,j,a,b) * f_ov(k,c) * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+    enddo
+    !$OMP MASTER
+    call wall_time(tb)
+    write(*,'(A1,F6.2,A5,F10.2,A2)') ' ', dble(c)/dble(nV)*100d0, '% in ', tb-ta, ' s'
+    !$OMP END MASTER
+  enddo
+  !$OMP END PARALLEL
+
+  do a = 2, nV
+    e_t(1) = e_t(1) + e_t(a)
+  enddo
+  
+  do a = 2, nV
+    e_dt(1) = e_dt(1) + e_dt(a)
+  enddo
+
+  e_t = e_t / 36d0
+  
+  do a = 1, nV
+    do i = 1, nO
+      e_st = e_st + s(i,a) * t1(i,a)
+    enddo
+  enddo
+  e_st = e_st * 0.25d0
+
+  e_dt = e_dt * 0.25d0
+
+  ! (T)
+  !print*,e_t(1),e_st,e_dt(1)
+  energy = e_t(1) + e_st + e_dt(1)
+
+  deallocate(t3_bc,s)
+  
+end

src/ccsd/org/ccsd_t_space_orb.org

@@ -0,0 +1,428 @@
+Ref:
+Integral-Direct and Parallel Implementation of the CCSD(T) Method:
+Algorithmic Developments and Large-Scale Applications
+László Gyevi-Nagy, Mihály Kállay, and Péter R. Nagy
+J. Chem. Theory Comput. 2020, 16, 1, 366–384
+https://doi.org/10.1021/acs.jctc.9b00957
+
+* Dumb way
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine ccsd_par_t_space(nO,nV,t1,t2,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO, nV)
+  double precision, intent(in)  :: t2(nO, nO, nV, nV)
+  double precision, intent(out) :: energy
+  
+  double precision, allocatable :: W(:,:,:,:,:,:)
+  double precision, allocatable :: V(:,:,:,:,:,:)
+  integer :: i,j,k,a,b,c
+  
+  allocate(W(nO,nO,nO,nV,nV,nV))
+  allocate(V(nO,nO,nO,nV,nV,nV))
+
+  call form_w(nO,nV,t2,W)   
+  call form_v(nO,nV,t1,W,V)
+
+  energy = 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
+              energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c))  !delta_ooovvv(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  
+  energy = energy / 3d0
+  
+  deallocate(V,W)
+end
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine form_w(nO,nV,t2,W)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t2(nO, nO, nV, nV)
+  double precision, intent(out) :: W(nO, nO, nO, nV, nV, nV)
+  
+  integer :: i,j,k,l,a,b,c,d
+
+  W = 0d0
+  do c = 1, nV
+    print*,'W:',c,'/',nV
+    do b = 1, nV
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+
+              do d = 1, nV
+                W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+                ! chem (bd|ai)
+                ! phys <ba|di>
+                + cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+                + cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+                + cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+                + cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+                + cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+                + cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+              enddo
+
+              do l = 1, nO
+                W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+                ! chem (ck|jl)
+                ! phys <cj|kl>
+                - cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
+                - cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
+                - cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
+                - cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
+                - cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
+                - cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
+              enddo
+
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine form_v(nO,nV,t1,w,v)
+
+implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO, nV)
+  double precision, intent(in)  :: W(nO, nO, nO, nV, nV, nV)
+  double precision, intent(out) :: V(nO, nO, nO, nV, nV, nV)
+
+  integer :: i,j,k,a,b,c
+
+  V = 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
+              V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+              + cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+              + cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+              + cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end
+#+END_SRC
+
+* Better way
+** Main
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine ccsd_par_t_space_v2(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV
+  double precision, intent(in)  :: t1(nO,nV), f_o(nO), f_v(nV)
+  double precision, intent(in)  :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO)
+  double precision, intent(out) :: energy
+  
+  double precision, allocatable :: W(:,:,:,:,:,:)
+  double precision, allocatable :: V(:,:,:,:,:,:)
+  double precision, allocatable :: W_ijk(:,:,:), V_ijk(:,:,:)
+  double precision, allocatable :: X_vvvo(:,:,:,:), X_ovoo(:,:,:,:), X_vvoo(:,:,:,:)
+  double precision, allocatable :: T_vvoo(:,:,:,:), T_ovvo(:,:,:,:), T_vo(:,:)
+  integer                       :: i,j,k,l,a,b,c,d
+  double precision              :: e,ta,tb, delta, delta_ijk
+ 
+  !allocate(W(nV,nV,nV,nO,nO,nO))
+  !allocate(V(nV,nV,nV,nO,nO,nO))
+  allocate(W_ijk(nV,nV,nV), V_ijk(nV,nV,nV))
+  allocate(X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO), X_vvoo(nV,nV,nO,nO))
+  allocate(T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO), T_vo(nV,nO))
+
+  ! Temporary arrays
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,T_vvoo,T_ovvo,T_vo,X_vvvo,X_ovoo,X_vvoo, &
+  !$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
+  !$OMP PRIVATE(a,b,c,d,i,j,k,l) &
+  !$OMP DEFAULT(NONE)
+  
+  !v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+  !X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j)
+  
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do a = 1, nV
+      do b = 1, nV
+        do d = 1, nV
+          X_vvvo(d,b,a,i) = v_vvvo(b,a,d,i)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do d = 1, nV
+          T_vvoo(d,c,k,j) = t2(k,j,c,d)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+ 
+  !v_vooo(c,j,k,l) * t2(i,l,a,b) &
+  !X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
+
+  !$OMP DO collapse(3)
+  do k = 1, nO
+    do j = 1, nO
+      do c = 1, nV
+        do l = 1, nO
+           X_ovoo(l,c,j,k) = v_vooo(c,j,k,l)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(3)
+  do i = 1, nO
+    do b = 1, nV
+      do a = 1, nV
+        do l = 1, nO
+          T_ovvo(l,a,b,i) = t2(i,l,a,b)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+                     
+  !v_vvoo(b,c,j,k) * t1(i,a) &
+  !X_vvoo(b,c,k,j) * T1_vo(a,i) &
+  
+  !$OMP DO collapse(3)
+  do j = 1, nO
+    do k = 1, nO
+      do c = 1, nV
+        do b = 1, nV
+          X_vvoo(b,c,k,j) = v_vvoo(b,c,j,k)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(1)
+  do i = 1, nO
+    do a = 1, nV
+      T_vo(a,i) = t1(i,a)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call wall_time(ta)
+  energy = 0d0
+  do i = 1, nO
+    do j = 1, nO
+      do k = 1, nO
+        delta_ijk = f_o(i) + f_o(j) + f_o(k)
+        call form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_ijk)
+        call form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,W_ijk,V_ijk)
+        !$OMP PARALLEL &
+        !$OMP SHARED(energy,nV,i,j,k,W_ijk,V_ijk,f_o,f_v,delta_ijk) &
+        !$OMP PRIVATE(a,b,c,e,delta) &
+        !$OMP DEFAULT(NONE)
+        e = 0d0
+        !$OMP DO
+        do c = 1, nV
+          do b = 1, nV
+            do a = 1, nV
+              delta = 1d0 / (delta_ijk - f_v(a) - f_v(b) - f_v(c))
+              !energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c))  !delta_ooovvv(i,j,k,a,b,c)
+              e = e + (4d0 * W_ijk(a,b,c) + W_ijk(b,c,a) + W_ijk(c,a,b)) &
+                       * (V_ijk(a,b,c) - V_ijk(c,b,a)) * delta
+            enddo
+          enddo
+        enddo
+        !$OMP END DO
+        !$OMP CRITICAL
+        energy = energy + e
+        !$OMP END CRITICAL
+        !$OMP END PARALLEL
+      enddo
+    enddo
+    call wall_time(tb)
+    write(*,'(F12.2,A5,F12.2,A2)') dble(i)/dble(nO)*100d0, '% in ', tb - ta, ' s'
+  enddo
+  
+  energy = energy / 3d0
+
+  deallocate(W_ijk,V_ijk,X_vvvo,X_ovoo,T_vvoo,T_ovvo,T_vo)
+  !deallocate(V,W)
+end
+#+END_SRC
+
+** W_ijk
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine form_w_ijk(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W)
+
+  implicit none
+
+  integer, intent(in)           :: nO,nV,i,j,k
+  !double precision, intent(in) :: t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO)
+  double precision, intent(in)  :: X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO)
+  double precision, intent(out) :: W(nV,nV,nV)!,nO,nO,nO)
+  
+  integer :: l,a,b,c,d
+
+  !W = 0d0
+  !do i = 1, nO
+  !  do j = 1, nO
+  !    do k = 1, nO
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,i,j,k,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W) &
+  !$OMP PRIVATE(a,b,c,d,l) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        W(a,b,c) = 0d0
+
+        do d = 1, nV
+          !W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+          W(a,b,c) = W(a,b,c) &
+          ! chem (bd|ai)
+          ! phys <ba|di>
+          !+ cc_space_v_vvvo(b,a,d,i) * t2(k,j,c,d) &
+          !+ cc_space_v_vvvo(c,a,d,i) * t2(j,k,b,d) & ! bc kj
+          !+ cc_space_v_vvvo(a,c,d,k) * t2(j,i,b,d) & ! prev ac ik
+          !+ cc_space_v_vvvo(b,c,d,k) * t2(i,j,a,d) & ! prev ab ij
+          !+ cc_space_v_vvvo(c,b,d,j) * t2(i,k,a,d) & ! prev bc kj
+          !+ cc_space_v_vvvo(a,b,d,j) * t2(k,i,c,d) ! prev ac ik
+          + X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) &
+          + X_vvvo(d,c,a,i) * T_vvoo(d,b,j,k) & ! bc kj
+          + X_vvvo(d,a,c,k) * T_vvoo(d,b,j,i) & ! prev ac ik
+          + X_vvvo(d,b,c,k) * T_vvoo(d,a,i,j) & ! prev ab ij
+          + X_vvvo(d,c,b,j) * T_vvoo(d,a,i,k) & ! prev bc kj
+          + X_vvvo(d,a,b,j) * T_vvoo(d,c,k,i) ! prev ac ik
+        enddo
+        
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+         
+        do l = 1, nO
+          !W(i,j,k,a,b,c) = W(i,j,k,a,b,c) &
+          W(a,b,c) = W(a,b,c) &
+          ! chem (ck|jl)
+          ! phys <cj|kl>
+          !- cc_space_v_vooo(c,j,k,l) * t2(i,l,a,b) &
+          !- cc_space_v_vooo(b,k,j,l) * t2(i,l,a,c) & ! bc kj
+          !- cc_space_v_vooo(b,i,j,l) * t2(k,l,c,a) & ! prev ac ik
+          !- cc_space_v_vooo(a,j,i,l) * t2(k,l,c,b) & ! prev ab ij
+          !- cc_space_v_vooo(a,k,i,l) * t2(j,l,b,c) & ! prev bc kj
+          !- cc_space_v_vooo(c,i,k,l) * t2(j,l,b,a) ! prev ac ik
+          - X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
+          - X_ovoo(l,b,k,j) * T_ovvo(l,a,c,i) & ! bc kj
+          - X_ovoo(l,b,i,j) * T_ovvo(l,c,a,k) & ! prev ac ik
+          - X_ovoo(l,a,j,i) * T_ovvo(l,c,b,k) & ! prev ab ij
+          - X_ovoo(l,a,k,i) * T_ovvo(l,b,c,j) & ! prev bc kj
+          - X_ovoo(l,c,i,k) * T_ovvo(l,b,a,j) ! prev ac ik
+        enddo
+
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+  !    enddo
+  !  enddo
+  !enddo
+
+end
+#+END_SRC
+
+** V_ijk
+#+BEGIN_SRC f90 :comments org :tangle ccsd_t_space_orb.irp.f
+subroutine form_v_ijk(nO,nV,i,j,k,T_vo,X_vvoo,w,v)
+
+implicit none
+
+  integer, intent(in)           :: nO,nV,i,j,k
+  !double precision, intent(in)  :: t1(nO,nV)
+  double precision, intent(in)  :: T_vo(nV,nO)
+  double precision, intent(in)  :: X_vvoo(nV,nV,nO,nO)
+  double precision, intent(in)  :: W(nV,nV,nV)!,nO,nO,nO)
+  double precision, intent(out) :: V(nV,nV,nV)!,nO,nO,nO)
+
+  integer :: a,b,c
+
+  !V = 0d0
+  !do i = 1, nO
+  !  do j = 1, nO
+  !    do k = 1, nO
+  
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,i,j,k,T_vo,X_vvoo,W,V) &
+  !$OMP PRIVATE(a,b,c) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(2)
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        !V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
+        V(a,b,c) = W(a,b,c) &
+        !+ cc_space_v_vvoo(b,c,j,k) * t1(i,a) &
+        !+ cc_space_v_vvoo(a,c,i,k) * t1(j,b) &
+        !+ cc_space_v_vvoo(a,b,i,j) * t1(k,c)
+        + X_vvoo(b,c,k,j) * T_vo(a,i) &
+        + X_vvoo(a,c,k,i) * T_vo(b,j) &
+        + X_vvoo(a,b,j,i) * T_vo(c,k)
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+  !    enddo
+  !  enddo
+  !enddo
+
+end
+#+END_SRC

src/ccsd/org/ccsd_t_spin_orb.org

@@ -0,0 +1,385 @@
+* CCSD(T) spin orb
+
+Ref:
+John D. Watts, Jürgen Gauss, and Rodney J. Bartlett
+J. Chem. Phys. 98, 8718 (1993)
+http://dx.doi.org/10.1063/1.464480
+
+** v1
+#+begin_src f90 :comments org :tangle ccsd_t_spin_orb.irp.f
+subroutine ccsd_par_t_spin(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,v_vvvo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+  double precision, intent(in)  :: v_ooov(nO,nO,nO,nV)
+  double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO), v_vvvo(nV,nV,nV,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: t3(:,:,:,:,:,:), s(:,:)
+  double precision              :: e_t, e_st, e_dt, delta_abc, delta
+  integer                       :: i,j,k,l,m,a,b,c,d,e
+
+  allocate(t3(nO,nO,nO,nV,nV,nV), s(nO,nV))
+
+  t3 = 0d0
+  
+  ! T3
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              do e = 1, nV
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+                  + t2(j,k,a,e) * v_vvvo(b,c,e,i) &
+                  - t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
+                  - t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
+                  - t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
+                  - t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+                  + t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+                  + t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+                  + t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+                  + t2(j,i,c,e) * v_vvvo(b,a,e,k)   ! + P(ik) P(ac)
+              enddo
+              do m = 1, nO
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) &
+                  + t2(m,i,b,c) * v_ooov(j,k,m,a) &
+                  - t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
+                  - t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
+                  - t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
+                  - t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+                  + t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+                  + t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+                  + t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+                  + t2(m,k,b,a) * v_ooov(j,i,m,c)   ! + P(ik) P(ac)
+              enddo
+                t3(i,j,k,a,b,c) = t3(i,j,k,a,b,c) * (1d0 / delta)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  
+
+  ! E_T
+  e_t = 0d0
+  do c = 1, nV
+    do b = 1, nV
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              e_t = e_t + t3(i,j,k,a,b,c) * delta * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  e_t = e_t / 36d0
+
+  ! E_ST
+  s = 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
+              s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  e_st = 0d0
+  do a = 1, nV
+    do i = 1, nO
+      e_st = e_st + s(i,a) * t1(i,a)
+    enddo
+  enddo
+  e_st = e_st * 0.25d0
+
+  ! E_DT
+  e_dt = 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
+              e_dt = e_dt + t2(i,j,a,b) * f_ov(k,c) * t3(i,j,k,a,b,c)
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+  e_dt = e_dt * 0.25d0
+
+  ! (T)
+  !print*,e_t,e_st,e_dt
+  energy = e_t + e_st + e_dt
+
+  deallocate(t3,s)
+  
+end
+#+end_src
+
+** v2
+#+begin_src f90 :comments org :tangle ccsd_t_spin_orb.irp.f
+subroutine ccsd_par_t_spin_v2(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,energy)
+
+  implicit none
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+  double precision, intent(in)  :: v_ooov(nO,nO,nO,nV)
+  double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO)
+  double precision, intent(out) :: energy
+
+  double precision, allocatable :: t3_bc(:,:,:,:), s(:,:), e_t(:), e_dt(:)
+  double precision, allocatable :: A_vovv(:,:,:,:), v_vvvo(:,:,:,:)
+  double precision, allocatable :: T_voov(:,:,:,:), B_ooov(:,:,:,:)
+  double precision              :: e_st, delta_abc, delta, ta, tb
+  integer                       :: i,j,k,l,m,a,b,c,d,e
+
+  allocate(t3_bc(nO,nO,nO,nV), s(nO,nV), e_t(nV), e_dt(nV))
+  allocate(A_vovv(nV,nO,nV,nV),v_vvvo(nV,nV,nV,nO),T_voov(nV,nO,nO,nV),B_ooov(nO,nO,nO,nV))
+
+  call gen_v_spin(cc_nV_m,cc_nV_m,cc_nV_m,cc_nO_m, &
+       cc_nV_S,cc_nV_S,cc_nV_S,cc_nO_S, &
+       cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, &
+       nV,nV,nV,nO, v_vvvo)
+  
+  ! Init 
+  s = 0d0
+  e_t = 0d0
+  e_st = 0d0
+  e_dt = 0d0
+
+  call wall_time(ta)
+  !$OMP PARALLEL &
+  !$OMP PRIVATE(i,j,k,m,a,b,c,e) &
+  !$OMP SHARED(A_vovv,ta,tb,t3_bc,s,e_t,e_st,e_dt,t2,v_vvvo,v_ooov, &
+  !$OMP v_vvoo,f_o,f_v,f_ov,delta,delta_abc,nO,nV,T_voov,B_ooov) &
+  !$OMP DEFAULT(NONE)
+
+  !$OMP DO collapse(3)
+  do c = 1, nV
+    do b = 1, nV
+      do i = 1, nO
+        do e = 1, nV
+          A_vovv(e,i,b,c) = v_vvvo(b,c,e,i)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO nowait
+
+  !$omp do collapse(3)
+  do a = 1, nV
+    do k = 1, nO
+      do j = 1, nO
+        do e = 1, nV
+          T_voov(e,j,k,a) = t2(j,k,a,e)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$omp end do nowait
+
+  !$omp do collapse(3)
+  do a = 1, nV
+    do k = 1, nO
+      do j = 1, nO
+        do m = 1, nO
+          B_ooov(m,j,k,a) = v_ooov(j,k,m,a)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$omp end do
+
+  do c = 1, nV
+    do b = 1, nV
+
+      ! T3(:,:,:,:,b,c)
+      ! Init
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              t3_bc(i,j,k,a) = 0d0
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              do e = 1, nV
+                t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
+                   !+ t2(j,k,a,e) * v_vvvo(b,c,e,i) &
+                   !- t2(i,k,a,e) * v_vvvo(b,c,e,j) & ! - P(ij)
+                   !- t2(j,i,a,e) * v_vvvo(b,c,e,k) & ! - P(ik)
+                   !- t2(j,k,b,e) * v_vvvo(a,c,e,i) & ! - P(ab)
+                   !- t2(j,k,c,e) * v_vvvo(b,a,e,i) & ! - P(ac)
+                   !+ t2(i,k,b,e) * v_vvvo(a,c,e,j) & ! + P(ij) P(ab)
+                   !+ t2(i,k,c,e) * v_vvvo(b,a,e,j) & ! + P(ij) P(ac)
+                   !+ t2(j,i,b,e) * v_vvvo(a,c,e,k) & ! + P(ik) P(ab)
+                   !+ t2(j,i,c,e) * v_vvvo(b,a,e,k)   ! + P(ik) P(ac)
+                   + T_voov(e,j,k,a) * A_vovv(e,i,b,c) &
+                   - T_voov(e,i,k,a) * A_vovv(e,j,b,c) & ! - P(ij)
+                   - T_voov(e,j,i,a) * A_vovv(e,k,b,c) & ! - P(ik)
+                   - T_voov(e,j,k,b) * A_vovv(e,i,a,c) & ! - P(ab)
+                   - T_voov(e,j,k,c) * A_vovv(e,i,b,a) & ! - P(ac)
+                   + T_voov(e,i,k,b) * A_vovv(e,j,a,c) & ! + P(ij) P(ab)
+                   + T_voov(e,i,k,c) * A_vovv(e,j,b,a) & ! + P(ij) P(ac)
+                   + T_voov(e,j,i,b) * A_vovv(e,k,a,c) & ! + P(ik) P(ab)
+                   + T_voov(e,j,i,c) * A_vovv(e,k,b,a)   ! + P(ik) P(ac)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      
+      !$OMP DO collapse(3)
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              do m = 1, nO
+                t3_bc(i,j,k,a) = t3_bc(i,j,k,a) &
+                   !+ t2(m,i,b,c) * v_ooov(j,k,m,a) &
+                   !- t2(m,j,b,c) * v_ooov(i,k,m,a) & ! - P(ij)
+                   !- t2(m,k,b,c) * v_ooov(j,i,m,a) & ! - P(ik)
+                   !- t2(m,i,a,c) * v_ooov(j,k,m,b) & ! - P(ab)
+                   !- t2(m,i,b,a) * v_ooov(j,k,m,c) & ! - P(ac)
+                   !+ t2(m,j,a,c) * v_ooov(i,k,m,b) & ! + P(ij) P(ab)
+                   !+ t2(m,j,b,a) * v_ooov(i,k,m,c) & ! + P(ij) P(ac)
+                   !+ t2(m,k,a,c) * v_ooov(j,i,m,b) & ! + P(ik) P(ab)
+                   !+ t2(m,k,b,a) * v_ooov(j,i,m,c)   ! + P(ik) P(ac)
+                   + t2(m,i,b,c) * B_ooov(m,j,k,a) &
+                   - t2(m,j,b,c) * B_ooov(m,i,k,a) & ! - P(ij)
+                   - t2(m,k,b,c) * B_ooov(m,j,i,a) & ! - P(ik)
+                   - t2(m,i,a,c) * B_ooov(m,j,k,b) & ! - P(ab)
+                   - t2(m,i,b,a) * B_ooov(m,j,k,c) & ! - P(ac)
+                   + t2(m,j,a,c) * B_ooov(m,i,k,b) & ! + P(ij) P(ab)
+                   + t2(m,j,b,a) * B_ooov(m,i,k,c) & ! + P(ij) P(ac)
+                   + t2(m,k,a,c) * B_ooov(m,j,i,b) & ! + P(ik) P(ab)
+                   + t2(m,k,b,a) * B_ooov(m,j,i,c)   ! + P(ik) P(ac)
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+
+      !$OMP DO
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+               delta = f_o(i) + f_o(j) + f_o(k) - delta_abc 
+               t3_bc(i,j,k,a) = t3_bc(i,j,k,a) * (1d0 / delta)
+             enddo
+           enddo
+         enddo
+       enddo
+      !$OMP END DO
+
+      ! E_T
+      !$OMP DO
+      do a = 1, nV
+        delta_abc = f_v(a) + f_v(b) + f_v(c)
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              delta = f_o(i) + f_o(j) + f_o(k) - delta_abc
+              e_t(a) = e_t(a) + t3_bc(i,j,k,a) * delta * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO nowait
+
+      ! E_ST
+      !$OMP DO
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              s(i,a) = s(i,a) + v_vvoo(b,c,j,k) * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO nowait
+
+      ! E_DT
+      !$OMP DO
+      do a = 1, nV
+        do k = 1, nO
+          do j = 1, nO
+            do i = 1, nO
+              e_dt(a) = e_dt(a) + t2(i,j,a,b) * f_ov(k,c) * t3_bc(i,j,k,a)
+            enddo
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+    enddo
+    !$OMP MASTER
+    call wall_time(tb)
+    write(*,'(A1,F6.2,A5,F10.2,A2)') ' ', dble(c)/dble(nV)*100d0, '% in ', tb-ta, ' s'
+    !$OMP END MASTER
+  enddo
+  !$OMP END PARALLEL
+
+  do a = 2, nV
+    e_t(1) = e_t(1) + e_t(a)
+  enddo
+  
+  do a = 2, nV
+    e_dt(1) = e_dt(1) + e_dt(a)
+  enddo
+
+  e_t = e_t / 36d0
+  
+  do a = 1, nV
+    do i = 1, nO
+      e_st = e_st + s(i,a) * t1(i,a)
+    enddo
+  enddo
+  e_st = e_st * 0.25d0
+
+  e_dt = e_dt * 0.25d0
+
+  ! (T)
+  !print*,e_t(1),e_st,e_dt(1)
+  energy = e_t(1) + e_st + e_dt(1)
+
+  deallocate(t3_bc,s)
+  
+end
+#+end_src

src/utils_cc/EZFIO.cfg

@@ -0,0 +1,77 @@
+[cc_thresh_conv]
+type: double precision
+doc: Threshold for the convergence of the residual equations.
+interface: ezfio,ocaml,provider
+default: 1e-6
+
+[cc_max_iter]
+type: integer
+doc: Maximum number of iterations.
+interface: ezfio,ocaml,provider
+default: 100
+
+[cc_diis_depth]
+type: integer
+doc: Maximum depth of the DIIS, i.e., maximum number of iterations that the DIIS keeps in memory. Warning, we allocate matrices with the diis depth at the beginning without update. If you don't have enough memory it should crash in memory.
+interface: ezfio,ocaml,provider
+default: 8
+
+[cc_level_shift]
+type: double precision
+doc: Level shift for the CC
+interface: ezfio,ocaml,provider
+default: 0.0
+
+[cc_level_shift_guess]
+type: double precision
+doc: Level shift for the guess of the CC amplitudes
+interface: ezfio,ocaml,provider
+default: 0.0
+
+[cc_update_method]
+type: character*(32)
+doc: Method used to update the CC amplitudes. none -> normal, diis -> with diis.
+interface: ezfio,ocaml,provider
+default: diis
+
+[cc_guess_t1]
+type: character*(32)
+doc: Guess used to initialize the T1 amplitudes. none -> 0, MP -> perturbation theory, read -> read from disk.
+interface: ezfio,ocaml,provider
+default: MP
+
+[cc_guess_t2]
+type: character*(32)
+doc: Guess used to initialize the T2 amplitudes. none -> 0, MP -> perturbation theory, read -> read from disk.
+interface: ezfio,ocaml,provider
+default: MP
+
+[cc_write_t1]
+type: logical
+doc: If true, it will write on disk the T1 amplitudes at the end of the calculation.
+interface: ezfio,ocaml,provider
+default: False
+
+[cc_write_t2]
+type: logical
+doc: If true, it will write on disk the T2 amplitudes at the end of the calculation.
+interface: ezfio,ocaml,provider
+default: False
+
+[cc_par_t]
+type: logical
+doc: If true, the CCSD(T) will be computed.
+interface: ezfio,ocaml,provider
+default: False
+
+[cc_dev]
+type: logical
+doc: Only for dev purposes.
+interface: ezfio,ocaml,provider
+default: False
+
+[cc_ref]
+type: integer
+doc: Index of the reference determinant in psi_det for CC calculation.
+interface: ezfio,ocaml,provider
+default: 1

src/utils_cc/NEED

@@ -0,0 +1,4 @@
+hartree_fock
+two_body_rdm
+bitmask
+determinants

src/utils_cc/README.md

@@ -0,0 +1,34 @@
+# Utils for CC
+
+Utils for the CC modules.
+
+## Contents
+- Providers related to reference occupancy
+- Integrals related to the reference
+- Diis for CC (but can be used for something else if you provide your own error vector)
+- Guess for CC amplitudes
+- Routines to update the CC amplitudes
+- Phase between to arbitrary determinants
+- print of the qp edit wf
+
+## Keywords
+- cc_thresh_conv: Threshold for the convergence of the residual equations. Default: 1e-6.
+- cc_max_iter: Maximum number of iterations. Default: 100.
+- cc_diis_depth: Diis depth. Default: 8.
+- cc_level_shift: Level shift for the CC. Default: 0.0.
+- cc_level_shift_guess: Level shift for the MP guess of the amplitudes. Default: 0.0.
+- cc_update_method: Method used to update the CC amplitudes. none -> normal, diis -> with diis. Default: diis.
+- cc_guess_t1: Guess used to initialize the T1 amplitudes. none -> 0, MP -> perturbation theory, read -> read from disk. Default: MP.
+- cc_guess_t2: Guess used to initialize the T2 amplitudes. none -> 0, MP -> perturbation theory, read -> read from disk. Default: MP.
+- cc_write_t1: If true, it will write on disk the T1 amplitudes at the end of the calculation. Default: False.
+- cc_write_t2: If true, it will write on disk the T2 amplitudes at the end of the calculation. Default: False.
+- cc_par_t: If true, the CCSD(T) will be computed.
+- cc_ref: Index of the reference determinant in psi_det for CC calculation. Default: 1.
+
+## Org files
+The org files are stored in the directory org in order to avoid overwriting on user changes.
+The org files can be modified, to export the change to the source code, run
+```
+./TANGLE_org_mode.sh and 
+mv *.irp.f ../.
+```

src/utils_cc/diis.irp.f

@@ -0,0 +1,529 @@
+! Code
+
+subroutine diis_cc(all_err,all_t,sze,m,iter,t)
+
+  implicit none
+
+  BEGIN_DOC
+  ! DIIS. Take the error vectors and the amplitudes of the previous
+  ! iterations to compute the new amplitudes
+  END_DOC
+  
+  ! {err_i}_{i=1}^{m_it} -> B -> c
+  ! {t_i}_{i=1}^{m_it}, c, {err_i}_{i=1}^{m_it} -> t_{m_it+1}
+
+  integer, intent(in)             :: m,iter,sze
+  double precision, intent(in)    :: all_err(sze,m)
+  double precision, intent(in)    :: all_t(sze,m)
+  
+  double precision, intent(out)   :: t(sze)
+  
+  double precision, allocatable   :: B(:,:), c(:), zero(:)
+  integer                         :: m_iter
+  integer                         :: i,j,k
+  integer                         :: info
+  integer, allocatable            :: ipiv(:)
+  double precision                :: accu
+  
+  m_iter = min(m,iter)
+  !print*,'m_iter',m_iter
+  allocate(B(m_iter+1,m_iter+1), c(m_iter), zero(m_iter+1))
+  allocate(ipiv(m+1))
+
+  ! B(i,j) =  < err(iter-m_iter+j),err(iter-m_iter+i) > ! iter-m_iter will be zero for us
+  B = 0d0
+  !$OMP PARALLEL &
+  !$OMP SHARED(B,m,m_iter,sze,all_err) &
+  !$OMP PRIVATE(i,j,k,accu) &
+  !$OMP DEFAULT(NONE)
+  do j = 1, m_iter
+    do i = 1, m_iter
+      accu = 0d0
+      !$OMP DO
+      do k = 1, sze
+        ! the errors of the ith iteration are in all_err(:,m+1-i)
+        accu = accu + all_err(k,m+1-i) * all_err(k,m+1-j)
+      enddo
+      !$OMP END DO NOWAIT
+      !$OMP CRITICAL
+      B(i,j) = B(i,j) + accu
+      !$OMP END CRITICAL
+    enddo
+  enddo
+  !$OMP END PARALLEL
+  
+  do i = 1, m_iter
+    B(i,m_iter+1) = -1
+  enddo
+  do j = 1, m_iter
+    B(m_iter+1,j) = -1
+  enddo
+  ! Debug
+  !print*,'B'
+  !do i = 1, m_iter+1
+  !  write(*,'(100(F10.6))') B(i,:)
+  !enddo
+
+  ! (0 0 .... 0 -1)
+  zero = 0d0
+  zero(m_iter+1) = -1d0
+
+  ! Solve B.c = zero
+  call dgesv(m_iter+1, 1, B, size(B,1), ipiv, zero, size(zero,1), info)
+  if (info /= 0) then
+    print*,'DIIS error in dgesv:', info
+    call abort
+  endif
+  ! c corresponds to the m_iter first solutions
+  c = zero(1:m_iter)
+  ! Debug
+  !print*,'c',c
+  !print*,'all_t' 
+  !do i = 1, m
+  !  write(*,'(100(F10.6))') all_t(:,i)
+  !enddo
+  !print*,'all_err' 
+  !do i = 1, m
+  !  write(*,'(100(F10.6))') all_err(:,i)
+  !enddo
+
+  ! update T
+  !$OMP PARALLEL &
+  !$OMP SHARED(t,c,m,all_err,all_t,sze,m_iter) &
+  !$OMP PRIVATE(i,j,accu) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do i = 1, sze
+    t(i) = 0d0
+  enddo
+  !$OMP END DO
+  do i = 1, m_iter
+    !$OMP DO
+    do j = 1, sze
+      t(j) = t(j) + c(i) * (all_t(j,m+1-i) + all_err(j,m+1-i))
+    enddo
+    !$OMP END DO
+  enddo
+  !$OMP END PARALLEL
+
+  !print*,'new t',t
+
+  deallocate(ipiv,B,c,zero)
+
+end
+
+! Update all err
+
+subroutine update_all_err(err,all_err,sze,m,iter)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Shift all the err vectors of the previous iterations to add the new one
+  ! The last err vector is placed in the last position and all the others are
+  ! moved toward the first one.
+  END_DOC
+
+  integer, intent(in)             :: m, iter, sze
+  double precision, intent(in)    :: err(sze)
+  double precision, intent(inout) :: all_err(sze,m)
+  integer                         :: i,j
+  integer                         :: m_iter
+
+  m_iter = min(m,iter)
+
+  ! Shift
+  !$OMP PARALLEL &
+  !$OMP SHARED(m,all_err,err,sze) &
+  !$OMP PRIVATE(i,j) &
+  !$OMP DEFAULT(NONE)
+  do i = 1, m-1
+    !$OMP DO
+    do j = 1, sze
+      all_err(j,i) = all_err(j,i+1)
+    enddo
+    !$OMP END DO
+  enddo
+  
+  ! Debug
+  !print*,'shift err'
+  !do i = 1, m
+  !  print*,i, all_err(:,i)
+  !enddo
+
+  ! New
+  !$OMP DO
+  do i = 1, sze
+    all_err(i,m) = err(i)
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Debug
+  !print*,'Updated err'
+  !do i = 1, m
+  !  print*,i, all_err(:,i)
+  !enddo
+
+end
+
+! Update all t
+
+subroutine update_all_t(t,all_t,sze,m,iter)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Shift all the t vectors of the previous iterations to add the new one
+  ! The last t vector is placed in the last position and all the others are
+  ! moved toward the first one.
+  END_DOC
+
+  integer, intent(in)             :: m, iter, sze
+  double precision, intent(in)    :: t(sze)
+  double precision, intent(inout) :: all_t(sze,m)
+  integer                         :: i,j
+  integer                         :: m_iter
+
+  m_iter = min(m,iter)
+
+  ! Shift
+  !$OMP PARALLEL &
+  !$OMP SHARED(m,all_t,t,sze) &
+  !$OMP PRIVATE(i,j) &
+  !$OMP DEFAULT(NONE)
+  do i = 1, m-1
+    !$OMP DO
+    do j = 1, sze
+      all_t(j,i) = all_t(j,i+1)
+    enddo
+    !$OMP END DO
+  enddo
+
+  ! New
+  !$OMP DO
+  do i = 1, sze
+    all_t(i,m) = t(i)
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Debug
+  !print*,'Updated t'
+  !do i = 1, m
+  !  print*,i, all_t(:,i)
+  !enddo
+
+end
+
+! Err1
+
+subroutine compute_err1(nO,nV,f_o,f_v,r1,err1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute the error vector for the t1
+  END_DOC
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: f_o(nO), f_v(nV), r1(nO,nV)
+  
+  double precision, intent(out) :: err1(nO,nV)
+
+  integer                       :: i,a
+  
+  !$OMP PARALLEL &
+  !$OMP SHARED(err1,r1,f_o,f_v,nO,nV,cc_level_shift) &
+  !$OMP PRIVATE(i,a) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do a = 1, nV
+    do i = 1, nO
+      err1(i,a) = - r1(i,a) / (f_o(i) - f_v(a) - cc_level_shift)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+end
+
+! Err2
+
+subroutine compute_err2(nO,nV,f_o,f_v,r2,err2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute the error vector for the t2
+  END_DOC
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: f_o(nO), f_v(nV), r2(nO,nO,nV,nV)
+  
+  double precision, intent(out) :: err2(nO,nO,nV,nV)
+
+  integer                       :: i,j,a,b
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(err2,r2,f_o,f_v,nO,nV,cc_level_shift) &
+  !$OMP PRIVATE(i,j,a,b) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(3)
+  do b = 1, nV
+    do a = 1, nV
+      do j = 1, nO
+        do i = 1, nO       
+          err2(i,j,a,b) = - r2(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+end
+
+! Update t
+
+subroutine update_t_ccsd(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
+
+  implicit none
+
+  integer, intent(in)             :: nO,nV,nb_iter
+  double precision, intent(in)    :: f_o(nO), f_v(nV)
+  double precision, intent(in)    :: r1(nO,nV), r2(nO,nO,nV,nV)
+  
+  double precision, intent(inout) :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(inout) :: all_err1(nO*nV, cc_diis_depth), all_err2(nO*nO*nV*nV, cc_diis_depth)
+  double precision, intent(inout) :: all_t1(nO*nV, cc_diis_depth), all_t2(nO*nO*nV*nV, cc_diis_depth)
+
+  double precision, allocatable   :: err1(:,:), err2(:,:,:,:)
+  double precision, allocatable   :: tmp_err1(:), tmp_err2(:)
+  double precision, allocatable   :: tmp_t1(:), tmp_t2(:)
+  
+  if (cc_update_method == 'diis') then
+
+    allocate(err1(nO,nV), err2(nO,nO,nV,nV))
+    allocate(tmp_err1(nO*nV), tmp_err2(nO*nO*nV*nV))
+    allocate(tmp_t1(nO*nV), tmp_t2(nO*nO*nV*nV))
+
+    ! DIIS T1, it is not always good since the t1 can be small
+    ! That's why there is a call to update the t1 in the standard way
+    ! T1 error tensor
+    !call compute_err1(nO,nV,f_o,f_v,r1,err1)
+    ! Transfo errors and parameters in vectors
+    !tmp_err1 = reshape(err1,(/nO*nV/))
+    !tmp_t1   = reshape(t1  ,(/nO*nV/))
+    ! Add the error and parameter vectors with those of the previous iterations
+    !call update_all_err(tmp_err1,all_err1,nO*nV,cc_diis_depth,nb_iter+1)
+    !call update_all_t  (tmp_t1  ,all_t1  ,nO*nV,cc_diis_depth,nb_iter+1)
+    ! Diis and reshape T as a tensor
+    !call diis_cc(all_err1,all_t1,nO*nV,cc_diis_depth,nb_iter+1,tmp_t1)
+    !t1 = reshape(tmp_t1  ,(/nO,nV/))
+    call update_t1(nO,nV,f_o,f_v,r1,t1)
+
+    ! DIIS T2
+    ! T2 error tensor
+    call compute_err2(nO,nV,f_o,f_v,r2,err2)
+    ! Transfo errors and parameters in vectors
+    tmp_err2 = reshape(err2,(/nO*nO*nV*nV/))
+    tmp_t2   = reshape(t2  ,(/nO*nO*nV*nV/))
+    ! Add the error and parameter vectors with those of the previous iterations
+    call update_all_err(tmp_err2,all_err2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+    call update_all_t  (tmp_t2  ,all_t2  ,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+    ! Diis and reshape T as a tensor
+    call diis_cc(all_err2,all_t2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp_t2)
+    t2 = reshape(tmp_t2  ,(/nO,nO,nV,nV/))
+
+    deallocate(tmp_t1,tmp_t2,tmp_err1,tmp_err2,err1,err2)
+
+  ! Standard update as T = T - Delta
+  elseif (cc_update_method == 'none') then
+     
+    call update_t1(nO,nV,f_o,f_v,r1,t1)
+    call update_t2(nO,nV,f_o,f_v,r2,t2)
+    
+  else
+    print*,'Unkonw cc_method_method: '//cc_update_method
+  endif
+  
+end
+
+! Update t v2
+
+subroutine update_t_ccsd_diis(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
+
+  implicit none
+
+  integer, intent(in)             :: nO,nV,nb_iter
+  double precision, intent(in)    :: f_o(nO), f_v(nV)
+  double precision, intent(in)    :: r1(nO,nV), r2(nO,nO,nV,nV)
+  
+  double precision, intent(inout) :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(inout) :: all_err1(nO*nV, cc_diis_depth), all_err2(nO*nO*nV*nV, cc_diis_depth)
+  double precision, intent(inout) :: all_t1(nO*nV, cc_diis_depth), all_t2(nO*nO*nV*nV, cc_diis_depth)
+
+  double precision, allocatable   :: all_t(:,:), all_err(:,:), tmp_t(:)
+  double precision, allocatable   :: err1(:,:), err2(:,:,:,:)
+  double precision, allocatable   :: tmp_err1(:), tmp_err2(:)
+  double precision, allocatable   :: tmp_t1(:), tmp_t2(:)
+
+  integer                         :: i,j
+  
+  ! Allocate
+  allocate(all_err(nO*nV+nO*nO*nV*nV,cc_diis_depth), all_t(nO*nV+nO*nO*nV*nV,cc_diis_depth))
+  allocate(tmp_t(nO*nV+nO*nO*nV*nV))
+  allocate(err1(nO,nV), err2(nO,nO,nV,nV))
+  allocate(tmp_err1(nO*nV), tmp_err2(nO*nO*nV*nV))
+  allocate(tmp_t1(nO*nV), tmp_t2(nO*nO*nV*nV))
+
+  ! Compute the errors and reshape them as vector
+  call compute_err1(nO,nV,f_o,f_v,r1,err1)
+  call compute_err2(nO,nV,f_o,f_v,r2,err2)
+  tmp_err1 = reshape(err1,(/nO*nV/))
+  tmp_err2 = reshape(err2,(/nO*nO*nV*nV/))
+  tmp_t1   = reshape(t1  ,(/nO*nV/))
+  tmp_t2   = reshape(t2  ,(/nO*nO*nV*nV/))
+  
+  ! Update the errors and parameters for the diis
+  call update_all_err(tmp_err1,all_err1,nO*nV,cc_diis_depth,nb_iter+1)
+  call update_all_t  (tmp_t1  ,all_t1  ,nO*nV,cc_diis_depth,nb_iter+1)
+  call update_all_err(tmp_err2,all_err2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+  call update_all_t  (tmp_t2  ,all_t2  ,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+
+  ! Gather the different parameters and errors
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,all_err,all_err1,all_err2,cc_diis_depth,&
+  !$OMP all_t,all_t1,all_t2) &
+  !$OMP PRIVATE(i,j) &
+  !$OMP DEFAULT(NONE)
+  do j = 1, cc_diis_depth
+    !$OMP DO 
+    do i = 1, nO*nV
+      all_err(i,j) = all_err1(i,j)
+    enddo
+    !$OMP END DO NOWAIT
+  enddo
+  do j = 1, cc_diis_depth
+    !$OMP DO
+    do i = 1, nO*nO*nV*nV
+      all_err(i+nO*nV,j) = all_err2(i,j)
+    enddo
+    !$OMP END DO NOWAIT
+  enddo
+  do j = 1, cc_diis_depth
+    !$OMP DO 
+    do i = 1, nO*nV
+      all_t(i,j) = all_t1(i,j)
+    enddo
+    !$OMP END DO NOWAIT
+  enddo
+  do j = 1, cc_diis_depth
+    !$OMP DO 
+    do i = 1, nO*nO*nV*nV
+      all_t(i+nO*nV,j) = all_t2(i,j)
+    enddo
+    !$OMP END DO
+  enddo
+  !$OMP END PARALLEL
+  
+  ! Diis
+  call diis_cc(all_err,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp_t)
+
+  ! Split the resulting vector
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,tmp_t,tmp_t1,tmp_t2) &
+  !$OMP PRIVATE(i) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do i = 1, nO*nV
+    tmp_t1(i) = tmp_t(i)
+  enddo
+  !$OMP END DO NOWAIT
+  !$OMP DO
+  do i = 1, nO*nO*nV*nV
+    tmp_t2(i) = tmp_t(i+nO*nV) 
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Reshape as tensors
+  t1 = reshape(tmp_t1 ,(/nO,nV/))
+  t2 = reshape(tmp_t2 ,(/nO,nO,nV,nV/))
+
+  ! Deallocate
+  deallocate(tmp_t1,tmp_t2,tmp_err1,tmp_err2,err1,err2,all_t,all_err)
+
+end
+
+! Update t v3
+
+subroutine update_t_ccsd_diis_v3(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err,all_t)
+
+  implicit none
+
+  integer, intent(in)             :: nO,nV,nb_iter
+  double precision, intent(in)    :: f_o(nO), f_v(nV)
+  double precision, intent(in)    :: r1(nO,nV), r2(nO,nO,nV,nV)
+  
+  double precision, intent(inout) :: t1(nO*nV), t2(nO*nO*nV*nV)
+  double precision, intent(inout) :: all_err(nO*nV+nO*nO*nV*nV, cc_diis_depth)
+  double precision, intent(inout) :: all_t(nO*nV+nO*nO*nV*nV, cc_diis_depth)
+
+  double precision, allocatable   :: tmp(:)
+
+  integer                         :: i,j
+  
+  ! Allocate
+  allocate(tmp(nO*nV+nO*nO*nV*nV))
+
+  ! Compute the errors
+  call compute_err1(nO,nV,f_o,f_v,r1,tmp(1:nO*nV))
+  call compute_err2(nO,nV,f_o,f_v,r2,tmp(nO*nV+1:nO*nV+nO*nO*nV*nV))
+  
+  ! Update the errors and parameters for the diis
+  call update_all_err(tmp,all_err,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,tmp,t1,t2) &
+  !$OMP PRIVATE(i) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do i = 1, nO*nV
+    tmp(i) = t1(i)
+  enddo
+  !$OMP END DO NOWAIT
+  !$OMP DO
+  do i = 1, nO*nO*nV*nV
+    tmp(i+nO*nV) = t2(i)
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+   
+  call update_all_t(tmp,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+
+  ! Diis
+  call diis_cc(all_err,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp)
+
+  ! Split the resulting vector
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,tmp,t1,t2) &
+  !$OMP PRIVATE(i) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do i = 1, nO*nV
+    t1(i) = tmp(i)
+  enddo
+  !$OMP END DO NOWAIT
+  !$OMP DO
+  do i = 1, nO*nO*nV*nV
+    t2(i) = tmp(i+nO*nV) 
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Deallocate
+  deallocate(tmp)
+
+end

src/utils_cc/energy.irp.f

@@ -0,0 +1,13 @@
+subroutine det_energy(det,energy)
+
+  implicit none
+
+  integer(bit_kind), intent(in) :: det
+
+  double precision, intent(out) :: energy
+
+  call i_H_j(det,det,N_int,energy)
+
+  energy = energy + nuclear_repulsion
+  
+end

src/utils_cc/guess_t.irp.f

@@ -0,0 +1,213 @@
+! T1
+
+subroutine guess_t1(nO,nV,f_o,f_v,f_ov,t1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Update the T1 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+
+  ! inout
+  double precision, intent(out) :: t1(nO, nV)
+
+  ! internal
+  integer                       :: i,a
+
+  if (trim(cc_guess_t1) == 'none') then
+     t1 = 0d0
+  else if (trim(cc_guess_t1) == 'MP') then
+    do a = 1, nV
+      do i = 1, nO
+        t1(i,a) = f_ov(i,a) / (f_o(i) - f_v(a) - cc_level_shift_guess)
+      enddo
+    enddo
+  else if (trim(cc_guess_t1) == 'read') then
+    call read_t1(nO,nV,t1)
+  else
+    print*, 'Unknown cc_guess_t1 type: '//trim(cc_guess_t1)
+    call abort
+  endif
+  
+end
+
+! T2
+
+subroutine guess_t2(nO,nV,f_o,f_v,v_oovv,t2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Update the T2 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: f_o(nO), f_v(nV), v_oovv(nO, nO, nV, nV)
+
+  ! inout
+  double precision, intent(out) :: t2(nO, nO, nV, nV)
+
+  ! internal
+  integer                       :: i,j,a,b
+
+  if (trim(cc_guess_t2) == 'none') then
+    t2 = 0d0
+  else if (trim(cc_guess_t2) == 'MP') then
+    do b = 1, nV
+      do a = 1, nV
+        do j = 1, nO
+          do i = 1, nO
+            t2(i,j,a,b) = v_oovv(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift_guess)
+          enddo
+        enddo
+      enddo
+    enddo
+  else if (trim(cc_guess_t2) == 'read') then
+    call read_t2(nO,nV,t2)
+  else
+    print*, 'Unknown cc_guess_t1 type: '//trim(cc_guess_t2)
+    call abort
+  endif
+  
+end
+
+! T1
+
+subroutine write_t1(nO,nV,t1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Write the T1 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)          :: nO, nV
+  double precision, intent(in) :: t1(nO, nV)
+
+  ! internal
+  integer                      :: i,a
+
+  if (cc_write_t1) then
+    open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
+    do a = 1, nV
+      do i = 1, nO
+         write(11,'(F20.12)') t1(i,a)
+      enddo
+    enddo
+    close(11)
+  endif
+  
+end
+
+! T2
+
+subroutine write_t2(nO,nV,t2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Write the T2 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)          :: nO, nV
+  double precision, intent(in) :: t2(nO, nO, nV, nV)
+
+  ! internal
+  integer                      :: i,j,a,b
+
+  if (cc_write_t2) then
+    open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
+    do b = 1, nV
+      do a = 1, nV
+        do j = 1, nO
+          do i = 1, nO
+             write(11,'(F20.12)') t2(i,j,a,b)
+          enddo
+        enddo
+      enddo
+    enddo
+    close(11)
+  endif
+  
+end
+
+! T1
+
+subroutine read_t1(nO,nV,t1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Read the T1 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)           :: nO, nV
+  double precision, intent(out) :: t1(nO, nV)
+
+  ! internal
+  integer                       :: i,a
+  logical                       :: ok
+
+  inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
+  if (.not. ok) then
+     print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
+     print*, 'Do a first calculation with cc_write_t1 = True'
+     print*, 'and cc_guess_t1 /= read before setting cc_guess_t1 = read'
+     call abort
+  endif
+  open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
+  do a = 1, nV
+    do i = 1, nO
+       read(11,'(F20.12)') t1(i,a)
+    enddo
+  enddo
+  close(11)
+  
+end
+
+! T2
+
+subroutine read_t2(nO,nV,t2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Read the T2 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)           :: nO, nV
+  double precision, intent(out) :: t2(nO, nO, nV, nV)
+
+  ! internal
+  integer                       :: i,j,a,b
+  logical                       :: ok
+
+  inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
+  if (.not. ok) then
+     print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
+     print*, 'Do a first calculation with cc_write_t2 = True'
+     print*, 'and cc_guess_t2 /= read before setting cc_guess_t2 = read'
+     call abort
+  endif
+  open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
+  do b = 1, nV
+    do a = 1, nV
+      do j = 1, nO
+        do i = 1, nO
+           read(11,'(F20.12)') t2(i,j,a,b)
+        enddo
+      enddo
+    enddo
+  enddo
+  close(11)
+  
+end

src/utils_cc/occupancy.irp.f

@@ -0,0 +1,328 @@
+! N spin orb
+
+subroutine extract_n_spin(det,n)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Returns the number of occupied alpha, occupied beta, virtual alpha, virtual beta spin orbitals
+  ! in det without counting the core and deleted orbitals in the format n(nOa,nOb,nVa,nVb)
+  END_DOC
+
+  integer(bit_kind), intent(in) :: det(N_int,2)
+  
+  integer, intent(out)          :: n(4)
+  
+  integer(bit_kind)             :: res(N_int,2)
+  integer                       :: i, si
+  logical                       :: ok, is_core, is_del
+
+  ! Init
+  n = 0
+
+  ! Loop over the spin
+  do si = 1, 2
+    do i = 1, mo_num
+      call apply_hole(det, si, i, res, ok, N_int)
+      
+      ! in core ?
+      if (is_core(i)) cycle
+      ! in del ?
+      if (is_del(i)) cycle
+      
+      if (ok) then
+        ! particle
+        n(si) = n(si) + 1
+      else
+        ! hole
+        n(si+2) = n(si+2) + 1
+      endif
+    enddo
+  enddo
+
+  !print*,n(1),n(2),n(3),n(4)
+
+end
+
+! Spin
+
+subroutine extract_list_orb_spin(det,nO_m,nV_m,list_occ,list_vir)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Returns the the list of occupied alpha/beta, virtual alpha/beta spin orbitals
+  ! size(nO_m,1) must be max(nOa,nOb) and size(nV_m,1) must be max(nVa,nVb)
+  END_DOC
+  
+  integer, intent(in)           :: nO_m, nV_m
+  integer(bit_kind), intent(in) :: det(N_int,2)
+  
+  integer, intent(out)          :: list_occ(nO_m,2), list_vir(nV_m,2)
+  
+  integer(bit_kind)             :: res(N_int,2)
+  integer                       :: i, si, idx_o, idx_v, idx_i, idx_b
+  logical                       :: ok, is_core, is_del
+
+  list_occ = 0
+  list_vir = 0
+
+  ! List of occ/vir alpha/beta
+
+  ! occ alpha -> list_occ(:,1)
+  ! occ beta -> list_occ(:,2)
+  ! vir alpha -> list_vir(:,1)
+  ! vir beta -> list_vir(:,2)
+
+  ! Loop over the spin 
+  do si = 1, 2
+    ! tmp idx
+    idx_o = 1
+    idx_v = 1
+    do i = 1, mo_num
+      call apply_hole(det, si, i, res, ok, N_int)
+
+      ! in core ?
+      if (is_core(i)) cycle
+      ! in del ?
+      if (is_del(i)) cycle
+            
+      if (ok) then
+        ! particle
+        list_occ(idx_o,si) = i
+        idx_o = idx_o + 1
+      else
+        ! hole
+        list_vir(idx_v,si) = i
+        idx_v = idx_v + 1
+      endif
+    enddo
+  enddo
+
+end
+
+! Space
+
+subroutine extract_list_orb_space(det,nO,nV,list_occ,list_vir)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Returns the the list of occupied and virtual alpha spin orbitals
+  END_DOC
+  
+  integer, intent(in)           :: nO, nV
+  integer(bit_kind), intent(in) :: det(N_int,2)
+  
+  integer, intent(out)          :: list_occ(nO), list_vir(nV)
+  
+  integer(bit_kind)             :: res(N_int,2)
+  integer                       :: i, si, idx_o, idx_v, idx_i, idx_b
+  logical                       :: ok, is_core, is_del
+  
+  if (elec_alpha_num /= elec_beta_num) then
+    print*,'Error elec_alpha_num /= elec_beta_num, impossible to create cc_list_occ and cc_list_vir, abort'
+    call abort
+  endif
+
+  list_occ = 0
+  list_vir = 0
+
+  ! List of occ/vir alpha
+
+  ! occ alpha -> list_occ(:,1)
+  ! vir alpha -> list_vir(:,1)
+
+  ! tmp idx
+  idx_o = 1
+  idx_v = 1
+  do i = 1, mo_num
+    call apply_hole(det, 1, i, res, ok, N_int)
+
+    ! in core ?
+    if (is_core(i)) cycle
+    ! in del ?
+    if (is_del(i)) cycle
+
+    if (ok) then
+      ! particle
+      list_occ(idx_o) = i
+      idx_o = idx_o + 1
+    else
+      ! hole
+      list_vir(idx_v) = i
+      idx_v = idx_v + 1
+    endif
+  enddo
+
+end
+
+! is_core
+
+function is_core(i)
+
+  implicit none
+
+  BEGIN_DOC
+  ! True if the orbital i is a core orbital
+  END_DOC
+
+  integer, intent(in) :: i
+  logical             :: is_core
+
+  integer             :: j
+
+  ! Init
+  is_core = .False.
+
+  ! Search
+  do j = 1, dim_list_core_orb
+    if (list_core(j) == i) then
+      is_core = .True.
+      exit
+    endif
+  enddo
+
+end
+
+! is_del
+
+function is_del(i)
+
+  implicit none
+
+  BEGIN_DOC
+  ! True if the orbital i is a deleted orbital
+  END_DOC
+
+  integer, intent(in) :: i
+  logical             :: is_del
+
+  integer             :: j
+
+  ! Init
+  is_del = .False.
+
+  ! Search
+  do j = 1, dim_list_core_orb
+    if (list_core(j) == i) then
+      is_del = .True.
+      exit
+    endif
+  enddo
+
+end
+
+! N orb
+
+BEGIN_PROVIDER [integer, cc_nO_m]
+&BEGIN_PROVIDER [integer, cc_nOa]
+&BEGIN_PROVIDER [integer, cc_nOb]
+&BEGIN_PROVIDER [integer, cc_nOab]
+&BEGIN_PROVIDER [integer, cc_nV_m]
+&BEGIN_PROVIDER [integer, cc_nVa]
+&BEGIN_PROVIDER [integer, cc_nVb]
+&BEGIN_PROVIDER [integer, cc_nVab]
+&BEGIN_PROVIDER [integer, cc_n_mo]
+&BEGIN_PROVIDER [integer, cc_nO_S, (2)]
+&BEGIN_PROVIDER [integer, cc_nV_S, (2)]
+
+  implicit none
+
+  BEGIN_DOC
+  ! Number of orbitals without core and deleted ones of the cc_ref det in psi_det
+  ! a: alpha, b: beta
+  ! nO_m: max(a,b) occupied 
+  ! nOa: nb a occupied 
+  ! nOb: nb b occupied 
+  ! nOab: nb a+b occupied 
+  ! nV_m: max(a,b) virtual 
+  ! nVa: nb a virtual 
+  ! nVb: nb b virtual 
+  ! nVab: nb a+b virtual 
+  END_DOC
+
+  integer :: n_spin(4)
+
+  ! Extract number of occ/vir alpha/beta spin orbitals
+  call extract_n_spin(psi_det(1,1,cc_ref),n_spin)
+
+  cc_nOa  = n_spin(1)
+  cc_nOb  = n_spin(2)
+  cc_nOab = cc_nOa + cc_nOb    !n_spin(1) + n_spin(2)
+  cc_nO_m = max(cc_nOa,cc_nOb) !max(n_spin(1), n_spin(2))
+  cc_nVa  = n_spin(3)
+  cc_nVb  = n_spin(4)
+  cc_nVab = cc_nVa + cc_nVb    !n_spin(3) + n_spin(4)
+  cc_nV_m = max(cc_nVa,cc_nVb) !max(n_spin(3), n_spin(4))
+  cc_n_mo = cc_nVa + cc_nVb    !n_spin(1) + n_spin(3)
+  cc_nO_S = (/cc_nOa,cc_nOb/)
+  cc_nV_S = (/cc_nVa,cc_nVb/)
+
+END_PROVIDER
+
+! General
+
+BEGIN_PROVIDER [integer, cc_list_gen, (cc_n_mo)]
+
+  implicit none
+
+  BEGIN_DOC
+  ! List of general orbitals without core and deleted ones
+  END_DOC
+
+  integer :: i,j
+  logical :: is_core, is_del
+  
+  j = 1
+  do i = 1, mo_num
+    ! in core ?
+    if (is_core(i)) cycle
+    ! in del ?
+    if (is_del(i)) cycle
+    cc_list_gen(j) = i
+    j = j+1
+  enddo
+
+END_PROVIDER
+
+! Space
+
+BEGIN_PROVIDER [integer, cc_list_occ, (cc_nOa)]
+&BEGIN_PROVIDER [integer, cc_list_vir, (cc_nVa)]
+
+  implicit none
+
+  BEGIN_DOC
+  ! List of occupied and virtual spatial orbitals without core and deleted ones
+  END_DOC
+
+  call extract_list_orb_space(psi_det(1,1,cc_ref),cc_nOa,cc_nVa,cc_list_occ,cc_list_vir)
+
+END_PROVIDER
+
+! Spin
+
+BEGIN_PROVIDER [integer, cc_list_occ_spin, (cc_nO_m,2)]
+&BEGIN_PROVIDER [integer, cc_list_vir_spin, (cc_nV_m,2)]
+&BEGIN_PROVIDER [logical, cc_ref_is_open_shell]
+
+  implicit none
+
+  BEGIN_DOC
+  ! List of occupied and virtual spin orbitals without core and deleted ones
+  END_DOC
+
+  integer :: i
+
+  call extract_list_orb_spin(psi_det(1,1,cc_ref),cc_nO_m,cc_nV_m,cc_list_occ_spin,cc_list_vir_spin)
+
+  cc_ref_is_open_shell = .False.
+  do i = 1, cc_nO_m
+    if (cc_list_occ_spin(i,1) /= cc_list_occ_spin(i,2)) then
+       cc_ref_is_open_shell = .True.
+    endif
+  enddo
+
+
+END_PROVIDER

src/utils_cc/org/TANGLE_org_mode.sh

@@ -0,0 +1,7 @@
+#!/bin/sh   
+
+list='ls *.org'
+for element in $list    
+do   
+		emacs --batch $element -f org-babel-tangle
+done

src/utils_cc/org/diis.org

@@ -0,0 +1,574 @@
+* DIIS
+https://hal.archives-ouvertes.fr/hal-02492983/document
+Maxime Chupin, Mi-Song Dupuy, Guillaume Legendre, Eric Séré. Convergence analysis of adaptive
+DIIS algorithms witerh application to electronic ground state calculations. 
+ESAIM: Mathematical Modelling and Numerical Analysis, EDP Sciences, 2021, 55 (6), pp.2785 - 2825. 10.1051/m2an/2021069ff.ffhal-02492983v5
+
+t_{k+1} = g(t_k)
+err_k = f(t_k) = t_{k+1} - t_k
+
+m_k = min(m,k)
+m maximal depth
+t_{k+1} = \sum_{i=0}^{m_k} c_i^k g(t_{k-m_k+i})
+\sum_{i=0}^{m_k} c_i^k = 1
+
+b_{ij}^k = < err^{k-m_k+j}, err^{k-m_k+i} >
+
+(b   -1) ( c^k    ) = (  0 )
+(-1   0) ( \lambda)   ( -1 )
+
+lambda is used to put the constraint \sum_{i=0}^{m_k} c_i^k = 1
+
+In: t_0, err_0, m
+err_0 = g(t_0)
+k = 0
+m_k = 0
+while ||err_k|| > CC
+  A.x=b
+  t_{k+1} = \sum_{i=0}^{m_k} c_i^k g(t_{k-m_k+i})
+  err_{k+1} = f(t_{k+1})
+  m_{k+1} = min(m_k+1,m)
+  k = k +1
+end
+
+* Code
+#+begin_src f90 :comments org :tangle diis.irp.f 
+subroutine diis_cc(all_err,all_t,sze,m,iter,t)
+
+  implicit none
+
+  BEGIN_DOC
+  ! DIIS. Take the error vectors and the amplitudes of the previous
+  ! iterations to compute the new amplitudes
+  END_DOC
+  
+  ! {err_i}_{i=1}^{m_it} -> B -> c
+  ! {t_i}_{i=1}^{m_it}, c, {err_i}_{i=1}^{m_it} -> t_{m_it+1}
+
+  integer, intent(in)             :: m,iter,sze
+  double precision, intent(in)    :: all_err(sze,m)
+  double precision, intent(in)    :: all_t(sze,m)
+  
+  double precision, intent(out)   :: t(sze)
+  
+  double precision, allocatable   :: B(:,:), c(:), zero(:)
+  integer                         :: m_iter
+  integer                         :: i,j,k
+  integer                         :: info
+  integer, allocatable            :: ipiv(:)
+  double precision                :: accu
+  
+  m_iter = min(m,iter)
+  !print*,'m_iter',m_iter
+  allocate(B(m_iter+1,m_iter+1), c(m_iter), zero(m_iter+1))
+  allocate(ipiv(m+1))
+
+  ! B(i,j) =  < err(iter-m_iter+j),err(iter-m_iter+i) > ! iter-m_iter will be zero for us
+  B = 0d0
+  !$OMP PARALLEL &
+  !$OMP SHARED(B,m,m_iter,sze,all_err) &
+  !$OMP PRIVATE(i,j,k,accu) &
+  !$OMP DEFAULT(NONE)
+  do j = 1, m_iter
+    do i = 1, m_iter
+      accu = 0d0
+      !$OMP DO
+      do k = 1, sze
+        ! the errors of the ith iteration are in all_err(:,m+1-i)
+        accu = accu + all_err(k,m+1-i) * all_err(k,m+1-j)
+      enddo
+      !$OMP END DO NOWAIT
+      !$OMP CRITICAL
+      B(i,j) = B(i,j) + accu
+      !$OMP END CRITICAL
+    enddo
+  enddo
+  !$OMP END PARALLEL
+  
+  do i = 1, m_iter
+    B(i,m_iter+1) = -1
+  enddo
+  do j = 1, m_iter
+    B(m_iter+1,j) = -1
+  enddo
+  ! Debug
+  !print*,'B'
+  !do i = 1, m_iter+1
+  !  write(*,'(100(F10.6))') B(i,:)
+  !enddo
+
+  ! (0 0 .... 0 -1)
+  zero = 0d0
+  zero(m_iter+1) = -1d0
+
+  ! Solve B.c = zero
+  call dgesv(m_iter+1, 1, B, size(B,1), ipiv, zero, size(zero,1), info)
+  if (info /= 0) then
+    print*,'DIIS error in dgesv:', info
+    call abort
+  endif
+  ! c corresponds to the m_iter first solutions
+  c = zero(1:m_iter)
+  ! Debug
+  !print*,'c',c
+  !print*,'all_t' 
+  !do i = 1, m
+  !  write(*,'(100(F10.6))') all_t(:,i)
+  !enddo
+  !print*,'all_err' 
+  !do i = 1, m
+  !  write(*,'(100(F10.6))') all_err(:,i)
+  !enddo
+
+  ! update T
+  !$OMP PARALLEL &
+  !$OMP SHARED(t,c,m,all_err,all_t,sze,m_iter) &
+  !$OMP PRIVATE(i,j,accu) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do i = 1, sze
+    t(i) = 0d0
+  enddo
+  !$OMP END DO
+  do i = 1, m_iter
+    !$OMP DO
+    do j = 1, sze
+      t(j) = t(j) + c(i) * (all_t(j,m+1-i) + all_err(j,m+1-i))
+    enddo
+    !$OMP END DO
+  enddo
+  !$OMP END PARALLEL
+
+  !print*,'new t',t
+
+  deallocate(ipiv,B,c,zero)
+
+end
+#+end_src
+
+** Update all err
+#+begin_src f90 :comments org :tangle diis.irp.f
+subroutine update_all_err(err,all_err,sze,m,iter)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Shift all the err vectors of the previous iterations to add the new one
+  ! The last err vector is placed in the last position and all the others are
+  ! moved toward the first one.
+  END_DOC
+
+  integer, intent(in)             :: m, iter, sze
+  double precision, intent(in)    :: err(sze)
+  double precision, intent(inout) :: all_err(sze,m)
+  integer                         :: i,j
+  integer                         :: m_iter
+
+  m_iter = min(m,iter)
+
+  ! Shift
+  !$OMP PARALLEL &
+  !$OMP SHARED(m,all_err,err,sze) &
+  !$OMP PRIVATE(i,j) &
+  !$OMP DEFAULT(NONE)
+  do i = 1, m-1
+    !$OMP DO
+    do j = 1, sze
+      all_err(j,i) = all_err(j,i+1)
+    enddo
+    !$OMP END DO
+  enddo
+  
+  ! Debug
+  !print*,'shift err'
+  !do i = 1, m
+  !  print*,i, all_err(:,i)
+  !enddo
+
+  ! New
+  !$OMP DO
+  do i = 1, sze
+    all_err(i,m) = err(i)
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Debug
+  !print*,'Updated err'
+  !do i = 1, m
+  !  print*,i, all_err(:,i)
+  !enddo
+
+end
+#+end_src
+
+** Update all t
+#+begin_src f90 :comments org :tangle diis.irp.f
+subroutine update_all_t(t,all_t,sze,m,iter)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Shift all the t vectors of the previous iterations to add the new one
+  ! The last t vector is placed in the last position and all the others are
+  ! moved toward the first one.
+  END_DOC
+
+  integer, intent(in)             :: m, iter, sze
+  double precision, intent(in)    :: t(sze)
+  double precision, intent(inout) :: all_t(sze,m)
+  integer                         :: i,j
+  integer                         :: m_iter
+
+  m_iter = min(m,iter)
+
+  ! Shift
+  !$OMP PARALLEL &
+  !$OMP SHARED(m,all_t,t,sze) &
+  !$OMP PRIVATE(i,j) &
+  !$OMP DEFAULT(NONE)
+  do i = 1, m-1
+    !$OMP DO
+    do j = 1, sze
+      all_t(j,i) = all_t(j,i+1)
+    enddo
+    !$OMP END DO
+  enddo
+
+  ! New
+  !$OMP DO
+  do i = 1, sze
+    all_t(i,m) = t(i)
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Debug
+  !print*,'Updated t'
+  !do i = 1, m
+  !  print*,i, all_t(:,i)
+  !enddo
+
+end
+#+end_src
+
+** Err
+*** Err1
+#+begin_src f90 :comments org :tangle diis.irp.f
+subroutine compute_err1(nO,nV,f_o,f_v,r1,err1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute the error vector for the t1
+  END_DOC
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: f_o(nO), f_v(nV), r1(nO,nV)
+  
+  double precision, intent(out) :: err1(nO,nV)
+
+  integer                       :: i,a
+  
+  !$OMP PARALLEL &
+  !$OMP SHARED(err1,r1,f_o,f_v,nO,nV,cc_level_shift) &
+  !$OMP PRIVATE(i,a) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do a = 1, nV
+    do i = 1, nO
+      err1(i,a) = - r1(i,a) / (f_o(i) - f_v(a) - cc_level_shift)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+end
+#+end_src
+
+*** Err2
+#+begin_src f90 :comments org :tangle diis.irp.f
+subroutine compute_err2(nO,nV,f_o,f_v,r2,err2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute the error vector for the t2
+  END_DOC
+
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: f_o(nO), f_v(nV), r2(nO,nO,nV,nV)
+  
+  double precision, intent(out) :: err2(nO,nO,nV,nV)
+
+  integer                       :: i,j,a,b
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(err2,r2,f_o,f_v,nO,nV,cc_level_shift) &
+  !$OMP PRIVATE(i,j,a,b) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(3)
+  do b = 1, nV
+    do a = 1, nV
+      do j = 1, nO
+        do i = 1, nO       
+          err2(i,j,a,b) = - r2(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+end
+#+end_src
+
+* Gather call diis
+** Update t
+#+begin_src f90 :comments org :tangle diis.irp.f
+subroutine update_t_ccsd(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
+
+  implicit none
+
+  integer, intent(in)             :: nO,nV,nb_iter
+  double precision, intent(in)    :: f_o(nO), f_v(nV)
+  double precision, intent(in)    :: r1(nO,nV), r2(nO,nO,nV,nV)
+  
+  double precision, intent(inout) :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(inout) :: all_err1(nO*nV, cc_diis_depth), all_err2(nO*nO*nV*nV, cc_diis_depth)
+  double precision, intent(inout) :: all_t1(nO*nV, cc_diis_depth), all_t2(nO*nO*nV*nV, cc_diis_depth)
+
+  double precision, allocatable   :: err1(:,:), err2(:,:,:,:)
+  double precision, allocatable   :: tmp_err1(:), tmp_err2(:)
+  double precision, allocatable   :: tmp_t1(:), tmp_t2(:)
+  
+  if (cc_update_method == 'diis') then
+
+    allocate(err1(nO,nV), err2(nO,nO,nV,nV))
+    allocate(tmp_err1(nO*nV), tmp_err2(nO*nO*nV*nV))
+    allocate(tmp_t1(nO*nV), tmp_t2(nO*nO*nV*nV))
+
+    ! DIIS T1, it is not always good since the t1 can be small
+    ! That's why there is a call to update the t1 in the standard way
+    ! T1 error tensor
+    !call compute_err1(nO,nV,f_o,f_v,r1,err1)
+    ! Transfo errors and parameters in vectors
+    !tmp_err1 = reshape(err1,(/nO*nV/))
+    !tmp_t1   = reshape(t1  ,(/nO*nV/))
+    ! Add the error and parameter vectors with those of the previous iterations
+    !call update_all_err(tmp_err1,all_err1,nO*nV,cc_diis_depth,nb_iter+1)
+    !call update_all_t  (tmp_t1  ,all_t1  ,nO*nV,cc_diis_depth,nb_iter+1)
+    ! Diis and reshape T as a tensor
+    !call diis_cc(all_err1,all_t1,nO*nV,cc_diis_depth,nb_iter+1,tmp_t1)
+    !t1 = reshape(tmp_t1  ,(/nO,nV/))
+    call update_t1(nO,nV,f_o,f_v,r1,t1)
+
+    ! DIIS T2
+    ! T2 error tensor
+    call compute_err2(nO,nV,f_o,f_v,r2,err2)
+    ! Transfo errors and parameters in vectors
+    tmp_err2 = reshape(err2,(/nO*nO*nV*nV/))
+    tmp_t2   = reshape(t2  ,(/nO*nO*nV*nV/))
+    ! Add the error and parameter vectors with those of the previous iterations
+    call update_all_err(tmp_err2,all_err2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+    call update_all_t  (tmp_t2  ,all_t2  ,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+    ! Diis and reshape T as a tensor
+    call diis_cc(all_err2,all_t2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp_t2)
+    t2 = reshape(tmp_t2  ,(/nO,nO,nV,nV/))
+
+    deallocate(tmp_t1,tmp_t2,tmp_err1,tmp_err2,err1,err2)
+
+  ! Standard update as T = T - Delta
+  elseif (cc_update_method == 'none') then
+     
+    call update_t1(nO,nV,f_o,f_v,r1,t1)
+    call update_t2(nO,nV,f_o,f_v,r2,t2)
+    
+  else
+    print*,'Unkonw cc_method_method: '//cc_update_method
+  endif
+  
+end
+ #+end_src
+ 
+** Update t v2
+#+begin_src f90 :comments org :tangle diis.irp.f
+subroutine update_t_ccsd_diis(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
+
+  implicit none
+
+  integer, intent(in)             :: nO,nV,nb_iter
+  double precision, intent(in)    :: f_o(nO), f_v(nV)
+  double precision, intent(in)    :: r1(nO,nV), r2(nO,nO,nV,nV)
+  
+  double precision, intent(inout) :: t1(nO,nV), t2(nO,nO,nV,nV)
+  double precision, intent(inout) :: all_err1(nO*nV, cc_diis_depth), all_err2(nO*nO*nV*nV, cc_diis_depth)
+  double precision, intent(inout) :: all_t1(nO*nV, cc_diis_depth), all_t2(nO*nO*nV*nV, cc_diis_depth)
+
+  double precision, allocatable   :: all_t(:,:), all_err(:,:), tmp_t(:)
+  double precision, allocatable   :: err1(:,:), err2(:,:,:,:)
+  double precision, allocatable   :: tmp_err1(:), tmp_err2(:)
+  double precision, allocatable   :: tmp_t1(:), tmp_t2(:)
+
+  integer                         :: i,j
+  
+  ! Allocate
+  allocate(all_err(nO*nV+nO*nO*nV*nV,cc_diis_depth), all_t(nO*nV+nO*nO*nV*nV,cc_diis_depth))
+  allocate(tmp_t(nO*nV+nO*nO*nV*nV))
+  allocate(err1(nO,nV), err2(nO,nO,nV,nV))
+  allocate(tmp_err1(nO*nV), tmp_err2(nO*nO*nV*nV))
+  allocate(tmp_t1(nO*nV), tmp_t2(nO*nO*nV*nV))
+
+  ! Compute the errors and reshape them as vector
+  call compute_err1(nO,nV,f_o,f_v,r1,err1)
+  call compute_err2(nO,nV,f_o,f_v,r2,err2)
+  tmp_err1 = reshape(err1,(/nO*nV/))
+  tmp_err2 = reshape(err2,(/nO*nO*nV*nV/))
+  tmp_t1   = reshape(t1  ,(/nO*nV/))
+  tmp_t2   = reshape(t2  ,(/nO*nO*nV*nV/))
+  
+  ! Update the errors and parameters for the diis
+  call update_all_err(tmp_err1,all_err1,nO*nV,cc_diis_depth,nb_iter+1)
+  call update_all_t  (tmp_t1  ,all_t1  ,nO*nV,cc_diis_depth,nb_iter+1)
+  call update_all_err(tmp_err2,all_err2,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+  call update_all_t  (tmp_t2  ,all_t2  ,nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+
+  ! Gather the different parameters and errors
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,all_err,all_err1,all_err2,cc_diis_depth,&
+  !$OMP all_t,all_t1,all_t2) &
+  !$OMP PRIVATE(i,j) &
+  !$OMP DEFAULT(NONE)
+  do j = 1, cc_diis_depth
+    !$OMP DO 
+    do i = 1, nO*nV
+      all_err(i,j) = all_err1(i,j)
+    enddo
+    !$OMP END DO NOWAIT
+  enddo
+  do j = 1, cc_diis_depth
+    !$OMP DO
+    do i = 1, nO*nO*nV*nV
+      all_err(i+nO*nV,j) = all_err2(i,j)
+    enddo
+    !$OMP END DO NOWAIT
+  enddo
+  do j = 1, cc_diis_depth
+    !$OMP DO 
+    do i = 1, nO*nV
+      all_t(i,j) = all_t1(i,j)
+    enddo
+    !$OMP END DO NOWAIT
+  enddo
+  do j = 1, cc_diis_depth
+    !$OMP DO 
+    do i = 1, nO*nO*nV*nV
+      all_t(i+nO*nV,j) = all_t2(i,j)
+    enddo
+    !$OMP END DO
+  enddo
+  !$OMP END PARALLEL
+  
+  ! Diis
+  call diis_cc(all_err,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp_t)
+
+  ! Split the resulting vector
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,tmp_t,tmp_t1,tmp_t2) &
+  !$OMP PRIVATE(i) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do i = 1, nO*nV
+    tmp_t1(i) = tmp_t(i)
+  enddo
+  !$OMP END DO NOWAIT
+  !$OMP DO
+  do i = 1, nO*nO*nV*nV
+    tmp_t2(i) = tmp_t(i+nO*nV) 
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Reshape as tensors
+  t1 = reshape(tmp_t1 ,(/nO,nV/))
+  t2 = reshape(tmp_t2 ,(/nO,nO,nV,nV/))
+
+  ! Deallocate
+  deallocate(tmp_t1,tmp_t2,tmp_err1,tmp_err2,err1,err2,all_t,all_err)
+
+end
+ #+end_src
+ 
+
+** Update t v3
+#+begin_src f90 :comments org :tangle diis.irp.f
+subroutine update_t_ccsd_diis_v3(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err,all_t)
+
+  implicit none
+
+  integer, intent(in)             :: nO,nV,nb_iter
+  double precision, intent(in)    :: f_o(nO), f_v(nV)
+  double precision, intent(in)    :: r1(nO,nV), r2(nO,nO,nV,nV)
+  
+  double precision, intent(inout) :: t1(nO*nV), t2(nO*nO*nV*nV)
+  double precision, intent(inout) :: all_err(nO*nV+nO*nO*nV*nV, cc_diis_depth)
+  double precision, intent(inout) :: all_t(nO*nV+nO*nO*nV*nV, cc_diis_depth)
+
+  double precision, allocatable   :: tmp(:)
+
+  integer                         :: i,j
+  
+  ! Allocate
+  allocate(tmp(nO*nV+nO*nO*nV*nV))
+
+  ! Compute the errors
+  call compute_err1(nO,nV,f_o,f_v,r1,tmp(1:nO*nV))
+  call compute_err2(nO,nV,f_o,f_v,r2,tmp(nO*nV+1:nO*nV+nO*nO*nV*nV))
+  
+  ! Update the errors and parameters for the diis
+  call update_all_err(tmp,all_err,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,tmp,t1,t2) &
+  !$OMP PRIVATE(i) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do i = 1, nO*nV
+    tmp(i) = t1(i)
+  enddo
+  !$OMP END DO NOWAIT
+  !$OMP DO
+  do i = 1, nO*nO*nV*nV
+    tmp(i+nO*nV) = t2(i)
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+   
+  call update_all_t(tmp,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1)
+
+  ! Diis
+  call diis_cc(all_err,all_t,nO*nV+nO*nO*nV*nV,cc_diis_depth,nb_iter+1,tmp)
+
+  ! Split the resulting vector
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,tmp,t1,t2) &
+  !$OMP PRIVATE(i) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do i = 1, nO*nV
+    t1(i) = tmp(i)
+  enddo
+  !$OMP END DO NOWAIT
+  !$OMP DO
+  do i = 1, nO*nO*nV*nV
+    t2(i) = tmp(i+nO*nV) 
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Deallocate
+  deallocate(tmp)
+
+end
+ #+end_src
+ 

src/utils_cc/org/energy.org

@@ -0,0 +1,15 @@
+#+begin_src f90 :comments org :tangle energy.irp.f
+subroutine det_energy(det,energy)
+
+  implicit none
+
+  integer(bit_kind), intent(in) :: det
+
+  double precision, intent(out) :: energy
+
+  call i_H_j(det,det,N_int,energy)
+
+  energy = energy + nuclear_repulsion
+  
+end  
+#+end_src

src/utils_cc/org/guess_t.org

@@ -0,0 +1,222 @@
+* Guess
+** T1
+#+begin_src f90 :comments org :tangle guess_t.irp.f
+subroutine guess_t1(nO,nV,f_o,f_v,f_ov,t1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Update the T1 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: f_o(nO), f_v(nV), f_ov(nO,nV)
+
+  ! inout
+  double precision, intent(out) :: t1(nO, nV)
+
+  ! internal
+  integer                       :: i,a
+
+  if (trim(cc_guess_t1) == 'none') then
+     t1 = 0d0
+  else if (trim(cc_guess_t1) == 'MP') then
+    do a = 1, nV
+      do i = 1, nO
+        t1(i,a) = f_ov(i,a) / (f_o(i) - f_v(a) - cc_level_shift_guess)
+      enddo
+    enddo
+  else if (trim(cc_guess_t1) == 'read') then
+    call read_t1(nO,nV,t1)
+  else
+    print*, 'Unknown cc_guess_t1 type: '//trim(cc_guess_t1)
+    call abort
+  endif
+  
+end  
+#+end_src
+
+** T2
+#+begin_src f90 :comments org :tangle guess_t.irp.f
+subroutine guess_t2(nO,nV,f_o,f_v,v_oovv,t2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Update the T2 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)           :: nO, nV
+  double precision, intent(in)  :: f_o(nO), f_v(nV), v_oovv(nO, nO, nV, nV)
+
+  ! inout
+  double precision, intent(out) :: t2(nO, nO, nV, nV)
+
+  ! internal
+  integer                       :: i,j,a,b
+
+  if (trim(cc_guess_t2) == 'none') then
+    t2 = 0d0
+  else if (trim(cc_guess_t2) == 'MP') then
+    do b = 1, nV
+      do a = 1, nV
+        do j = 1, nO
+          do i = 1, nO
+            t2(i,j,a,b) = v_oovv(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift_guess)
+          enddo
+        enddo
+      enddo
+    enddo
+  else if (trim(cc_guess_t2) == 'read') then
+    call read_t2(nO,nV,t2)
+  else
+    print*, 'Unknown cc_guess_t1 type: '//trim(cc_guess_t2)
+    call abort
+  endif
+  
+end  
+#+end_src
+
+* Write
+** T1
+#+begin_src f90 :comments org :tangle guess_t.irp.f
+subroutine write_t1(nO,nV,t1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Write the T1 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)          :: nO, nV
+  double precision, intent(in) :: t1(nO, nV)
+
+  ! internal
+  integer                      :: i,a
+
+  if (cc_write_t1) then
+    open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
+    do a = 1, nV
+      do i = 1, nO
+         write(11,'(F20.12)') t1(i,a)
+      enddo
+    enddo
+    close(11)
+  endif
+  
+end  
+#+end_src
+
+** T2
+#+begin_src f90 :comments org :tangle guess_t.irp.f
+subroutine write_t2(nO,nV,t2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Write the T2 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)          :: nO, nV
+  double precision, intent(in) :: t2(nO, nO, nV, nV)
+
+  ! internal
+  integer                      :: i,j,a,b
+
+  if (cc_write_t2) then
+    open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
+    do b = 1, nV
+      do a = 1, nV
+        do j = 1, nO
+          do i = 1, nO
+             write(11,'(F20.12)') t2(i,j,a,b)
+          enddo
+        enddo
+      enddo
+    enddo
+    close(11)
+  endif
+  
+end  
+#+end_src
+
+* Read
+** T1
+#+begin_src f90 :comments org :tangle guess_t.irp.f
+subroutine read_t1(nO,nV,t1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Read the T1 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)           :: nO, nV
+  double precision, intent(out) :: t1(nO, nV)
+
+  ! internal
+  integer                       :: i,a
+  logical                       :: ok
+
+  inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
+  if (.not. ok) then
+     print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
+     print*, 'Do a first calculation with cc_write_t1 = True'
+     print*, 'and cc_guess_t1 /= read before setting cc_guess_t1 = read'
+     call abort
+  endif
+  open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
+  do a = 1, nV
+    do i = 1, nO
+       read(11,'(F20.12)') t1(i,a)
+    enddo
+  enddo
+  close(11)
+  
+end  
+#+end_src
+
+** T2
+#+begin_src f90 :comments org :tangle guess_t.irp.f
+subroutine read_t2(nO,nV,t2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Read the T2 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)           :: nO, nV
+  double precision, intent(out) :: t2(nO, nO, nV, nV)
+
+  ! internal
+  integer                       :: i,j,a,b
+  logical                       :: ok
+
+  inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
+  if (.not. ok) then
+     print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
+     print*, 'Do a first calculation with cc_write_t2 = True'
+     print*, 'and cc_guess_t2 /= read before setting cc_guess_t2 = read'
+     call abort
+  endif
+  open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
+  do b = 1, nV
+    do a = 1, nV
+      do j = 1, nO
+        do i = 1, nO
+           read(11,'(F20.12)') t2(i,j,a,b)
+        enddo
+      enddo
+    enddo
+  enddo
+  close(11)
+  
+end  
+#+end_src

src/utils_cc/org/occupancy.org

@@ -0,0 +1,341 @@
+* N spin orb
+#+begin_src f90 :comments org :tangle occupancy.irp.f
+subroutine extract_n_spin(det,n)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Returns the number of occupied alpha, occupied beta, virtual alpha, virtual beta spin orbitals
+  ! in det without counting the core and deleted orbitals in the format n(nOa,nOb,nVa,nVb)
+  END_DOC
+
+  integer(bit_kind), intent(in) :: det(N_int,2)
+  
+  integer, intent(out)          :: n(4)
+  
+  integer(bit_kind)             :: res(N_int,2)
+  integer                       :: i, si
+  logical                       :: ok, is_core, is_del
+
+  ! Init
+  n = 0
+
+  ! Loop over the spin
+  do si = 1, 2
+    do i = 1, mo_num
+      call apply_hole(det, si, i, res, ok, N_int)
+      
+      ! in core ?
+      if (is_core(i)) cycle
+      ! in del ?
+      if (is_del(i)) cycle
+      
+      if (ok) then
+        ! particle
+        n(si) = n(si) + 1
+      else
+        ! hole
+        n(si+2) = n(si+2) + 1
+      endif
+    enddo
+  enddo
+
+  !print*,n(1),n(2),n(3),n(4)
+
+end  
+#+end_src
+
+* List_orb
+** Spin
+#+begin_src f90 :comments org :tangle occupancy.irp.f
+subroutine extract_list_orb_spin(det,nO_m,nV_m,list_occ,list_vir)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Returns the the list of occupied alpha/beta, virtual alpha/beta spin orbitals
+  ! size(nO_m,1) must be max(nOa,nOb) and size(nV_m,1) must be max(nVa,nVb)
+  END_DOC
+  
+  integer, intent(in)           :: nO_m, nV_m
+  integer(bit_kind), intent(in) :: det(N_int,2)
+  
+  integer, intent(out)          :: list_occ(nO_m,2), list_vir(nV_m,2)
+  
+  integer(bit_kind)             :: res(N_int,2)
+  integer                       :: i, si, idx_o, idx_v, idx_i, idx_b
+  logical                       :: ok, is_core, is_del
+
+  list_occ = 0
+  list_vir = 0
+
+  ! List of occ/vir alpha/beta
+
+  ! occ alpha -> list_occ(:,1)
+  ! occ beta -> list_occ(:,2)
+  ! vir alpha -> list_vir(:,1)
+  ! vir beta -> list_vir(:,2)
+
+  ! Loop over the spin 
+  do si = 1, 2
+    ! tmp idx
+    idx_o = 1
+    idx_v = 1
+    do i = 1, mo_num
+      call apply_hole(det, si, i, res, ok, N_int)
+
+      ! in core ?
+      if (is_core(i)) cycle
+      ! in del ?
+      if (is_del(i)) cycle
+            
+      if (ok) then
+        ! particle
+        list_occ(idx_o,si) = i
+        idx_o = idx_o + 1
+      else
+        ! hole
+        list_vir(idx_v,si) = i
+        idx_v = idx_v + 1
+      endif
+    enddo
+  enddo
+
+end
+#+end_src
+
+** Space
+#+begin_src f90 :comments org :tangle occupancy.irp.f
+subroutine extract_list_orb_space(det,nO,nV,list_occ,list_vir)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Returns the the list of occupied and virtual alpha spin orbitals
+  END_DOC
+  
+  integer, intent(in)           :: nO, nV
+  integer(bit_kind), intent(in) :: det(N_int,2)
+  
+  integer, intent(out)          :: list_occ(nO), list_vir(nV)
+  
+  integer(bit_kind)             :: res(N_int,2)
+  integer                       :: i, si, idx_o, idx_v, idx_i, idx_b
+  logical                       :: ok, is_core, is_del
+  
+  if (elec_alpha_num /= elec_beta_num) then
+    print*,'Error elec_alpha_num /= elec_beta_num, impossible to create cc_list_occ and cc_list_vir, abort'
+    call abort
+  endif
+
+  list_occ = 0
+  list_vir = 0
+
+  ! List of occ/vir alpha
+
+  ! occ alpha -> list_occ(:,1)
+  ! vir alpha -> list_vir(:,1)
+
+  ! tmp idx
+  idx_o = 1
+  idx_v = 1
+  do i = 1, mo_num
+    call apply_hole(det, 1, i, res, ok, N_int)
+
+    ! in core ?
+    if (is_core(i)) cycle
+    ! in del ?
+    if (is_del(i)) cycle
+
+    if (ok) then
+      ! particle
+      list_occ(idx_o) = i
+      idx_o = idx_o + 1
+    else
+      ! hole
+      list_vir(idx_v) = i
+      idx_v = idx_v + 1
+    endif
+  enddo
+
+end
+#+end_src
+
+** is_core
+#+begin_src f90 :comments org :tangle occupancy.irp.f
+function is_core(i)
+
+  implicit none
+
+  BEGIN_DOC
+  ! True if the orbital i is a core orbital
+  END_DOC
+
+  integer, intent(in) :: i
+  logical             :: is_core
+
+  integer             :: j
+
+  ! Init
+  is_core = .False.
+
+  ! Search
+  do j = 1, dim_list_core_orb
+    if (list_core(j) == i) then
+      is_core = .True.
+      exit
+    endif
+  enddo
+
+end
+#+end_src
+
+** is_del
+#+begin_src f90 :comments org :tangle occupancy.irp.f
+function is_del(i)
+
+  implicit none
+
+  BEGIN_DOC
+  ! True if the orbital i is a deleted orbital
+  END_DOC
+
+  integer, intent(in) :: i
+  logical             :: is_del
+
+  integer             :: j
+
+  ! Init
+  is_del = .False.
+
+  ! Search
+  do j = 1, dim_list_core_orb
+    if (list_core(j) == i) then
+      is_del = .True.
+      exit
+    endif
+  enddo
+
+end
+#+end_src
+
+* Providers
+** N orb
+#+BEGIN_SRC f90 :comments org :tangle occupancy.irp.f
+ BEGIN_PROVIDER [integer, cc_nO_m]
+&BEGIN_PROVIDER [integer, cc_nOa]
+&BEGIN_PROVIDER [integer, cc_nOb]
+&BEGIN_PROVIDER [integer, cc_nOab]
+&BEGIN_PROVIDER [integer, cc_nV_m]
+&BEGIN_PROVIDER [integer, cc_nVa]
+&BEGIN_PROVIDER [integer, cc_nVb]
+&BEGIN_PROVIDER [integer, cc_nVab]
+&BEGIN_PROVIDER [integer, cc_n_mo]
+&BEGIN_PROVIDER [integer, cc_nO_S, (2)]
+&BEGIN_PROVIDER [integer, cc_nV_S, (2)]
+
+  implicit none
+
+  BEGIN_DOC
+  ! Number of orbitals without core and deleted ones of the cc_ref det in psi_det
+  ! a: alpha, b: beta
+  ! nO_m: max(a,b) occupied 
+  ! nOa: nb a occupied 
+  ! nOb: nb b occupied 
+  ! nOab: nb a+b occupied 
+  ! nV_m: max(a,b) virtual 
+  ! nVa: nb a virtual 
+  ! nVb: nb b virtual 
+  ! nVab: nb a+b virtual 
+  END_DOC
+
+  integer :: n_spin(4)
+
+  ! Extract number of occ/vir alpha/beta spin orbitals
+  call extract_n_spin(psi_det(1,1,cc_ref),n_spin)
+
+  cc_nOa  = n_spin(1)
+  cc_nOb  = n_spin(2)
+  cc_nOab = cc_nOa + cc_nOb    !n_spin(1) + n_spin(2)
+  cc_nO_m = max(cc_nOa,cc_nOb) !max(n_spin(1), n_spin(2))
+  cc_nVa  = n_spin(3)
+  cc_nVb  = n_spin(4)
+  cc_nVab = cc_nVa + cc_nVb    !n_spin(3) + n_spin(4)
+  cc_nV_m = max(cc_nVa,cc_nVb) !max(n_spin(3), n_spin(4))
+  cc_n_mo = cc_nVa + cc_nVb    !n_spin(1) + n_spin(3)
+  cc_nO_S = (/cc_nOa,cc_nOb/)
+  cc_nV_S = (/cc_nVa,cc_nVb/)
+
+END_PROVIDER
+#+end_src
+
+** List orb
+
+*** General
+#+BEGIN_SRC f90 :comments org :tangle occupancy.irp.f
+ BEGIN_PROVIDER [integer, cc_list_gen, (cc_n_mo)]
+
+  implicit none
+
+  BEGIN_DOC
+  ! List of general orbitals without core and deleted ones
+  END_DOC
+
+  integer :: i,j
+  logical :: is_core, is_del
+  
+  j = 1
+  do i = 1, mo_num
+    ! in core ?
+    if (is_core(i)) cycle
+    ! in del ?
+    if (is_del(i)) cycle
+    cc_list_gen(j) = i
+    j = j+1
+  enddo
+
+END_PROVIDER
+#+end_src
+
+*** Space
+#+BEGIN_SRC f90 :comments org :tangle occupancy.irp.f
+ BEGIN_PROVIDER [integer, cc_list_occ, (cc_nOa)]
+&BEGIN_PROVIDER [integer, cc_list_vir, (cc_nVa)]
+
+  implicit none
+
+  BEGIN_DOC
+  ! List of occupied and virtual spatial orbitals without core and deleted ones
+  END_DOC
+
+  call extract_list_orb_space(psi_det(1,1,cc_ref),cc_nOa,cc_nVa,cc_list_occ,cc_list_vir)
+
+END_PROVIDER
+#+end_src
+
+*** Spin
+#+BEGIN_SRC f90 :comments org :tangle occupancy.irp.f
+ BEGIN_PROVIDER [integer, cc_list_occ_spin, (cc_nO_m,2)]
+&BEGIN_PROVIDER [integer, cc_list_vir_spin, (cc_nV_m,2)]
+&BEGIN_PROVIDER [logical, cc_ref_is_open_shell]
+
+  implicit none
+
+  BEGIN_DOC
+  ! List of occupied and virtual spin orbitals without core and deleted ones
+  END_DOC
+
+  integer :: i
+
+  call extract_list_orb_spin(psi_det(1,1,cc_ref),cc_nO_m,cc_nV_m,cc_list_occ_spin,cc_list_vir_spin)
+
+  cc_ref_is_open_shell = .False.
+  do i = 1, cc_nO_m
+    if (cc_list_occ_spin(i,1) /= cc_list_occ_spin(i,2)) then
+       cc_ref_is_open_shell = .True.
+    endif
+  enddo
+
+
+END_PROVIDER
+#+end_src

src/utils_cc/org/phase.org

@@ -0,0 +1,178 @@
+#+begin_src f90 :comments org :notangle phase.irp.f
+program run
+  implicit none
+
+  integer :: n(2), degree1, degree2, exc(0:2,2,2)
+  integer, allocatable :: list_anni(:,:), list_crea(:,:)
+  double precision :: phase1, phase2
+  integer :: h1,h2,p1,p2,s1,s2,i,j
+
+  allocate(list_anni(N_int*bit_kind_size,2))
+  allocate(list_crea(N_int*bit_kind_size,2))
+
+  do i = 1, N_det-1
+    do j = i+1, N_det
+      !call print_det(psi_det(1,1,j),N_int)
+      call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree1,phase1,N_int)
+      call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2)
+      !print*,'old',degree1,phase1
+      !print*,'h1:',h1,'h2:',h2,'s1:',s1,'s2:',s2
+      !print*,'p1:',p1,'p2:',p2
+      call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree1,N_int)
+      call get_excitation_general(psi_det(1,1,i),psi_det(1,1,j),degree2,n,list_anni,list_crea,phase2,N_int)
+      !print*,'new',degree2,phase2
+      !print*,'ha:',list_anni(1:n(1),1),'hb',list_anni(1:n(2),2)
+      !print*,'pa:',list_crea(1:n(1),1),'pb',list_crea(1:n(2),2)
+      !print*,''
+      if (degree1 /= degree2) then
+        print*,'Error degree:',degree1,degree2
+        call abort
+      endif
+      if (degree1 <= 2 .and. phase1 /= phase2) then
+        print*,'Error phase',phase1,phase2
+        call abort
+      endif
+    enddo
+  enddo
+  
+end
+#+end_src
+
+** phase
+#+begin_src f90 :comments org :tangle phase.irp.f
+subroutine get_phase_general(det1,det2,phase,degree,Nint)
+  implicit none
+
+  integer, intent(in)           :: Nint
+  integer(bit_kind), intent(in) :: det1(Nint,2), det2(Nint,2)
+  double precision, intent(out) :: phase
+  integer, intent(out)          :: degree
+  integer :: n(2)
+  integer, allocatable :: list_anni(:,:), list_crea(:,:)
+
+  allocate(list_anni(N_int*bit_kind_size,2))
+  allocate(list_crea(N_int*bit_kind_size,2))
+
+  call get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,Nint)
+end
+#+end_src
+
+** Get excitation general
+#+begin_src f90 :comments org :tangle phase.irp.f
+subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,Nint)
+
+  use bitmasks
+
+  implicit none
+
+  integer, intent(in) :: Nint
+  integer(bit_kind), intent(in)  :: det1(Nint,2), det2(Nint,2)
+  double precision, intent(out)  :: phase
+  integer, intent(out)           :: list_crea(Nint*bit_kind_size,2)
+  integer, intent(out)           :: list_anni(Nint*bit_kind_size,2)
+  integer, intent(out)           :: degree, n(2)
+  
+  integer, allocatable           :: l1(:,:), l2(:,:) 
+  integer(bit_kind), allocatable :: det_crea(:,:), det_anni(:,:)
+  integer, allocatable           :: pos_anni(:,:), pos_crea(:,:)
+
+  integer :: n1(2),n2(2),n_crea(2),n_anni(2),i,j,k,d
+
+  allocate(l1(Nint*bit_kind_size,2))
+  allocate(l2(Nint*bit_kind_size,2))
+  allocate(det_crea(Nint,2),det_anni(Nint,2))
+
+  ! 1      111010
+  ! 2      110101
+  !
+  !not 1-> 000101
+  !    2   110101
+  !and     000101 -> crea
+  !
+  !    1   111010
+  !not 2-> 001010
+  !        001010 -> anni
+
+  do j = 1, 2
+    do i = 1, Nint
+      det_crea(i,j) = iand(not(det1(i,j)),det2(i,j))
+    enddo
+  enddo
+  
+  do j = 1, 2
+    do i = 1, Nint
+      det_anni(i,j) = iand(det1(i,j),not(det2(i,j)))
+    enddo
+  enddo
+
+  call bitstring_to_list_ab(det1,l1,n1,Nint)
+  call bitstring_to_list_ab(det2,l2,n2,Nint)
+  call bitstring_to_list_ab(det_crea,list_crea,n_crea,Nint)
+  call bitstring_to_list_ab(det_anni,list_anni,n_anni,Nint)
+
+  do i = 1, 2
+    if (n_crea(i) /= n_anni(i)) then
+      print*,'Well, it seems we have a problem here...'
+      call abort
+    endif
+  enddo
+
+  !1    11110011001  1 2 3 4 7 8  11
+  !pos               1 2 3 4 5 6  7 
+  !2    11100101011  1 2 3 6 8 10 11
+  !anni 00010010000  4 7
+  !pos               4 5
+  !crea 00000100010  6 10
+  !pos               4 6
+  !4 -> 6  pos(4 -> 4)
+  !7 -> 10 pos(5 -> 6)
+
+  n = n_anni
+  degree = n_anni(1) + n_anni(2)
+
+  allocate(pos_anni(max(n(1),n(2)),2))
+  allocate(pos_crea(max(n(1),n(2)),2))
+  
+  ! Search pos anni
+  do j = 1, 2
+    k = 1
+    do i = 1, n1(j)
+       if (l1(i,j) /= list_anni(k,j)) cycle
+       pos_anni(k,j) = i
+       k = k + 1
+    enddo
+  enddo
+
+  ! Search pos crea
+  do j = 1, 2
+    k = 1
+    do i = 1, n2(j)
+       if (l2(i,j) /= list_crea(k,j)) cycle
+       pos_crea(k,j) = i
+       k = k + 1
+    enddo
+  enddo
+
+  ! Distance between the ith anni and the ith crea op
+  ! By doing so there is no crossing between the different pairs of anni/crea
+  ! and the phase is determined by the sum of the distances
+  ! -> (-1)^{sum of the distances}
+  d = 0
+  do j = 1, 2
+    do i = 1, n(j)
+      d = d + abs(pos_anni(i,j) - pos_crea(i,j))
+    enddo
+  enddo
+  
+  phase = dble((-1)**d)
+
+  ! Debug
+  !print*,l2(1:n2(1),1)
+  !print*,l2(1:n2(2),2)
+  !!call print_det(det1,Nint)
+  !!call print_det(det2,Nint)
+  !print*,phase
+  !print*,''
+end
+#+end_src
+

src/utils_cc/org/print_wf_qp_edit.org

@@ -0,0 +1,33 @@
+#+begin_src f90 :comments org :tangle print_wf_qp_edit.irp.f
+program run
+
+  implicit none
+
+  read_wf = .true.
+  touch read_wf
+
+  call print_wf_qp_edit()
+
+end 
+#+end_src
+
+#+begin_src f90 :comments org :tangle print_wf_qp_edit.irp.f
+subroutine print_wf_qp_edit()
+  
+  implicit none
+
+  BEGIN_DOC
+  ! Print the psi_det wave function up to n_det_qp_edit
+  END_DOC
+  
+  integer :: i
+
+  do i = 1, n_det_qp_edit
+    print*,i
+    write(*,'(100(1pE12.4))') psi_coef(i,:)
+    call print_det(psi_det(1,1,i),N_int)
+    print*,''
+  enddo
+
+end
+#+end_src

src/utils_cc/org/update_t.org

@@ -0,0 +1,76 @@
+* T1
+#+begin_src f90 :comments org :tangle update_t.irp.f
+subroutine update_t1(nO,nV,f_o,f_v,r1,t1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Update the T1 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)             :: nO, nV
+  double precision, intent(in)    :: f_o(nO), f_v(nV), r1(nO, nV)
+
+  ! inout
+  double precision, intent(inout) :: t1(nO, nV)
+
+  ! internal
+  integer                         :: i,a
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,t1,r1,cc_level_shift,f_o,f_v) &
+  !$OMP PRIVATE(i,a) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(1)
+  do a = 1, nV
+    do i = 1, nO
+      t1(i,a) = t1(i,a) - r1(i,a) / (f_o(i) - f_v(a) - cc_level_shift)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+end  
+#+end_src
+
+* T2
+#+begin_src f90 :comments org :tangle update_t.irp.f
+subroutine update_t2(nO,nV,f_o,f_v,r2,t2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Update the T2 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)             :: nO, nV
+  double precision, intent(in)    :: f_o(nO), f_v(nV), r2(nO, nO, nV, nV)
+
+  ! inout
+  double precision, intent(inout) :: t2(nO, nO, nV, nV)
+
+  ! internal
+  integer                         :: i,j,a,b
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,t2,r2,cc_level_shift,f_o,f_v) &
+  !$OMP PRIVATE(i,j,a,b) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(3)
+  do b = 1, nV
+    do a = 1, nV
+      do j = 1, nO
+        do i = 1, nO
+          t2(i,j,a,b) = t2(i,j,a,b) - r2(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+end  
+#+end_src
+

src/utils_cc/phase.irp.f

@@ -0,0 +1,135 @@
+! phase
+
+subroutine get_phase_general(det1,det2,phase,degree,Nint)
+  implicit none
+
+  integer, intent(in)           :: Nint
+  integer(bit_kind), intent(in) :: det1(Nint,2), det2(Nint,2)
+  double precision, intent(out) :: phase
+  integer, intent(out)          :: degree
+  integer :: n(2)
+  integer, allocatable :: list_anni(:,:), list_crea(:,:)
+
+  allocate(list_anni(N_int*bit_kind_size,2))
+  allocate(list_crea(N_int*bit_kind_size,2))
+
+  call get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,Nint)
+end
+
+! Get excitation general
+
+subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,Nint)
+
+  use bitmasks
+
+  implicit none
+
+  integer, intent(in) :: Nint
+  integer(bit_kind), intent(in)  :: det1(Nint,2), det2(Nint,2)
+  double precision, intent(out)  :: phase
+  integer, intent(out)           :: list_crea(Nint*bit_kind_size,2)
+  integer, intent(out)           :: list_anni(Nint*bit_kind_size,2)
+  integer, intent(out)           :: degree, n(2)
+  
+  integer, allocatable           :: l1(:,:), l2(:,:) 
+  integer(bit_kind), allocatable :: det_crea(:,:), det_anni(:,:)
+  integer, allocatable           :: pos_anni(:,:), pos_crea(:,:)
+
+  integer :: n1(2),n2(2),n_crea(2),n_anni(2),i,j,k,d
+
+  allocate(l1(Nint*bit_kind_size,2))
+  allocate(l2(Nint*bit_kind_size,2))
+  allocate(det_crea(Nint,2),det_anni(Nint,2))
+
+  ! 1      111010
+  ! 2      110101
+  !
+  !not 1-> 000101
+  !    2   110101
+  !and     000101 -> crea
+  !
+  !    1   111010
+  !not 2-> 001010
+  !        001010 -> anni
+
+  do j = 1, 2
+    do i = 1, Nint
+      det_crea(i,j) = iand(not(det1(i,j)),det2(i,j))
+    enddo
+  enddo
+  
+  do j = 1, 2
+    do i = 1, Nint
+      det_anni(i,j) = iand(det1(i,j),not(det2(i,j)))
+    enddo
+  enddo
+
+  call bitstring_to_list_ab(det1,l1,n1,Nint)
+  call bitstring_to_list_ab(det2,l2,n2,Nint)
+  call bitstring_to_list_ab(det_crea,list_crea,n_crea,Nint)
+  call bitstring_to_list_ab(det_anni,list_anni,n_anni,Nint)
+
+  do i = 1, 2
+    if (n_crea(i) /= n_anni(i)) then
+      print*,'Well, it seems we have a problem here...'
+      call abort
+    endif
+  enddo
+
+  !1    11110011001  1 2 3 4 7 8  11
+  !pos               1 2 3 4 5 6  7 
+  !2    11100101011  1 2 3 6 8 10 11
+  !anni 00010010000  4 7
+  !pos               4 5
+  !crea 00000100010  6 10
+  !pos               4 6
+  !4 -> 6  pos(4 -> 4)
+  !7 -> 10 pos(5 -> 6)
+
+  n = n_anni
+  degree = n_anni(1) + n_anni(2)
+
+  allocate(pos_anni(max(n(1),n(2)),2))
+  allocate(pos_crea(max(n(1),n(2)),2))
+  
+  ! Search pos anni
+  do j = 1, 2
+    k = 1
+    do i = 1, n1(j)
+       if (l1(i,j) /= list_anni(k,j)) cycle
+       pos_anni(k,j) = i
+       k = k + 1
+    enddo
+  enddo
+
+  ! Search pos crea
+  do j = 1, 2
+    k = 1
+    do i = 1, n2(j)
+       if (l2(i,j) /= list_crea(k,j)) cycle
+       pos_crea(k,j) = i
+       k = k + 1
+    enddo
+  enddo
+
+  ! Distance between the ith anni and the ith crea op
+  ! By doing so there is no crossing between the different pairs of anni/crea
+  ! and the phase is determined by the sum of the distances
+  ! -> (-1)^{sum of the distances}
+  d = 0
+  do j = 1, 2
+    do i = 1, n(j)
+      d = d + abs(pos_anni(i,j) - pos_crea(i,j))
+    enddo
+  enddo
+  
+  phase = dble((-1)**d)
+
+  ! Debug
+  !print*,l2(1:n2(1),1)
+  !print*,l2(1:n2(2),2)
+  !!call print_det(det1,Nint)
+  !!call print_det(det2,Nint)
+  !print*,phase
+  !print*,''
+end

src/utils_cc/print_wf_qp_edit.irp.f

@@ -0,0 +1,29 @@
+program run
+
+  implicit none
+
+  read_wf = .true.
+  touch read_wf
+
+  call print_wf_qp_edit()
+
+end
+
+subroutine print_wf_qp_edit()
+  
+  implicit none
+
+  BEGIN_DOC
+  ! Print the psi_det wave function up to n_det_qp_edit
+  END_DOC
+  
+  integer :: i
+
+  do i = 1, n_det_qp_edit
+    print*,i
+    write(*,'(100(1pE12.4))') psi_coef(i,:)
+    call print_det(psi_det(1,1,i),N_int)
+    print*,''
+  enddo
+
+end

src/utils_cc/update_t.irp.f

@@ -0,0 +1,73 @@
+! T1
+
+subroutine update_t1(nO,nV,f_o,f_v,r1,t1)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Update the T1 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)             :: nO, nV
+  double precision, intent(in)    :: f_o(nO), f_v(nV), r1(nO, nV)
+
+  ! inout
+  double precision, intent(inout) :: t1(nO, nV)
+
+  ! internal
+  integer                         :: i,a
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,t1,r1,cc_level_shift,f_o,f_v) &
+  !$OMP PRIVATE(i,a) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(1)
+  do a = 1, nV
+    do i = 1, nO
+      t1(i,a) = t1(i,a) - r1(i,a) / (f_o(i) - f_v(a) - cc_level_shift)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+end
+
+! T2
+
+subroutine update_t2(nO,nV,f_o,f_v,r2,t2)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Update the T2 amplitudes for CC
+  END_DOC
+
+  ! in
+  integer, intent(in)             :: nO, nV
+  double precision, intent(in)    :: f_o(nO), f_v(nV), r2(nO, nO, nV, nV)
+
+  ! inout
+  double precision, intent(inout) :: t2(nO, nO, nV, nV)
+
+  ! internal
+  integer                         :: i,j,a,b
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(nO,nV,t2,r2,cc_level_shift,f_o,f_v) &
+  !$OMP PRIVATE(i,j,a,b) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO collapse(3)
+  do b = 1, nV
+    do a = 1, nV
+      do j = 1, nO
+        do i = 1, nO
+          t2(i,j,a,b) = t2(i,j,a,b) - r2(i,j,a,b) / (f_o(i) + f_o(j) - f_v(a) - f_v(b) - cc_level_shift)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  
+end