diff --git a/src/casscf/50.casscf.bats b/src/casscf/50.casscf.bats
deleted file mode 100644
index a0db725d..00000000
--- a/src/casscf/50.casscf.bats
+++ /dev/null
@@ -1,49 +0,0 @@
-#!/usr/bin/env bats
-
-source $QP_ROOT/tests/bats/common.bats.sh
-source $QP_ROOT/quantum_package.rc
-
-
-function run_stoch() {
-  thresh=$2
-  test_exe casscf || skip
-  qp set perturbation do_pt2 True
-  qp set determinants n_det_max $3
-  qp set davidson threshold_davidson 1.e-10
-  qp set davidson n_states_diag 4
-  qp run casscf | tee casscf.out 
-  energy1="$(ezfio get casscf energy_pt2 | tr '[]' ' ' | cut -d ',' -f 1)"
-  eq $energy1 $1 $thresh
-}
-
-@test "F2" { # 18.0198s
-  rm -rf f2_casscf
-  qp_create_ezfio -b aug-cc-pvdz ../input/f2.zmt  -o f2_casscf
-  qp set_file f2_casscf
-  qp run scf 
-  qp set_mo_class --core="[1-6,8-9]" --act="[7,10]" --virt="[11-46]"
-  run_stoch  -198.773366970 1.e-4  100000
-}
-
-@test "N2" { # 18.0198s
-  rm -rf n2_casscf
-  qp_create_ezfio -b aug-cc-pvdz ../input/n2.xyz  -o n2_casscf
-  qp set_file n2_casscf
-  qp run scf 
-  qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]"
-  run_stoch  -109.0961643162 1.e-4  100000
-}
-
-@test "N2_stretched" {
-  rm -rf n2_stretched_casscf
-  qp_create_ezfio -b aug-cc-pvdz -m 7 ../input/n2_stretched.xyz  -o n2_stretched_casscf
-  qp set_file n2_stretched_casscf 
-  qp run scf  | tee scf.out 
-  qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]"
-  qp set electrons elec_alpha_num 7 
-  qp set electrons elec_beta_num 7 
-  run_stoch  -108.7860471300 1.e-4  100000
-#   
-
-}
-
diff --git a/src/casscf/EZFIO.cfg b/src/casscf/EZFIO.cfg
deleted file mode 100644
index 4e4d3d3a..00000000
--- a/src/casscf/EZFIO.cfg
+++ /dev/null
@@ -1,31 +0,0 @@
-[energy]
-type: double precision
-doc: Calculated Selected |FCI| energy
-interface: ezfio
-size: (determinants.n_states)
-
-[energy_pt2]
-type: double precision
-doc: Calculated |FCI| energy + |PT2|
-interface: ezfio
-size: (determinants.n_states)
-
-[cisd_guess]
-type: logical
-doc: If true, the CASSCF starts with a CISD wave function
-interface: ezfio,provider,ocaml
-default: True 
-
-[state_following_casscf]
-type: logical                                                                                                                                
-doc: If |true|, the CASSCF will try to follow the guess CI vector and orbitals
-interface: ezfio,provider,ocaml
-default: False
-
-
-[level_shift_casscf]
-type: Positive_float
-doc: Energy shift on the virtual MOs to improve SCF convergence
-interface: ezfio,provider,ocaml
-default: 0.005
-
diff --git a/src/casscf/MORALITY b/src/casscf/MORALITY
deleted file mode 100644
index 9701a647..00000000
--- a/src/casscf/MORALITY
+++ /dev/null
@@ -1 +0,0 @@
-the CASCF can be obtained if a proper guess is given to the WF part
diff --git a/src/casscf/NEED b/src/casscf/NEED
deleted file mode 100644
index d9da718e..00000000
--- a/src/casscf/NEED
+++ /dev/null
@@ -1,4 +0,0 @@
-cipsi
-selectors_full
-generators_cas
-two_body_rdm
diff --git a/src/casscf/README.rst b/src/casscf/README.rst
deleted file mode 100644
index 08bfd95b..00000000
--- a/src/casscf/README.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-======
-casscf
-======
-
-|CASSCF| program with the CIPSI algorithm.
diff --git a/src/casscf/bavard.irp.f b/src/casscf/bavard.irp.f
deleted file mode 100644
index 463c3ea4..00000000
--- a/src/casscf/bavard.irp.f
+++ /dev/null
@@ -1,6 +0,0 @@
-! -*- F90 -*- 
-BEGIN_PROVIDER [logical, bavard]
-!  bavard=.true.
-   bavard=.false.
-END_PROVIDER
-
diff --git a/src/casscf/bielec.irp.f b/src/casscf/bielec.irp.f
deleted file mode 100644
index 0a44f994..00000000
--- a/src/casscf/bielec.irp.f
+++ /dev/null
@@ -1,155 +0,0 @@
-BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
-  BEGIN_DOC
-  ! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active
-  ! indices are unshifted orbital numbers
-  END_DOC
-  implicit none
-  integer                        :: i,j,ii,jj,p,q,i3,j3,t3,v3
-  real*8                         :: mo_two_e_integral
-  
-  bielec_PQxx(:,:,:,:) = 0.d0
-  PROVIDE mo_two_e_integrals_in_map
-  
-  !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP PRIVATE(i,ii,j,jj,i3,j3) &
-  !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx, &
-  !$OMP  n_act_orb,mo_integrals_map,list_act)
-
-  !$OMP DO
-  do i=1,n_core_inact_orb
-    ii=list_core_inact(i)
-    do j=i,n_core_inact_orb
-      jj=list_core_inact(j)
-      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j),mo_integrals_map)
-      bielec_PQxx(:,:,j,i)=bielec_PQxx(:,:,i,j)
-    end do
-    do j=1,n_act_orb
-      jj=list_act(j)
-      j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j3),mo_integrals_map)
-      bielec_PQxx(:,:,j3,i)=bielec_PQxx(:,:,i,j3)
-    end do
-  end do
-  !$OMP END DO
-
-
-  !$OMP DO
-  do i=1,n_act_orb
-    ii=list_act(i)
-    i3=i+n_core_inact_orb
-    do j=i,n_act_orb
-      jj=list_act(j)
-      j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i3,j3),mo_integrals_map)
-      bielec_PQxx(:,:,j3,i3)=bielec_PQxx(:,:,i3,j3)
-    end do
-  end do
-  !$OMP END DO
-
-  !$OMP END PARALLEL
-
-END_PROVIDER
-
-
-
-BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
-  BEGIN_DOC
-  ! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
-  ! indices are unshifted orbital numbers
-  END_DOC
-  implicit none
-  integer                        :: i,j,ii,jj,p,q,i3,j3,t3,v3
-  double precision, allocatable  :: integrals_array(:,:)
-  real*8                         :: mo_two_e_integral
-  
-  PROVIDE mo_two_e_integrals_in_map
-  bielec_PxxQ = 0.d0
-
-  !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP PRIVATE(i,ii,j,jj,i3,j3,integrals_array) &
-  !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ, &
-  !$OMP  n_act_orb,mo_integrals_map,list_act)
-
-  allocate(integrals_array(mo_num,mo_num))
-  
-  !$OMP DO
-  do i=1,n_core_inact_orb
-    ii=list_core_inact(i)
-    do j=i,n_core_inact_orb
-      jj=list_core_inact(j)
-      call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
-      do q=1,mo_num
-        do p=1,mo_num
-          bielec_PxxQ(p,i,j,q)=integrals_array(p,q)
-          bielec_PxxQ(p,j,i,q)=integrals_array(q,p)
-        end do
-      end do
-    end do
-    do j=1,n_act_orb
-      jj=list_act(j)
-      j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
-      do q=1,mo_num
-        do p=1,mo_num
-          bielec_PxxQ(p,i,j3,q)=integrals_array(p,q)
-          bielec_PxxQ(p,j3,i,q)=integrals_array(q,p)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO
-
-
-  ! (ip|qj)
-  !$OMP DO
-  do i=1,n_act_orb
-    ii=list_act(i)
-    i3=i+n_core_inact_orb
-    do j=i,n_act_orb
-      jj=list_act(j)
-      j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
-      do q=1,mo_num
-        do p=1,mo_num
-          bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q)
-          bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO
-
-  deallocate(integrals_array)
-  !$OMP END PARALLEL
-
-END_PROVIDER
-
-
-BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
-  BEGIN_DOC
-  ! bielecCI : integrals (tu|vp) with p arbitrary, tuv active
-  ! index p runs over the whole basis, t,u,v only over the active orbitals
-  END_DOC
-  implicit none
-  integer                        :: i,j,k,p,t,u,v
-  double precision, external     :: mo_two_e_integral
-  PROVIDE mo_two_e_integrals_in_map
-  
-  !$OMP PARALLEL DO DEFAULT(NONE) &
-  !$OMP PRIVATE(i,j,k,p,t,u,v) &
-  !$OMP SHARED(mo_num,n_act_orb,list_act,bielecCI)
-  do p=1,mo_num
-    do j=1,n_act_orb
-      u=list_act(j)
-      do k=1,n_act_orb
-        v=list_act(k)
-        do i=1,n_act_orb
-          t=list_act(i)
-          bielecCI(i,k,j,p) = mo_two_e_integral(t,u,v,p)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END PARALLEL DO
-
-END_PROVIDER
diff --git a/src/casscf/bielec_natorb.irp.f b/src/casscf/bielec_natorb.irp.f
deleted file mode 100644
index 9968530c..00000000
--- a/src/casscf/bielec_natorb.irp.f
+++ /dev/null
@@ -1,369 +0,0 @@
- BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
-  BEGIN_DOC
-  ! integral (pq|xx) in the basis of natural MOs
-  ! indices are unshifted orbital numbers
-  END_DOC
-  implicit none
-  integer                        :: i,j,k,l,t,u,p,q
-  double precision, allocatable  :: f(:,:,:), d(:,:,:)
-
-
-
-  !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP PRIVATE(j,k,l,p,d,f) &
-  !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
-  !$OMP   bielec_PQxx_no,bielec_PQxx,list_act,natorbsCI)
-
-  allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), &
-      d(n_act_orb,mo_num,n_core_inact_act_orb))
-
-  !$OMP DO
-  do l=1,n_core_inact_act_orb
-    bielec_PQxx_no(:,:,:,l) = bielec_PQxx(:,:,:,l)
-
-    do k=1,n_core_inact_act_orb
-      do j=1,mo_num
-        do p=1,n_act_orb
-          f(p,j,k)=bielec_PQxx_no(list_act(p),j,k,l)
-        end do
-      end do
-    end do
-    call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0,  &
-          natorbsCI, size(natorbsCI,1),                              &
-          f, n_act_orb,                                              &
-          0.d0,                                                      &
-          d, n_act_orb)
-    do k=1,n_core_inact_act_orb
-      do j=1,mo_num
-        do p=1,n_act_orb
-          bielec_PQxx_no(list_act(p),j,k,l)=d(p,j,k)
-        end do
-      end do
-
-      do j=1,mo_num
-        do p=1,n_act_orb
-          f(p,j,k)=bielec_PQxx_no(j,list_act(p),k,l)
-        end do
-      end do
-    end do
-    call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0,  &
-          natorbsCI, n_act_orb,                                      &
-          f, n_act_orb,                                              &
-          0.d0,                                                      &
-          d, n_act_orb)
-    do k=1,n_core_inact_act_orb
-      do p=1,n_act_orb
-        do j=1,mo_num
-          bielec_PQxx_no(j,list_act(p),k,l)=d(p,j,k)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO NOWAIT
-
-  deallocate (f,d)
-
-  allocate (f(mo_num,mo_num,n_act_orb),d(mo_num,mo_num,n_act_orb))
-
-  !$OMP DO 
-  do l=1,n_core_inact_act_orb
-
-    do p=1,n_act_orb
-      do k=1,mo_num
-        do j=1,mo_num
-          f(j,k,p) = bielec_PQxx_no(j,k,n_core_inact_orb+p,l)
-        end do
-      end do
-    end do
-    call dgemm('N','N',mo_num*mo_num,n_act_orb,n_act_orb,1.d0,       &
-          f, mo_num*mo_num,                                          &
-          natorbsCI, n_act_orb,                                      &
-          0.d0,                                                      &
-          d, mo_num*mo_num)
-    do p=1,n_act_orb
-      do k=1,mo_num
-        do j=1,mo_num
-          bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,p)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO NOWAIT
-
-  !$OMP BARRIER 
-
-  !$OMP DO
-  do l=1,n_core_inact_act_orb
-    do p=1,n_act_orb
-      do k=1,mo_num
-        do j=1,mo_num
-          f(j,k,p) = bielec_PQxx_no(j,k,l,n_core_inact_orb+p)
-        end do
-      end do
-    end do
-    call dgemm('N','N',mo_num*mo_num,n_act_orb,n_act_orb,1.d0,       &
-          f, mo_num*mo_num,                                          &
-          natorbsCI, n_act_orb,                                      &
-          0.d0,                                                      &
-          d, mo_num*mo_num)
-    do p=1,n_act_orb
-      do k=1,mo_num
-        do j=1,mo_num
-          bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,p)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO
-
-  deallocate (f,d)
-  !$OMP END PARALLEL
-
-END_PROVIDER
-
-
-
-BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
-  BEGIN_DOC
-  ! integral (px|xq) in the basis of natural MOs
-  ! indices are unshifted orbital numbers
-  END_DOC
-  implicit none
-  integer                        :: i,j,k,l,t,u,p,q
-  double precision, allocatable  :: f(:,:,:), d(:,:,:)
-
-  !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP PRIVATE(j,k,l,p,d,f) &
-  !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
-  !$OMP   bielec_PxxQ_no,bielec_PxxQ,list_act,natorbsCI)
-
-
-  allocate (f(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb), &
-      d(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb))
-
-  !$OMP DO
-  do j=1,mo_num
-    bielec_PxxQ_no(:,:,:,j) = bielec_PxxQ(:,:,:,j)
-    do l=1,n_core_inact_act_orb
-      do k=1,n_core_inact_act_orb
-        do p=1,n_act_orb
-            f(p,k,l) = bielec_PxxQ_no(list_act(p),k,l,j)
-        end do
-      end do
-    end do
-    call dgemm('T','N',n_act_orb,n_core_inact_act_orb**2,n_act_orb,1.d0,  &
-          natorbsCI, size(natorbsCI,1),                              &
-          f, n_act_orb,                                              &
-          0.d0,                                                      &
-          d, n_act_orb)
-    do l=1,n_core_inact_act_orb
-      do k=1,n_core_inact_act_orb
-        do p=1,n_act_orb
-          bielec_PxxQ_no(list_act(p),k,l,j)=d(p,k,l)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO NOWAIT
-
-  deallocate (f,d)
-
-  allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), &
-    d(n_act_orb,mo_num,n_core_inact_act_orb))
-
-  !$OMP DO
-  do k=1,mo_num
-    do l=1,n_core_inact_act_orb
-      do j=1,mo_num
-        do p=1,n_act_orb
-          f(p,j,l) = bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)
-        end do
-      end do
-    end do
-    call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0,  &
-          natorbsCI, size(natorbsCI,1),                              &
-          f, n_act_orb,                                              &
-          0.d0,                                                      &
-          d, n_act_orb)
-    do l=1,n_core_inact_act_orb
-      do j=1,mo_num
-        do p=1,n_act_orb
-          bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(p,j,l)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO NOWAIT
-
-  deallocate(f,d)
-
-  allocate(f(mo_num,n_core_inact_act_orb,n_act_orb), &
-    d(mo_num,n_core_inact_act_orb,n_act_orb) )
-
-  !$OMP DO
-  do k=1,mo_num
-    do p=1,n_act_orb
-      do l=1,n_core_inact_act_orb
-        do j=1,mo_num
-          f(j,l,p) = bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)
-        end do
-      end do
-    end do
-    call dgemm('N','N',mo_num*n_core_inact_act_orb,n_act_orb,n_act_orb,1.d0,  &
-          f, mo_num*n_core_inact_act_orb,                                       &
-          natorbsCI, size(natorbsCI,1),                              &
-          0.d0,                                                      &
-          d, mo_num*n_core_inact_act_orb)
-    do p=1,n_act_orb
-      do l=1,n_core_inact_act_orb
-        do j=1,mo_num
-          bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(j,l,p)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO NOWAIT
-
-  !$OMP BARRIER
-
-  !$OMP DO 
-  do l=1,n_core_inact_act_orb
-    do p=1,n_act_orb
-      do k=1,n_core_inact_act_orb
-        do j=1,mo_num
-          f(j,k,p) = bielec_PxxQ_no(j,k,l,n_core_inact_orb+p)
-        end do
-      end do
-    end do
-    call dgemm('N','N',mo_num*n_core_inact_act_orb,n_act_orb,n_act_orb,1.d0,  &
-          f, mo_num*n_core_inact_act_orb,                                       &
-          natorbsCI, size(natorbsCI,1),                              &
-          0.d0,                                                      &
-          d, mo_num*n_core_inact_act_orb)
-    do p=1,n_act_orb
-      do k=1,n_core_inact_act_orb
-        do j=1,mo_num
-          bielec_PxxQ_no(j,k,l,n_core_inact_orb+p)=d(j,k,p)
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO NOWAIT
-  deallocate(f,d)
-  !$OMP END PARALLEL
-
-END_PROVIDER
-
-
-BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
-  BEGIN_DOC
-  ! integrals (tu|vp) in the basis of natural MOs
-  ! index p runs over the whole basis, t,u,v only over the active orbitals
-  END_DOC
-  implicit none
-  integer                        :: i,j,k,l,t,u,p,q
-  double precision, allocatable  :: f(:,:,:), d(:,:,:)
-  
-  !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP PRIVATE(j,k,l,p,d,f) &
-  !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
-  !$OMP   bielecCI_no,bielecCI,list_act,natorbsCI)
-
-  allocate (f(n_act_orb,n_act_orb,mo_num), &
-      d(n_act_orb,n_act_orb,mo_num))
-
-  !$OMP DO
-  do l=1,mo_num
-    bielecCI_no(:,:,:,l) = bielecCI(:,:,:,l)
-    do k=1,n_act_orb
-      do j=1,n_act_orb
-        do p=1,n_act_orb
-          f(p,j,k)=bielecCI_no(p,j,k,l)
-        end do
-      end do
-    end do
-    call dgemm('T','N',n_act_orb,n_act_orb*n_act_orb,n_act_orb,1.d0,  &
-          natorbsCI, size(natorbsCI,1),                              &
-          f, n_act_orb,                                              &
-          0.d0,                                                      &
-          d, n_act_orb)
-    do k=1,n_act_orb
-      do j=1,n_act_orb
-        do p=1,n_act_orb
-          bielecCI_no(p,j,k,l)=d(p,j,k)
-        end do
-      end do
-
-      do j=1,n_act_orb
-        do p=1,n_act_orb
-          f(p,j,k)=bielecCI_no(j,p,k,l)
-        end do
-      end do
-    end do
-    call dgemm('T','N',n_act_orb,n_act_orb*n_act_orb,n_act_orb,1.d0,  &
-          natorbsCI, n_act_orb,                                      &
-          f, n_act_orb,                                              &
-          0.d0,                                                      &
-          d, n_act_orb)
-    do k=1,n_act_orb
-      do p=1,n_act_orb
-        do j=1,n_act_orb
-          bielecCI_no(j,p,k,l)=d(p,j,k)
-        end do
-      end do
-    end do
-
-    do p=1,n_act_orb
-      do k=1,n_act_orb
-        do j=1,n_act_orb
-          f(j,k,p)=bielecCI_no(j,k,p,l)
-        end do
-      end do
-    end do
-    call dgemm('N','N',n_act_orb*n_act_orb,n_act_orb,n_act_orb,1.d0,       &
-          f, n_act_orb*n_act_orb,                                          &
-          natorbsCI, n_act_orb,                                      &
-          0.d0,                                                      &
-          d, n_act_orb*n_act_orb)
-
-    do p=1,n_act_orb
-      do k=1,n_act_orb
-        do j=1,n_act_orb
-          bielecCI_no(j,k,p,l)=d(j,k,p)
-        end do
-      end do 
-    end do  
-  end do  
-  !$OMP END DO
-
-  !$OMP DO
-  do l=1,n_act_orb
-    do p=1,n_act_orb
-      do k=1,n_act_orb
-        do j=1,n_act_orb
-          f(j,k,p)=bielecCI_no(j,k,l,list_act(p))
-        end do
-      end do
-    end do
-    call dgemm('N','N',n_act_orb*n_act_orb,n_act_orb,n_act_orb,1.d0,       &
-          f, n_act_orb*n_act_orb,                                          &
-          natorbsCI, n_act_orb,                                      &
-          0.d0,                                                      &
-          d, n_act_orb*n_act_orb)
-
-    do p=1,n_act_orb
-      do k=1,n_act_orb
-        do j=1,n_act_orb
-          bielecCI_no(j,k,l,list_act(p))=d(j,k,p)
-        end do
-      end do 
-    end do  
-  end do  
-  !$OMP END DO
-  
-  deallocate(d,f)
-  !$OMP END PARALLEL
-
-
-END_PROVIDER
-
diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f
deleted file mode 100644
index 950cfd55..00000000
--- a/src/casscf/casscf.irp.f
+++ /dev/null
@@ -1,58 +0,0 @@
-program casscf
-  implicit none
-  BEGIN_DOC
-! TODO : Put the documentation of the program here
-  END_DOC
-  call reorder_orbitals_for_casscf
-  no_vvvv_integrals = .True.
-  touch no_vvvv_integrals
-  pt2_max = 0.02
-  SOFT_TOUCH no_vvvv_integrals pt2_max
-  call run_stochastic_cipsi
-  call run
-end
-
-subroutine run
-  implicit none
-  double precision               :: energy_old, energy
-  logical                        :: converged,state_following_casscf_save
-  integer                        :: iteration
-  converged = .False.
-
-  energy = 0.d0
-  mo_label = "MCSCF"
-  iteration = 1
-  state_following_casscf_save = state_following_casscf
-  state_following_casscf = .True.
-  touch state_following_casscf
-  do while (.not.converged)
-    call run_stochastic_cipsi
-    energy_old = energy
-    energy = eone+etwo+ecore
-
-    call write_time(6)
-    call write_int(6,iteration,'CAS-SCF iteration')
-    call write_double(6,energy,'CAS-SCF energy')
-    call write_double(6,energy_improvement, 'Predicted energy improvement')
-
-    converged = dabs(energy_improvement) < thresh_scf
-    pt2_max = dabs(energy_improvement / pt2_relative_error)
-
-    mo_coef = NewOrbs
-    mo_occ  = occnum  
-    call save_mos
-    iteration += 1
-    N_det = max(N_det/2 ,N_states)
-    psi_det = psi_det_sorted
-    psi_coef = psi_coef_sorted
-    read_wf = .True.
-    call clear_mo_map
-    SOFT_TOUCH mo_coef N_det pt2_max  psi_det psi_coef 
-    if(iteration .gt. 3)then
-     state_following_casscf = state_following_casscf_save 
-     touch state_following_casscf
-    endif
-
-  enddo
-
-end
diff --git a/src/casscf/class.irp.f b/src/casscf/class.irp.f
deleted file mode 100644
index 7360a661..00000000
--- a/src/casscf/class.irp.f
+++ /dev/null
@@ -1,12 +0,0 @@
- BEGIN_PROVIDER [ logical, do_only_1h1p ]
-&BEGIN_PROVIDER [ logical, do_only_cas  ]
-&BEGIN_PROVIDER [ logical, do_ddci ]
- implicit none
- BEGIN_DOC
- ! In the CAS case, all those are always false except do_only_cas
- END_DOC
- do_only_cas  = .True.
- do_only_1h1p = .False.
- do_ddci = .False.
-END_PROVIDER
-
diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f
deleted file mode 100644
index d181d732..00000000
--- a/src/casscf/densities.irp.f
+++ /dev/null
@@ -1,63 +0,0 @@
-use bitmasks
-
-BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
-  implicit none
-  BEGIN_DOC
-  ! the first-order density matrix in the basis of the starting MOs.
-  ! matrix is state averaged.
-  END_DOC
-  integer                        :: t,u
-  
-  do u=1,n_act_orb
-    do t=1,n_act_orb
-      D0tu(t,u) = one_e_dm_mo_alpha_average( list_act(t), list_act(u) ) + &
-                  one_e_dm_mo_beta_average ( list_act(t), list_act(u) ) 
-    enddo
-  enddo
-  
-END_PROVIDER
-
-BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
-   BEGIN_DOC
-   ! The second-order density matrix in the basis of the starting MOs ONLY IN THE RANGE OF ACTIVE MOS
-   ! The values are state averaged
-   !
-   ! We use the spin-free generators of mono-excitations
-   ! E_pq destroys q and creates p
-   ! D_pq   =     <0|E_pq|0> = D_qp
-   ! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0>
-   !
-   ! P0tuvx(p,q,r,s) = chemist notation : 1/2 <0|E_pq E_rs - delta_qr E_ps|0>
-   END_DOC
-   implicit none
-   integer                        :: t,u,v,x
-   integer                        :: tt,uu,vv,xx
-   integer                        :: mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
-   integer                        :: ierr
-   real*8                         :: phase1,phase11,phase12,phase2,phase21,phase22
-   integer                        :: nu1,nu2,nu11,nu12,nu21,nu22
-   integer                        :: ierr1,ierr2,ierr11,ierr12,ierr21,ierr22
-   real*8                         :: cI_mu(N_states),term
-   integer(bit_kind), dimension(N_int,2) :: det_mu, det_mu_ex
-   integer(bit_kind), dimension(N_int,2) :: det_mu_ex1, det_mu_ex11, det_mu_ex12
-   integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22
-   
-  if (bavard) then
-    write(6,*) ' providing the 2 body RDM on the active part'
-  endif
-
-  P0tuvx= 0.d0
-  do istate=1,N_states
-   do x = 1, n_act_orb
-    do v = 1, n_act_orb
-     do u = 1, n_act_orb
-      do t = 1, n_act_orb
-       !      1 1 2 2                                     1 2 1 2
-       P0tuvx(t,u,v,x) = state_av_act_2_rdm_spin_trace_mo(t,v,u,x)
-      enddo
-     enddo 
-    enddo
-   enddo
-  enddo
-
-END_PROVIDER
diff --git a/src/casscf/det_manip.irp.f b/src/casscf/det_manip.irp.f
deleted file mode 100644
index d8c309a4..00000000
--- a/src/casscf/det_manip.irp.f
+++ /dev/null
@@ -1,125 +0,0 @@
-use bitmasks
-
-subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart,       &
-      ispin,phase,ierr)
-  BEGIN_DOC
-  ! we create the mono-excitation, and determine, if possible,
-  ! the phase and the number in the list of determinants
-  END_DOC
-  implicit none
-  integer(bit_kind)              :: key1(N_int,2),key2(N_int,2)
-  integer(bit_kind), allocatable :: keytmp(:,:)
-  integer                        :: exc(0:2,2,2),ihole,ipart,ierr,nu,ispin
-  real*8                         :: phase
-  logical                        :: found
-  allocate(keytmp(N_int,2))
-  
-  nu=-1
-  phase=1.D0
-  ierr=0
-  call det_copy(key1,key2,N_int)
-  !        write(6,*) ' key2 before excitation ',ihole,' -> ',ipart,' spin = ',ispin
-  !        call print_det(key2,N_int)
-  call do_single_excitation(key2,ihole,ipart,ispin,ierr)
-  !        write(6,*) ' key2 after ',ihole,' -> ',ipart,' spin = ',ispin
-  !        call print_det(key2,N_int)
-  !        write(6,*) ' excitation ',ihole,' -> ',ipart,' gives ierr = ',ierr
-  if (ierr.eq.1) then
-    ! excitation is possible
-    ! get the phase
-    call get_single_excitation(key1,key2,exc,phase,N_int)
-    ! get the number in the list
-    found=.false.
-    nu=0
-
-    !TODO BOTTLENECK
-    do while (.not.found)
-      nu+=1
-      if (nu.gt.N_det) then
-        ! the determinant is possible, but not in the list
-        found=.true.
-        nu=-1
-      else
-        call det_extract(keytmp,nu,N_int)
-        integer                        :: i,ii
-        found=.true.
-        do ii=1,2
-          do i=1,N_int
-            if (keytmp(i,ii).ne.key2(i,ii)) then
-              found=.false.
-            end if
-          end do
-        end do
-      end if
-    end do
-  end if
-  !
-  ! we found the new string, the phase, and possibly the number in the list
-  !
-end subroutine do_signed_mono_excitation
-
-subroutine det_extract(key,nu,Nint)
-  BEGIN_DOC
-  ! extract a determinant from the list of determinants
-  END_DOC
-  implicit none
-  integer                        :: ispin,i,nu,Nint
-  integer(bit_kind)              :: key(Nint,2)
-  do ispin=1,2
-    do i=1,Nint
-      key(i,ispin)=psi_det(i,ispin,nu)
-    end do
-  end do
-end subroutine det_extract
-
-subroutine det_copy(key1,key2,Nint)
-  use bitmasks ! you need to include the bitmasks_module.f90 features
-  BEGIN_DOC
-  ! copy a determinant from key1 to key2
-  END_DOC
-  implicit none
-  integer                        :: ispin,i,Nint
-  integer(bit_kind)              :: key1(Nint,2),key2(Nint,2)
-  do ispin=1,2
-    do i=1,Nint
-      key2(i,ispin)=key1(i,ispin)
-    end do
-  end do
-end subroutine det_copy
-
-subroutine do_spinfree_mono_excitation(key_in,key_out1,key_out2      &
-      ,nu1,nu2,ihole,ipart,phase1,phase2,ierr,jerr)
-  BEGIN_DOC
-  ! we create the spin-free mono-excitation E_pq=(a^+_p a_q + a^+_P a_Q)
-  ! we may create two determinants as result
-  !
-  END_DOC
-  implicit none
-  integer(bit_kind)              :: key_in(N_int,2),key_out1(N_int,2)
-  integer(bit_kind)              :: key_out2(N_int,2)
-  integer                        :: ihole,ipart,ierr,jerr,nu1,nu2
-  integer                        :: ispin
-  real*8                         :: phase1,phase2
-  
-  !        write(6,*) ' applying E_',ipart,ihole,' on determinant '
-  !        call print_det(key_in,N_int)
-  
-  ! spin alpha
-  ispin=1
-  call do_signed_mono_excitation(key_in,key_out1,nu1,ihole           &
-      ,ipart,ispin,phase1,ierr)
-  !        if (ierr.eq.1) then
-  !         write(6,*) ' 1 result is ',nu1,phase1
-  !         call print_det(key_out1,N_int)
-  !        end if
-  ! spin beta
-  ispin=2
-  call do_signed_mono_excitation(key_in,key_out2,nu2,ihole           &
-      ,ipart,ispin,phase2,jerr)
-  !        if (jerr.eq.1) then
-  !         write(6,*) ' 2 result is ',nu2,phase2
-  !         call print_det(key_out2,N_int)
-  !        end if
-  
-end subroutine do_spinfree_mono_excitation
-
diff --git a/src/casscf/driver_optorb.irp.f b/src/casscf/driver_optorb.irp.f
deleted file mode 100644
index 2e3e02dc..00000000
--- a/src/casscf/driver_optorb.irp.f
+++ /dev/null
@@ -1,3 +0,0 @@
-subroutine driver_optorb
-  implicit none
-end 
diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f
deleted file mode 100644
index cfb26b59..00000000
--- a/src/casscf/get_energy.irp.f
+++ /dev/null
@@ -1,51 +0,0 @@
-program print_2rdm
- implicit none
- BEGIN_DOC
- ! get the active part of the bielectronic energy on a given wave function.
- !
- ! useful to test the active part of the spin trace 2 rdms
- END_DOC
-!no_vvvv_integrals = .True.
- read_wf = .True.
-!touch read_wf no_vvvv_integrals
-!call routine
-!call routine_bis
- call print_grad
-end
-
-subroutine print_grad
- implicit none
- integer :: i 
- do i = 1, nMonoEx
-  if(dabs(gradvec2(i)).gt.1.d-5)then
-   print*,''
-   print*,i,gradvec2(i),excit(:,i)
-  endif
- enddo
-end
-
-subroutine routine
- integer :: i,j,k,l
- integer :: ii,jj,kk,ll
- double precision :: accu(4),twodm,thr,act_twodm2,integral,get_two_e_integral
- thr = 1.d-10
-
-
- accu = 0.d0
- do ll = 1, n_act_orb
-  l = list_act(ll)
-  do kk = 1, n_act_orb
-   k = list_act(kk)
-   do jj = 1, n_act_orb
-    j = list_act(jj)
-    do ii = 1, n_act_orb
-     i = list_act(ii)
-     integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
-     accu(1) += state_av_act_2_rdm_spin_trace_mo(ii,jj,kk,ll) * integral
-    enddo
-   enddo
-  enddo
- enddo
- print*,'accu             = ',accu(1)
-
-end
diff --git a/src/casscf/grad_old.irp.f b/src/casscf/grad_old.irp.f
deleted file mode 100644
index d60a60c8..00000000
--- a/src/casscf/grad_old.irp.f
+++ /dev/null
@@ -1,74 +0,0 @@
-
-BEGIN_PROVIDER [real*8, gradvec_old, (nMonoEx)]
-  BEGIN_DOC
-  ! calculate the orbital gradient <Psi| H E_pq |Psi> by hand, i.e. for
-  ! each determinant I we determine the string E_pq |I> (alpha and beta
-  ! separately) and generate <Psi|H E_pq |I>
-  ! sum_I c_I <Psi|H E_pq |I> is then the pq component of the orbital
-  ! gradient
-  ! E_pq = a^+_pa_q + a^+_Pa_Q
-  END_DOC
-  implicit none
-  integer                        :: ii,tt,aa,indx,ihole,ipart,istate
-  real*8                         :: res
-  
-  do indx=1,nMonoEx
-    ihole=excit(1,indx)
-    ipart=excit(2,indx)
-    call calc_grad_elem(ihole,ipart,res)
-    gradvec_old(indx)=res
-  end do
-  
-  real*8                         :: norm_grad
-  norm_grad=0.d0
-  do indx=1,nMonoEx
-    norm_grad+=gradvec_old(indx)*gradvec_old(indx)
-  end do
-  norm_grad=sqrt(norm_grad)
-  if (bavard) then
-    write(6,*)
-    write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
-    write(6,*)
-  endif
-  
-  
-END_PROVIDER
-
-subroutine calc_grad_elem(ihole,ipart,res)
-  BEGIN_DOC
-  ! eq 18 of Siegbahn et al, Physica Scripta 1980
-  ! we calculate 2 <Psi| H E_pq | Psi>, q=hole, p=particle
-  END_DOC
-  implicit none
-  integer                        :: ihole,ipart,mu,iii,ispin,ierr,nu,istate
-  real*8                         :: res
-  integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:)
-  real*8                         :: i_H_psi_array(N_states),phase
-  allocate(det_mu(N_int,2))
-  allocate(det_mu_ex(N_int,2))
-  
-  res=0.D0
-  
-  do mu=1,n_det
-    ! get the string of the determinant
-    call det_extract(det_mu,mu,N_int)
-    do ispin=1,2
-      ! do the monoexcitation on it
-      call det_copy(det_mu,det_mu_ex,N_int)
-      call do_signed_mono_excitation(det_mu,det_mu_ex,nu             &
-          ,ihole,ipart,ispin,phase,ierr)
-      if (ierr.eq.1) then
-        call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int                &
-            ,N_det,N_det,N_states,i_H_psi_array)
-        do istate=1,N_states
-          res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase
-        end do
-      end if
-    end do
-  end do
-  
-  ! state-averaged gradient
-  res*=2.D0/dble(N_states)
-  
-end subroutine calc_grad_elem
-
diff --git a/src/casscf/gradient.irp.f b/src/casscf/gradient.irp.f
deleted file mode 100644
index e717e822..00000000
--- a/src/casscf/gradient.irp.f
+++ /dev/null
@@ -1,171 +0,0 @@
-use bitmasks
-
-BEGIN_PROVIDER [ integer, nMonoEx ]
-  BEGIN_DOC
-  ! Number of single excitations
-  END_DOC
-  implicit none
-  nMonoEx=n_core_inact_orb*n_act_orb+n_core_inact_orb*n_virt_orb+n_act_orb*n_virt_orb
-END_PROVIDER
-
- BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
-&BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)]
-  BEGIN_DOC
-  ! a list of the orbitals involved in the excitation
-  END_DOC
-  
-  implicit none
-  integer                        :: i,t,a,ii,tt,aa,indx
-  indx=0
-  do ii=1,n_core_inact_orb
-    i=list_core_inact(ii)
-    do tt=1,n_act_orb
-      t=list_act(tt)
-      indx+=1
-      excit(1,indx)=i
-      excit(2,indx)=t
-      excit_class(indx)='c-a'
-    end do
-  end do
-  
-  do ii=1,n_core_inact_orb
-    i=list_core_inact(ii)
-    do aa=1,n_virt_orb
-      a=list_virt(aa)
-      indx+=1
-      excit(1,indx)=i
-      excit(2,indx)=a
-      excit_class(indx)='c-v'
-    end do
-  end do
-  
-  do tt=1,n_act_orb
-    t=list_act(tt)
-    do aa=1,n_virt_orb
-      a=list_virt(aa)
-      indx+=1
-      excit(1,indx)=t
-      excit(2,indx)=a
-      excit_class(indx)='a-v'
-    end do
-  end do
-  
-  if (bavard) then
-    write(6,*) ' Filled the table of the Monoexcitations '
-    do indx=1,nMonoEx
-      write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> '              &
-          ,excit(2,indx),'  ',excit_class(indx)
-    end do
-  end if
-  
-END_PROVIDER
-
-BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
-  BEGIN_DOC
-  ! calculate the orbital gradient <Psi| H E_pq |Psi> from density
-  ! matrices and integrals; Siegbahn et al, Phys Scr 1980
-  ! eqs 14 a,b,c
-  END_DOC
-  implicit none
-  integer                        :: i,t,a,indx
-  real*8                         :: gradvec_it,gradvec_ia,gradvec_ta
-  real*8                         :: norm_grad
-  
-  indx=0
-  do i=1,n_core_inact_orb
-    do t=1,n_act_orb
-      indx+=1
-      gradvec2(indx)=gradvec_it(i,t)
-    end do
-  end do
-  
-  do i=1,n_core_inact_orb
-    do a=1,n_virt_orb
-      indx+=1
-      gradvec2(indx)=gradvec_ia(i,a)
-    end do
-  end do
-  
-  do t=1,n_act_orb
-    do a=1,n_virt_orb
-      indx+=1
-      gradvec2(indx)=gradvec_ta(t,a)
-    end do
-  end do
-  
-  norm_grad=0.d0
-  do indx=1,nMonoEx
-    norm_grad+=gradvec2(indx)*gradvec2(indx)
-  end do
-  norm_grad=sqrt(norm_grad)
-    write(6,*)
-    write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
-    write(6,*)
-  
-END_PROVIDER
-
-real*8 function gradvec_it(i,t)
-  BEGIN_DOC
-  ! the orbital gradient core/inactive -> active
-  ! we assume natural orbitals
-  END_DOC
-  implicit none
-  integer                        :: i,t
-  
-  integer                        :: ii,tt,v,vv,x,y
-  integer                        :: x3,y3
-  
-  ii=list_core_inact(i)
-  tt=list_act(t)
-  gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii))
-  gradvec_it-=occnum(tt)*Fipq(ii,tt)
-  do v=1,n_act_orb
-    vv=list_act(v)
-    do x=1,n_act_orb
-      x3=x+n_core_inact_orb
-      do y=1,n_act_orb
-        y3=y+n_core_inact_orb
-        gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3)
-      end do
-    end do
-  end do
-  gradvec_it*=2.D0
-end function gradvec_it
-
-real*8 function gradvec_ia(i,a)
-  BEGIN_DOC
-  ! the orbital gradient core/inactive -> virtual
-  END_DOC
-  implicit none
-  integer                        :: i,a,ii,aa
-  
-  ii=list_core_inact(i)
-  aa=list_virt(a)
-  gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii))
-  gradvec_ia*=2.D0
-  
-end function gradvec_ia
-
-real*8 function gradvec_ta(t,a)
-  BEGIN_DOC
-  ! the orbital gradient active -> virtual
-  ! we assume natural orbitals
-  END_DOC
-  implicit none
-  integer                        :: t,a,tt,aa,v,vv,x,y
-  
-  tt=list_act(t)
-  aa=list_virt(a)
-  gradvec_ta=0.D0
-  gradvec_ta+=occnum(tt)*Fipq(aa,tt)
-  do v=1,n_act_orb
-    do x=1,n_act_orb
-      do y=1,n_act_orb
-        gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa)
-      end do
-    end do
-  end do
-  gradvec_ta*=2.D0
-  
-end function gradvec_ta
-
diff --git a/src/casscf/hessian.irp.f b/src/casscf/hessian.irp.f
deleted file mode 100644
index 52be1b76..00000000
--- a/src/casscf/hessian.irp.f
+++ /dev/null
@@ -1,656 +0,0 @@
-use bitmasks
-
-BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
-  BEGIN_DOC
-  ! calculate the orbital hessian 2 <Psi| E_pq H E_rs |Psi>
-  ! + <Psi| E_pq E_rs H |Psi> + <Psi| E_rs E_pq H |Psi> by hand,
-  ! determinant per determinant, as for the gradient
-  !
-  ! we assume that we have natural active orbitals
-  END_DOC
-  implicit none
-  integer                        :: indx,ihole,ipart
-  integer                        :: jndx,jhole,jpart
-  character*3                    :: iexc,jexc
-  real*8                         :: res
-  
-  if (bavard) then
-    write(6,*) ' providing Hessian matrix hessmat '
-    write(6,*) '  nMonoEx = ',nMonoEx
-  endif
-  
-  do indx=1,nMonoEx
-    do jndx=1,nMonoEx
-      hessmat(indx,jndx)=0.D0
-    end do
-  end do
-  
-  do indx=1,nMonoEx
-    ihole=excit(1,indx)
-    ipart=excit(2,indx)
-    iexc=excit_class(indx)
-    do jndx=indx,nMonoEx
-      jhole=excit(1,jndx)
-      jpart=excit(2,jndx)
-      jexc=excit_class(jndx)
-      call calc_hess_elem(ihole,ipart,jhole,jpart,res)
-      hessmat(indx,jndx)=res
-      hessmat(jndx,indx)=res
-    end do
-  end do
-  
-END_PROVIDER
-
-subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res)
-  BEGIN_DOC
-  ! eq 19 of Siegbahn et al, Physica Scripta 1980
-  ! we calculate  2 <Psi| E_pq H E_rs |Psi>
-  !  + <Psi| E_pq E_rs H |Psi> + <Psi| E_rs E_pq H |Psi>
-  ! average over all states is performed.
-  ! no transition between states.
-  END_DOC
-  implicit none
-  integer                        :: ihole,ipart,ispin,mu,istate
-  integer                        :: jhole,jpart,jspin
-  integer                        :: mu_pq, mu_pqrs, mu_rs, mu_rspq, nu_rs,nu
-  real*8                         :: res
-  integer(bit_kind), allocatable :: det_mu(:,:)
-  integer(bit_kind), allocatable :: det_nu(:,:)
-  integer(bit_kind), allocatable :: det_mu_pq(:,:)
-  integer(bit_kind), allocatable :: det_mu_rs(:,:)
-  integer(bit_kind), allocatable :: det_nu_rs(:,:)
-  integer(bit_kind), allocatable :: det_mu_pqrs(:,:)
-  integer(bit_kind), allocatable :: det_mu_rspq(:,:)
-  real*8                         :: i_H_psi_array(N_states),phase,phase2,phase3
-  real*8                         :: i_H_j_element
-  allocate(det_mu(N_int,2))
-  allocate(det_nu(N_int,2))
-  allocate(det_mu_pq(N_int,2))
-  allocate(det_mu_rs(N_int,2))
-  allocate(det_nu_rs(N_int,2))
-  allocate(det_mu_pqrs(N_int,2))
-  allocate(det_mu_rspq(N_int,2))
-  integer                        :: mu_pq_possible
-  integer                        :: mu_rs_possible
-  integer                        :: nu_rs_possible
-  integer                        :: mu_pqrs_possible
-  integer                        :: mu_rspq_possible
-  
-  res=0.D0
-  
-  ! the terms <0|E E H |0>
-  do mu=1,n_det
-    ! get the string of the determinant
-    call det_extract(det_mu,mu,N_int)
-    do ispin=1,2
-      ! do the monoexcitation pq on it
-      call det_copy(det_mu,det_mu_pq,N_int)
-      call do_signed_mono_excitation(det_mu,det_mu_pq,mu_pq          &
-          ,ihole,ipart,ispin,phase,mu_pq_possible)
-      if (mu_pq_possible.eq.1) then
-        ! possible, but not necessarily in the list
-        ! do the second excitation
-        do jspin=1,2
-          call det_copy(det_mu_pq,det_mu_pqrs,N_int)
-          call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs&
-              ,jhole,jpart,jspin,phase2,mu_pqrs_possible)
-          ! excitation possible
-          if (mu_pqrs_possible.eq.1) then
-            call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int          &
-                ,N_det,N_det,N_states,i_H_psi_array)
-            do istate=1,N_states
-              res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2
-            end do
-          end if
-          ! try the de-excitation with opposite sign
-          call det_copy(det_mu_pq,det_mu_pqrs,N_int)
-          call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs&
-              ,jpart,jhole,jspin,phase2,mu_pqrs_possible)
-          phase2=-phase2
-          ! excitation possible
-          if (mu_pqrs_possible.eq.1) then
-            call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int          &
-                ,N_det,N_det,N_states,i_H_psi_array)
-            do istate=1,N_states
-              res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2
-            end do
-          end if
-        end do
-      end if
-      ! exchange the notion of pq and rs
-      ! do the monoexcitation rs on the initial determinant
-      call det_copy(det_mu,det_mu_rs,N_int)
-      call do_signed_mono_excitation(det_mu,det_mu_rs,mu_rs          &
-          ,jhole,jpart,ispin,phase2,mu_rs_possible)
-      if (mu_rs_possible.eq.1) then
-        ! do the second excitation
-        do jspin=1,2
-          call det_copy(det_mu_rs,det_mu_rspq,N_int)
-          call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq&
-              ,ihole,ipart,jspin,phase3,mu_rspq_possible)
-          ! excitation possible (of course, the result is outside the CAS)
-          if (mu_rspq_possible.eq.1) then
-            call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int          &
-                ,N_det,N_det,N_states,i_H_psi_array)
-            do istate=1,N_states
-              res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3
-            end do
-          end if
-          ! we may try the de-excitation, with opposite sign
-          call det_copy(det_mu_rs,det_mu_rspq,N_int)
-          call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq&
-              ,ipart,ihole,jspin,phase3,mu_rspq_possible)
-          phase3=-phase3
-          ! excitation possible (of course, the result is outside the CAS)
-          if (mu_rspq_possible.eq.1) then
-            call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int          &
-                ,N_det,N_det,N_states,i_H_psi_array)
-            do istate=1,N_states
-              res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3
-            end do
-          end if
-        end do
-      end if
-      !
-      ! the operator E H E, we have to do a double loop over the determinants
-      !  we still have the determinant mu_pq and the phase in memory
-      if (mu_pq_possible.eq.1) then
-        do nu=1,N_det
-          call det_extract(det_nu,nu,N_int)
-          do jspin=1,2
-            call det_copy(det_nu,det_nu_rs,N_int)
-            call do_signed_mono_excitation(det_nu,det_nu_rs,nu_rs    &
-                ,jhole,jpart,jspin,phase2,nu_rs_possible)
-            ! excitation possible ?
-            if (nu_rs_possible.eq.1) then
-              call i_H_j(det_mu_pq,det_nu_rs,N_int,i_H_j_element)
-              do istate=1,N_states
-                res+=2.D0*i_H_j_element*psi_coef(mu,istate)          &
-                    *psi_coef(nu,istate)*phase*phase2
-              end do
-            end if
-          end do
-        end do
-      end if
-    end do
-  end do
-  
-  ! state-averaged Hessian
-  res*=1.D0/dble(N_states)
-  
-end subroutine calc_hess_elem
-
-BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
-  BEGIN_DOC
-  ! explicit hessian matrix from density matrices and integrals
-  ! of course, this will be used for a direct Davidson procedure later
-  ! we will not store the matrix in real life
-  ! formulas are broken down as functions for the 6 classes of matrix elements
-  !
-  END_DOC
-  implicit none
-  integer                        :: i,j,t,u,a,b,indx,jndx,bstart,ustart,indx_shift
-  
-  real*8                         :: hessmat_itju
-  real*8                         :: hessmat_itja
-  real*8                         :: hessmat_itua
-  real*8                         :: hessmat_iajb
-  real*8                         :: hessmat_iatb
-  real*8                         :: hessmat_taub
-  
-  if (bavard) then
-    write(6,*) ' providing Hessian matrix hessmat2 '
-    write(6,*) '  nMonoEx = ',nMonoEx
-  endif
-  
-  !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP SHARED(hessmat2,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) &
-  !$OMP PRIVATE(i,indx,jndx,j,ustart,t,u,a,bstart,indx_shift)
-
-  !$OMP DO
-  do i=1,n_core_inact_orb
-    do t=1,n_act_orb
-      indx = t + (i-1)*n_act_orb
-      jndx=indx
-      do j=i,n_core_inact_orb
-        if (i.eq.j) then
-          ustart=t
-        else
-          ustart=1
-        end if
-        do u=ustart,n_act_orb
-          hessmat2(jndx,indx)=hessmat_itju(i,t,j,u)
-          jndx+=1
-        end do
-      end do
-      do j=1,n_core_inact_orb
-        do a=1,n_virt_orb
-          hessmat2(jndx,indx)=hessmat_itja(i,t,j,a)
-          jndx+=1
-        end do
-      end do
-      do u=1,n_act_orb
-        do a=1,n_virt_orb
-          hessmat2(jndx,indx)=hessmat_itua(i,t,u,a)
-          jndx+=1
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO NOWAIT
-  
-  indx_shift = n_core_inact_orb*n_act_orb
-  !$OMP DO
-  do a=1,n_virt_orb
-    do i=1,n_core_inact_orb
-      indx = a + (i-1)*n_virt_orb + indx_shift
-      jndx=indx
-      do j=i,n_core_inact_orb
-        if (i.eq.j) then
-          bstart=a
-        else
-          bstart=1
-        end if
-        do b=bstart,n_virt_orb
-          hessmat2(jndx,indx)=hessmat_iajb(i,a,j,b)
-          jndx+=1
-        end do
-      end do
-      do t=1,n_act_orb
-        do b=1,n_virt_orb
-          hessmat2(jndx,indx)=hessmat_iatb(i,a,t,b)
-          jndx+=1
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO NOWAIT
-  
-  indx_shift += n_core_inact_orb*n_virt_orb
-  !$OMP DO 
-  do a=1,n_virt_orb
-    do t=1,n_act_orb
-      indx = a + (t-1)*n_virt_orb + indx_shift
-      jndx=indx
-      do u=t,n_act_orb
-        if (t.eq.u) then
-          bstart=a
-        else
-          bstart=1
-        end if
-        do b=bstart,n_virt_orb
-          hessmat2(jndx,indx)=hessmat_taub(t,a,u,b)
-          jndx+=1
-        end do
-      end do
-    end do
-  end do
-  !$OMP END DO 
-
-  !$OMP END PARALLEL
-  
-  do jndx=1,nMonoEx
-    do indx=1,jndx-1
-      hessmat2(indx,jndx) = hessmat2(jndx,indx)
-    enddo
-  enddo
-
-  
-END_PROVIDER
-
-real*8 function hessmat_itju(i,t,j,u)
-  BEGIN_DOC
-  ! the orbital hessian for core/inactive -> active, core/inactive -> active
-  ! i, t, j, u are list indices, the corresponding orbitals are ii,tt,jj,uu
-  !
-  ! we assume natural orbitals
-  END_DOC
-  implicit none
-  integer                        :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj
-  real*8                         :: term,t2
-  
-  ii=list_core_inact(i)
-  tt=list_act(t)
-  if (i.eq.j) then
-    if (t.eq.u) then
-      ! diagonal element
-      term=occnum(tt)*Fipq(ii,ii)+2.D0*(Fipq(tt,tt)+Fapq(tt,tt))     &
-          -2.D0*(Fipq(ii,ii)+Fapq(ii,ii))
-      term+=2.D0*(3.D0*bielec_pxxq_no(tt,i,i,tt)-bielec_pqxx_no(tt,tt,i,i))
-      term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq_no(tt,i,i,tt)             &
-          -bielec_pqxx_no(tt,tt,i,i))
-      term-=occnum(tt)*Fipq(tt,tt)
-      do v=1,n_act_orb
-        vv=list_act(v)
-        do x=1,n_act_orb
-          xx=list_act(x)
-          term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(vv,xx,i,i)      &
-              +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))*              &
-              bielec_pxxq_no(vv,i,i,xx))
-          do y=1,n_act_orb
-            term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx)
-          end do
-        end do
-      end do
-    else
-      ! it/iu, t != u
-      uu=list_act(u)
-      term=2.D0*(Fipq(tt,uu)+Fapq(tt,uu))
-      term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) &
-          -bielec_PQxx_no(tt,uu,i,j))
-      term-=occnum(tt)*Fipq(uu,tt)
-      term-=(occnum(tt)+occnum(uu))                                  &
-          *(3.D0*bielec_PxxQ_no(tt,i,i,uu)-bielec_PQxx_no(uu,tt,i,i))
-      do v=1,n_act_orb
-        vv=list_act(v)
-        !         term-=D0tu(u,v)*Fipq(tt,vv)   ! published, but inverting t and u seems more correct
-        do x=1,n_act_orb
-          xx=list_act(x)
-          term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,i)      &
-              +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v))               &
-              *bielec_pxxq_no(vv,i,i,xx))
-          do y=1,n_act_orb
-            term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(u,v,y,xx)
-          end do
-        end do
-      end do
-    end if
-  else
-    ! it/ju
-    jj=list_core_inact(j)
-    uu=list_act(u)
-    if (t.eq.u) then
-      term=occnum(tt)*Fipq(ii,jj)
-      term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj))
-    else
-      term=0.D0
-    end if
-    term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt)   &
-        -bielec_PQxx_no(tt,uu,i,j))
-    term-=(occnum(tt)+occnum(uu))*                                   &
-        (4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt)          &
-        -bielec_PQxx_no(uu,tt,i,j))
-    do v=1,n_act_orb
-      vv=list_act(v)
-      do x=1,n_act_orb
-        xx=list_act(x)
-        term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,j)        &
-            +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v))                 &
-            *bielec_pxxq_no(vv,i,j,xx))
-      end do
-    end do
-  end if
-  
-  term*=2.D0
-  hessmat_itju=term
-  
-end function hessmat_itju
-
-real*8 function hessmat_itja(i,t,j,a)
-  BEGIN_DOC
-  ! the orbital hessian for core/inactive -> active, core/inactive -> virtual
-  END_DOC
-  implicit none
-  integer                        :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y
-  real*8                         :: term
-  
-  ! it/ja
-  ii=list_core_inact(i)
-  tt=list_act(t)
-  jj=list_core_inact(j)
-  aa=list_virt(a)
-  term=2.D0*(4.D0*bielec_pxxq_no(aa,j,i,tt)                             &
-      -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt))
-  term-=occnum(tt)*(4.D0*bielec_pxxq_no(aa,j,i,tt)                      &
-      -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt))
-  if (i.eq.j) then
-    term+=2.D0*(Fipq(aa,tt)+Fapq(aa,tt))
-    term-=0.5D0*occnum(tt)*Fipq(aa,tt)
-    do v=1,n_act_orb
-      do x=1,n_act_orb
-        do y=1,n_act_orb
-          term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa)
-        end do
-      end do
-    end do
-  end if
-  term*=2.D0
-  hessmat_itja=term
-  
-end function hessmat_itja
-
-real*8 function hessmat_itua(i,t,u,a)
-  BEGIN_DOC
-  ! the orbital hessian for core/inactive -> active, active -> virtual
-  END_DOC
-  implicit none
-  integer                        :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3
-  real*8                         :: term
-  
-  ii=list_core_inact(i)
-  tt=list_act(t)
-  t3=t+n_core_inact_orb
-  uu=list_act(u)
-  u3=u+n_core_inact_orb
-  aa=list_virt(a)
-  if (t.eq.u) then
-    term=-occnum(tt)*Fipq(aa,ii)
-  else
-    term=0.D0
-  end if
-  term-=occnum(uu)*(bielec_pqxx_no(aa,ii,t3,u3)-4.D0*bielec_pqxx_no(aa,uu,t3,i)&
-      +bielec_pxxq_no(aa,t3,u3,ii))
-  do v=1,n_act_orb
-    vv=list_act(v)
-    v3=v+n_core_inact_orb
-    do x=1,n_act_orb
-      integer                        :: x3
-      xx=list_act(x)
-      x3=x+n_core_inact_orb
-      term-=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,ii,v3,x3)        &
-          +(P0tuvx_no(t,v,u,x)+P0tuvx_no(t,v,x,u))                   &
-          *bielec_pqxx_no(aa,xx,v3,i))
-    end do
-  end do
-  if (t.eq.u) then
-    term+=Fipq(aa,ii)+Fapq(aa,ii)
-  end if
-  term*=2.D0
-  hessmat_itua=term
-  
-end function hessmat_itua
-
-real*8 function hessmat_iajb(i,a,j,b)
-  BEGIN_DOC
-  ! the orbital hessian for core/inactive -> virtual, core/inactive -> virtual
-  END_DOC
-  implicit none
-  integer                        :: i,a,j,b,ii,aa,jj,bb
-  real*8                         :: term
-  
-  ii=list_core_inact(i)
-  aa=list_virt(a)
-  if (i.eq.j) then
-    if (a.eq.b) then
-      ! ia/ia
-      term=2.D0*(Fipq(aa,aa)+Fapq(aa,aa)-Fipq(ii,ii)-Fapq(ii,ii))
-      term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,aa)-bielec_pqxx_no(aa,aa,i,i))
-    else
-      bb=list_virt(b)
-      ! ia/ib
-      term=2.D0*(Fipq(aa,bb)+Fapq(aa,bb))
-      term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,bb)-bielec_pqxx_no(aa,bb,i,i))
-    end if
-  else
-    ! ia/jb
-    jj=list_core_inact(j)
-    bb=list_virt(b)
-    term=2.D0*(4.D0*bielec_pxxq_no(aa,i,j,bb)-bielec_pqxx_no(aa,bb,i,j)    &
-        -bielec_pxxq_no(aa,j,i,bb))
-    if (a.eq.b) then
-      term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj))
-    end if
-  end if
-  term*=2.D0
-  hessmat_iajb=term
-  
-end function hessmat_iajb
-
-real*8 function hessmat_iatb(i,a,t,b)
-  BEGIN_DOC
-  ! the orbital hessian for core/inactive -> virtual, active -> virtual
-  END_DOC
-  implicit none
-  integer                        :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3
-  real*8                         :: term
-  
-  ii=list_core_inact(i)
-  aa=list_virt(a)
-  tt=list_act(t)
-  bb=list_virt(b)
-  t3=t+n_core_inact_orb
-  term=occnum(tt)*(4.D0*bielec_pxxq_no(aa,i,t3,bb)-bielec_pxxq_no(aa,t3,i,bb)&
-      -bielec_pqxx_no(aa,bb,i,t3))
-  if (a.eq.b) then
-    term-=Fipq(tt,ii)+Fapq(tt,ii)
-    term-=0.5D0*occnum(tt)*Fipq(tt,ii)
-    do v=1,n_act_orb
-      do x=1,n_act_orb
-        do y=1,n_act_orb
-          term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,ii)
-        end do
-      end do
-    end do
-  end if
-  term*=2.D0
-  hessmat_iatb=term
-  
-end function hessmat_iatb
-
-real*8 function hessmat_taub(t,a,u,b)
-  BEGIN_DOC
-  ! the orbital hessian for act->virt,act->virt
-  END_DOC
-  implicit none
-  integer                        :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y
-  integer                        :: v3,x3
-  real*8                         :: term,t1,t2,t3
-  
-  tt=list_act(t)
-  aa=list_virt(a)
-  if (t == u) then
-    if (a == b) then
-      ! ta/ta
-      t1=occnum(tt)*Fipq(aa,aa)
-      t2=0.D0
-      t3=0.D0
-      t1-=occnum(tt)*Fipq(tt,tt)
-      do v=1,n_act_orb
-        vv=list_act(v)
-        v3=v+n_core_inact_orb
-        do x=1,n_act_orb
-          xx=list_act(x)
-          x3=x+n_core_inact_orb
-          t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3)      &
-              +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))*              &
-              bielec_pxxq_no(aa,x3,v3,aa))
-          do y=1,n_act_orb
-            t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx)
-          end do
-        end do
-      end do
-      term=t1+t2+t3
-    else
-      bb=list_virt(b)
-      ! ta/tb b/=a
-      term=occnum(tt)*Fipq(aa,bb)
-      do v=1,n_act_orb
-        vv=list_act(v)
-        v3=v+n_core_inact_orb
-        do x=1,n_act_orb
-          xx=list_act(x)
-          x3=x+n_core_inact_orb
-          term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3)    &
-              +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))               &
-              *bielec_pxxq_no(aa,x3,v3,bb))
-        end do
-      end do
-    end if
-  else
-    ! ta/ub t/=u
-    uu=list_act(u)
-    bb=list_virt(b)
-    term=0.D0
-    do v=1,n_act_orb
-      vv=list_act(v)
-      v3=v+n_core_inact_orb
-      do x=1,n_act_orb
-        xx=list_act(x)
-        x3=x+n_core_inact_orb
-        term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3)      &
-            +(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v))                 &
-            *bielec_pxxq_no(aa,x3,v3,bb))
-      end do
-    end do
-    if (a.eq.b) then
-      term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu))
-      do v=1,n_act_orb
-        do y=1,n_act_orb
-          do x=1,n_act_orb
-            term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu)
-            term-=P0tuvx_no(u,v,x,y)*bielecCI_no(x,y,v,tt)
-          end do
-        end do
-      end do
-    end if
-    
-  end if
-  
-  term*=2.D0
-  hessmat_taub=term
-  
-end function hessmat_taub
-
-BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)]
-  BEGIN_DOC
-  ! the diagonal of the Hessian, needed for the Davidson procedure
-  END_DOC
-  implicit none
-  integer                        :: i,t,a,indx,indx_shift
-  real*8                         :: hessmat_itju,hessmat_iajb,hessmat_taub
-  
-  !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP SHARED(hessdiag,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) &
-  !$OMP PRIVATE(i,indx,t,a,indx_shift)
-
-  !$OMP DO
-  do i=1,n_core_inact_orb
-    do t=1,n_act_orb
-      indx = t + (i-1)*n_act_orb
-      hessdiag(indx)=hessmat_itju(i,t,i,t)
-    end do
-  end do
-  !$OMP END DO NOWAIT
-  
-  indx_shift = n_core_inact_orb*n_act_orb
-  !$OMP DO 
-  do a=1,n_virt_orb
-    do i=1,n_core_inact_orb
-      indx = a + (i-1)*n_virt_orb + indx_shift
-      hessdiag(indx)=hessmat_iajb(i,a,i,a)
-    end do
-  end do
-  !$OMP END DO NOWAIT
-  
-  indx_shift += n_core_inact_orb*n_virt_orb
-  !$OMP DO 
-  do a=1,n_virt_orb
-    do t=1,n_act_orb
-      indx = a + (t-1)*n_virt_orb + indx_shift
-      hessdiag(indx)=hessmat_taub(t,a,t,a)
-    end do
-  end do
-  !$OMP END DO
-  !$OMP END PARALLEL
-  
-END_PROVIDER
diff --git a/src/casscf/mcscf_fock.irp.f b/src/casscf/mcscf_fock.irp.f
deleted file mode 100644
index e4568405..00000000
--- a/src/casscf/mcscf_fock.irp.f
+++ /dev/null
@@ -1,80 +0,0 @@
-BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ]
-   BEGIN_DOC
-   ! the inactive Fock matrix, in molecular orbitals
-   END_DOC
-   implicit none
-   integer                        :: p,q,k,kk,t,tt,u,uu
-   
-   do q=1,mo_num
-     do p=1,mo_num
-       Fipq(p,q)=one_ints_no(p,q)
-     end do
-   end do
-   
-   ! the inactive Fock matrix
-   do k=1,n_core_inact_orb
-     kk=list_core_inact(k)
-     do q=1,mo_num
-       do p=1,mo_num
-         Fipq(p,q)+=2.D0*bielec_pqxx_no(p,q,k,k) -bielec_pxxq_no(p,k,k,q)
-       end do
-     end do
-   end do
-   
-   if (bavard) then
-     integer                        :: i
-     write(6,*)
-     write(6,*) ' the diagonal of the inactive effective Fock matrix '
-     write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
-     write(6,*)
-   end if
-   
-   
-END_PROVIDER
- 
- 
-BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
-   BEGIN_DOC
-   ! the active active Fock matrix, in molecular orbitals
-   ! we create them in MOs, quite expensive
-   !
-   ! for an implementation in AOs we need first the natural orbitals
-   ! for forming an active density matrix in AOs
-   !
-   END_DOC
-   implicit none
-   integer                        :: p,q,k,kk,t,tt,u,uu
-   
-   Fapq = 0.d0
-   
-   ! the active Fock matrix, D0tu is diagonal
-   do t=1,n_act_orb
-     tt=list_act(t)
-     do q=1,mo_num
-       do p=1,mo_num
-         Fapq(p,q)+=occnum(tt)                                       &
-             *(bielec_pqxx_no(p,q,tt,tt)-0.5D0*bielec_pxxq_no(p,tt,tt,q))
-       end do
-     end do
-   end do
-   
-   if (bavard) then
-     integer                        :: i
-     write(6,*)
-     write(6,*) ' the effective Fock matrix over MOs'
-     write(6,*)
-     
-     write(6,*)
-     write(6,*) ' the diagonal of the inactive effective Fock matrix '
-     write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
-     write(6,*)
-     write(6,*)
-     write(6,*) ' the diagonal of the active Fock matrix '
-     write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num)
-     write(6,*)
-   end if
-   
-   
-END_PROVIDER
- 
- 
diff --git a/src/casscf/natorb.irp.f b/src/casscf/natorb.irp.f
deleted file mode 100644
index 9ce90304..00000000
--- a/src/casscf/natorb.irp.f
+++ /dev/null
@@ -1,231 +0,0 @@
- BEGIN_PROVIDER [real*8, occnum, (mo_num)]
-  implicit none
-  BEGIN_DOC
-  ! MO occupation numbers
-  END_DOC
-  
-  integer :: i
-  occnum=0.D0
-  do i=1,n_core_inact_orb
-    occnum(list_core_inact(i))=2.D0
-  end do
-  
-  do i=1,n_act_orb
-    occnum(list_act(i))=occ_act(i)
-  end do
-
-  if (bavard) then
-    write(6,*) ' occupation numbers '
-    do i=1,mo_num
-      write(6,*) i,occnum(i)
-    end do
-  endif
-
-END_PROVIDER
-
-
- BEGIN_PROVIDER [ real*8, natorbsCI, (n_act_orb,n_act_orb) ]
-&BEGIN_PROVIDER [ real*8, occ_act, (n_act_orb) ]
- implicit none
- BEGIN_DOC
- ! Natural orbitals of CI
- END_DOC
- integer                        :: i, j
- double precision               :: Vt(n_act_orb,n_act_orb)
- 
-! call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb)
- call svd(D0tu, size(D0tu,1), natorbsCI,size(natorbsCI,1), occ_act, Vt, size(Vt,1),n_act_orb,n_act_orb) 
- 
- if (bavard) then
-   write(6,*) ' found occupation numbers as '
-   do i=1,n_act_orb
-     write(6,*) i,occ_act(i)
-   end do
- 
-   integer                        :: nmx
-   real*8                         :: xmx
-   do i=1,n_act_orb
-     ! largest element of the eigenvector should be positive
-     xmx=0.D0
-     nmx=0
-     do j=1,n_act_orb
-       if (abs(natOrbsCI(j,i)).gt.xmx) then
-         nmx=j
-         xmx=abs(natOrbsCI(j,i))
-       end if
-     end do
-     xmx=sign(1.D0,natOrbsCI(nmx,i))
-     do j=1,n_act_orb
-       natOrbsCI(j,i)*=xmx
-     end do
-     
-     write(6,*) ' Eigenvector No ',i
-     write(6,'(5(I3,F12.5))') (j,natOrbsCI(j,i),j=1,n_act_orb)
-   end do
- end if
-  
-END_PROVIDER
-
-
-BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
-  implicit none
-  BEGIN_DOC
-  ! 4-index transformation of 2part matrices
-  END_DOC
-  integer                        :: i,j,k,l,p,q
-  real*8                         :: d(n_act_orb)
-
-  ! index per index
-  ! first quarter
-  P0tuvx_no(:,:,:,:) = P0tuvx(:,:,:,:)
-
-  do j=1,n_act_orb
-    do k=1,n_act_orb
-      do l=1,n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          P0tuvx_no(p,j,k,l)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 2nd quarter
-  do j=1,n_act_orb
-    do k=1,n_act_orb
-      do l=1,n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          P0tuvx_no(j,p,k,l)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 3rd quarter
-  do j=1,n_act_orb
-    do k=1,n_act_orb
-      do l=1,n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          P0tuvx_no(j,k,p,l)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 4th quarter
-  do j=1,n_act_orb
-    do k=1,n_act_orb
-      do l=1,n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          P0tuvx_no(j,k,l,p)=d(p)
-        end do
-      end do
-    end do
-  end do
-
-END_PROVIDER
-
-
-
-BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
-  implicit none
-  BEGIN_DOC
-  ! Transformed one-e integrals
-  END_DOC
-  integer :: i,j, p, q
-  real*8 :: d(n_act_orb)
-  one_ints_no(:,:)=mo_one_e_integrals(:,:)
-
-  ! 1st half-trf
-  do j=1,mo_num
-    do p=1,n_act_orb
-      d(p)=0.D0
-    end do
-    do p=1,n_act_orb
-      do q=1,n_act_orb
-        d(p)+=one_ints_no(list_act(q),j)*natorbsCI(q,p)
-      end do
-    end do
-    do p=1,n_act_orb
-      one_ints_no(list_act(p),j)=d(p)
-    end do
-  end do
-
-  ! 2nd half-trf
-  do j=1,mo_num
-    do p=1,n_act_orb
-      d(p)=0.D0
-    end do
-    do p=1,n_act_orb
-      do q=1,n_act_orb
-        d(p)+=one_ints_no(j,list_act(q))*natorbsCI(q,p)
-      end do
-    end do
-    do p=1,n_act_orb
-      one_ints_no(j,list_act(p))=d(p)
-    end do
-  end do
-END_PROVIDER
-
-
-BEGIN_PROVIDER [ double precision, NatOrbsCI_mos, (mo_num, mo_num) ]
-  implicit none
-  BEGIN_DOC
-  ! Rotation matrix from current MOs to the CI natural MOs
-  END_DOC
-  integer :: p,q
-
-  NatOrbsCI_mos(:,:) = 0.d0
-
-  do q = 1,mo_num
-   NatOrbsCI_mos(q,q) = 1.d0
-  enddo
-
-  do q = 1,n_act_orb
-    do p = 1,n_act_orb
-      NatOrbsCI_mos(list_act(p),list_act(q)) = natorbsCI(p,q)
-    enddo
-  enddo
-END_PROVIDER
-
-
-BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
-  implicit none
-  BEGIN_DOC
-! FCI natural orbitals
-  END_DOC
-
-  call dgemm('N','N', ao_num,mo_num,mo_num,1.d0,                     &
-      mo_coef, size(mo_coef,1),                                      &
-      NatOrbsCI_mos, size(NatOrbsCI_mos,1), 0.d0,                    &
-      NatOrbsFCI, size(NatOrbsFCI,1))
-END_PROVIDER
-
diff --git a/src/casscf/neworbs.irp.f b/src/casscf/neworbs.irp.f
deleted file mode 100644
index 06a89318..00000000
--- a/src/casscf/neworbs.irp.f
+++ /dev/null
@@ -1,221 +0,0 @@
-BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
-  implicit none
-  BEGIN_DOC
-  ! Single-excitation matrix
-  END_DOC
-  
-  integer                        :: i,j
-  
-  do i=1,nMonoEx+1
-    do j=1,nMonoEx+1
-      SXmatrix(i,j)=0.D0
-    end do
-  end do
-  
-  do i=1,nMonoEx
-    SXmatrix(1,i+1)=gradvec2(i)
-    SXmatrix(1+i,1)=gradvec2(i)
-  end do
-  
-  do i=1,nMonoEx
-    do j=1,nMonoEx
-      SXmatrix(i+1,j+1)=hessmat2(i,j)
-      SXmatrix(j+1,i+1)=hessmat2(i,j)
-    end do
-  end do
-  
-  do i = 1, nMonoEx
-   SXmatrix(i+1,i+1) += level_shift_casscf
-  enddo
-  if (bavard) then
-    do i=2,nMonoEx
-      write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i)
-    end do
-  end if
-  
-  
-END_PROVIDER
-
- BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,nMonoEx+1)]
-&BEGIN_PROVIDER [real*8, SXeigenval, (nMonoEx+1)]
-  implicit none
-  BEGIN_DOC
-  ! Eigenvectors/eigenvalues of the single-excitation matrix
-  END_DOC
-  call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
-  if (bavard) then
-    write(6,*) ' SXdiag : lowest 5 eigenvalues '
-    write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1)
-    if(nmonoex.gt.0)then
-     write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2)
-     write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3)
-     write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4)
-     write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5)
-    endif
-    write(6,*)
-    write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1)
-  endif
-END_PROVIDER
-
- BEGIN_PROVIDER [real*8, energy_improvement]
- implicit none
- if(state_following_casscf)then
-  energy_improvement = SXeigenval(best_vector_ovrlp_casscf)
- else 
-  energy_improvement = SXeigenval(1)
- endif
- END_PROVIDER 
-
-
-
- BEGIN_PROVIDER [ integer, best_vector_ovrlp_casscf ]
-&BEGIN_PROVIDER [ double precision,  best_overlap_casscf ]
-  implicit none
-  integer :: i
-  double precision :: c0
-  best_overlap_casscf = 0.D0
-  best_vector_ovrlp_casscf = -1000
-  do i=1,nMonoEx+1
-    if (SXeigenval(i).lt.0.D0) then
-      if (abs(SXeigenvec(1,i)).gt.best_overlap_casscf) then
-        best_overlap_casscf=abs(SXeigenvec(1,i))
-        best_vector_ovrlp_casscf = i
-      end if
-    end if
-  end do
-  if(best_vector_ovrlp_casscf.lt.0)then 
-   best_vector_ovrlp_casscf = minloc(SXeigenval,nMonoEx+1) 
-  endif
-  c0=SXeigenvec(1,best_vector_ovrlp_casscf)
-  if (bavard) then
-    write(6,*) ' SXdiag : eigenvalue for best overlap with '
-    write(6,*) ' previous orbitals = ',SXeigenval(best_vector_ovrlp_casscf)
-    write(6,*) ' weight of the 1st element ',c0
-  endif
- END_PROVIDER 
-
- BEGIN_PROVIDER [double precision, SXvector, (nMonoEx+1)]
-  implicit none
-  BEGIN_DOC
-  ! Best eigenvector of the single-excitation matrix
-  END_DOC
-  integer           :: i
-  double precision  :: c0
-  c0=SXeigenvec(1,best_vector_ovrlp_casscf)
-  do i=1,nMonoEx+1
-    SXvector(i)=SXeigenvec(i,best_vector_ovrlp_casscf)/c0
-  end do
- END_PROVIDER
-
-
-BEGIN_PROVIDER [double precision, NewOrbs, (ao_num,mo_num) ]
-  implicit none
-  BEGIN_DOC
-  ! Updated orbitals
-  END_DOC
-  integer                        :: i,j,ialph
-  
-  if(state_following_casscf)then
-   print*,'Using the state following casscf '
-   call dgemm('N','T', ao_num,mo_num,mo_num,1.d0,                     &
-       NatOrbsFCI, size(NatOrbsFCI,1),                                &
-       Umat, size(Umat,1), 0.d0,                                      &
-       NewOrbs, size(NewOrbs,1))
-
-    level_shift_casscf *= 0.5D0
-    level_shift_casscf = max(level_shift_casscf,0.002d0)
-   !touch level_shift_casscf
-  else
-   if(best_vector_ovrlp_casscf.ne.1.and.n_orb_swap.ne.0)then
-     print*,'Taking the lowest root for the CASSCF'
-     print*,'!!! SWAPPING MOS !!!!!!'
-     level_shift_casscf *= 2.D0
-    level_shift_casscf = min(level_shift_casscf,0.5d0)
-     print*,'level_shift_casscf = ',level_shift_casscf
-     NewOrbs = switch_mo_coef
-    !mo_coef = switch_mo_coef
-    !soft_touch mo_coef
-    !call save_mos_no_occ
-    !stop
-   else 
-    level_shift_casscf *= 0.5D0
-    level_shift_casscf = max(level_shift_casscf,0.002d0)
-   !touch level_shift_casscf
-    call dgemm('N','T', ao_num,mo_num,mo_num,1.d0,                     &
-        NatOrbsFCI, size(NatOrbsFCI,1),                                &
-        Umat, size(Umat,1), 0.d0,                                      &
-        NewOrbs, size(NewOrbs,1))
-   endif
-  endif
-  
-END_PROVIDER
-
-BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
-  implicit none
-  BEGIN_DOC
-  ! Orbital rotation matrix
-  END_DOC
-  integer                        :: i,j,indx,k,iter,t,a,ii,tt,aa
-  logical                        :: converged
-  
-  real*8 :: Tpotmat (mo_num,mo_num), Tpotmat2 (mo_num,mo_num) 
-  real*8 :: Tmat(mo_num,mo_num) 
-  real*8 :: f
-  
-  ! the orbital rotation matrix T
-  Tmat(:,:)=0.D0
-  indx=1
-  do i=1,n_core_inact_orb
-    ii=list_core_inact(i)
-    do t=1,n_act_orb
-      tt=list_act(t)
-      indx+=1
-      Tmat(ii,tt)= SXvector(indx)
-      Tmat(tt,ii)=-SXvector(indx)
-    end do
-  end do
-  do i=1,n_core_inact_orb
-    ii=list_core_inact(i)
-    do a=1,n_virt_orb
-      aa=list_virt(a)
-      indx+=1
-      Tmat(ii,aa)= SXvector(indx)
-      Tmat(aa,ii)=-SXvector(indx)
-    end do
-  end do
-  do t=1,n_act_orb
-    tt=list_act(t)
-    do a=1,n_virt_orb
-      aa=list_virt(a)
-      indx+=1
-      Tmat(tt,aa)= SXvector(indx)
-      Tmat(aa,tt)=-SXvector(indx)
-    end do
-  end do
-  
-  ! Form the exponential
-
-  Tpotmat(:,:)=0.D0
-  Umat(:,:)   =0.D0
-  do i=1,mo_num
-    Tpotmat(i,i)=1.D0
-    Umat(i,i)   =1.d0
-  end do
-  iter=0
-  converged=.false.
-  do while (.not.converged)
-    iter+=1
-    f = 1.d0 / dble(iter)
-    Tpotmat2(:,:) = Tpotmat(:,:) * f
-    call dgemm('N','N', mo_num,mo_num,mo_num,1.d0,                   &
-        Tpotmat2, size(Tpotmat2,1),                                  &
-        Tmat, size(Tmat,1), 0.d0,                                    &
-        Tpotmat, size(Tpotmat,1))
-    Umat(:,:) = Umat(:,:) + Tpotmat(:,:)
-    
-    converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
-  end do
-END_PROVIDER
-  
-
-
diff --git a/src/casscf/reorder_orb.irp.f b/src/casscf/reorder_orb.irp.f
deleted file mode 100644
index 3cb90522..00000000
--- a/src/casscf/reorder_orb.irp.f
+++ /dev/null
@@ -1,70 +0,0 @@
-subroutine reorder_orbitals_for_casscf
- implicit none
- BEGIN_DOC
-! routine that reorders the orbitals of the CASSCF in terms block of core, active and virtual
- END_DOC
- integer :: i,j,iorb
- integer, allocatable :: iorder(:),array(:)
- allocate(iorder(mo_num),array(mo_num))
- do i = 1, n_core_orb
-  iorb = list_core(i)
-  array(iorb) = i
- enddo
-
- do i = 1, n_inact_orb
-  iorb = list_inact(i)
-  array(iorb) = mo_num + i
- enddo
-
- do i = 1, n_act_orb
-  iorb = list_act(i)
-  array(iorb) = 2 * mo_num + i
- enddo
-
- do i = 1, n_virt_orb
-  iorb = list_virt(i)
-  array(iorb) = 3 * mo_num + i
- enddo
-
- do i = 1, mo_num
-  iorder(i) = i
- enddo
- call isort(array,iorder,mo_num)
- double precision, allocatable :: mo_coef_new(:,:)
- allocate(mo_coef_new(ao_num,mo_num))
- do i = 1, mo_num
-  mo_coef_new(:,i) = mo_coef(:,iorder(i))
- enddo
- mo_coef = mo_coef_new 
- touch mo_coef
-
- list_core_reverse = 0
- do i = 1, n_core_orb
-  list_core(i) = i
-  list_core_reverse(i) = i
-  mo_class(i) = "Core"
- enddo
-
- list_inact_reverse = 0
- do i = 1, n_inact_orb
-  list_inact(i) = i + n_core_orb
-  list_inact_reverse(i+n_core_orb) = i
-  mo_class(i+n_core_orb) = "Inactive"
- enddo
-
- list_act_reverse = 0
- do i = 1, n_act_orb
-  list_act(i) = n_core_inact_orb + i
-  list_act_reverse(n_core_inact_orb + i) = i
-  mo_class(n_core_inact_orb + i) = "Active"
- enddo
-
- list_virt_reverse = 0
- do i = 1, n_virt_orb
-  list_virt(i) = n_core_inact_orb + n_act_orb + i
-  list_virt_reverse(n_core_inact_orb + n_act_orb + i) = i
-  mo_class(n_core_inact_orb + n_act_orb + i) = "Virtual"
- enddo
- touch list_core_reverse list_core list_inact list_inact_reverse list_act list_act_reverse list_virt list_virt_reverse 
-
-end
diff --git a/src/casscf/save_energy.irp.f b/src/casscf/save_energy.irp.f
deleted file mode 100644
index 8729c5af..00000000
--- a/src/casscf/save_energy.irp.f
+++ /dev/null
@@ -1,9 +0,0 @@
-subroutine save_energy(E,pt2)
-  implicit none
-  BEGIN_DOC
-! Saves the energy in |EZFIO|.
-  END_DOC
-  double precision, intent(in) :: E(N_states), pt2(N_states)
-  call ezfio_set_casscf_energy(E(1:N_states))
-  call ezfio_set_casscf_energy_pt2(E(1:N_states)+pt2(1:N_states))
-end
diff --git a/src/casscf/superci_dm.irp.f b/src/casscf/superci_dm.irp.f
deleted file mode 100644
index ee831c35..00000000
--- a/src/casscf/superci_dm.irp.f
+++ /dev/null
@@ -1,207 +0,0 @@
- BEGIN_PROVIDER [double precision, super_ci_dm, (mo_num,mo_num)]
- implicit none 
- BEGIN_DOC
-! density matrix of the super CI matrix, in the basis of NATURAL ORBITALS OF THE CASCI WF 
-! 
-! This is obtained from annex B of Roos et. al. Chemical Physics 48 (1980) 157-173
-!
-! WARNING ::: in the equation B3.d there is a TYPO with a forgotten MINUS SIGN (see variable mat_tmp_dm_super_ci ) 
- END_DOC
- super_ci_dm = 0.d0
- integer :: i,j,iorb,jorb
- integer :: a,aorb,b,borb
- integer :: t,torb,v,vorb,u,uorb,x,xorb
- double precision :: c0,ci
- c0 = SXeigenvec(1,1)
- ! equation B3.a of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173
- ! loop over the core/inact 
- do i = 1, n_core_inact_orb
-  iorb = list_core_inact(i)
-  super_ci_dm(iorb,iorb) = 2.d0 ! first term of B3.a
-  ! loop over the core/inact 
-  do j = 1, n_core_inact_orb
-   jorb = list_core_inact(j)
-   ! loop over the virtual 
-   do a = 1, n_virt_orb
-    aorb = list_virt(a)
-    super_ci_dm(jorb,iorb) += -2.d0 * lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,jorb) ! second term in B3.a
-   enddo
-   do t = 1, n_act_orb
-    torb = list_act(t)
-    ! thrid term of the B3.a 
-    super_ci_dm(jorb,iorb) += - lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(jorb,torb) * (2.d0 - occ_act(t))
-   enddo
-  enddo
- enddo
-
- ! equation B3.b of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173
- do i = 1, n_core_inact_orb
-  iorb = list_core_inact(i)
-  do t = 1, n_act_orb
-   torb = list_act(t)
-   super_ci_dm(iorb,torb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t))
-   super_ci_dm(torb,iorb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t))
-   do a = 1, n_virt_orb
-    aorb = list_virt(a)
-    super_ci_dm(iorb,torb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t)
-    super_ci_dm(torb,iorb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t)
-   enddo
-  enddo
- enddo
-
- ! equation B3.c of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173
- do i = 1, n_core_inact_orb
-  iorb = list_core_inact(i)
-  do a = 1, n_virt_orb
-   aorb = list_virt(a)
-   super_ci_dm(aorb,iorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb)
-   super_ci_dm(iorb,aorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb)
-  enddo
- enddo
-
- ! equation B3.d of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173
- do t = 1, n_act_orb
-  torb = list_act(t)
-  super_ci_dm(torb,torb) = occ_act(t) ! first term of equation B3.d 
-  do x = 1, n_act_orb
-   xorb = list_act(x) 
-   super_ci_dm(torb,torb) += - occ_act(x) * occ_act(t)* mat_tmp_dm_super_ci(x,x) ! second term involving the ONE-rdm 
-  enddo
-  do u = 1, n_act_orb
-   uorb = list_act(u)
-
-   ! second term of equation B3.d 
-   do x = 1, n_act_orb
-    xorb = list_act(x) 
-    do v = 1, n_act_orb
-     vorb = list_act(v) 
-     super_ci_dm(torb,uorb) +=  2.d0 * P0tuvx_no(v,x,t,u) * mat_tmp_dm_super_ci(v,x) ! second term involving the TWO-rdm 
-    enddo
-   enddo
-   
-   ! third term of equation B3.d 
-   do i = 1, n_core_inact_orb
-    iorb = list_core_inact(i)
-    super_ci_dm(torb,uorb) += lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(iorb,uorb) * (2.d0 - occ_act(t) - occ_act(u)) 
-   enddo
-
-  enddo
- enddo
-
- ! equation B3.e of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173
- do t = 1, n_act_orb
-  torb = list_act(t)
-  do a = 1, n_virt_orb
-   aorb = list_virt(a)
-   super_ci_dm(aorb,torb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t)
-   super_ci_dm(torb,aorb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t)
-   do i = 1, n_core_inact_orb
-    iorb = list_core_inact(i)
-    super_ci_dm(aorb,torb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t))
-    super_ci_dm(torb,aorb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t))
-   enddo
-  enddo
- enddo
-
- ! equation B3.f of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173
- do a = 1, n_virt_orb
-  aorb = list_virt(a)
-  do b = 1, n_virt_orb
-   borb= list_virt(b)
-
-   ! First term of equation B3.f 
-   do i = 1, n_core_inact_orb
-    iorb = list_core_inact(i)
-    super_ci_dm(borb,aorb) += 2.d0 * lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,borb) 
-   enddo
-
-   ! Second term of equation B3.f
-   do t = 1, n_act_orb
-    torb = list_act(t)
-    super_ci_dm(borb,aorb) += lowest_super_ci_coef_mo(torb,aorb) * lowest_super_ci_coef_mo(torb,borb) * occ_act(t)
-   enddo
-  enddo
- enddo
-
- END_PROVIDER 
-
- BEGIN_PROVIDER [double precision, superci_natorb, (ao_num,mo_num)
-&BEGIN_PROVIDER [double precision, superci_nat_occ, (mo_num)
- implicit none
- call general_mo_coef_new_as_svd_vectors_of_mo_matrix_eig(super_ci_dm,mo_num,mo_num,mo_num,NatOrbsFCI,superci_nat_occ,superci_natorb)
-
-END_PROVIDER 
-
- BEGIN_PROVIDER [double precision, mat_tmp_dm_super_ci, (n_act_orb,n_act_orb)]
- implicit none
- BEGIN_DOC
- ! computation of the term in [ ] in the equation B3.d of Roos et. al. Chemical Physics 48 (1980) 157-173
- !
- ! !!!!! WARNING !!!!!! there is a TYPO: a MINUS SIGN SHOULD APPEAR in that term 
- END_DOC
- integer :: a,aorb,i,iorb
- integer :: x,xorb,v,vorb
- mat_tmp_dm_super_ci = 0.d0 
- do v = 1, n_act_orb
-  vorb = list_act(v)
-  do x = 1, n_act_orb
-   xorb = list_act(x)
-   do a = 1, n_virt_orb
-    aorb = list_virt(a)
-    mat_tmp_dm_super_ci(x,v) += lowest_super_ci_coef_mo(aorb,vorb) * lowest_super_ci_coef_mo(aorb,xorb)
-   enddo
-   do i = 1, n_core_inact_orb
-    iorb = list_core_inact(i)
-    ! MARK THE MINUS SIGN HERE !!!!!!!!!!! BECAUSE OF TYPO IN THE ORIGINAL PAPER 
-    mat_tmp_dm_super_ci(x,v) -= lowest_super_ci_coef_mo(iorb,vorb) * lowest_super_ci_coef_mo(iorb,xorb)
-   enddo
-  enddo
- enddo
- END_PROVIDER 
- 
- BEGIN_PROVIDER [double precision, lowest_super_ci_coef_mo, (mo_num,mo_num)]
- implicit none
- integer :: i,j,iorb,jorb
- integer :: a, aorb,t, torb
- double precision :: sqrt2
-
- sqrt2 = 1.d0/dsqrt(2.d0)
- do i = 1, nMonoEx
-  iorb = excit(1,i)
-  jorb = excit(2,i)
-  lowest_super_ci_coef_mo(iorb,jorb) = SXeigenvec(i+1,1)
-  lowest_super_ci_coef_mo(jorb,iorb) = SXeigenvec(i+1,1)
- enddo
-
- ! a_{it} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173
- do i = 1, n_core_inact_orb
-  iorb = list_core_inact(i)
-  do t = 1, n_act_orb
-   torb = list_act(t)
-   lowest_super_ci_coef_mo(torb,iorb) *= (2.d0 - occ_act(t))**(-0.5d0)
-   lowest_super_ci_coef_mo(iorb,torb) *= (2.d0 - occ_act(t))**(-0.5d0)
-  enddo
- enddo
-
- ! a_{ia} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173
- do i = 1, n_core_inact_orb
-  iorb = list_core_inact(i)
-  do a = 1, n_virt_orb
-   aorb = list_virt(a)
-   lowest_super_ci_coef_mo(aorb,iorb) *= sqrt2
-   lowest_super_ci_coef_mo(iorb,aorb) *= sqrt2
-  enddo
- enddo
- 
- ! a_{ta} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173
- do a = 1, n_virt_orb
-  aorb = list_virt(a)
-  do t = 1, n_act_orb
-   torb = list_act(t)
-   lowest_super_ci_coef_mo(torb,aorb) *= occ_act(t)**(-0.5d0)
-   lowest_super_ci_coef_mo(aorb,torb) *= occ_act(t)**(-0.5d0)
-  enddo
- enddo
-
- END_PROVIDER 
-
diff --git a/src/casscf/swap_orb.irp.f b/src/casscf/swap_orb.irp.f
deleted file mode 100644
index 5d442157..00000000
--- a/src/casscf/swap_orb.irp.f
+++ /dev/null
@@ -1,132 +0,0 @@
- BEGIN_PROVIDER [double precision, SXvector_lowest, (nMonoEx)]
- implicit none
- integer :: i 
-  do i=2,nMonoEx+1
-   SXvector_lowest(i-1)=SXeigenvec(i,1)
-  enddo
- END_PROVIDER 
-
- BEGIN_PROVIDER [double precision, thresh_overlap_switch]
- implicit none
- thresh_overlap_switch = 0.5d0
- END_PROVIDER 
-
- BEGIN_PROVIDER [integer, max_overlap, (nMonoEx)]
-&BEGIN_PROVIDER [integer, n_max_overlap]
-&BEGIN_PROVIDER [integer, dim_n_max_overlap]
- implicit none
- double precision, allocatable :: vec_tmp(:)
- integer, allocatable :: iorder(:)
- allocate(vec_tmp(nMonoEx),iorder(nMonoEx))
- integer :: i
- do i = 1, nMonoEx
-  iorder(i)  = i
-  vec_tmp(i) = -dabs(SXvector_lowest(i))
- enddo
- call dsort(vec_tmp,iorder,nMonoEx)
- n_max_overlap = 0
- do i = 1, nMonoEx
-  if(dabs(vec_tmp(i)).gt.thresh_overlap_switch)then
-   n_max_overlap += 1
-   max_overlap(n_max_overlap) = iorder(i)
-  endif
- enddo
- dim_n_max_overlap = max(1,n_max_overlap)
- END_PROVIDER 
-
- BEGIN_PROVIDER [integer, orb_swap, (2,dim_n_max_overlap)]
-&BEGIN_PROVIDER [integer, index_orb_swap, (dim_n_max_overlap)]
-&BEGIN_PROVIDER [integer, n_orb_swap ]
- implicit none
- use bitmasks ! you need to include the bitmasks_module.f90 features
- integer :: i,imono,iorb,jorb,j
- n_orb_swap = 0
- do i = 1, n_max_overlap
-  imono = max_overlap(i)
-  iorb = excit(1,imono)
-  jorb = excit(2,imono)
-  if (excit_class(imono) == "c-a" .and.hessmat2(imono,imono).gt.0.d0)then ! core --> active rotation
-   n_orb_swap += 1
-   orb_swap(1,n_orb_swap) = iorb ! core 
-   orb_swap(2,n_orb_swap) = jorb ! active
-   index_orb_swap(n_orb_swap) = imono
-  else if (excit_class(imono) == "a-v" .and.hessmat2(imono,imono).gt.0.d0)then ! active --> virtual rotation
-   n_orb_swap += 1
-   orb_swap(1,n_orb_swap) = jorb ! virtual 
-   orb_swap(2,n_orb_swap) = iorb ! active
-   index_orb_swap(n_orb_swap) = imono
-  endif
- enddo
-
- integer,allocatable :: orb_swap_tmp(:,:)
- allocate(orb_swap_tmp(2,dim_n_max_overlap))
- do i = 1, n_orb_swap
-  orb_swap_tmp(1,i) = orb_swap(1,i)
-  orb_swap_tmp(2,i) = orb_swap(2,i)
- enddo
- 
- integer(bit_kind), allocatable :: det_i(:),det_j(:)
- allocate(det_i(N_int),det_j(N_int))
- logical, allocatable :: good_orb_rot(:)
- allocate(good_orb_rot(n_orb_swap))
- integer, allocatable ::  index_orb_swap_tmp(:) 
- allocate(index_orb_swap_tmp(dim_n_max_overlap))
- index_orb_swap_tmp = index_orb_swap
- good_orb_rot = .True.
- integer :: icount,k
- do i = 1, n_orb_swap
-  if(.not.good_orb_rot(i))cycle
-  det_i = 0_bit_kind 
-  call set_bit_to_integer(orb_swap(1,i),det_i,N_int)
-  call set_bit_to_integer(orb_swap(2,i),det_i,N_int)
-  do j = i+1, n_orb_swap
-   det_j = 0_bit_kind 
-   call set_bit_to_integer(orb_swap(1,j),det_j,N_int)
-   call set_bit_to_integer(orb_swap(2,j),det_j,N_int)
-   icount = 0
-   do k = 1, N_int
-    icount += popcnt(ior(det_i(k),det_j(k)))
-   enddo
-   if (icount.ne.4)then
-    good_orb_rot(i) = .False.
-    good_orb_rot(j) = .False.
-    exit
-   endif
-  enddo
- enddo
- icount = n_orb_swap
- n_orb_swap = 0
- do i = 1, icount
-  if(good_orb_rot(i))then
-   n_orb_swap += 1
-   index_orb_swap(n_orb_swap) = index_orb_swap_tmp(i)
-   orb_swap(1,n_orb_swap) = orb_swap_tmp(1,i)
-   orb_swap(2,n_orb_swap) = orb_swap_tmp(2,i)
-  endif
- enddo
-
- if(n_orb_swap.gt.0)then
-  print*,'n_orb_swap = ',n_orb_swap
- endif
- do i = 1, n_orb_swap
-  print*,'imono = ',index_orb_swap(i)
-  print*,orb_swap(1,i),'-->',orb_swap(2,i)
- enddo
- END_PROVIDER 
-
- BEGIN_PROVIDER [double precision, switch_mo_coef, (ao_num,mo_num)]
-  implicit none
-  integer :: i,j,iorb,jorb
-  switch_mo_coef = NatOrbsFCI
-  do i = 1, n_orb_swap
-   iorb = orb_swap(1,i)
-   jorb = orb_swap(2,i)
-   do j = 1, ao_num
-    switch_mo_coef(j,jorb) = NatOrbsFCI(j,iorb)
-   enddo
-   do j = 1, ao_num
-    switch_mo_coef(j,iorb) = NatOrbsFCI(j,jorb)
-   enddo
-  enddo
-
- END_PROVIDER 
diff --git a/src/casscf/test_pert_2rdm.irp.f b/src/casscf/test_pert_2rdm.irp.f
deleted file mode 100644
index 7c40de0f..00000000
--- a/src/casscf/test_pert_2rdm.irp.f
+++ /dev/null
@@ -1,29 +0,0 @@
-program test_pert_2rdm
- implicit none
- read_wf = .True.
- touch read_wf
-!call get_pert_2rdm
- integer :: i,j,k,l,ii,jj,kk,ll
- double precision :: accu , get_two_e_integral, integral
- accu = 0.d0
- print*,'n_orb_pert_rdm = ',n_orb_pert_rdm
- do ii = 1, n_orb_pert_rdm
-  i = list_orb_pert_rdm(ii) 
-  do jj = 1, n_orb_pert_rdm
-   j = list_orb_pert_rdm(jj) 
-   do kk = 1, n_orb_pert_rdm
-    k= list_orb_pert_rdm(kk) 
-    do ll = 1, n_orb_pert_rdm
-     l = list_orb_pert_rdm(ll) 
-     integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
-!    if(dabs(pert_2rdm_provider(ii,jj,kk,ll) * integral).gt.1.d-12)then
-!     print*,i,j,k,l
-!     print*,pert_2rdm_provider(ii,jj,kk,ll) * integral,pert_2rdm_provider(ii,jj,kk,ll), pert_2rdm_provider(ii,jj,kk,ll), integral
-!    endif
-     accu    += pert_2rdm_provider(ii,jj,kk,ll) * integral
-    enddo
-   enddo
-  enddo
- enddo
- print*,'accu = ',accu
-end
diff --git a/src/casscf/tot_en.irp.f b/src/casscf/tot_en.irp.f
deleted file mode 100644
index 1d70e087..00000000
--- a/src/casscf/tot_en.irp.f
+++ /dev/null
@@ -1,101 +0,0 @@
- BEGIN_PROVIDER [real*8, etwo]
-&BEGIN_PROVIDER [real*8, eone]
-&BEGIN_PROVIDER [real*8, eone_bis]
-&BEGIN_PROVIDER [real*8, etwo_bis]
-&BEGIN_PROVIDER [real*8, etwo_ter]
-&BEGIN_PROVIDER [real*8, ecore]
-&BEGIN_PROVIDER [real*8, ecore_bis]
-   implicit none
-   integer                        :: t,u,v,x,i,ii,tt,uu,vv,xx,j,jj,t3,u3,v3,x3
-   real*8                         :: e_one_all,e_two_all
-   e_one_all=0.D0
-   e_two_all=0.D0
-   do i=1,n_core_inact_orb
-     ii=list_core_inact(i)
-     e_one_all+=2.D0*mo_one_e_integrals(ii,ii)
-     do j=1,n_core_inact_orb
-       jj=list_core_inact(j)
-       e_two_all+=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i)
-     end do
-     do t=1,n_act_orb
-       tt=list_act(t)
-       t3=t+n_core_inact_orb
-       do u=1,n_act_orb
-         uu=list_act(u)
-         u3=u+n_core_inact_orb
-         e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i)           &
-             -bielec_PQxx(tt,ii,i,u3))
-       end do
-     end do
-   end do
-   do t=1,n_act_orb
-     tt=list_act(t)
-     do u=1,n_act_orb
-       uu=list_act(u)
-       e_one_all+=D0tu(t,u)*mo_one_e_integrals(tt,uu)
-       do v=1,n_act_orb
-         v3=v+n_core_inact_orb
-         do x=1,n_act_orb
-           x3=x+n_core_inact_orb
-           e_two_all  +=P0tuvx(t,u,v,x)*bielec_PQxx(tt,uu,v3,x3)
-         end do
-       end do
-     end do
-   end do
-   ecore    =nuclear_repulsion
-   ecore_bis=nuclear_repulsion
-   do i=1,n_core_inact_orb
-     ii=list_core_inact(i)
-     ecore    +=2.D0*mo_one_e_integrals(ii,ii)
-     ecore_bis+=2.D0*mo_one_e_integrals(ii,ii)
-     do j=1,n_core_inact_orb
-       jj=list_core_inact(j)
-       ecore    +=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i)
-       ecore_bis+=2.D0*bielec_PxxQ(ii,i,j,jj)-bielec_PxxQ(ii,j,j,ii)
-     end do
-   end do
-   eone    =0.D0
-   eone_bis=0.D0
-   etwo    =0.D0
-   etwo_bis=0.D0
-   etwo_ter=0.D0
-   do t=1,n_act_orb
-     tt=list_act(t)
-     t3=t+n_core_inact_orb
-     do u=1,n_act_orb
-       uu=list_act(u)
-       u3=u+n_core_inact_orb
-       eone    +=D0tu(t,u)*mo_one_e_integrals(tt,uu)
-       eone_bis+=D0tu(t,u)*mo_one_e_integrals(tt,uu)
-       do i=1,n_core_inact_orb
-         ii=list_core_inact(i)
-         eone    +=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i)            &
-             -bielec_PQxx(tt,ii,i,u3))
-         eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQ(tt,u3,i,ii)           &
-             -bielec_PxxQ(tt,i,i,uu))
-       end do
-       do v=1,n_act_orb
-         vv=list_act(v)
-         v3=v+n_core_inact_orb
-         do x=1,n_act_orb
-           xx=list_act(x)
-           x3=x+n_core_inact_orb
-           real*8                         :: h1,h2,h3
-           h1=bielec_PQxx(tt,uu,v3,x3)
-           h2=bielec_PxxQ(tt,u3,v3,xx)
-           h3=bielecCI(t,u,v,xx)
-           etwo    +=P0tuvx(t,u,v,x)*h1
-           etwo_bis+=P0tuvx(t,u,v,x)*h2
-           etwo_ter+=P0tuvx(t,u,v,x)*h3
-           if ((h1.ne.h2).or.(h1.ne.h3)) then
-             write(6,9901) t,u,v,x,h1,h2,h3
-             9901 format('aie: ',4I4,3E20.12)
-           end if
-         end do
-       end do
-     end do
-   end do
-   
-END_PROVIDER
- 
-