added casscf_cipsi

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_save
+  integer                        :: iteration
+  converged = .False.
+
+  energy = 0.d0
+  mo_label = "MCSCF"
+  iteration = 1
+  state_following_casscf_save = state_following_casscf
+  state_following_casscf = .True.
+  touch state_following_casscf
+  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_energy_pt2(E_PT2)
+     call ezfio_get_casscf_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_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_energy(E(1:N_states))
+  call ezfio_set_casscf_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
+ 
+