Merge branch 'dev-stable' of https://github.com/QuantumPackage/qp2 into dev-stable

.github/workflows/configuration.yml

@@ -22,7 +22,7 @@ jobs:
        - uses: actions/checkout@v3
        - name: Install dependencies
          run: |
-           sudo apt install gfortran gcc liblapack-dev libblas-dev wget python3 make m4 pkg-config hdf5
+           sudo apt install gfortran gcc liblapack-dev libblas-dev wget python3 make m4 pkg-config libhdf5-dev
        - name: zlib
          run: |
               ./configure -i zlib || echo OK

src/casscf_cipsi/50.casscf.bats

@@ -0,0 +1,49 @@
+#!/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
+#   
+
+}
+

src/casscf_cipsi/EZFIO.cfg

@@ -0,0 +1,75 @@
+[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)
+
+[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
+
+
+[diag_hess_cas]
+type: logical
+doc: If |true|, only the DIAGONAL part of the hessian is retained for the CASSCF
+interface: ezfio,provider,ocaml
+default: False
+
+[hess_cv_cv]
+type: logical
+doc: If |true|, the core-virtual - core-virtual part of the hessian is computed
+interface: ezfio,provider,ocaml
+default: True
+
+
+[level_shift_casscf]
+type: Positive_float
+doc: Energy shift on the virtual MOs to improve SCF convergence
+interface: ezfio,provider,ocaml
+default: 0.005
+
+
+[fast_2rdm]
+type: logical
+doc: If true, the two-rdm are computed with a fast algo
+interface: ezfio,provider,ocaml
+default: True
+
+[criterion_casscf]
+type: character*(32)
+doc: choice of the criterion for the convergence of the casscf: can be  energy  or  gradients  or e_pt2
+interface: ezfio, provider, ocaml
+default: e_pt2
+
+[thresh_casscf]
+type: Threshold
+doc: Threshold on the convergence of the CASCF energy.
+interface: ezfio,provider,ocaml
+default: 1.e-06
+
+
+[pt2_min_casscf]
+type: Threshold
+doc: Minimum value of the pt2_max parameter for the CIPSI in the CASSCF iterations.
+interface: ezfio,provider,ocaml
+default: 1.e-04
+
+[n_big_act_orb]
+type: integer
+doc: Number of active orbitals from which the active space is considered as large, and therefore pt2_min_casscf is activated.
+interface: ezfio,provider,ocaml
+default: 16
+
+[adaptive_pt2_max]
+type: logical
+doc: If |true|, the pt2_max value in the CIPSI iterations will automatically adapt, otherwise it is fixed at the value given in the EZFIO folder
+interface: ezfio,provider,ocaml
+default: True

src/casscf_cipsi/NEED

@@ -0,0 +1,5 @@
+cipsi
+selectors_full
+generators_cas
+two_body_rdm
+dav_general_mat

src/casscf_cipsi/README.rst

@@ -0,0 +1,5 @@
+======
+casscf
+======
+
+|CASSCF| program with the CIPSI algorithm.

src/casscf_cipsi/bavard.irp.f

@@ -0,0 +1,6 @@
+! -*- F90 -*- 
+BEGIN_PROVIDER [logical, bavard]
+!  bavard=.true.
+   bavard=.false.
+END_PROVIDER
+

src/casscf_cipsi/bielec.irp.f

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

src/casscf_cipsi/bielec_natorb.irp.f

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

src/casscf_cipsi/casscf.irp.f

@@ -0,0 +1,110 @@
+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
+  n_det_max_full = 500
+  touch n_det_max_full
+  pt2_relative_error = 0.04
+  touch pt2_relative_error 
+!  call run_stochastic_cipsi
+  call run
+end
+
+subroutine run
+  implicit none
+  double precision               :: energy_old, energy, pt2_max_before, ept2_before,delta_E
+  logical                        :: converged,state_following_casscf_cipsi_save
+  integer                        :: iteration
+  converged = .False.
+
+  energy = 0.d0
+  mo_label = "MCSCF"
+  iteration = 1
+  state_following_casscf_cipsi_save = state_following_casscf
+  state_following_casscf = .True.
+  touch state_following_casscf
+  ept2_before = 0.d0
+  if(adaptive_pt2_max)then
+   pt2_max = 0.005
+   SOFT_TOUCH pt2_max
+  endif
+  do while (.not.converged)
+    print*,'pt2_max = ',pt2_max
+    call run_stochastic_cipsi
+    energy_old = energy
+    energy = eone+etwo+ecore
+    pt2_max_before = pt2_max
+
+    call write_time(6)
+    call write_int(6,iteration,'CAS-SCF iteration = ')
+    call write_double(6,energy,'CAS-SCF energy = ')
+    if(n_states == 1)then
+     double precision :: E_PT2, PT2
+     call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
+     call ezfio_get_casscf_cipsi_energy(PT2)
+     PT2 -= E_PT2
+     call write_double(6,E_PT2,'E + PT2 energy = ')
+     call write_double(6,PT2,'  PT2          = ')
+     call write_double(6,pt2_max,' PT2_MAX       = ')
+    endif
+
+    print*,''
+    call write_double(6,norm_grad_vec2,'Norm of gradients = ')
+    call write_double(6,norm_grad_vec2_tab(1), '     Core-active  gradients = ')
+    call write_double(6,norm_grad_vec2_tab(2), '     Core-virtual gradients = ')
+    call write_double(6,norm_grad_vec2_tab(3), '   Active-virtual gradients = ')
+    print*,''
+    call write_double(6,energy_improvement, 'Predicted energy improvement = ')
+
+    if(criterion_casscf == "energy")then
+     converged = dabs(energy_improvement) < thresh_scf
+    else if (criterion_casscf == "gradients")then
+     converged = norm_grad_vec2 < thresh_scf
+    else if (criterion_casscf == "e_pt2")then
+     delta_E = dabs(E_PT2 - ept2_before)
+     converged = dabs(delta_E) < thresh_casscf
+    endif
+    ept2_before = E_PT2
+    if(adaptive_pt2_max)then
+     pt2_max = dabs(energy_improvement / (pt2_relative_error))
+     pt2_max = min(pt2_max, pt2_max_before)
+     if(n_act_orb.ge.n_big_act_orb)then
+      pt2_max = max(pt2_max,pt2_min_casscf)
+     endif
+    endif
+    print*,''
+    call write_double(6,pt2_max, 'PT2_MAX for next iteration = ')
+
+    mo_coef = NewOrbs
+    mo_occ  = occnum
+    call save_mos
+    if(.not.converged)then
+     iteration += 1
+     if(norm_grad_vec2.gt.0.01d0)then
+      N_det = N_states
+     else
+      N_det = max(N_det/8 ,N_states)
+     endif
+     psi_det = psi_det_sorted
+     psi_coef = psi_coef_sorted
+     read_wf = .True.
+     call clear_mo_map
+     SOFT_TOUCH mo_coef N_det psi_det psi_coef
+     if(adaptive_pt2_max)then
+       SOFT_TOUCH pt2_max  
+     endif
+     if(iteration .gt. 3)then
+      state_following_casscf = state_following_casscf_cipsi_save
+      soft_touch state_following_casscf
+     endif
+    endif
+
+  enddo
+
+end
+
+

src/casscf_cipsi/class.irp.f

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

src/casscf_cipsi/dav_sx_mat.irp.f

@@ -0,0 +1,45 @@
+
+
+subroutine davidson_diag_sx_mat(N_st, u_in, energies)
+ implicit none
+ integer, intent(in) :: N_st
+ double precision, intent(out) :: u_in(nMonoEx+1,n_states_diag), energies(N_st)
+ integer :: i,j,N_st_tmp, dim_in, sze, N_st_diag_in
+ integer, allocatable :: list_guess(:)
+ double precision, allocatable :: H_jj(:)
+ logical    :: converged                                                                                                
+ N_st_diag_in = n_states_diag
+ provide SXmatrix
+ sze = nMonoEx+1
+ dim_in = sze
+ allocate(H_jj(sze), list_guess(sze))
+ H_jj(1) = 0.d0
+ N_st_tmp = 1
+ list_guess(1) = 1
+ do j = 2, nMonoEx+1
+  H_jj(j) = SXmatrix(j,j)
+  if(H_jj(j).lt.0.d0)then 
+   list_guess(N_st_tmp) = j
+   N_st_tmp += 1
+  endif
+ enddo
+ if(N_st_tmp .ne. N_st)then
+  print*,'Pb in davidson_diag_sx_mat'
+  print*,'N_st_tmp .ne. N_st'
+  print*,N_st_tmp, N_st
+  stop
+ endif
+ print*,'Number of possibly interesting states = ',N_st
+ print*,'Corresponding diagonal elements of the SX matrix '
+ u_in = 0.d0
+ do i = 1, min(N_st, N_st_diag_in)
+! do i = 1, N_st
+  j = list_guess(i) 
+  print*,'i,j',i,j
+  print*,'SX(i,i) = ',H_jj(j)
+  u_in(j,i) = 1.d0
+ enddo
+ call davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,SXmatrix)
+ print*,'energies = ',energies
+
+end

src/casscf_cipsi/densities.irp.f

@@ -0,0 +1,67 @@
+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
+  if(fast_2rdm)then
+   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) = 0.5d0 * state_av_act_2_rdm_spin_trace_mo(t,v,u,x)
+       enddo
+      enddo 
+     enddo
+    enddo
+   enddo
+  else
+   P0tuvx = P0tuvx_peter
+  endif
+
+END_PROVIDER

src/casscf_cipsi/densities_peter.irp.f

@@ -0,0 +1,150 @@
+use bitmasks
+
+BEGIN_PROVIDER [real*8, P0tuvx_peter, (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
+   ! matrices 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>
+   !
+   END_DOC
+   implicit none
+   integer                        :: t,u,v,x,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 density matrix P0'
+  endif
+
+   P0tuvx_peter = 0.d0
+   
+   ! first loop: we apply E_tu, once for D_tu, once for -P_tvvu
+   do mu=1,n_det
+     call det_extract(det_mu,mu,N_int)
+     do istate=1,n_states
+       cI_mu(istate)=psi_coef(mu,istate)
+     end do
+     do t=1,n_act_orb
+       ipart=list_act(t)
+       do u=1,n_act_orb
+         ihole=list_act(u)
+         ! apply E_tu
+         call det_copy(det_mu,det_mu_ex1,N_int)
+         call det_copy(det_mu,det_mu_ex2,N_int)
+         call do_spinfree_mono_excitation(det_mu,det_mu_ex1          &
+             ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
+         ! det_mu_ex1 is in the list
+         if (nu1.ne.-1) then
+           do istate=1,n_states
+             term=cI_mu(istate)*psi_coef(nu1,istate)*phase1
+             ! and we fill P0_tvvu
+             do v=1,n_act_orb
+               P0tuvx_peter(t,v,v,u)-=term
+             end do
+           end do
+         end if
+         ! det_mu_ex2 is in the list
+         if (nu2.ne.-1) then
+           do istate=1,n_states
+             term=cI_mu(istate)*psi_coef(nu2,istate)*phase2
+             do v=1,n_act_orb
+               P0tuvx_peter(t,v,v,u)-=term
+             end do
+           end do
+         end if
+       end do
+     end do
+   end do
+   ! now we do the double excitation E_tu E_vx |0>
+   do mu=1,n_det
+     call det_extract(det_mu,mu,N_int)
+     do istate=1,n_states
+       cI_mu(istate)=psi_coef(mu,istate)
+     end do
+     do v=1,n_act_orb
+       ipart=list_act(v)
+       do x=1,n_act_orb
+         ihole=list_act(x)
+         ! apply E_vx
+         call det_copy(det_mu,det_mu_ex1,N_int)
+         call det_copy(det_mu,det_mu_ex2,N_int)
+         call do_spinfree_mono_excitation(det_mu,det_mu_ex1          &
+             ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
+         ! we apply E_tu to the first resultant determinant, thus E_tu E_vx |0>
+         if (ierr1.eq.1) then
+           do t=1,n_act_orb
+             jpart=list_act(t)
+             do u=1,n_act_orb
+               jhole=list_act(u)
+               call det_copy(det_mu_ex1,det_mu_ex11,N_int)
+               call det_copy(det_mu_ex1,det_mu_ex12,N_int)
+               call do_spinfree_mono_excitation(det_mu_ex1,det_mu_ex11&
+                   ,det_mu_ex12,nu11,nu12,jhole,jpart,phase11,phase12,ierr11,ierr12)
+               if (nu11.ne.-1) then
+                 do istate=1,n_states
+                   P0tuvx_peter(t,u,v,x)+=cI_mu(istate)*psi_coef(nu11,istate)&
+                       *phase11*phase1
+                 end do
+               end if
+               if (nu12.ne.-1) then
+                 do istate=1,n_states
+                   P0tuvx_peter(t,u,v,x)+=cI_mu(istate)*psi_coef(nu12,istate)&
+                       *phase12*phase1
+                 end do
+               end if
+             end do
+           end do
+         end if
+         
+         ! we apply E_tu to the second resultant determinant
+         if (ierr2.eq.1) then
+           do t=1,n_act_orb
+             jpart=list_act(t)
+             do u=1,n_act_orb
+               jhole=list_act(u)
+               call det_copy(det_mu_ex2,det_mu_ex21,N_int)
+               call det_copy(det_mu_ex2,det_mu_ex22,N_int)
+               call do_spinfree_mono_excitation(det_mu_ex2,det_mu_ex21&
+                   ,det_mu_ex22,nu21,nu22,jhole,jpart,phase21,phase22,ierr21,ierr22)
+               if (nu21.ne.-1) then
+                 do istate=1,n_states
+                   P0tuvx_peter(t,u,v,x)+=cI_mu(istate)*psi_coef(nu21,istate)&
+                       *phase21*phase2
+                 end do
+               end if
+               if (nu22.ne.-1) then
+                 do istate=1,n_states
+                   P0tuvx_peter(t,u,v,x)+=cI_mu(istate)*psi_coef(nu22,istate)&
+                       *phase22*phase2
+                 end do
+               end if
+             end do
+           end do
+         end if
+         
+       end do
+     end do
+   end do
+   
+   ! we average by just dividing by the number of 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
+           P0tuvx_peter(t,u,v,x)*=0.5D0/dble(N_states)
+         end do
+       end do
+     end do
+   end do
+   
+END_PROVIDER

src/casscf_cipsi/det_manip.irp.f

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

src/casscf_cipsi/driver_optorb.irp.f

@@ -0,0 +1,3 @@
+subroutine driver_optorb
+  implicit none
+end 

src/casscf_cipsi/get_energy.irp.f

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

src/casscf_cipsi/grad_old.irp.f

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

src/casscf_cipsi/gradient.irp.f

@@ -0,0 +1,215 @@
+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, n_c_a_prov]
+&BEGIN_PROVIDER [integer, n_c_v_prov]
+&BEGIN_PROVIDER [integer, n_a_v_prov]
+  implicit none
+  n_c_a_prov = n_core_inact_orb * n_act_orb
+  n_c_v_prov = n_core_inact_orb * n_virt_orb
+  n_a_v_prov = n_act_orb * n_virt_orb 
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
+&BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)]
+&BEGIN_PROVIDER [integer, list_idx_c_a, (3,n_c_a_prov) ]
+&BEGIN_PROVIDER [integer, list_idx_c_v, (3,n_c_v_prov) ]
+&BEGIN_PROVIDER [integer, list_idx_a_v, (3,n_a_v_prov) ]
+&BEGIN_PROVIDER [integer, mat_idx_c_a, (n_core_inact_orb,n_act_orb)
+&BEGIN_PROVIDER [integer, mat_idx_c_v, (n_core_inact_orb,n_virt_orb)
+&BEGIN_PROVIDER [integer, mat_idx_a_v, (n_act_orb,n_virt_orb)
+  BEGIN_DOC
+  ! a list of the orbitals involved in the excitation
+  END_DOC
+  
+  implicit none
+  integer                        :: i,t,a,ii,tt,aa,indx,indx_tmp
+  indx=0
+  indx_tmp = 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'
+      indx_tmp += 1
+      list_idx_c_a(1,indx_tmp) = indx
+      list_idx_c_a(2,indx_tmp) = ii
+      list_idx_c_a(3,indx_tmp) = tt
+      mat_idx_c_a(ii,tt) = indx
+    end do
+  end do
+  
+  indx_tmp = 0
+  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'
+      indx_tmp += 1
+      list_idx_c_v(1,indx_tmp) = indx
+      list_idx_c_v(2,indx_tmp) = ii
+      list_idx_c_v(3,indx_tmp) = aa 
+      mat_idx_c_v(ii,aa) = indx
+    end do
+  end do
+  
+  indx_tmp = 0
+  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'
+      indx_tmp += 1
+      list_idx_a_v(1,indx_tmp) = indx
+      list_idx_a_v(2,indx_tmp) = tt
+      list_idx_a_v(3,indx_tmp) = aa
+      mat_idx_a_v(tt,aa) = indx
+    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_PROVIDER [real*8, norm_grad_vec2]
+&BEGIN_PROVIDER [real*8, norm_grad_vec2_tab, (3)]
+  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
+  
+  indx=0
+  norm_grad_vec2_tab = 0.d0
+  do i=1,n_core_inact_orb
+    do t=1,n_act_orb
+      indx+=1
+      gradvec2(indx)=gradvec_it(i,t)
+      norm_grad_vec2_tab(1) += gradvec2(indx)*gradvec2(indx)
+    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)
+      norm_grad_vec2_tab(2) += gradvec2(indx)*gradvec2(indx)
+    end do
+  end do
+  
+  do t=1,n_act_orb
+    do a=1,n_virt_orb
+      indx+=1
+      gradvec2(indx)=gradvec_ta(t,a)
+      norm_grad_vec2_tab(3) += gradvec2(indx)*gradvec2(indx)
+    end do
+  end do
+  
+  norm_grad_vec2=0.d0
+  do indx=1,nMonoEx
+    norm_grad_vec2+=gradvec2(indx)*gradvec2(indx)
+  end do
+  do i = 1, 3
+   norm_grad_vec2_tab(i) = dsqrt(norm_grad_vec2_tab(i))
+  enddo
+  norm_grad_vec2=sqrt(norm_grad_vec2)
+  if(bavard)then
+   write(6,*)
+   write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad_vec2
+   write(6,*)
+  endif
+  
+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 ! active 
+    vv=list_act(v)
+    do x=1,n_act_orb ! active 
+      x3=x+n_core_inact_orb ! list_act(x) 
+      do y=1,n_act_orb ! active
+        y3=y+n_core_inact_orb ! list_act(y)
+        !                Gamma(2)  a a a a      1/r12      i  a  a  a
+        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
+

src/casscf_cipsi/hessian.irp.f

@@ -0,0 +1,539 @@
+use bitmasks
+
+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
+
+
+BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
+ implicit none
+ integer                        :: i,j,t,u,a,b
+ integer                        :: indx,indx_tmp, jndx, jndx_tmp
+ integer                        :: ustart,bstart
+ real*8                         :: hessmat_itju
+ real*8                         :: hessmat_itja
+ real*8                         :: hessmat_itua
+ real*8                         :: hessmat_iajb
+ real*8                         :: hessmat_iatb
+ real*8                         :: hessmat_taub
+ !       c-a c-v a-v
+ !  c-a | X   X  X
+ !  c-v |     X  X 
+ !  a-v |        X
+
+  provide mo_two_e_integrals_in_map
+
+ hessmat = 0.d0
+
+ !$OMP PARALLEL DEFAULT(NONE) &
+ !$OMP SHARED(hessmat,n_c_a_prov,list_idx_c_a,n_core_inact_orb,n_act_orb,mat_idx_c_a) &
+ !$OMP PRIVATE(indx_tmp,indx,i,t,j,u,ustart,jndx)
+
+ !$OMP DO
+!!!! < Core-active| H |Core-active >
+ ! Core-active excitations 
+ do indx_tmp = 1, n_c_a_prov
+  indx = list_idx_c_a(1,indx_tmp)
+  i    = list_idx_c_a(2,indx_tmp)
+  t    = list_idx_c_a(3,indx_tmp)
+  ! Core-active excitations 
+  do j = 1, n_core_inact_orb
+   if (i.eq.j) then
+     ustart=t
+   else
+     ustart=1
+   end if
+   do u=ustart,n_act_orb
+    jndx = mat_idx_c_a(j,u)
+    hessmat(jndx,indx) = hessmat_itju(i,t,j,u)
+    hessmat(indx,jndx) = hessmat(jndx,indx)
+   enddo
+  enddo
+ enddo
+  !$OMP END DO NOWAIT
+  !$OMP END PARALLEL
+
+ !$OMP PARALLEL DEFAULT(NONE) &
+ !$OMP SHARED(hessmat,n_c_a_prov,n_c_v_prov,list_idx_c_a,list_idx_c_v) &
+ !$OMP PRIVATE(indx_tmp,jndx_tmp,indx,i,t,j,a,jndx)
+
+ !$OMP DO
+!!!! < Core-active| H |Core-VIRTUAL >
+ ! Core-active excitations 
+ do indx_tmp = 1, n_c_a_prov
+  indx = list_idx_c_a(1,indx_tmp)
+  i    = list_idx_c_a(2,indx_tmp)
+  t    = list_idx_c_a(3,indx_tmp)
+  ! Core-VIRTUAL excitations 
+  do jndx_tmp = 1, n_c_v_prov
+   jndx = list_idx_c_v(1,jndx_tmp)
+   j    = list_idx_c_v(2,jndx_tmp)
+   a    = list_idx_c_v(3,jndx_tmp)
+   hessmat(jndx,indx) = hessmat_itja(i,t,j,a)
+   hessmat(indx,jndx) = hessmat(jndx,indx)
+  enddo
+ enddo
+  !$OMP END DO NOWAIT
+  !$OMP END PARALLEL
+
+ !$OMP PARALLEL DEFAULT(NONE) &
+ !$OMP SHARED(hessmat,n_c_a_prov,n_a_v_prov,list_idx_c_a,list_idx_a_v) &
+ !$OMP PRIVATE(indx_tmp,jndx_tmp,indx,i,t,u,a,jndx)
+
+ !$OMP DO
+!!!! < Core-active| H |ACTIVE-VIRTUAL >
+ ! Core-active excitations 
+ do indx_tmp = 1, n_c_a_prov
+  indx = list_idx_c_a(1,indx_tmp)
+  i    = list_idx_c_a(2,indx_tmp)
+  t    = list_idx_c_a(3,indx_tmp)
+  ! ACTIVE-VIRTUAL excitations 
+  do jndx_tmp = 1, n_a_v_prov
+   jndx = list_idx_a_v(1,jndx_tmp)
+   u    = list_idx_a_v(2,jndx_tmp)
+   a    = list_idx_a_v(3,jndx_tmp)
+   hessmat(jndx,indx) = hessmat_itua(i,t,u,a)
+   hessmat(indx,jndx) = hessmat(jndx,indx)
+  enddo
+ enddo
+
+  !$OMP END DO NOWAIT
+  !$OMP END PARALLEL
+
+
+ if(hess_cv_cv)then
+ !$OMP PARALLEL DEFAULT(NONE) &
+ !$OMP SHARED(hessmat,n_c_v_prov,list_idx_c_v,n_core_inact_orb,n_virt_orb,mat_idx_c_v) &
+ !$OMP PRIVATE(indx_tmp,indx,i,a,j,b,bstart,jndx)
+  !$OMP DO
+!!!!! < Core-VIRTUAL | H |Core-VIRTUAL >
+  ! Core-VIRTUAL excitations 
+  do indx_tmp = 1, n_c_v_prov
+   indx = list_idx_c_v(1,indx_tmp)
+   i    = list_idx_c_v(2,indx_tmp)
+   a    = list_idx_c_v(3,indx_tmp)
+   ! Core-VIRTUAL excitations 
+   do j = 1, n_core_inact_orb
+    if (i.eq.j) then
+      bstart=a
+    else
+      bstart=1
+    end if
+    do b=bstart,n_virt_orb
+     jndx = mat_idx_c_v(j,b)
+     hessmat(jndx,indx) = hessmat_iajb(i,a,j,b)
+     hessmat(indx,jndx) = hessmat(jndx,indx)
+    enddo
+   enddo
+  enddo
+ 
+  !$OMP END DO NOWAIT
+  !$OMP END PARALLEL
+ endif
+
+ !$OMP PARALLEL DEFAULT(NONE) &
+ !$OMP SHARED(hessmat,n_c_v_prov,n_a_v_prov,list_idx_c_v,list_idx_a_v) &
+ !$OMP PRIVATE(indx_tmp,jndx_tmp,indx,i,a,t,b,jndx)
+
+ !$OMP DO
+!!!! < Core-VIRTUAL | H |Active-VIRTUAL >
+ ! Core-VIRTUAL excitations 
+ do indx_tmp = 1, n_c_v_prov
+  indx = list_idx_c_v(1,indx_tmp)
+  i    = list_idx_c_v(2,indx_tmp)
+  a    = list_idx_c_v(3,indx_tmp)
+  ! Active-VIRTUAL excitations 
+  do jndx_tmp = 1, n_a_v_prov
+   jndx = list_idx_a_v(1,jndx_tmp)
+   t    = list_idx_a_v(2,jndx_tmp)
+   b    = list_idx_a_v(3,jndx_tmp)
+   hessmat(jndx,indx) = hessmat_iatb(i,a,t,b)
+   hessmat(indx,jndx) = hessmat(jndx,indx)
+  enddo
+ enddo
+  !$OMP END DO NOWAIT
+  !$OMP END PARALLEL
+
+
+ !$OMP PARALLEL DEFAULT(NONE) &
+ !$OMP SHARED(hessmat,n_a_v_prov,list_idx_a_v,n_act_orb,n_virt_orb,mat_idx_a_v) &
+ !$OMP PRIVATE(indx_tmp,indx,t,a,u,b,bstart,jndx)
+
+ !$OMP DO
+!!!! < Active-VIRTUAL | H |Active-VIRTUAL >
+ ! Active-VIRTUAL excitations 
+ do indx_tmp = 1, n_a_v_prov
+  indx = list_idx_a_v(1,indx_tmp)
+  t    = list_idx_a_v(2,indx_tmp)
+  a    = list_idx_a_v(3,indx_tmp)
+  ! Active-VIRTUAL excitations 
+  do u=t,n_act_orb
+   if (t.eq.u) then
+     bstart=a
+   else
+     bstart=1
+   end if
+   do b=bstart,n_virt_orb
+    jndx = mat_idx_a_v(u,b)
+    hessmat(jndx,indx) = hessmat_taub(t,a,u,b)
+    hessmat(indx,jndx) = hessmat(jndx,indx)
+   enddo
+  enddo
+ enddo
+  !$OMP END DO NOWAIT
+  !$OMP END PARALLEL
+
+END_PROVIDER 

src/casscf_cipsi/hessian_old.irp.f

@@ -0,0 +1,310 @@
+
+use bitmasks
+BEGIN_PROVIDER [real*8, hessmat_old, (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_old '
+    write(6,*) '  nMonoEx = ',nMonoEx
+  endif
+  
+  do indx=1,nMonoEx
+    do jndx=1,nMonoEx
+      hessmat_old(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_old(indx,jndx)=res
+      hessmat_old(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, hessmat_peter, (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 hessmat_peter '
+    write(6,*) '  nMonoEx = ',nMonoEx
+  endif
+  provide mo_two_e_integrals_in_map
+  
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP SHARED(hessmat_peter,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
+  ! (DOUBLY OCCUPIED ---> ACT ) 
+  do i=1,n_core_inact_orb
+    do t=1,n_act_orb
+      indx = t + (i-1)*n_act_orb
+      jndx=indx
+      ! (DOUBLY OCCUPIED ---> ACT )
+      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
+          hessmat_peter(jndx,indx)=hessmat_itju(i,t,j,u)
+          jndx+=1
+        end do
+      end do
+      ! (DOUBLY OCCUPIED ---> VIRTUAL) 
+      do j=1,n_core_inact_orb
+        do a=1,n_virt_orb
+          hessmat_peter(jndx,indx)=hessmat_itja(i,t,j,a)
+          jndx+=1
+        end do
+      end do
+      ! (ACTIVE ---> VIRTUAL) 
+      do u=1,n_act_orb
+        do a=1,n_virt_orb
+          hessmat_peter(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
+  ! (DOUBLY OCCUPIED ---> VIRTUAL) 
+  do a=1,n_virt_orb
+    do i=1,n_core_inact_orb
+      indx = a + (i-1)*n_virt_orb + indx_shift
+      jndx=indx
+      ! (DOUBLY OCCUPIED ---> VIRTUAL) 
+      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
+          hessmat_peter(jndx,indx)=hessmat_iajb(i,a,j,b)
+          jndx+=1
+        end do
+      end do
+      ! (ACT ---> VIRTUAL) 
+      do t=1,n_act_orb
+        do b=1,n_virt_orb
+          hessmat_peter(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 
+  ! (ACT ---> VIRTUAL) 
+  do a=1,n_virt_orb
+    do t=1,n_act_orb
+      indx = a + (t-1)*n_virt_orb + indx_shift
+      jndx=indx
+      ! (ACT ---> VIRTUAL) 
+      do u=t,n_act_orb
+        if (t.eq.u) then
+          bstart=a
+        else
+          bstart=1
+        end if
+        do b=bstart,n_virt_orb
+          hessmat_peter(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
+      hessmat_peter(indx,jndx) = hessmat_peter(jndx,indx)
+    enddo
+  enddo
+
+  
+END_PROVIDER
+

src/casscf_cipsi/mcscf_fock.irp.f

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

src/casscf_cipsi/natorb.irp.f

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

src/casscf_cipsi/neworbs.irp.f

@@ -0,0 +1,253 @@
+ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
+&BEGIN_PROVIDER [integer, n_guess_sx_mat ]
+  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
+  if(diag_hess_cas)then  
+   do i = 1, nMonoEx
+    SXmatrix(i+1,i+1) = hessdiag(i)
+   enddo
+  else
+   do i=1,nMonoEx
+     do j=1,nMonoEx
+       SXmatrix(i+1,j+1)=hessmat(i,j)
+       SXmatrix(j+1,i+1)=hessmat(i,j)
+     end do
+   end do
+  endif
+  
+  do i = 1, nMonoEx
+   SXmatrix(i+1,i+1) += level_shift_casscf
+  enddo
+  n_guess_sx_mat = 1
+  do i = 1, nMonoEx
+   if(SXmatrix(i+1,i+1).lt.0.d0 )then
+    n_guess_sx_mat += 1
+   endif
+  enddo
+  if (bavard) then
+    do i=2,nMonoEx
+      write(6,*) ' diagonal of the Hessian : ',i,hessmat(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
+  if(nMonoEx+1.gt.n_det_max_full)then
+   if(bavard)then
+    print*,'Using the Davidson algorithm to diagonalize the SXmatrix'
+   endif
+   double precision, allocatable :: u_in(:,:),energies(:)
+   allocate(u_in(nMonoEx+1,n_states_diag),energies(n_guess_sx_mat))
+   call davidson_diag_sx_mat(n_guess_sx_mat, u_in, energies)
+   integer :: i,j
+   SXeigenvec = 0.d0
+   SXeigenval = 0.d0
+   do i = 1, n_guess_sx_mat
+    SXeigenval(i) = energies(i)
+    do j = 1, nMonoEx+1
+     SXeigenvec(j,i) = u_in(j,i)
+    enddo
+   enddo
+  else
+   if(bavard)then
+    print*,'Diagonalize the SXmatrix with Jacobi'
+   endif
+   call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
+  endif
+  if (bavard) then
+    write(6,*) ' SXdiag : lowest eigenvalues '
+    write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1)
+    if(n_guess_sx_mat.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 (dabs(SXeigenvec(1,i)).gt.best_overlap_casscf) then
+        best_overlap_casscf=dabs(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
+  
+
+

src/casscf_cipsi/reorder_orb.irp.f

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

src/casscf_cipsi/save_energy.irp.f

@@ -0,0 +1,9 @@
+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_cipsi_energy(E(1:N_states))
+  call ezfio_set_casscf_cipsi_energy_pt2(E(1:N_states)+pt2(1:N_states))
+end

src/casscf_cipsi/superci_dm.irp.f

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

src/casscf_cipsi/swap_orb.irp.f

@@ -0,0 +1,132 @@
+ 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.hessmat(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.hessmat(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 

src/casscf_cipsi/tot_en.irp.f

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

src/casscf_tc_bi/grad_dm.irp.f

@@ -0,0 +1,132 @@
+ BEGIN_PROVIDER [real*8, gradvec_tc_r, (0:3,nMonoEx)]
+&BEGIN_PROVIDER [real*8, gradvec_tc_l,  (0:3,nMonoEx)]
+ BEGIN_DOC
+! gradvec_tc_r(0:3,i) = <Chi_0| H E_q^p(i) |Phi_0> 
+! 
+! gradvec_tc_l(0:3,i) = <Chi_0| E_p^q(i) H |Phi_0> 
+!
+! where the indices "i" corresponds to E_q^p(i) 
+! 
+! i = mat_idx_c_a(q,p) 
+!
+! and gradvec_tc_r/l(0) = full matrix element
+!
+!     gradvec_tc_r/l(1) = one-body part 
+
+!     gradvec_tc_r/l(2) = two-body part 
+
+!     gradvec_tc_r/l(3) = three-body part 
+ END_DOC
+ implicit none
+ integer                        :: ii,tt,aa,indx
+ integer                        :: i,t,a,fff
+ double precision :: res_l(0:3), res_r(0:3)
+ gradvec_tc_l = 0.d0
+ gradvec_tc_r = 0.d0
+ ! computing the core/inactive --> virtual orbitals gradients
+  do i=1,n_core_inact_orb
+    ii=list_core_inact(i)
+    do t=1,n_act_orb
+      tt=list_act(t)
+      indx = mat_idx_c_a(i,t) 
+      call gradvec_tc_it(ii,tt,res_l,res_r)
+      do fff = 0,3
+       gradvec_tc_l(fff,indx)=res_l(fff)
+       gradvec_tc_r(fff,indx)=res_r(fff)
+      enddo
+    end do
+  end do
+  
+  do i=1,n_core_inact_orb
+    ii=list_core_inact(i)
+    do a=1,n_virt_orb
+      indx = mat_idx_c_v(i,a) 
+      aa=list_virt(a)
+      call gradvec_tc_ia(ii,aa,res_l,res_r)
+      do fff = 0,3
+       gradvec_tc_l(fff,indx)=res_l(fff)
+       gradvec_tc_r(fff,indx)=res_r(fff)
+      enddo
+    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 = mat_idx_a_v(t,a) 
+      call gradvec_tc_ta(tt,aa,res_l, res_r)
+      do fff = 0,3
+       gradvec_tc_l(fff,indx)=res_l(fff)
+       gradvec_tc_r(fff,indx)=res_r(fff)
+      enddo
+    end do
+  end do
+END_PROVIDER 
+
+subroutine gradvec_tc_ia(i,a,res_l, res_r)
+ implicit none
+ BEGIN_DOC
+! doubly occupied --> virtual TC gradient 
+!
+! Corresponds to res_r = <X0|H E_i^a|Phi_0>, 
+!
+!                res_l = <X0|E_a^i H|Phi_0>
+ END_DOC
+ integer, intent(in) :: i,a
+ double precision, intent(out) :: res_l(0:3), res_r(0:3)
+ res_l = 0.d0
+ res_r = 0.d0
+ res_l(1) = -2 * mo_bi_ortho_tc_one_e(a,i)
+ res_r(1) = -2 * mo_bi_ortho_tc_one_e(i,a)
+ 
+end
+
+subroutine gradvec_tc_it(i,t,res_l, res_r)
+ implicit none
+ BEGIN_DOC
+! doubly occupied --> active TC gradient 
+!
+! Corresponds to res_r = <X0|H E_i^t|Phi_0>
+!
+!                res_l = <X0|E_t^i H |Phi_0>
+ END_DOC
+ integer, intent(in) :: i,t
+ double precision, intent(out) :: res_l(0:3),res_r(0:3)
+ integer :: rr,r,ss,s,m,mm
+ double precision :: dm
+ res_r = 0.d0
+ res_l = 0.d0
+ res_r(1)  += -2.d0 * mo_bi_ortho_tc_one_e(i,t) 
+ res_l(1)  +=  2.D0 * mo_bi_ortho_tc_one_e(t,i) 
+
+ do rr = 1, n_act_orb
+  r = list_act(rr)
+  res_r(1) +=  mo_bi_ortho_tc_one_e(i,r) * tc_transition_matrix_mo(t,r,1,1)   
+  res_l(1) += -mo_bi_ortho_tc_one_e(r,i) * tc_transition_matrix_mo(r,t,1,1)
+ enddo
+ 
+end
+
+subroutine gradvec_tc_ta(t,a,res_l, res_r)
+ implicit none
+ BEGIN_DOC
+! active --> virtual TC gradient 
+!
+! Corresponds to res_r = <X0|H E_t^a|Phi_0>
+!
+!                res_l = <X0|E_a^t H |Phi_0>
+ END_DOC
+ integer, intent(in) :: t,a
+ double precision, intent(out) :: res_l(0:3),res_r(0:3)
+ integer :: rr,r,m
+ double precision :: dm
+ res_r = 0.d0
+ res_l = 0.d0
+ do rr = 1, n_act_orb
+  r = list_act(rr) 
+  res_l(1) +=  mo_bi_ortho_tc_one_e(a,r) * tc_transition_matrix_mo(t,r,1,1)
+  res_r(1) += -mo_bi_ortho_tc_one_e(r,a) * tc_transition_matrix_mo(r,t,1,1)
+ enddo
+ 
+end

src/casscf_tc_bi/grad_old.irp.f

@@ -0,0 +1,132 @@
+
+ BEGIN_PROVIDER [real*8, gradvec_detail_right_old, (0:3,nMonoEx)]
+&BEGIN_PROVIDER [real*8, gradvec_detail_left_old,  (0:3,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,ll
+  real*8                         :: res_l(0:3), res_r(0:3)
+  
+ do ii = 1, n_core_inact_orb
+  ihole = list_core_inact(ii)
+  do aa = 1, n_virt_orb
+   ipart = list_virt(aa)
+   indx = mat_idx_c_v(ii,aa) 
+   call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
+   do ll = 0, 3
+    gradvec_detail_left_old (ll,indx)=res_l(ll)
+    gradvec_detail_right_old(ll,indx)=res_r(ll)
+   enddo
+  enddo
+ enddo
+
+ do ii = 1, n_core_inact_orb
+  ihole = list_core_inact(ii)
+  do tt = 1, n_act_orb
+   ipart = list_act(tt)
+   indx = mat_idx_c_a(ii,tt) 
+   call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
+   do ll = 0, 3
+    gradvec_detail_left_old (ll,indx)=res_l(ll)
+    gradvec_detail_right_old(ll,indx)=res_r(ll)
+   enddo
+  enddo
+ enddo
+
+ do tt = 1, n_act_orb
+  ihole = list_act(tt)
+  do aa = 1, n_virt_orb
+   ipart = list_virt(aa)
+   indx = mat_idx_a_v(tt,aa) 
+   call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
+   do ll = 0, 3
+    gradvec_detail_left_old (ll,indx)=res_l(ll)
+    gradvec_detail_right_old(ll,indx)=res_r(ll)
+   enddo
+  enddo
+ enddo
+  
+  real*8                         :: norm_grad_left, norm_grad_right
+  norm_grad_left=0.d0
+  norm_grad_right=0.d0
+  do indx=1,nMonoEx
+    norm_grad_left+=gradvec_detail_left_old(0,indx)*gradvec_detail_left_old(0,indx)
+    norm_grad_right+=gradvec_detail_right_old(0,indx)*gradvec_detail_right_old(0,indx)
+  end do
+  norm_grad_left=sqrt(norm_grad_left)
+  norm_grad_right=sqrt(norm_grad_right)
+!  if (bavard) then
+    write(6,*)
+    write(6,*) ' Norm of the LEFT  orbital gradient (via <0|EH|0>) : ', norm_grad_left
+    write(6,*) ' Norm of the RIGHT orbital gradient (via <0|HE|0>) : ', norm_grad_right
+    write(6,*)
+!  endif
+  
+  
+END_PROVIDER
+
+subroutine calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
+  BEGIN_DOC
+  ! Computes the gradient with respect to orbital rotation BRUT FORCE
+  ! 
+  ! res_l = <Chi| E_qp H^tc | Phi>
+  !
+  ! res_r = <Chi| H^tc E_pq | Phi>
+  !
+  ! q=hole, p=particle. NOTE that on res_l it is E_qp and on res_r it is E_pq 
+  ! 
+  ! res_l(0) =   total matrix element, res_l(1) =   one-electron part,
+  ! 
+  ! res_l(2) =   two-electron part,  res_l(3) = three-electron part 
+  ! 
+  END_DOC
+  implicit none
+  integer, intent(in)            :: ihole,ipart
+  double precision, intent(out)  :: res_l(0:3), res_r(0:3)
+  integer                        :: mu,iii,ispin,ierr,nu,istate,ll
+  integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:)
+  real*8                         :: chi_H_mu_ex_array(0:3,N_states),mu_ex_H_phi_array(0:3,N_states),phase
+  allocate(det_mu(N_int,2))
+  allocate(det_mu_ex(N_int,2))
+  
+  res_l=0.D0
+  res_r=0.D0
+  
+  do mu=1,n_det
+    ! get the string of the determinant |mu>
+    call det_extract(det_mu,mu,N_int)
+    do ispin=1,2
+      ! do the monoexcitation on it: |det_mu_ex> = a^dagger_{p,ispin} a_{q,ispin} |mu> 
+      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)
+      ! |det_mu_ex> = a^dagger_{p,ispin} a_{q,ispin} |mu> 
+      if (ierr.eq.1) then
+        call i_H_tc_psi_phi(det_mu_ex,psi_det,psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,N_int & 
+            ,N_det,psi_det_size,N_states,chi_H_mu_ex_array,mu_ex_H_phi_array)
+        ! chi_H_mu_ex_array = <Chi|H E_qp |mu >
+        ! mu_ex_H_phi_array = <mu |E_qp H |Phi>
+        do istate=1,N_states
+         do ll = 0,3 ! loop over the body components (1e,2e,3e)
+          !res_l =  \sum_mu c_mu^l <mu|E_qp H |Phi> = <Chi|E_qp H |Phi>
+          res_l(ll)+= mu_ex_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase 
+          !res_r =  \sum_mu c_mu^r <Chi|H E_qp |mu> = <Chi|H E_qp |Phi>
+          res_r(ll)+= chi_H_mu_ex_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase
+         enddo
+        end do
+      end if
+    end do
+  end do
+  
+  ! state-averaged gradient
+  res_l*=1.d0/dble(N_states)
+  res_r*=1.d0/dble(N_states)
+  
+end 
+

src/casscf_tc_bi/gradient.irp.f

@@ -0,0 +1,94 @@
+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, n_c_a_prov]
+&BEGIN_PROVIDER [integer, n_c_v_prov]
+&BEGIN_PROVIDER [integer, n_a_v_prov]
+  implicit none
+  n_c_a_prov = n_core_inact_orb * n_act_orb
+  n_c_v_prov = n_core_inact_orb * n_virt_orb
+  n_a_v_prov = n_act_orb * n_virt_orb 
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
+&BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)]
+&BEGIN_PROVIDER [integer, list_idx_c_a, (3,n_c_a_prov) ]
+&BEGIN_PROVIDER [integer, list_idx_c_v, (3,n_c_v_prov) ]
+&BEGIN_PROVIDER [integer, list_idx_a_v, (3,n_a_v_prov) ]
+&BEGIN_PROVIDER [integer, mat_idx_c_a, (n_core_inact_orb,n_act_orb)
+&BEGIN_PROVIDER [integer, mat_idx_c_v, (n_core_inact_orb,n_virt_orb)
+&BEGIN_PROVIDER [integer, mat_idx_a_v, (n_act_orb,n_virt_orb)
+  BEGIN_DOC
+  ! a list of the orbitals involved in the excitation
+  END_DOC
+  
+  implicit none
+  integer                        :: i,t,a,ii,tt,aa,indx,indx_tmp
+  indx=0
+  indx_tmp = 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'
+      indx_tmp += 1
+      list_idx_c_a(1,indx_tmp) = indx
+      list_idx_c_a(2,indx_tmp) = ii
+      list_idx_c_a(3,indx_tmp) = tt
+      mat_idx_c_a(ii,tt) = indx
+    end do
+  end do
+  
+  indx_tmp = 0
+  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'
+      indx_tmp += 1
+      list_idx_c_v(1,indx_tmp) = indx
+      list_idx_c_v(2,indx_tmp) = ii
+      list_idx_c_v(3,indx_tmp) = aa 
+      mat_idx_c_v(ii,aa) = indx
+    end do
+  end do
+  
+  indx_tmp = 0
+  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'
+      indx_tmp += 1
+      list_idx_a_v(1,indx_tmp) = indx
+      list_idx_a_v(2,indx_tmp) = tt
+      list_idx_a_v(3,indx_tmp) = aa
+      mat_idx_a_v(tt,aa) = indx
+    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

src/ccsd/ccsd_space_orb_sub_chol.irp.f

@@ -1024,56 +1024,26 @@ subroutine compute_A1_chol(nO,nV,t1,t2,tau,A1)
 
   integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta
 
-  double precision, allocatable :: X_vooo(:,:,:,:), Y_oooo(:,:,:,:)
-  allocate(X_vooo(nV,nO,nO,nO), Y_oooo(nO,nO,nO,nO))
+  double precision, allocatable :: Y_oooo(:,:,:,:)
+  allocate(Y_oooo(nO,nO,nO,nO))
 
   ! A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j)
-  !$omp parallel &
-  !$omp shared(nO,nV,A1,cc_space_v_oooo,cc_space_v_ovoo,X_vooo) &
-  !$omp private(u,v,i,j) &
-  !$omp default(none)
-  !$omp do collapse(2)
-  do j = 1, nO
-    do i = 1, nO
-      do v = 1, nO
-        do u = 1, nO
-          A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do nowait
-
   ! A1(u,v,i,j) += cc_space_v_ovoo(u,a,i,j) * t1(v,a) &
 
-  !$omp do collapse(2)
-  do j = 1, nO
-    do i = 1, nO
-      do u = 1, nO
-        do a = 1, nV
-          X_vooo(a,u,i,j) = cc_space_v_ovoo(u,a,i,j)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-
   call dgemm('N','N', nO, nO*nO*nO, nV, &
              1d0, t1    , size(t1,1), &
-                  X_vooo, size(X_vooo,1), &
+                  cc_space_v_vooo, size(cc_space_v_vooo,1), &
              0d0, Y_oooo, size(Y_oooo,1))
 
   !$omp parallel &
-  !$omp shared(nO,nV,A1,Y_oooo) &
   !$omp private(u,v,i,j) &
-  !$omp default(none)
+  !$omp default(shared)
   !$omp do collapse(2)
   do j = 1, nO
     do i = 1, nO
       do v = 1, nO
         do u = 1, nO
-          A1(u,v,i,j) = A1(u,v,i,j) + Y_oooo(v,u,i,j)
+          A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j) + Y_oooo(v,u,j,i) + Y_oooo(u,v,i,j)
         enddo
       enddo
     enddo
@@ -1081,13 +1051,7 @@ subroutine compute_A1_chol(nO,nV,t1,t2,tau,A1)
   !$omp end do
   !$omp end parallel
 
-  deallocate(X_vooo,Y_oooo)
-
-  ! A1(u,v,i,j) += cc_space_v_vooo(a,v,i,j) * t1(u,a)
-  call dgemm('N','N', nO, nO*nO*nO, nV, &
-             1d0, t1      , size(t1,1), &
-                  cc_space_v_vooo, size(cc_space_v_vooo,1), &
-             1d0, A1      , size(A1,1))
+  deallocate(Y_oooo)
 
   ! A1(u,v,i,j) += cc_space_v_vvoo(a,b,i,j) * tau(u,v,a,b)
   call dgemm('N','N', nO*nO, nO*nO, nV*nV, &

src/cipsi_tc_bi_ortho/get_d0_good.irp.f

@@ -53,7 +53,7 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         if (hij == (0.d0,0.d0)) cycle
         !DIR$ LOOP COUNT AVG(4)
         do k=1,N_states
-          mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * hij  ! HOTSPOT
+          mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,2) * hij  ! HOTSPOT
         enddo
       end do
       !!!!!!!!!! <phi|H|alpha>
@@ -72,7 +72,7 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         if (hji == (0.d0,0.d0)) cycle
         !DIR$ LOOP COUNT AVG(4)
         do k=1,N_states
-          mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,2) * hji  ! HOTSPOT
+          mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,1) * hji  ! HOTSPOT
         enddo
       end do
     end do
@@ -109,7 +109,7 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         end if
         !DIR$ LOOP COUNT AVG(4)
         do k=1,N_states
-          mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
+          mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
         enddo
       end do
 
@@ -128,7 +128,7 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         end if
         !DIR$ LOOP COUNT AVG(4)
         do k=1,N_states
-          mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
+          mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
         enddo
       end do
     end do

src/cipsi_tc_bi_ortho/get_d1_good.irp.f

@@ -76,7 +76,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,1)
+            tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,2)
           enddo
         endif
       end do
@@ -88,7 +88,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,1)
+            tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,2)
           enddo
         endif
       end do
@@ -114,7 +114,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hji = hji * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,2)
+            tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,1)
           enddo
         endif
       end do
@@ -126,7 +126,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hji = hji * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,2)
+            tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,1)
           enddo
         endif
       end do
@@ -169,7 +169,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,1)
+            tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,2)
           enddo
         endif
       end if
@@ -180,7 +180,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         if (hij /= 0.d0) then
           hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
           do k=1,N_states
-            tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,1)
+            tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,2)
           enddo
         endif
       end if
@@ -211,7 +211,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hji = hji * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,2)
+            tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,1)
           enddo
         endif
       end if
@@ -222,7 +222,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         if (hji /= 0.d0) then
           hji = hji * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
           do k=1,N_states
-            tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,2)
+            tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,1)
           enddo
         endif
       end if
@@ -265,7 +265,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hij = hij_cache(putj,1) - hij_cache(putj,2)
           if (hij /= 0.d0) then
             hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
-            tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,1)
+            tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,2)
           endif
         end do
         do putj=hfix+1,mo_num
@@ -274,7 +274,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hij = hij_cache(putj,2) - hij_cache(putj,1)
           if (hij /= 0.d0) then
             hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
-            tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,1)
+            tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,2)
           endif
         end do
 
@@ -293,7 +293,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hji = hji_cache(putj,1) - hji_cache(putj,2)
           if (hji /= 0.d0) then
             hji = hji * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
-            tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,2)
+            tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,1)
           endif
         end do
         do putj=hfix+1,mo_num
@@ -302,7 +302,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hji = hji_cache(putj,2) - hji_cache(putj,1)
           if (hji /= 0.d0) then
             hji = hji * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
-            tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,2)
+            tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,1)
           endif
         end do
 
@@ -342,7 +342,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
             hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
             !DIR$ LOOP COUNT AVG(4)
             do k=1,N_states
-              tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,1)
+              tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,2)
             enddo
           endif
         end if
@@ -353,7 +353,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           if (hij /= 0.d0) then
             hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
             do k=1,N_states
-              tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,1)
+              tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,2)
             enddo
           endif
         end if
@@ -385,7 +385,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
             hji = hji * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
             !DIR$ LOOP COUNT AVG(4)
             do k=1,N_states
-              tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,2)
+              tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,1)
             enddo
           endif
         end if
@@ -396,7 +396,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           if (hji /= 0.d0) then
             hji = hji * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
             do k=1,N_states
-              tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,2)
+              tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,1)
             enddo
           endif
         end if
@@ -445,8 +445,8 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         do k=1,N_states
           ! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
 !          mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * dconjg(hij)
-          mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * hij
-          mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,2) * hji
+          mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,2) * hij
+          mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,1) * hji
         enddo
       end do
     end do

src/cipsi_tc_bi_ortho/get_d2_good.irp.f

@@ -79,12 +79,12 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         if(ma == 1) then ! if particle spins are (alpha,alpha,alpha,beta), then puti is beta and putj is alpha
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,1) * hij
+            mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij
           enddo
         else            ! if particle spins are (beta,beta,beta,alpha), then puti is alpha and putj is beta
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
+            mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
           enddo
         end if
       end do
@@ -103,12 +103,12 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         if(ma == 1) then ! if particle spins are (alpha,alpha,alpha,beta), then puti is beta and putj is alpha
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,2) * hji
+            mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji
           enddo
         else            ! if particle spins are (beta,beta,beta,alpha), then puti is alpha and putj is beta
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
+            mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
           enddo
         end if
       end do
@@ -135,7 +135,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
             hij = hij * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
             !DIR$ LOOP COUNT AVG(4)
             do k=1,N_states
-              mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
+              mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
             enddo
           endif
         end do
@@ -154,7 +154,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
             hji = hji * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
             !DIR$ LOOP COUNT AVG(4)
             do k=1,N_states
-              mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
+              mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
             enddo
           endif
         end do
@@ -189,7 +189,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,1) * hij
+            mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,2) * hij
           enddo
         end do
       end do
@@ -210,7 +210,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hji = hji * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,2) * hji
+            mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,1) * hji
           enddo
         end do
       end do
@@ -239,12 +239,12 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         if (puti < putj) then
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
+            mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
           enddo
         else
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,1) * hij
+            mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij
           enddo
         endif
       end do
@@ -262,12 +262,12 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
         if (puti < putj) then
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
+            mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
           enddo
         else
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,2) * hji
+            mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji
           enddo
         endif
       end do
@@ -290,7 +290,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hij = hij * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
+            mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
           enddo
         end if
       !! <phi|H|alpha>
@@ -299,7 +299,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
           hji = hji * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
           do k=1,N_states
-            mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
+            mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
           enddo
         end if
       end if

src/cipsi_tc_bi_ortho/selection.irp.f

@@ -893,20 +893,45 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
               alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function
             enddo
           else if(debug_tc_pt2 == 2)then !! debugging the new version
+!            psi_h_alpha_tmp = 0.d0
+!            alpha_h_psi_tmp = 0.d0
+!            do iii = 1, N_det_selectors ! old version
+!              call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
+!              call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
+!              psi_h_alpha_tmp += i_h_alpha * psi_selectors_coef_tc(iii,1,1) ! left function
+!              alpha_h_psi_tmp += alpha_h_i * psi_selectors_coef_tc(iii,2,1) ! right function
+!            enddo
             psi_h_alpha_tmp = mat_l(istate, p1, p2) ! new version
             alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version
             psi_h_alpha = 0.d0
             alpha_h_psi = 0.d0
-            do iii = 1, N_det_selectors ! old version
-              call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
-              call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
-              psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function
-              alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function
+            do iii = 1, N_det ! old version
+              call htilde_mu_mat_opt_bi_ortho_no_3e(psi_det(1,1,iii), det, N_int, i_h_alpha)
+              call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_det(1,1,iii), N_int, alpha_h_i)
+              psi_h_alpha += i_h_alpha * psi_l_coef_bi_ortho(iii,1) ! left function
+              alpha_h_psi += alpha_h_i * psi_r_coef_bi_ortho(iii,1) ! right function
             enddo
             if(dabs(psi_h_alpha*alpha_h_psi/delta_E).gt.1.d-10)then
               error = dabs(psi_h_alpha * alpha_h_psi - psi_h_alpha_tmp * alpha_h_psi_tmp)/dabs(psi_h_alpha * alpha_h_psi)
               if(error.gt.1.d-2)then
+                call debug_det(det, N_int)
                 print*,'error =',error,psi_h_alpha * alpha_h_psi/delta_E,psi_h_alpha_tmp * alpha_h_psi_tmp/delta_E
+                print*,psi_h_alpha , alpha_h_psi
+                print*,psi_h_alpha_tmp , alpha_h_psi_tmp
+                print*,'selectors '
+                do iii = 1, N_det_selectors ! old version
+                 print*,'iii',iii,psi_selectors_coef_tc(iii,1,1),psi_selectors_coef_tc(iii,2,1)
+                 call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
+                 call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
+                 print*,i_h_alpha,alpha_h_i
+                 call debug_det(psi_selectors(1,1,iii),N_int) 
+                enddo 
+!                print*,'psi_det '
+!                do iii = 1, N_det! old version
+!                 print*,'iii',iii,psi_l_coef_bi_ortho(iii,1),psi_r_coef_bi_ortho(iii,1)
+!                 call debug_det(psi_det(1,1,iii),N_int) 
+!                enddo 
+                stop
               endif
             endif
           else

src/cisd/EZFIO.cfg

@@ -5,3 +5,11 @@ interface: ezfio
 size: (determinants.n_states)
 
 
+
+[lcc_energy]
+type: double precision
+doc: lccsd energy 
+interface: ezfio
+size: (determinants.n_states)
+
+

src/cisd/NEED

@@ -1,3 +1,4 @@
 selectors_full
 single_ref_method
 davidson_undressed
+dav_general_mat

src/cisd/lccsd.irp.f

@@ -0,0 +1,95 @@
+program lccsd
+  implicit none
+  BEGIN_DOC
+! Linerarized CCSD 
+!
+  ! This program takes a reference Slater determinant of ROHF-like occupancy,
+  !
+  ! and performs all single and double excitations on top of it, disregarding
+  ! spatial symmetry and compute the "n_states" lowest eigenstates of that CI
+  ! matrix (see :option:`determinants n_states`).
+  !
+  ! This program can be useful in many cases:
+  !
+  ! * **Ground state calculation**: if even after a :c:func:`cis` calculation, natural
+  !   orbitals (see :c:func:`save_natorb`) and then :c:func:`scf` optimization, you are not sure to have the lowest scf
+  !   solution,
+  !   do the same strategy with the :c:func:`cisd` executable instead of the :c:func:`cis` exectuable to generate the natural
+  !   orbitals as a guess for the :c:func:`scf`.
+  !
+  !
+  !
+  ! * **Excited states calculations**: the lowest excited states are much likely to
+  !   be dominanted by single- or double-excitations.
+  !   Therefore, running a :c:func:`cisd` will save the "n_states" lowest states within
+  !   the CISD space
+  !   in the |EZFIO| directory, which can afterward be used as guess wave functions
+  !   for a further multi-state fci calculation if you specify "read_wf" = True
+  !   before running the fci executable (see :option:`determinants read_wf`).
+  !   Also, if you specify "s2_eig" = True, the cisd will only retain states
+  !   having the good value :math:`S^2` value
+  !   (see :option:`determinants expected_s2` and :option:`determinants s2_eig`).
+  !   If "s2_eig" = False, it will take the lowest n_states, whatever
+  !   multiplicity they are.
+  !
+  !
+  !
+  !   Note: if you would like to discard some orbitals, use
+  !   :ref:`qp_set_mo_class` to specify:
+  !
+  !   * "core" orbitals which will be always doubly occupied
+  !
+  !   * "act" orbitals where an electron can be either excited from or to
+  !
+  !   * "del" orbitals which will be never occupied
+  !
+  END_DOC
+  PROVIDE N_states
+  read_wf = .False.
+  TOUCH read_wf
+  call run
+end
+
+subroutine run
+  implicit none
+
+  if(pseudo_sym)then
+   call H_apply_cisd_sym
+  else
+   call H_apply_cisd
+  endif
+  call get_lccsd_2
+end
+
+subroutine get_lccsd_2
+  implicit none
+  integer                        :: i,k
+  double precision               :: cisdq(N_states), delta_e
+  double precision,external      :: diag_h_mat_elem
+  psi_coef = lccsd_coef
+  SOFT_TOUCH psi_coef
+  call save_wavefunction_truncated(save_threshold)
+  call ezfio_set_cisd_lcc_energy(lccsd_energies)
+
+  print *,  'N_det = ', N_det
+  print*,''
+  print*,'******************************'
+  print *,  'LCCSD Energies'
+  do i = 1,N_states
+    print *,  i, lccsd_energies(i)
+  enddo
+  if (N_states > 1) then
+      print*,'******************************'
+      print*,'Excitation energies (au)    (LCCSD)'
+      do i = 2, N_states
+        print*, i ,lccsd_energies(i) - lccsd_energies(1)
+      enddo
+      print*,''
+      print*,'******************************'
+      print*,'Excitation energies (eV)    (LCCSD)'
+      do i = 2, N_states
+        print*, i ,(lccsd_energies(i) - lccsd_energies(1)) * ha_to_ev
+      enddo
+  endif
+
+end

src/cisd/lccsd_prov.irp.f

@@ -0,0 +1,46 @@
+ BEGIN_PROVIDER [ double precision, lccsd_coef, (N_det, N_states)]
+&BEGIN_PROVIDER [ double precision, lccsd_energies, (N_states)]
+ implicit none
+ double precision, allocatable :: Dress_jj(:), H_jj(:), u_in(:,:)
+ double precision :: ebefore, eafter, ecorr, thresh
+ integer :: i,it
+ logical :: converged
+ external H_u_0_nstates_openmp
+ allocate(Dress_jj(N_det),H_jj(N_det),u_in(N_det,N_states_diag))
+ thresh = 1.d-6
+ converged = .False.
+ Dress_jj = 0.d0
+ u_in = 0.d0
+ it = 0
+ ! initial guess
+ do i = 1, N_states_diag
+  u_in(i,i) = 1.d0
+ enddo
+ do i = 1,N_det
+  call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,H_jj(i))
+ enddo
+ ebefore = H_jj(1)
+ do while (.not.converged)
+  it += 1
+  print*,'N_det = ',N_det
+  call davidson_general_ext_rout_diag_dressed(u_in,H_jj,Dress_jj,lccsd_energies,& 
+                   N_det,N_states,N_states_diag,converged,H_u_0_nstates_openmp)
+  ecorr = lccsd_energies(1) - H_jj(1) 
+  print*,'---------------------'
+  print*,'it = ',it
+  print*,'ecorr = ',ecorr
+  Dress_jj(1) = 0.d0
+  do i = 2, N_det
+    if(ecorr + H_jj(i) .gt. H_jj(1))then
+     Dress_jj(i) = ecorr
+    endif
+  enddo
+  eafter = lccsd_energies(1)
+  converged = (dabs(eafter - ebefore).lt.thresh)
+  ebefore = eafter
+ enddo
+ do i = 1, N_states
+  lccsd_coef(1:N_det,i) = u_in(1:N_det,i)
+ enddo
+
+END_PROVIDER 

src/fci_tc_bi/selectors.irp.f

@@ -27,6 +27,8 @@ END_PROVIDER
   implicit none
   BEGIN_DOC
   ! Determinants on which we apply <i|H|psi> for perturbation.
+  ! psi_selectors_coef_tc(iii,1,istate) = left  coefficient of the iii determinant
+  ! psi_selectors_coef_tc(iii,2,istate) = right coefficient of the iii determinant
   END_DOC
   integer                        :: i,k
 

src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f

@@ -90,3 +90,96 @@ subroutine htcdag_bi_ortho_calc_tdav_slow(v, u, N_st, sze)
 
 end 
 
+subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,chi_H_i_array,i_H_phi_array)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Computes $\langle i|H|Phi \rangle  = \sum_J c^R_J \langle i | H | J \rangle$.
+!
+! AND      $\langle Chi|H|  i \rangle  = \sum_J c^L_J \langle J | H | i \rangle$.
+!
+! CONVENTION: i_H_phi_array(0) =   total matrix element,  
+!
+!             i_H_phi_array(1) =   one-electron matrix element,
+!
+!             i_H_phi_array(2) =   two-electron matrix element,
+!
+!             i_H_phi_array(3) = three-electron matrix element,
+!
+! Uses filter_connected_i_H_psi0 to get all the $|J \rangle$ to which $|i \rangle$
+! is connected.
+!
+! The i_H_psi_minilist is much faster but requires to build the
+! minilists.
+  END_DOC
+  integer, intent(in)            :: Nint, Ndet,Ndet_max,Nstate
+  integer(bit_kind), intent(in)  :: keys(Nint,2,Ndet)
+  integer(bit_kind), intent(in)  :: key(Nint,2)
+  double precision, intent(in)   :: coef_l(Ndet_max,Nstate),coef_r(Ndet_max,Nstate)
+  double precision, intent(out)  :: chi_H_i_array(0:3,Nstate),i_H_phi_array(0:3,Nstate)
+
+  integer                        :: i, ii,j
+  double precision               :: phase
+  integer                        :: exc(0:2,2,2)
+  double precision               :: hmono, htwoe, hthree, htot
+  integer, allocatable           :: idx(:)
+
+  ASSERT (Nint > 0)
+  ASSERT (N_int == Nint)
+  ASSERT (Nstate > 0)
+  ASSERT (Ndet > 0)
+  ASSERT (Ndet_max >= Ndet)
+  allocate(idx(0:Ndet))
+
+  chi_H_i_array = 0.d0
+  i_H_phi_array = 0.d0
+
+  call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
+  if (Nstate == 1) then
+
+    do ii=1,idx(0)
+      i = idx(ii)
+      ! computes <Chi|H_tc|i>
+      !DIR$ FORCEINLINE
+      call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot)
+      chi_H_i_array(0,1) = chi_H_i_array(0,1) + coef_l(i,1)*htot
+      chi_H_i_array(1,1) = chi_H_i_array(1,1) + coef_l(i,1)*hmono
+      chi_H_i_array(2,1) = chi_H_i_array(2,1) + coef_l(i,1)*htwoe
+      chi_H_i_array(3,1) = chi_H_i_array(3,1) + coef_l(i,1)*hthree
+      ! computes <i|H_tc|Phi>
+      !DIR$ FORCEINLINE
+      call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot)
+      i_H_phi_array(0,1) = i_H_phi_array(0,1) + coef_r(i,1)*htot
+      i_H_phi_array(1,1) = i_H_phi_array(1,1) + coef_r(i,1)*hmono
+      i_H_phi_array(2,1) = i_H_phi_array(2,1) + coef_r(i,1)*htwoe
+      i_H_phi_array(3,1) = i_H_phi_array(3,1) + coef_r(i,1)*hthree
+    enddo
+
+  else
+
+    do ii=1,idx(0)
+      i = idx(ii)
+      ! computes <Chi|H_tc|i>
+      !DIR$ FORCEINLINE
+      call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot)
+      do j = 1, Nstate
+        chi_H_i_array(0,j) = chi_H_i_array(0,j) + coef_l(i,j)*htot
+        chi_H_i_array(1,j) = chi_H_i_array(1,j) + coef_l(i,j)*hmono
+        chi_H_i_array(2,j) = chi_H_i_array(2,j) + coef_l(i,j)*htwoe
+        chi_H_i_array(3,j) = chi_H_i_array(3,j) + coef_l(i,j)*hthree
+      enddo
+      ! computes <i|H_tc|Phi>
+      !DIR$ FORCEINLINE
+      call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot)
+      do j = 1, Nstate
+        i_H_phi_array(0,j) = i_H_phi_array(0,j) + coef_r(i,j)*htot
+        i_H_phi_array(1,j) = i_H_phi_array(1,j) + coef_r(i,j)*hmono
+        i_H_phi_array(2,j) = i_H_phi_array(2,j) + coef_r(i,j)*htwoe
+        i_H_phi_array(3,j) = i_H_phi_array(3,j) + coef_r(i,j)*hthree
+      enddo
+    enddo
+
+  endif
+
+end
+

src/tc_bi_ortho/normal_ordered_old.irp.f

@@ -120,6 +120,13 @@ END_PROVIDER
 subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
 
   use bitmasks ! you need to include the bitmasks_module.f90 features
+  BEGIN_DOC
+!  give the contribution for a double excitation of opposite spin BUT averaged over spin
+!
+! it is the average of <p1_down p2_up |h1_down h2_up> and <p1_up p2_down |h1_up h2_down> 
+!
+! because the orbitals h1,h2,p1,p2 are spatial orbitals and therefore can be of different spins
+  END_DOC
 
   implicit none
   integer, intent(in)           :: Nint, h1, h2, p1, p2
@@ -158,7 +165,8 @@ subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
     call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
     int_exc_12 = -1.d0 * integral
 
-    hthree += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12)
+    hthree += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12) ! spin average
+!    hthree += 1.d0 * int_direct - 1.0d0 * (int_exc_13 + int_exc_12)
   enddo
 
   return

src/tc_bi_ortho/slater_tc_3e_slow.irp.f

@@ -184,7 +184,7 @@ subroutine single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
       ii = occ(i,s1)
       do j = i+1, Ne(s1)
        jj = occ(j,s1)
-!       ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1)
+!      !ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1)
        hthree += three_e_single_parrallel_spin(jj,ii,p1,h1) ! USES THE 4-IDX TENSOR
       enddo
      enddo

src/tc_bi_ortho/tc_h_eigvectors.irp.f

@@ -328,6 +328,11 @@ end
   TOUCH psi_r_coef_bi_ortho
   call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(buffer)
   deallocate(buffer)
+!  print*,'After diag'
+!  do i = 1, N_det! old version
+!   print*,'i',i,psi_l_coef_bi_ortho(i,1),psi_r_coef_bi_ortho(i,1)
+!   call debug_det(psi_det(1,1,i),N_int) 
+!  enddo 
 
 END_PROVIDER 
 

src/tc_bi_ortho/tc_natorb.irp.f

@@ -23,7 +23,7 @@
 
   dm_tmp(1:mo_num,1:mo_num) = -tc_transition_matrix_mo(1:mo_num,1:mo_num,1,1)
 
-  print *, ' dm_tmp'
+  print *, ' Transition density matrix '
   do i = 1, mo_num
     fock_diag(i) = fock_matrix_tc_mo_tot(i,i)
     write(*, '(100(F16.10,X))') -dm_tmp(:,i)
@@ -32,8 +32,15 @@
   thr_d   = 1.d-6
   thr_nd  = 1.d-6
   thr_deg = 1.d-3
-  call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & 
-                              , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
+  if(n_core_orb.ne.0)then
+!   print*,'core orbitals'
+!   pause
+   call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg & 
+                                    , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
+  else
+   call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & 
+                               , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
+  endif
 !   call non_hrmt_bieig( mo_num, dm_tmp&
 !                      , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& 
 !                      , mo_num, natorb_tc_eigval )

src/tc_bi_ortho/tc_prop.irp.f

@@ -29,7 +29,7 @@
          tc_transition_matrix_mo_alpha(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
         enddo
         do p = 1, n_occ_ab(2) ! browsing the beta electrons
-         m = occ(p,1)
+         m = occ(p,2)
          tc_transition_matrix_mo_beta(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
         enddo
        else
@@ -38,12 +38,14 @@
           ! Single alpha
           h = exc(1,1,1) ! hole in psi_det(1,1,j) 
           p = exc(1,2,1) ! particle in psi_det(1,1,j) 
-          tc_transition_matrix_mo_alpha(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
+          tc_transition_matrix_mo_alpha(p,h,istate,jstate)+= &
+          phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
         else
           ! Single beta
           h = exc(1,1,2) ! hole in psi_det(1,1,j) 
           p = exc(1,2,2) ! particle in psi_det(1,1,j) 
-          tc_transition_matrix_mo_beta(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
+          tc_transition_matrix_mo_beta(p,h,istate,jstate)+=  &
+          phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
         endif
        endif
      enddo

src/tc_bi_ortho/test_normal_order.irp.f

@@ -33,7 +33,7 @@ subroutine test
  integer :: h1,h2,p1,p2,s1,s2,i_ok,degree,Ne(2)
  integer                        :: exc(0:2,2,2)
  integer(bit_kind), allocatable :: det_i(:,:)
- double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal
+ double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal,hthree_tmp
  integer,           allocatable :: occ(:,:)
  allocate( occ(N_int*bit_kind_size,2) )
  call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
@@ -45,15 +45,44 @@ subroutine test
   do p1 = elec_alpha_num+1, mo_num
    do h2 = 1, elec_beta_num
     do p2 = elec_beta_num+1, mo_num
+     hthree = 0.d0
+
      det_i = ref_bitmask
+ s1 = 1 
+ s2 = 2
      call do_single_excitation(det_i,h1,p1,s1,i_ok)
+    if(i_ok.ne.1)cycle
      call do_single_excitation(det_i,h2,p2,s2,i_ok)
-     call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+    if(i_ok.ne.1)cycle
+     call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,htilde_ij)
      call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
      call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
-     hthree *= phase
-!    !normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1)
-     call three_comp_two_e_elem(det_i,h1,h2,p1,p2,s1,s2,normal)
+     hthree_tmp *= phase
+     hthree += 0.5d0 * hthree_tmp
+     det_i = ref_bitmask
+ s1 = 2 
+ s2 = 1
+     call do_single_excitation(det_i,h1,p1,s1,i_ok)
+    if(i_ok.ne.1)cycle
+     call do_single_excitation(det_i,h2,p2,s2,i_ok)
+    if(i_ok.ne.1)cycle
+     call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,htilde_ij)
+     call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
+     call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
+     hthree_tmp *= phase
+     hthree += 0.5d0 * hthree_tmp
+
+
+!     normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1)
+     call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, normal)
+     if(dabs(hthree).lt.1.d-10)cycle
+     if(dabs(hthree-normal).gt.1.d-10)then
+!      print*,pp2,pp1,hh2,hh1
+      print*,p2,p1,h2,h1
+      print*,hthree,normal,dabs(hthree-normal)
+      stop
+     endif
+!     call three_comp_two_e_elem(det_i,h1,h2,p1,p2,s1,s2,normal)
 !     normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1)
      accu += dabs(hthree-normal)
     enddo
@@ -86,8 +115,8 @@ do h1 = 1, elec_alpha_num
     integer :: hh1, pp1, hh2, pp2, ss1, ss2
     call decode_exc(exc, 2, hh1, pp1, hh2, pp2, ss1, ss2)
     hthree *= phase
-!    normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
-     normal = eff_2_e_from_3_e_aa(p2,p1,h2,h1)
+    normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
+!     normal = eff_2_e_from_3_e_aa(p2,p1,h2,h1)
     if(dabs(hthree).lt.1.d-10)cycle
     if(dabs(hthree-normal).gt.1.d-10)then
      print*,pp2,pp1,hh2,hh1

src/tc_bi_ortho/test_tc_fock.irp.f

@@ -152,9 +152,7 @@ subroutine routine_tot()
 !   do i = 1, elec_num_tab(s1)
 !     do a = elec_num_tab(s1)+1, mo_num ! virtual 
    do i = 1, elec_beta_num
-     do a = elec_beta_num+1, elec_alpha_num! virtual 
-!   do i = elec_beta_num+1, elec_alpha_num
-!     do a = elec_alpha_num+1, mo_num! virtual 
+     do a = elec_beta_num+1, mo_num! virtual 
        print*,i,a 
  
        det_i = ref_bitmask
@@ -167,7 +165,7 @@ subroutine routine_tot()
  
        call htilde_mu_mat_bi_ortho_slow(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
        print*,htilde_ij
-       if(dabs(htilde_ij).lt.1.d-10)cycle
+!       if(dabs(htilde_ij).lt.1.d-10)cycle
        print*, ' excited det'
        call debug_det(det_i, N_int)
 
@@ -184,9 +182,12 @@ subroutine routine_tot()
 !       endif
        err_ai = dabs(dabs(ref) - dabs(new))
        if(err_ai .gt. 1d-7) then
+         print*,'---------'
          print*,'s1 = ',s1
          print*, ' warning on', i, a
          print*, ref,new,err_ai
+         print*,hmono, htwoe, hthree
+         print*,'---------'
        endif
        print*, ref,new,err_ai
        err_tot += err_ai

src/tc_bi_ortho/test_tc_two_rdm.irp.f

@@ -0,0 +1,73 @@
+program test_tc_rdm
+
+  BEGIN_DOC
+  !
+  ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together 
+  !        with the energy. Saves the left-right wave functions at the end. 
+  !
+  END_DOC
+
+  my_grid_becke = .True.
+  PROVIDE tc_grid1_a tc_grid1_r
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  read_wf = .True.
+  touch read_wf
+
+  print*, ' nb of states = ', N_states
+  print*, ' nb of det    = ', N_det
+
+  call test()
+
+end
+
+subroutine  test
+ implicit none
+ integer :: h1,p1,h2,p2,i,j,istate
+ double precision :: rdm, integral, accu,ref
+ double precision :: hmono, htwoe, hthree, htot
+ accu = 0.d0
+ do h1 = 1, mo_num
+  do p1 = 1, mo_num
+   do h2 = 1, mo_num
+    do p2 = 1, mo_num
+     integral = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
+     rdm = tc_two_rdm(p1,h1,p2,h2)
+!     if(dabs(rdm).gt.1.d-10)then
+!      print*,h1,p1,h2,p2
+!      print*,rdm,integral,rdm*integral
+!     endif
+     accu += integral * rdm
+    enddo
+   enddo
+  enddo
+ enddo
+ accu *= 0.5d0
+ print*,'accu = ',accu
+! print*,mo_bi_ortho_tc_two_e(2,15,2,1)
+! print*,mo_bi_ortho_tc_two_e(15,2,2,1)
+! print*,mo_bi_ortho_tc_two_e(2,1,2,15)
+! print*,mo_bi_ortho_tc_two_e(2,1,15,2)
+ ref = 0.d0
+ do i = 1, N_det
+  do j = 1, N_det
+!   if(i.eq.j)cycle
+   call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
+   do istate = 1,N_states
+!  print*,'i,j',i,j
+!  print*,psi_l_coef_bi_ortho(i,istate) , psi_r_coef_bi_ortho(j,istate) , htwoe 
+!  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe 
+!   if(i.ne.j)then
+!  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) ,  htwoe 
+!  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe 
+!   endif
+    ref += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe 
+   enddo
+  enddo
+ enddo
+ print*,' ref = ',ref 
+ print*,'delta= ',ref-accu
+
+end

src/tc_bi_ortho/two_rdm_naive.irp.f

@@ -0,0 +1,124 @@
+BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! tc_two_rdm(p,s,q,r) = <Phi| a^dager_p
+ END_DOC
+ integer :: i,j,istate,m,mm,nn
+ integer                        :: exc(0:2,2,2)
+ double precision               :: phase
+ double precision               :: contrib
+ integer                        :: h1,p1,s1,h2,p2,s2,degree
+ integer, allocatable           :: occ(:,:)
+ integer                        :: n_occ_ab(2),other_spin(2)
+ other_spin(1) = 2
+ other_spin(2) = 1
+ allocate(occ(N_int*bit_kind_size,2))
+ tc_two_rdm = 0.d0
+
+ do i = 1, N_det ! psi_left 
+  do j = 1, N_det ! psi_right 
+   call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+   if(degree.gt.2)cycle
+   if(degree.gt.0)then
+    ! get excitation operators: from psi_det(j) --> psi_det(i)
+    call get_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,degree,phase,N_int)
+    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+    contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * phase * state_average_weight(1)
+    do istate = 2, N_states
+     contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * phase * state_average_weight(istate)
+    enddo
+    if(degree == 2)then
+     call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+    else if(degree==1)then
+     ! occupation of the determinant psi_det(j)
+     call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) 
+ 
+     ! run over the electrons of opposite spin than the excitation
+     s2 = other_spin(s1)
+     do mm = 1, n_occ_ab(s2) 
+      m = occ(mm,s2)
+      h2 = m 
+      p2 = m 
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+     enddo
+     ! run over the electrons of same spin than the excitation
+     s2 = s1
+     do mm = 1, n_occ_ab(s2) 
+      m = occ(mm,s2)
+      h2 = m 
+      p2 = m 
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+     enddo
+    endif
+   else if(degree == 0)then
+!    cycle
+    contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) *  state_average_weight(1)
+    do istate = 2, N_states
+     contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) *  state_average_weight(istate)
+    enddo
+    ! occupation of the determinant psi_det(j)
+    call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) 
+    s1 = 1 ! alpha electrons
+     do nn = 1, n_occ_ab(s1)
+      h1 = occ(nn,s1)
+      p1 = occ(nn,s1)
+      ! run over the couple of alpha-beta electrons 
+      s2 = other_spin(s1)
+      do mm = 1, n_occ_ab(s2) 
+       m = occ(mm,s2)
+       h2 = m 
+       p2 = m 
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      enddo
+      ! run over the couple of alpha-alpha electrons 
+      s2 = s1
+      do mm = 1, n_occ_ab(s2) 
+       m = occ(mm,s2)
+       h2 = m 
+       p2 = m 
+       if(h2.le.h1)cycle
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      enddo
+    enddo
+    s1 = 2
+     do nn = 1, n_occ_ab(s1)
+      h1 = occ(nn,s1)
+      p1 = occ(nn,s1)
+      ! run over the couple of beta-beta electrons 
+      s2 = s1
+      do mm = 1, n_occ_ab(s2) 
+       m = occ(mm,s2)
+       h2 = m 
+       p2 = m 
+       if(h2.le.h1)cycle
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      enddo
+    enddo
+   endif
+  enddo
+ enddo
+
+END_PROVIDER 
+
+subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,contrib)
+ implicit none
+ integer, intent(in) :: h1,p1,h2,p2,s1,s2,sze
+ double precision, intent(in) :: contrib
+ double precision, intent(inout) :: array(sze, sze, sze, sze)
+ integer :: istate
+ if(s1.ne.s2)then
+   array(p1,h1,p2,h2) += contrib
+   ! permutation for particle symmetry
+   array(p2,h2,p1,h1) += contrib
+ else ! same spin double excitation 
+   array(p1,h1,p2,h2) += contrib
+   ! exchange 
+   ! exchanging the particles 
+   array(p2,h1,p1,h2) -= contrib
+   ! exchanging the 
+   array(p1,h2,p2,h1) -= contrib
+   ! permutation for particle symmetry
+   array(p2,h2,p1,h1) += contrib
+ endif
+
+end

src/tc_keywords/EZFIO.cfg

@@ -214,6 +214,12 @@ doc: Threshold to determine if diagonal elements of the bi-orthogonal condition
 interface: ezfio,provider,ocaml
 default: 1.e-6
 
+[thresh_lr_angle]
+type: double precision
+doc: Maximum value of the angle between the couple of left and right orbital for the rotations
+interface: ezfio,provider,ocaml
+default: 20.0
+
 [thresh_biorthog_nondiag]
 type: Threshold
 doc: Threshold to determine if non-diagonal elements of L.T x R are close enouph to 0

src/tc_scf/fock_tc.irp.f

@@ -208,10 +208,10 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
 
     if(three_body_h_tc) then
       !call wall_time(tt0)
-      !PROVIDE fock_a_tot_3e_bi_orth
-      !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth
-      PROVIDE fock_3e_uhf_mo_a
-      Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a
+      PROVIDE fock_a_tot_3e_bi_orth
+      Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth
+!      PROVIDE fock_3e_uhf_mo_a
+!      Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a
       !call wall_time(tt1)
       !print*, ' 3-e term:', tt1-tt0
     endif
@@ -241,21 +241,13 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
 
   if(bi_ortho) then
 
-    !allocate(tmp(ao_num,ao_num))
-    !tmp = Fock_matrix_tc_ao_beta
-    !if(three_body_h_tc) then
-    !  tmp += fock_3e_uhf_ao_b
-    !endif
-    !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1))
-    !deallocate(tmp)
-
     call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
                           , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
     if(three_body_h_tc) then
-      !PROVIDE fock_b_tot_3e_bi_orth
-      !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth
-      PROVIDE fock_3e_uhf_mo_b
-      Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b
+       PROVIDE fock_b_tot_3e_bi_orth
+       Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth
+!      PROVIDE fock_3e_uhf_mo_b
+!      Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b
     endif
 
   else

src/tc_scf/routines_rotates.irp.f

@@ -140,7 +140,11 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
  ! compute the overlap between the left and rescaled right
   call build_s_matrix(ao_num, mo_num, mo_r_coef_new, mo_r_coef_new, ao_overlap, s_mat)
 ! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
-  call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list)
+  if(n_core_orb.ne.0)then
+   call give_degen_full_listcore(fock_diag, mo_num, list_core, n_core_orb, thr_deg, list_degen, n_degen_list)
+  else
+   call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list)
+  endif
   print *, ' fock_matrix_mo'
   do i = 1, mo_num
     print *, i, fock_diag(i), angle_left_right(i)
@@ -152,6 +156,8 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
 !  n_degen = ilast - ifirst +1
 
     n_degen = list_degen(i,0)
+    if(n_degen .ge. 1000)n_degen = 1 ! convention for core orbitals
+     
     if(n_degen .eq. 1) cycle
 
     allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen))
@@ -279,7 +285,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
   allocate(new_angles(mo_num))
   new_angles(1:mo_num) = dabs(angle_left_right(1:mo_num))
   max_angle = maxval(new_angles)
-  good_angles = max_angle.lt.45.d0
+  good_angles = max_angle.lt.thresh_lr_angle
   print *, ' max_angle = ', max_angle
   deallocate(new_angles)
 
@@ -397,11 +403,11 @@ subroutine print_energy_and_mos(good_angles)
   print *, ' TC SCF energy gradient = ', grad_non_hermit
   print *, ' Max angle Left/right   = ', max_angle_left_right
 
-  if(max_angle_left_right .lt. 45.d0) then
+  if(max_angle_left_right .lt. thresh_lr_angle) then
     print *, ' Maximum angle BELOW 45 degrees, everthing is OK !'
     good_angles = .true.
-  else if(max_angle_left_right .gt. 45.d0 .and. max_angle_left_right .lt. 75.d0) then
-    print *, ' Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...'
+  else if(max_angle_left_right .gt. thresh_lr_angle .and. max_angle_left_right .lt. 75.d0) then
+    print *, ' Maximum angle between thresh_lr_angle and 75 degrees, this is not the best for TC-CI calculations ...'
     good_angles = .false.
   else if(max_angle_left_right .gt. 75.d0) then
     print *, ' Maximum angle between ABOVE 75 degrees, YOU WILL CERTAINLY FIND TROUBLES IN TC-CI calculations ...'

src/utils/block_diag_degen_core.irp.f

@@ -0,0 +1,244 @@
+
+subroutine diag_mat_per_fock_degen_core(fock_diag, mat_ref, listcore,ncore, n, thr_d, thr_nd, thr_deg, leigvec, reigvec, eigval)
+
+
+  BEGIN_DOC
+  !
+  ! subroutine that diagonalizes a matrix mat_ref BY BLOCK
+  !
+  ! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag"
+  !
+  ! the elements of listcore are untouched 
+  !
+  ! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together
+  !
+  !          : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together
+  !
+  !          : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together
+  !
+  ! etc... the advantage is to guarentee no spurious mixing because of numerical problems. 
+  !
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: n,ncore, listcore(ncore)
+  double precision, intent(in)  :: fock_diag(n), mat_ref(n,n), thr_d, thr_nd, thr_deg
+  double precision, intent(out) :: leigvec(n,n), reigvec(n,n), eigval(n)
+
+  integer                       :: n_degen_list, n_degen,size_mat, i, j, k, icount, m, index_degen
+  integer                       :: ii, jj, i_good, j_good, n_real
+  integer                       :: icount_eigval
+  logical,          allocatable :: is_ok(:)
+  integer,          allocatable :: list_degen(:,:), list_same_degen(:)
+  integer,          allocatable :: iorder(:), list_degen_sorted(:)
+  double precision, allocatable :: leigvec_unsrtd(:,:), reigvec_unsrtd(:,:), eigval_unsrtd(:)
+  double precision, allocatable :: mat_tmp(:,:), eigval_tmp(:), leigvec_tmp(:,:), reigvec_tmp(:,:)
+
+  allocate(leigvec_unsrtd(n,n), reigvec_unsrtd(n,n), eigval_unsrtd(n))
+  leigvec_unsrtd = 0.d0
+  reigvec_unsrtd = 0.d0
+  eigval_unsrtd  = 0.d0
+
+  ! obtain degeneracies 
+  allocate(list_degen(n,0:n))
+  call give_degen_full_listcore(fock_diag, n, listcore, ncore, thr_deg, list_degen, n_degen_list)
+
+  allocate(iorder(n_degen_list), list_degen_sorted(n_degen_list))
+  do i = 1, n_degen_list
+    n_degen = list_degen(i,0) 
+    list_degen_sorted(i) = n_degen
+    iorder(i) = i
+  enddo
+
+  ! sort by number of degeneracies 
+  call isort(list_degen_sorted, iorder, n_degen_list)
+
+  allocate(is_ok(n_degen_list))
+  is_ok = .True.
+  icount_eigval = 0
+
+  ! loop over degeneracies 
+  do i = 1, n_degen_list
+    if(.not.is_ok(i)) cycle
+
+    is_ok(i) = .False.
+    n_degen  = list_degen_sorted(i)
+
+    
+    if(n_degen.ge.1000)then
+     print*,'core orbital '
+    else
+     print *, ' diagonalizing for n_degen = ', n_degen
+    endif
+
+    k = 1
+
+   ! group all the entries having the same degeneracies 
+!!  do while (list_degen_sorted(i+k)==n_degen)
+    do m = i+1, n_degen_list
+      if(list_degen_sorted(m)==n_degen) then
+        is_ok(i+k) = .False.
+        k += 1
+      endif
+    enddo
+
+    print *, ' number of identical degeneracies = ', k
+    if(n_degen.ge.1000)then
+     n_degen = 1
+    endif
+    size_mat = k*n_degen  
+    print *, ' size_mat = ', size_mat
+    allocate(mat_tmp(size_mat,size_mat), list_same_degen(size_mat))
+    allocate(eigval_tmp(size_mat), leigvec_tmp(size_mat,size_mat), reigvec_tmp(size_mat,size_mat))
+    ! group all the elements sharing the same degeneracy
+    icount = 0
+    do j = 1, k ! jth set of degeneracy
+      index_degen = iorder(i+j-1)
+      do m = 1, n_degen
+        icount += 1
+        list_same_degen(icount) = list_degen(index_degen,m)
+      enddo
+    enddo
+
+    print *, ' list of elements '
+    do icount = 1, size_mat
+      print *, icount, list_same_degen(icount)
+    enddo
+
+    ! you copy subset of matrix elements having all the same degeneracy in mat_tmp
+    do ii = 1, size_mat
+      i_good = list_same_degen(ii)
+      do jj = 1, size_mat
+        j_good = list_same_degen(jj)
+        mat_tmp(jj,ii) = mat_ref(j_good,i_good)
+      enddo
+    enddo
+
+    call non_hrmt_bieig( size_mat, mat_tmp, thr_d, thr_nd &
+                       , leigvec_tmp, reigvec_tmp         & 
+                       , n_real, eigval_tmp )
+
+    do ii = 1, size_mat
+      icount_eigval += 1
+      eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues 
+      do jj = 1, size_mat ! copy the eigenvectors 
+        j_good = list_same_degen(jj)
+        leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii)
+        reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii)
+      enddo
+    enddo
+
+    deallocate(mat_tmp, list_same_degen)
+    deallocate(eigval_tmp, leigvec_tmp, reigvec_tmp)
+  enddo
+
+  if(icount_eigval .ne. n) then
+    print *, ' pb !! (icount_eigval.ne.n)'
+    print *, ' icount_eigval,n', icount_eigval, n
+    stop
+  endif
+ 
+  deallocate(iorder)
+  allocate(iorder(n))
+  do i = 1, n
+    iorder(i) = i
+  enddo
+  call dsort(eigval_unsrtd, iorder, n)
+
+  do i = 1, n
+    print*,'sorted eigenvalues '
+    i_good = iorder(i)
+    eigval(i) = eigval_unsrtd(i)
+    print*,'i,eigval(i) = ',i,eigval(i)
+    do j = 1, n
+      leigvec(j,i) = leigvec_unsrtd(j,i_good)
+      reigvec(j,i) = reigvec_unsrtd(j,i_good)
+    enddo
+  enddo
+
+  deallocate(leigvec_unsrtd, reigvec_unsrtd, eigval_unsrtd)
+  deallocate(list_degen)
+  deallocate(iorder, list_degen_sorted)
+  deallocate(is_ok)
+
+end
+
+! ---
+
+subroutine give_degen_full_listcore(A, n, listcore, ncore, thr, list_degen, n_degen_list)
+
+  BEGIN_DOC
+  ! you enter with an array A(n) and spits out all the elements degenerated up to thr
+  !
+  ! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL 
+  !
+  ! list_degen(i,0) = number of degenerate entries 
+  !
+  ! list_degen(i,1) = index of the first degenerate entry
+  !
+  ! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries 
+  !
+  ! if list_degen(i,0) == 1 it means that there is no degeneracy for that element
+  !
+  ! if list_degen(i,0) >= 1000 it means that it is core orbitals
+  END_DOC
+
+  implicit none
+
+  double precision, intent(in)  :: A(n)
+  double precision, intent(in)  :: thr
+  integer,          intent(in)  :: n,ncore, listcore(ncore)
+  integer,          intent(out) :: list_degen(n,0:n), n_degen_list
+  integer                       :: i, j, icount, icheck,k
+  logical, allocatable          :: is_ok(:)
+
+
+  allocate(is_ok(n))
+  n_degen_list = 0
+  is_ok = .True.
+  ! you first exclude the "core" orbitals  
+  do i = 1, ncore
+   j=listcore(i)
+   is_ok(j) = .False.
+  enddo
+  do i = 1, n
+    if(.not.is_ok(i)) cycle
+    n_degen_list +=1
+    is_ok(i) = .False.
+    list_degen(n_degen_list,1) = i
+    icount = 1
+    do j = i+1, n
+      if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then
+        is_ok(j) = .False.
+        icount += 1
+        list_degen(n_degen_list,icount) = j
+      endif
+    enddo
+
+    list_degen(n_degen_list,0) = icount
+  enddo
+  ! you set all the core orbitals as separate entities 
+  icheck = 0
+  do i = 1, n_degen_list
+    icheck += list_degen(i,0)
+  enddo
+  if(icheck.ne.(n-ncore))then
+    print *, ' pb ! :: icheck.ne.n-ncore'
+    print *, icheck, n-ncore
+    stop
+  endif
+  k=1000
+  do i = 1, ncore
+   n_degen_list+= 1
+   j=listcore(i)
+   list_degen(n_degen_list,1) = i
+   list_degen(n_degen_list,0) = k
+   k+=1
+  enddo
+
+
+
+end
+
+! ---
+

src/utils/fortran_mmap.c

@@ -9,7 +9,6 @@
 
 void* mmap_fortran(char* filename, size_t bytes, int* file_descr, int read_only)
 {
-    int i;
     int fd;
     int result;
     void* map;

src/utils/linear_algebra.irp.f

@@ -1661,7 +1661,15 @@ subroutine restore_symmetry(m,n,A,LDA,thresh)
     ! Update i
     i = i + 1
   enddo
-  copy(i:) = 0.d0
+
+  ! To nullify the remaining elements that are below the threshold
+  if (i == sze) then
+    if (-copy(i) <= thresh) then
+      copy(i) = 0d0
+    endif
+  else
+    copy(i:) = 0.d0
+  endif
 
   !$OMP PARALLEL if (sze>10000) &
   !$OMP SHARED(m,sze,copy_sign,copy,key,A,ii,jj) &