Merge branch 'casscf' into dev-lct

INSTALL.rst

@@ -143,7 +143,7 @@ IRPF90
 to Parameters (IRP) method. 
 
 * Download the latest version of IRPF90
-  here : `<https://github.com/scemama/irpf90/releases/latest>`_ and move
+  here : `<https://gitlab.com/scemama/irpf90/-/archive/v1.7.2/irpf90-v1.7.2.tar.gz>`_ and move
   the downloaded archive in the :file:`${QP_ROOT}/external` directory
 
 * Extract the archive and go into the :file:`irpf90-*` directory to run

configure

@@ -297,12 +297,13 @@ EOF
 ${QP_ROOT}/bin
 
 
+
 EOF
 
                 rm ${QP_ROOT}/external/opam_installer.sh
                 source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true
 
-                ${QP_ROOT}/bin/opam init --verbose --yes 
+                ${QP_ROOT}/bin/opam init --verbose --yes --comp=4.07.1 --disable-sandboxing
 
                 eval $(${QP_ROOT}/bin/opam env)
                 opam install -y ${OCAML_PACKAGES}  || exit 1
@@ -310,13 +311,14 @@ EOF
                 # Conventional commands
                 execute << EOF
                   chmod +x "\${QP_ROOT}"/external/opam_installer.sh
+                  "\${QP_ROOT}"/external/opam_installer.sh --no-backup 
                   rm --force \${QP_ROOT}/bin/opam
                   export OPAMROOT=\${OPAMROOT:-\${QP_ROOT}/external/opam}
                   echo \${QP_ROOT}/bin \
                   | sh \${QP_ROOT}/external/opam_installer.sh 
                   rm \${QP_ROOT}/external/opam_installer.sh
                   source \${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true
-                  \${QP_ROOT}/bin/opam init --verbose --yes
+                  \${QP_ROOT}/bin/opam init --verbose --yes --comp=4.07.1 --disable-sandboxing
                   eval \$(\${QP_ROOT}/bin/opam env)
                   opam install -y \${OCAML_PACKAGES}  || exit 1
 EOF

src/bitmask/bitmasks.irp.f

@@ -443,15 +443,6 @@ BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ]
   
 END_PROVIDER
 
- BEGIN_PROVIDER [ integer, n_core_inact_orb ]
- implicit none
- integer :: i
- n_core_inact_orb = 0
- do i = 1, N_int
-  n_core_inact_orb += popcnt(reunion_of_core_inact_bitmask(i,1))
- enddo
- ENd_PROVIDER
-
 BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)]
   implicit none
   BEGIN_DOC
@@ -465,6 +456,20 @@ END_PROVIDER
 END_PROVIDER
 
 
+BEGIN_PROVIDER [integer(bit_kind), reunion_of_inact_act_bitmask, (N_int,2)]
+  implicit none
+  BEGIN_DOC
+  ! Reunion of the  inactive and active bitmasks
+  END_DOC
+  integer                        :: i,j
+  
+  do i = 1, N_int
+    reunion_of_inact_act_bitmask(i,1) = ior(inact_bitmask(i,1),act_bitmask(i,1))
+    reunion_of_inact_act_bitmask(i,2) = ior(inact_bitmask(i,2),act_bitmask(i,2))
+  enddo
+END_PROVIDER
+
+
 BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)]
   implicit none
   BEGIN_DOC
@@ -550,19 +555,3 @@ END_PROVIDER
   enddo
 END_PROVIDER
 
-
- BEGIN_PROVIDER [integer, n_core_orb_allocate]
- implicit none
- n_core_orb_allocate = max(n_core_orb,1)
- END_PROVIDER
-
- BEGIN_PROVIDER [integer, n_inact_orb_allocate]
- implicit none
- n_inact_orb_allocate = max(n_inact_orb,1)
- END_PROVIDER
-
- BEGIN_PROVIDER [integer, n_virt_orb_allocate]
- implicit none
- n_virt_orb_allocate = max(n_virt_orb,1)
- END_PROVIDER
-

src/bitmask/core_inact_act_virt.irp.f

@@ -1,250 +1,382 @@
 use bitmasks
 
-
 BEGIN_PROVIDER [ integer, n_core_orb]
- &BEGIN_PROVIDER [ integer, n_inact_orb ]
- &BEGIN_PROVIDER [ integer, n_act_orb]
- &BEGIN_PROVIDER [ integer, n_virt_orb ]
- &BEGIN_PROVIDER [ integer, n_del_orb ]
   implicit none
   BEGIN_DOC
- ! inact_bitmask : Bitmask of the inactive orbitals which are supposed to be doubly excited 
+  ! Number of core MOs
- ! in post CAS methods
- ! n_inact_orb   : Number of inactive orbitals
- ! virt_bitmask  : Bitmaks of vritual orbitals which are supposed to be recieve electrons 
- ! in post CAS methods
- ! n_virt_orb    : Number of virtual orbitals
- ! list_inact : List of the inactive orbitals which are supposed to be doubly excited 
- ! in post CAS methods
- ! list_virt  : List of vritual orbitals which are supposed to be recieve electrons 
- ! in post CAS methods
- ! list_inact_reverse : reverse list of inactive orbitals 
- ! list_inact_reverse(i) = 0 ::> not an inactive 
- ! list_inact_reverse(i) = k ::> IS the kth inactive 
- ! list_virt_reverse : reverse list of virtual orbitals 
- ! list_virt_reverse(i) = 0 ::> not an virtual 
- ! list_virt_reverse(i) = k ::> IS the kth virtual 
- ! list_act(i) = index of the ith active orbital
- ! 
- ! list_act_reverse : reverse list of active orbitals 
- ! list_act_reverse(i) = 0 ::> not an active 
- ! list_act_reverse(i) = k ::> IS the kth active orbital
   END_DOC
- logical                        :: exists
+  integer                        :: i
- integer                        :: j,i
   
   n_core_orb = 0
- n_inact_orb = 0
- n_act_orb = 0
- n_virt_orb = 0
- n_del_orb = 0
   do i = 1, mo_num
     if(mo_class(i) == 'Core')then
       n_core_orb += 1
-  else if (mo_class(i) == 'Inactive')then
+    endif
+  enddo
+  
+  call write_int(6,n_core_orb, 'Number of core     MOs')
+   
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_inact_orb ]
+  implicit none
+  BEGIN_DOC
+  ! Number of inactive MOs
+  END_DOC
+  integer                        :: i
+  
+  n_inact_orb = 0
+  do i = 1, mo_num
+    if (mo_class(i) == 'Inactive')then
       n_inact_orb += 1
-  else if (mo_class(i) == 'Active')then
+    endif
+  enddo
+  
+  call write_int(6,n_inact_orb,'Number of inactive MOs')
+  
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_act_orb]
+  implicit none
+  BEGIN_DOC
+  ! Number of active MOs
+  END_DOC
+  integer                        :: i
+  
+  n_act_orb = 0
+  do i = 1, mo_num
+    if (mo_class(i) == 'Active')then
       n_act_orb += 1
-  else if (mo_class(i) == 'Virtual')then
+    endif
+  enddo
+  
+  call write_int(6,n_act_orb, 'Number of active   MOs')
+   
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_virt_orb ]
+  implicit none
+  BEGIN_DOC
+  ! Number of virtual MOs
+  END_DOC
+  integer                        :: i
+  
+  n_virt_orb = 0
+  do i = 1, mo_num
+    if (mo_class(i) == 'Virtual')then
       n_virt_orb += 1
-  else if (mo_class(i) == 'Deleted')then
+    endif
+  enddo
+  
+  call write_int(6,n_virt_orb, 'Number of virtual  MOs')
+   
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_del_orb ]
+  implicit none
+  BEGIN_DOC
+  ! Number of deleted MOs
+  END_DOC
+  integer                        :: i
+  
+  n_del_orb = 0
+  do i = 1, mo_num
+    if (mo_class(i) == 'Deleted')then
       n_del_orb += 1
     endif
   enddo
   
-
- call write_int(6,n_core_orb, 'Number of core     MOs')
- call write_int(6,n_inact_orb,'Number of inactive MOs')
- call write_int(6,n_act_orb, 'Number of active   MOs')
- call write_int(6,n_virt_orb, 'Number of virtual  MOs')
   call write_int(6,n_del_orb, 'Number of deleted  MOs')
    
 END_PROVIDER
 
 
- BEGIN_PROVIDER [integer, dim_list_core_orb]
+BEGIN_PROVIDER [ integer, n_core_inact_orb ]
-&BEGIN_PROVIDER [integer, dim_list_inact_orb]
-&BEGIN_PROVIDER [integer, dim_list_virt_orb]
-&BEGIN_PROVIDER [integer, dim_list_act_orb]
-&BEGIN_PROVIDER [integer, dim_list_del_orb]
   implicit none
   BEGIN_DOC
-! dimensions for the allocation of list_inact, list_virt, list_core and list_act
+  ! n_core + n_inact
+  END_DOC
+  integer                        :: i
+  n_core_inact_orb = 0
+  do i = 1, N_int
+    n_core_inact_orb += popcnt(reunion_of_core_inact_bitmask(i,1))
+  enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [integer, n_inact_act_orb ]
+   implicit none
+  BEGIN_DOC
+  ! n_inact + n_act
+  END_DOC
+   n_inact_act_orb = (n_inact_orb+n_act_orb)
+END_PROVIDER
+ 
+BEGIN_PROVIDER [integer, dim_list_core_orb]
+  implicit none
+  BEGIN_DOC
+  ! dimensions for the allocation of list_core.
   ! it is at least 1
   END_DOC
    dim_list_core_orb = max(n_core_orb,1)
+END_PROVIDER
+
+BEGIN_PROVIDER [integer, dim_list_inact_orb]
+   implicit none
+   BEGIN_DOC
+   ! dimensions for the allocation of list_inact.
+   ! it is at least 1
+   END_DOC
    dim_list_inact_orb = max(n_inact_orb,1)
- dim_list_virt_orb = max(n_virt_orb,1)
+END_PROVIDER
+
+BEGIN_PROVIDER [integer, dim_list_act_orb]
+   implicit none
+   BEGIN_DOC
+   ! dimensions for the allocation of list_act.
+   ! it is at least 1
+   END_DOC
    dim_list_act_orb = max(n_act_orb,1)
+END_PROVIDER
+
+BEGIN_PROVIDER [integer, dim_list_virt_orb]
+   implicit none
+   BEGIN_DOC
+   ! dimensions for the allocation of list_virt.
+   ! it is at least 1
+   END_DOC
+   dim_list_virt_orb = max(n_virt_orb,1)
+END_PROVIDER
+
+BEGIN_PROVIDER [integer, dim_list_del_orb]
+   implicit none
+   BEGIN_DOC
+   ! dimensions for the allocation of list_del.
+   ! it is at least 1
+   END_DOC
    dim_list_del_orb = max(n_del_orb,1)
 END_PROVIDER
 
- BEGIN_PROVIDER [ integer, list_inact, (dim_list_inact_orb)]
+BEGIN_PROVIDER [integer, n_core_inact_act_orb ]
-&BEGIN_PROVIDER [ integer, list_virt, (dim_list_virt_orb)]
+  implicit none
-&BEGIN_PROVIDER [ integer, list_inact_reverse, (mo_num)]
+  BEGIN_DOC
-&BEGIN_PROVIDER [ integer, list_virt_reverse, (mo_num)]
+  !  Number of core inactive and active MOs
-&BEGIN_PROVIDER [ integer, list_del_reverse, (mo_num)]
+  END_DOC
-&BEGIN_PROVIDER [ integer, list_del, (mo_num)]
+  n_core_inact_act_orb = (n_core_orb + n_inact_orb + n_act_orb)
-&BEGIN_PROVIDER [integer, list_core, (dim_list_core_orb)]
+END_PROVIDER
-&BEGIN_PROVIDER [integer, list_core_reverse, (mo_num)]
+ 
-&BEGIN_PROVIDER [integer, list_act, (dim_list_act_orb)]
+
-&BEGIN_PROVIDER [integer, list_act_reverse, (mo_num)]
+
-&BEGIN_PROVIDER [ integer(bit_kind), core_bitmask,  (N_int,2)]
+ 
+ BEGIN_PROVIDER [ integer(bit_kind), core_bitmask , (N_int,2) ]
 &BEGIN_PROVIDER [ integer(bit_kind), inact_bitmask, (N_int,2) ]
 &BEGIN_PROVIDER [ integer(bit_kind), act_bitmask  , (N_int,2) ]
 &BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask , (N_int,2) ]
 &BEGIN_PROVIDER [ integer(bit_kind), del_bitmask  , (N_int,2) ]
    implicit none
    BEGIN_DOC
- ! inact_bitmask : Bitmask of the inactive orbitals which are supposed to be doubly excited 
+   ! Bitmask identifying the core/inactive/active/virtual/deleted MOs 
- ! in post CAS methods
- ! n_inact_orb   : Number of inactive orbitals
- ! virt_bitmask  : Bitmaks of vritual orbitals which are supposed to be recieve electrons 
- ! in post CAS methods
- ! n_virt_orb    : Number of virtual orbitals
- ! list_inact : List of the inactive orbitals which are supposed to be doubly excited 
- ! in post CAS methods
- ! list_virt  : List of vritual orbitals which are supposed to be recieve electrons 
- ! in post CAS methods
- ! list_inact_reverse : reverse list of inactive orbitals 
- ! list_inact_reverse(i) = 0 ::> not an inactive 
- ! list_inact_reverse(i) = k ::> IS the kth inactive 
- ! list_virt_reverse : reverse list of virtual orbitals 
- ! list_virt_reverse(i) = 0 ::> not an virtual 
- ! list_virt_reverse(i) = k ::> IS the kth virtual 
- ! list_act(i) = index of the ith active orbital
- ! 
- ! list_act_reverse : reverse list of active orbitals 
- ! list_act_reverse(i) = 0 ::> not an active 
- ! list_act_reverse(i) = k ::> IS the kth active orbital
    END_DOC
- logical                        :: exists
+
- integer                        :: j,i
- integer                        :: n_core_orb_tmp, n_inact_orb_tmp, n_act_orb_tmp, n_virt_orb_tmp,n_del_orb_tmp
- integer                        :: list_core_tmp(N_int*bit_kind_size)
- integer                        :: list_inact_tmp(N_int*bit_kind_size)
- integer                        :: list_act_tmp(N_int*bit_kind_size)
- integer                        :: list_virt_tmp(N_int*bit_kind_size)
- integer                        :: list_del_tmp(N_int*bit_kind_size)
- list_core = 0
- list_inact = 0
- list_act = 0
- list_virt = 0
- list_del = 0
- list_core_reverse = 0
- list_inact_reverse = 0
- list_act_reverse = 0
- list_virt_reverse = 0
- list_del_reverse = 0
- n_core_orb_tmp = 0
- n_inact_orb_tmp = 0
- n_act_orb_tmp = 0
- n_virt_orb_tmp = 0
- n_del_orb_tmp = 0
    core_bitmask  = 0_bit_kind
    inact_bitmask = 0_bit_kind
    act_bitmask   = 0_bit_kind
    virt_bitmask  = 0_bit_kind
- do i = 1, mo_num
+   del_bitmask   = 0_bit_kind
-  if(mo_class(i) == 'Core')then
-   n_core_orb_tmp += 1
-   list_core(n_core_orb_tmp) = i
-   list_core_tmp(n_core_orb_tmp) = i
-   list_core_reverse(i) = n_core_orb_tmp
-  else if (mo_class(i) == 'Inactive')then
-   n_inact_orb_tmp += 1
-   list_inact(n_inact_orb_tmp) = i
-   list_inact_tmp(n_inact_orb_tmp) = i
-   list_inact_reverse(i) = n_inact_orb_tmp
-  else if (mo_class(i) == 'Active')then
-   n_act_orb_tmp += 1
-   list_act(n_act_orb_tmp) = i
-   list_act_tmp(n_act_orb_tmp) = i
-   list_act_reverse(i) = n_act_orb_tmp
-  else if (mo_class(i) == 'Virtual')then
-   n_virt_orb_tmp += 1
-   list_virt(n_virt_orb_tmp) = i
-   list_virt_tmp(n_virt_orb_tmp) = i
-   list_virt_reverse(i) = n_virt_orb_tmp
-  else if (mo_class(i) == 'Deleted')then
-   n_del_orb_tmp += 1
-   list_del(n_del_orb_tmp) = i
-   list_del_tmp(n_del_orb_tmp) = i
-   list_del_reverse(i) = n_del_orb_tmp
-  endif
- enddo
    
- if(n_core_orb.ne.0)then
+   if(n_core_orb > 0)then
      call list_to_bitstring( core_bitmask(1,1), list_core, n_core_orb, N_int)
      call list_to_bitstring( core_bitmask(1,2), list_core, n_core_orb, N_int)
    endif
- if(n_inact_orb.ne.0)then
+   if(n_inact_orb > 0)then
      call list_to_bitstring( inact_bitmask(1,1), list_inact, n_inact_orb, N_int)
      call list_to_bitstring( inact_bitmask(1,2), list_inact, n_inact_orb, N_int)
    endif
- if(n_act_orb.ne.0)then
+   if(n_act_orb > 0)then
      call list_to_bitstring( act_bitmask(1,1), list_act, n_act_orb, N_int)
      call list_to_bitstring( act_bitmask(1,2), list_act, n_act_orb, N_int)
    endif
- if(n_virt_orb.ne.0)then
+   if(n_virt_orb > 0)then
      call list_to_bitstring( virt_bitmask(1,1), list_virt, n_virt_orb, N_int)
      call list_to_bitstring( virt_bitmask(1,2), list_virt, n_virt_orb, N_int)
    endif
- if(n_del_orb.ne.0)then
+   if(n_del_orb > 0)then
      call list_to_bitstring( del_bitmask(1,1), list_del, n_del_orb, N_int)
      call list_to_bitstring( del_bitmask(1,2), list_del, n_del_orb, N_int)
    endif
    
-
 END_PROVIDER
 
-BEGIN_PROVIDER [integer, n_inact_act_orb ]
+
+
+
+ 
+ BEGIN_PROVIDER [ integer, list_core        , (dim_list_core_orb) ]
+&BEGIN_PROVIDER [ integer, list_core_reverse, (mo_num) ]
    implicit none
- n_inact_act_orb = (n_inact_orb+n_act_orb)
+   BEGIN_DOC
+   ! List of MO indices which are in the core.
+   END_DOC
+   integer                        :: i, n
+   list_core = 0
+   list_core_reverse = 0
+
+   n=0
+   do i = 1, mo_num
+     if(mo_class(i) == 'Core')then
+       n += 1
+       list_core(n) = i
+       list_core_reverse(i) = n
+     endif
+   enddo
+   print *,  'Core MOs:'
+   print *,  list_core(1:n_core_orb)
    
 END_PROVIDER
  
-BEGIN_PROVIDER [integer, list_inact_act, (n_inact_act_orb)]
+ BEGIN_PROVIDER [ integer, list_inact        , (dim_list_inact_orb) ]
+&BEGIN_PROVIDER [ integer, list_inact_reverse, (mo_num) ]
+   implicit none
+   BEGIN_DOC
+   ! List of MO indices which are inactive.
+   END_DOC
+   integer                        :: i, n
+   list_inact = 0
+   list_inact_reverse = 0
+
+   n=0
+   do i = 1, mo_num
+     if (mo_class(i) == 'Inactive')then
+       n += 1
+       list_inact(n) = i
+       list_inact_reverse(i) = n
+     endif
+   enddo
+   print *,  'Inactive MOs:'
+   print *,  list_inact(1:n_inact_orb)
+   
+END_PROVIDER
+ 
+ BEGIN_PROVIDER [ integer, list_virt        , (dim_list_virt_orb) ]
+&BEGIN_PROVIDER [ integer, list_virt_reverse, (mo_num) ]
+   implicit none
+   BEGIN_DOC
+   ! List of MO indices which are virtual
+   END_DOC
+   integer                        :: i, n
+   list_virt = 0
+   list_virt_reverse = 0
+
+   n=0
+   do i = 1, mo_num
+     if (mo_class(i) == 'Virtual')then
+       n += 1
+       list_virt(n) = i
+       list_virt_reverse(i) = n
+     endif
+   enddo
+   print *,  'Virtual MOs:'
+   print *,  list_virt(1:n_virt_orb)
+   
+END_PROVIDER
+ 
+ BEGIN_PROVIDER [ integer, list_del        , (dim_list_del_orb) ]
+&BEGIN_PROVIDER [ integer, list_del_reverse, (mo_num) ]
+   implicit none
+   BEGIN_DOC
+   ! List of MO indices which are deleted.
+   END_DOC
+   integer                        :: i, n
+   list_del = 0
+   list_del_reverse = 0
+
+   n=0
+   do i = 1, mo_num
+     if (mo_class(i) == 'Deleted')then
+       n += 1
+       list_del(n) = i
+       list_del_reverse(i) = n
+     endif
+   enddo
+   print *,  'Deleted MOs:'
+   print *,  list_del(1:n_del_orb)
+   
+END_PROVIDER
+ 
+ BEGIN_PROVIDER [ integer, list_act        , (dim_list_act_orb) ]
+&BEGIN_PROVIDER [ integer, list_act_reverse, (mo_num) ]
+   implicit none
+   BEGIN_DOC
+   ! List of MO indices which are in the active.
+   END_DOC
+   integer                        :: i, n
+   list_act = 0
+   list_act_reverse = 0
+
+   n=0
+   do i = 1, mo_num
+     if (mo_class(i) == 'Active')then
+       n += 1
+       list_act(n) = i
+       list_act_reverse(i) = n
+     endif
+   enddo
+   print *,  'Active MOs:'
+   print *,  list_act(1:n_act_orb)
+   
+END_PROVIDER
+ 
+
+ 
+ BEGIN_PROVIDER [ integer, list_core_inact        , (n_core_inact_orb) ]
+&BEGIN_PROVIDER [ integer, list_core_inact_reverse, (mo_num) ]
+   implicit none
+   BEGIN_DOC
+   ! List of indices of the core and inactive MOs
+   END_DOC
    integer                        :: i,itmp
- itmp = 0
+   call bitstring_to_list(reunion_of_core_inact_bitmask(1,1), list_core_inact, itmp, N_int)
- do i = 1, n_inact_orb
+   list_core_inact_reverse = 0
-  itmp += 1
+   ASSERT (itmp == n_core_inact_orb)
-  list_inact_act(itmp) = list_inact(i) 
+   do i = 1, n_core_inact_orb
- enddo
+     list_core_inact_reverse(list_core_inact(i)) = i
- do i = 1, n_act_orb
-  itmp += 1
-  list_inact_act(itmp) = list_act(i) 
    enddo
+   print *,  'Core and Inactive MOs:'
+   print *,  list_core_inact(1:n_core_inact_orb)
 END_PROVIDER
  
-BEGIN_PROVIDER [integer, n_core_inact_act_orb ]
- implicit none
- n_core_inact_act_orb = (n_core_orb + n_inact_orb + n_act_orb)
-
-END_PROVIDER 
  
  BEGIN_PROVIDER [ integer, list_core_inact_act        , (n_core_inact_act_orb) ]
-&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (n_core_inact_act_orb)]
+&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (mo_num) ]
+   implicit none
+   BEGIN_DOC
+   ! List of indices of the core inactive and active MOs
+   END_DOC
    integer                        :: i,itmp
- itmp = 0
+   call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), list_core_inact_act, itmp, N_int)
- do i = 1, n_core_orb
+   list_core_inact_act_reverse = 0
-  itmp += 1
+   ASSERT (itmp == n_core_inact_act_orb)
-  list_core_inact_act(itmp) = list_core(i) 
- enddo
- do i = 1, n_inact_orb
-  itmp += 1
-  list_core_inact_act(itmp) = list_inact(i) 
- enddo
- do i = 1, n_act_orb
-  itmp += 1
-  list_core_inact_act(itmp) = list_act(i) 
- enddo
-
- integer :: occ_inact(N_int*bit_kind_size)
- occ_inact = 0
- call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), occ_inact(1), itest, N_int)
- list_inact_reverse = 0
    do i = 1, n_core_inact_act_orb
-  list_core_inact_act_reverse(occ_inact(i)) = i
+     list_core_inact_act_reverse(list_core_inact_act(i)) = i
    enddo
+   print *,  'Core, Inactive and Active MOs:'
+   print *,  list_core_inact_act(1:n_core_inact_act_orb)
 END_PROVIDER
+ 
+ 
+ BEGIN_PROVIDER [ integer, list_inact_act        , (n_inact_act_orb) ]
+&BEGIN_PROVIDER [ integer, list_inact_act_reverse, (mo_num) ]
+   implicit none
+   BEGIN_DOC
+   ! List of indices of the inactive and active MOs
+   END_DOC
+   integer                        :: i,itmp
+   call bitstring_to_list(reunion_of_inact_act_bitmask(1,1), list_inact_act, itmp, N_int)
+   list_inact_act_reverse = 0
+   ASSERT (itmp == n_inact_act_orb)
+   do i = 1, n_inact_act_orb
+     list_inact_act_reverse(list_inact_act(i)) = i
+   enddo
+   print *,  'Inactive and Active MOs:'
+   print *,  list_inact_act(1:n_inact_act_orb)
+END_PROVIDER
+ 

src/casscf/bielec.irp.f

@@ -1,65 +1,58 @@
- BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)]
+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
-  double precision, allocatable  :: integrals_array(:,:)
   real*8                         :: mo_two_e_integral
   
-  allocate(integrals_array(mo_num,mo_num))
+  bielec_PQxx(:,:,:,:) = 0.d0
+  PROVIDE mo_two_e_integrals_in_map
   
-  bielec_PQxx = 0.d0
+  !$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)
 
-  do i=1,n_core_orb
+  !$OMP DO
-    ii=list_core(i)
+  do i=1,n_core_inact_orb
-    do j=i,n_core_orb
+    ii=list_core_inact(i)
-      jj=list_core(j)
+    do j=i,n_core_inact_orb
-      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map)
+      jj=list_core_inact(j)
-      do p=1,mo_num
+      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j),mo_integrals_map)
-        do q=1,mo_num
+      bielec_PQxx(:,:,j,i)=bielec_PQxx(:,:,i,j)
-          bielec_PQxx(p,q,i,j)=integrals_array(p,q)
-          bielec_PQxx(p,q,j,i)=integrals_array(p,q)
-        end do
-      end do
     end do
     do j=1,n_act_orb
       jj=list_act(j)
-      j3=j+n_core_orb
+      j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map)
+      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j3),mo_integrals_map)
-      do p=1,mo_num
+      bielec_PQxx(:,:,j3,i)=bielec_PQxx(:,:,i,j3)
-        do q=1,mo_num
-          bielec_PQxx(p,q,i,j3)=integrals_array(p,q)
-          bielec_PQxx(p,q,j3,i)=integrals_array(p,q)
-        end do
-      end do
     end do
   end do
+  !$OMP END DO
 
 
-  ! (ij|pq)
+  !$OMP DO
   do i=1,n_act_orb
     ii=list_act(i)
-    i3=i+n_core_orb
+    i3=i+n_core_inact_orb
     do j=i,n_act_orb
       jj=list_act(j)
-      j3=j+n_core_orb
+      j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map)
+      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i3,j3),mo_integrals_map)
-      do p=1,mo_num
+      bielec_PQxx(:,:,j3,i3)=bielec_PQxx(:,:,i3,j3)
-        do q=1,mo_num
-          bielec_PQxx(p,q,i3,j3)=integrals_array(p,q)
-          bielec_PQxx(p,q,j3,i3)=integrals_array(p,q)
-        end do
-      end do
     end do
   end do
+  !$OMP END DO
+
+  !$OMP END PARALLEL
 
 END_PROVIDER
 
 
 
-BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
+BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_orb+n_act_orb,n_core_inact_orb+n_act_orb, mo_num)]
   BEGIN_DOC
   ! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
   ! indices are unshifted orbital numbers
@@ -69,17 +62,24 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_a
   double precision, allocatable  :: integrals_array(:,:)
   real*8                         :: mo_two_e_integral
   
-  allocate(integrals_array(mo_num,mo_num))
+  PROVIDE mo_two_e_integrals_in_map
-  
   bielec_PxxQ = 0.d0
 
-  do i=1,n_core_orb
+  !$OMP PARALLEL DEFAULT(NONE) &
-    ii=list_core(i)
+  !$OMP PRIVATE(i,ii,j,jj,i3,j3,integrals_array) &
-    do j=i,n_core_orb
+  !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ, &
-      jj=list_core(j)
+  !$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 p=1,mo_num
       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
@@ -87,34 +87,41 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_a
     end do
     do j=1,n_act_orb
       jj=list_act(j)
-      j3=j+n_core_orb
+      j3=j+n_core_inact_orb
       call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
-      do p=1,mo_num
       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_orb
+    i3=i+n_core_inact_orb
     do j=i,n_act_orb
       jj=list_act(j)
-      j3=j+n_core_orb
+      j3=j+n_core_inact_orb
       call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
-      do p=1,mo_num
       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
 
 
@@ -125,24 +132,25 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
   END_DOC
   implicit none
   integer                        :: i,j,k,p,t,u,v
-  double precision, allocatable  :: integrals_array(:)
+  double precision, external     :: mo_two_e_integral
-  real*8                         :: mo_two_e_integral
+  PROVIDE mo_two_e_integrals_in_map
   
-  allocate(integrals_array(mo_num))
+  !$OMP PARALLEL DO DEFAULT(NONE) &
-  
+  !$OMP PRIVATE(i,j,k,p,t,u,v) &
-  do i=1,n_act_orb
+  !$OMP SHARED(mo_num,n_act_orb,list_act,bielecCI)
-    t=list_act(i)
+  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)
-        ! (tu|vp)
+        do i=1,n_act_orb
-        call get_mo_two_e_integrals(t,u,v,mo_num,integrals_array,mo_integrals_map)
+          t=list_act(i)
-        do p=1,mo_num
+          bielecCI(i,k,j,p) = mo_two_e_integral(t,u,v,p)
-          bielecCI(i,k,j,p)=integrals_array(p)
         end do
       end do
     end do
   end do
+  !$OMP END PARALLEL DO
+
 END_PROVIDER
   

src/casscf/bielec_natorb.irp.f

@@ -1,180 +1,264 @@
- BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)]
+ 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,pp
-  real*8                         :: d(n_act_orb)
+  double precision, allocatable  :: f(:,:,:), d(:,:,:)
 
-  bielec_PQxx_no(:,:,:,:) = bielec_PQxx(:,:,:,:)
+
+
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP PRIVATE(j,k,l,p,pp,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
+          pp=n_act_orb-p+1
+          bielec_PQxx_no(list_act(p),j,k,l)=d(pp,j,k)
+        end do
+      end do
 
       do j=1,mo_num
-    do k=1,n_core_orb+n_act_orb
-      do l=1,n_core_orb+n_act_orb
         do p=1,n_act_orb
-          d(p)=0.D0
+          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
         pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PQxx_no(list_act(q),j,k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PQxx_no(list_act(p),j,k,l)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 2nd quarter
         do j=1,mo_num
-    do k=1,n_core_orb+n_act_orb
+          bielec_PQxx_no(j,list_act(p),k,l)=d(pp,j,k)
-      do l=1,n_core_orb+n_act_orb
+        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
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PQxx_no(j,list_act(q),k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PQxx_no(j,list_act(p),k,l)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 3rd quarter
-  do j=1,mo_num
       do k=1,mo_num
-      do l=1,n_core_orb+n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PQxx_no(j,k,n_core_orb+q,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PQxx_no(j,k,n_core_orb+p,l)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 4th quarter
         do j=1,mo_num
-    do k=1,mo_num
+          f(j,k,p) = bielec_PQxx_no(j,k,n_core_inact_orb+p,l)
-      do l=1,n_core_orb+n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
         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
       pp=n_act_orb-p+1
-          do q=1,n_act_orb
+      do k=1,mo_num
-            d(pp)+=bielec_PQxx_no(j,k,l,n_core_orb+q)*natorbsCI(q,p)
+        do j=1,mo_num
+          bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,pp)
         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
-          bielec_PQxx_no(j,k,l,n_core_orb+p)=d(p)
+      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
+      pp=n_act_orb-p+1
+      do k=1,mo_num
+        do j=1,mo_num
+          bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,pp)
         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_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
+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,pp
-  real*8                         :: d(n_act_orb)
+  double precision, allocatable  :: f(:,:,:), d(:,:,:)
 
-  bielec_PxxQ_no(:,:,:,:) = bielec_PxxQ(:,:,:,:)
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP PRIVATE(j,k,l,p,pp,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
-    do k=1,n_core_orb+n_act_orb
+    bielec_PxxQ_no(:,:,:,j) = bielec_PxxQ(:,:,:,j)
-      do l=1,n_core_orb+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
-          d(p)=0.D0             
+            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
           pp=n_act_orb-p+1
-          do q=1,n_act_orb
+          bielec_PxxQ_no(list_act(p),k,l,j)=d(pp,k,l)
-            d(pp)+=bielec_PxxQ_no(list_act(q),k,l,j)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PxxQ_no(list_act(p),k,l,j)=d(p)
         end do
       end do
     end do
   end do
-  ! 2nd quarter
+  !$OMP END DO NOWAIT
-  do j=1,mo_num
+
-    do k=1,n_core_orb+n_act_orb
+  deallocate (f,d)
-      do l=1,n_core_orb+n_act_orb
+
-        do p=1,n_act_orb
+  allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), &
-          d(p)=0.D0
+    d(n_act_orb,mo_num,n_core_inact_act_orb))
-        end do
+
-        do p=1,n_act_orb
+  !$OMP DO
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PxxQ_no(j,k,l,list_act(q))*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PxxQ_no(j,k,l,list_act(p))=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 3rd quarter
-  do j=1,mo_num
   do k=1,mo_num
-      do l=1,n_core_orb+n_act_orb
+    do l=1,n_core_inact_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PxxQ_no(j,n_core_orb+q,l,k)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PxxQ_no(j,n_core_orb+p,l,k)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 4th quarter
       do j=1,mo_num
-    do k=1,mo_num
-      do l=1,n_core_orb+n_act_orb
         do p=1,n_act_orb
-          d(p)=0.D0
+          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
           pp=n_act_orb-p+1
-          do q=1,n_act_orb
+          bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(pp,j,l)
-            d(pp)+=bielec_PxxQ_no(j,l,n_core_orb+q,k)*natorbsCI(q,p)
         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
-          bielec_PxxQ_no(j,l,n_core_orb+p,k)=d(p)
+      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
+      pp=n_act_orb-p+1
+      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,pp)
         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
+      pp=n_act_orb-p+1
+      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,pp)
+        end do
+      end do
+    end do
+  end do
+  !$OMP END DO NOWAIT
+  deallocate(f,d)
+  !$OMP END PARALLEL
 
 END_PROVIDER
 
@@ -186,85 +270,112 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
   END_DOC
   implicit none
   integer                        :: i,j,k,l,t,u,p,q,pp
-  real*8                         :: d(n_act_orb)
+  double precision, allocatable  :: f(:,:,:), d(:,:,:)
   
-  bielecCI_no(:,:,:,:) = bielecCI(:,:,:,:)
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP PRIVATE(j,k,l,p,pp,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
+          pp=n_act_orb-p+1
+          bielecCI_no(p,j,k,l)=d(pp,j,k)
+        end do
+      end do
 
       do j=1,n_act_orb
-    do k=1,n_act_orb
-      do l=1,mo_num
         do p=1,n_act_orb
-          d(p)=0.D0
+          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
         pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielecCI_no(q,j,k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielecCI_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
+          bielecCI_no(j,p,k,l)=d(pp,j,k)
-      do l=1,mo_num
-        do p=1,n_act_orb
-          d(p)=0.D0
         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
       pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielecCI_no(j,q,k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielecCI_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,mo_num
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielecCI_no(j,k,q,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielecCI_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
+          bielecCI_no(j,k,p,l)=d(j,k,pp)
+        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
-          d(p)=0.D0
+      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
       pp=n_act_orb-p+1
-          do q=1,n_act_orb
+      do k=1,n_act_orb
-            d(pp)+=bielecCI_no(j,k,l,list_act(q))*natorbsCI(q,p)
+        do j=1,n_act_orb
-          end do
+          bielecCI_no(j,k,l,list_act(p))=d(j,k,pp)
-        end do
-        do p=1,n_act_orb
-          bielecCI_no(j,k,l,list_act(p))=d(p)
         end do
       end do 
     end do  
   end do  
+  !$OMP END DO
+  
+  deallocate(d,f)
+  !$OMP END PARALLEL
+
 
 END_PROVIDER
 

src/casscf/gradient.irp.f

@@ -5,7 +5,7 @@ BEGIN_PROVIDER [ integer, nMonoEx ]
   ! Number of single excitations
   END_DOC
   implicit none
-  nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb
+  nMonoEx=n_core_inact_orb*n_act_orb+n_core_inact_orb*n_virt_orb+n_act_orb*n_virt_orb
 END_PROVIDER
 
  BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
@@ -17,8 +17,8 @@ END_PROVIDER
   implicit none
   integer                        :: i,t,a,ii,tt,aa,indx
   indx=0
-  do ii=1,n_core_orb
+  do ii=1,n_core_inact_orb
-    i=list_core(ii)
+    i=list_core_inact(ii)
     do tt=1,n_act_orb
       t=list_act(tt)
       indx+=1
@@ -28,8 +28,8 @@ END_PROVIDER
     end do
   end do
   
-  do ii=1,n_core_orb
+  do ii=1,n_core_inact_orb
-    i=list_core(ii)
+    i=list_core_inact(ii)
     do aa=1,n_virt_orb
       a=list_virt(aa)
       indx+=1
@@ -145,14 +145,14 @@ BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
   real*8                         :: norm_grad
   
   indx=0
-  do i=1,n_core_orb
+  do i=1,n_core_inact_orb
     do t=1,n_act_orb
       indx+=1
       gradvec2(indx)=gradvec_it(i,t)
     end do
   end do
   
-  do i=1,n_core_orb
+  do i=1,n_core_inact_orb
     do a=1,n_virt_orb
       indx+=1
       gradvec2(indx)=gradvec_ia(i,a)
@@ -181,7 +181,7 @@ END_PROVIDER
 
 real*8 function gradvec_it(i,t)
   BEGIN_DOC
-  ! the orbital gradient core -> active
+  ! the orbital gradient core/inactive -> active
   ! we assume natural orbitals
   END_DOC
   implicit none
@@ -190,16 +190,16 @@ real*8 function gradvec_it(i,t)
   integer                        :: ii,tt,v,vv,x,y
   integer                        :: x3,y3
   
-  ii=list_core(i)
+  ii=list_core_inact(i)
   tt=list_act(t)
   gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii))
   gradvec_it-=occnum(tt)*Fipq(ii,tt)
   do v=1,n_act_orb
     vv=list_act(v)
     do x=1,n_act_orb
-      x3=x+n_core_orb
+      x3=x+n_core_inact_orb
       do y=1,n_act_orb
-        y3=y+n_core_orb
+        y3=y+n_core_inact_orb
         gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3)
       end do
     end do
@@ -209,12 +209,12 @@ end function gradvec_it
 
 real*8 function gradvec_ia(i,a)
   BEGIN_DOC
-  ! the orbital gradient core -> virtual
+  ! the orbital gradient core/inactive -> virtual
   END_DOC
   implicit none
   integer                        :: i,a,ii,aa
   
-  ii=list_core(i)
+  ii=list_core_inact(i)
   aa=list_virt(a)
   gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii))
   gradvec_ia*=2.D0

src/casscf/hessian.irp.f

@@ -204,10 +204,10 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
   endif
   
   indx=1
-  do i=1,n_core_orb
+  do i=1,n_core_inact_orb
     do t=1,n_act_orb
       jndx=indx
-      do j=i,n_core_orb
+      do j=i,n_core_inact_orb
         if (i.eq.j) then
           ustart=t
         else
@@ -219,7 +219,7 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
           jndx+=1
         end do
       end do
-      do j=1,n_core_orb
+      do j=1,n_core_inact_orb
         do a=1,n_virt_orb
           hessmat2(indx,jndx)=hessmat_itja(i,t,j,a)
           hessmat2(jndx,indx)=hessmat2(indx,jndx)
@@ -237,10 +237,10 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
     end do
   end do
   
-  do i=1,n_core_orb
+  do i=1,n_core_inact_orb
     do a=1,n_virt_orb
       jndx=indx
-      do j=i,n_core_orb
+      do j=i,n_core_inact_orb
         if (i.eq.j) then
           bstart=a
         else
@@ -286,7 +286,7 @@ END_PROVIDER
 
 real*8 function hessmat_itju(i,t,j,u)
   BEGIN_DOC
-  ! the orbital hessian for core->act,core->act
+  ! 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
@@ -295,7 +295,7 @@ real*8 function hessmat_itju(i,t,j,u)
   integer                        :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj
   real*8                         :: term,t2
   
-  ii=list_core(i)
+  ii=list_core_inact(i)
   tt=list_act(t)
   if (i.eq.j) then
     if (t.eq.u) then
@@ -343,7 +343,7 @@ real*8 function hessmat_itju(i,t,j,u)
     end if
   else
     ! it/ju
-    jj=list_core(j)
+    jj=list_core_inact(j)
     uu=list_act(u)
     if (t.eq.u) then
       term=occnum(tt)*Fipq(ii,jj)
@@ -374,16 +374,16 @@ end function hessmat_itju
 
 real*8 function hessmat_itja(i,t,j,a)
   BEGIN_DOC
-  ! the orbital hessian for core->act,core->virt
+  ! 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(i)
+  ii=list_core_inact(i)
   tt=list_act(t)
-  jj=list_core(j)
+  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))
@@ -407,17 +407,17 @@ end function hessmat_itja
 
 real*8 function hessmat_itua(i,t,u,a)
   BEGIN_DOC
-  ! the orbital hessian for core->act,act->virt
+  ! 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(i)
+  ii=list_core_inact(i)
   tt=list_act(t)
-  t3=t+n_core_orb
+  t3=t+n_core_inact_orb
   uu=list_act(u)
-  u3=u+n_core_orb
+  u3=u+n_core_inact_orb
   aa=list_virt(a)
   if (t.eq.u) then
     term=-occnum(tt)*Fipq(aa,ii)
@@ -428,11 +428,11 @@ real*8 function hessmat_itua(i,t,u,a)
       +bielec_pxxq_no(aa,t3,u3,ii))
   do v=1,n_act_orb
     vv=list_act(v)
-    v3=v+n_core_orb
+    v3=v+n_core_inact_orb
     do x=1,n_act_orb
       integer                        :: x3
       xx=list_act(x)
-      x3=x+n_core_orb
+      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))
@@ -448,13 +448,13 @@ end function hessmat_itua
 
 real*8 function hessmat_iajb(i,a,j,b)
   BEGIN_DOC
-  ! the orbital hessian for core->virt,core->virt
+  ! 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(i)
+  ii=list_core_inact(i)
   aa=list_virt(a)
   if (i.eq.j) then
     if (a.eq.b) then
@@ -469,7 +469,7 @@ real*8 function hessmat_iajb(i,a,j,b)
     end if
   else
     ! ia/jb
-    jj=list_core(j)
+    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))
@@ -484,17 +484,17 @@ end function hessmat_iajb
 
 real*8 function hessmat_iatb(i,a,t,b)
   BEGIN_DOC
-  ! the orbital hessian for core->virt,act->virt
+  ! 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(i)
+  ii=list_core_inact(i)
   aa=list_virt(a)
   tt=list_act(t)
   bb=list_virt(b)
-  t3=t+n_core_orb
+  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
@@ -533,10 +533,10 @@ real*8 function hessmat_taub(t,a,u,b)
       t1-=occnum(tt)*Fipq(tt,tt)
       do v=1,n_act_orb
         vv=list_act(v)
-        v3=v+n_core_orb
+        v3=v+n_core_inact_orb
         do x=1,n_act_orb
           xx=list_act(x)
-          x3=x+n_core_orb
+          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))
@@ -552,10 +552,10 @@ real*8 function hessmat_taub(t,a,u,b)
       term=occnum(tt)*Fipq(aa,bb)
       do v=1,n_act_orb
         vv=list_act(v)
-        v3=v+n_core_orb
+        v3=v+n_core_inact_orb
         do x=1,n_act_orb
           xx=list_act(x)
-          x3=x+n_core_orb
+          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))
@@ -569,10 +569,10 @@ real*8 function hessmat_taub(t,a,u,b)
     term=0.D0
     do v=1,n_act_orb
       vv=list_act(v)
-      v3=v+n_core_orb
+      v3=v+n_core_inact_orb
       do x=1,n_act_orb
         xx=list_act(x)
-        x3=x+n_core_orb
+        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))
@@ -606,14 +606,14 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)]
   real*8                         :: hessmat_itju,hessmat_iajb,hessmat_taub
   
   indx=0
-  do i=1,n_core_orb
+  do i=1,n_core_inact_orb
     do t=1,n_act_orb
       indx+=1
       hessdiag(indx)=hessmat_itju(i,t,i,t)
     end do
   end do
   
-  do i=1,n_core_orb
+  do i=1,n_core_inact_orb
     do a=1,n_virt_orb
       indx+=1
       hessdiag(indx)=hessmat_iajb(i,a,i,a)

src/casscf/mcscf_fock.irp.f

@@ -12,8 +12,8 @@ BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ]
    end do
    
    ! the inactive Fock matrix
-   do k=1,n_core_orb
+   do k=1,n_core_inact_orb
-     kk=list_core(k)
+     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)

src/casscf/natorb.irp.f

@@ -6,8 +6,8 @@
   
   integer :: i
   occnum=0.D0
-  do i=1,n_core_orb
+  do i=1,n_core_inact_orb
-    occnum(list_core(i))=2.D0
+    occnum(list_core_inact(i))=2.D0
   end do
   
   do i=1,n_act_orb

src/casscf/neworbs.irp.f

@@ -125,8 +125,8 @@ BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
   ! the orbital rotation matrix T
   Tmat(:,:)=0.D0
   indx=1
-  do i=1,n_core_orb
+  do i=1,n_core_inact_orb
-    ii=list_core(i)
+    ii=list_core_inact(i)
     do t=1,n_act_orb
       tt=list_act(t)
       indx+=1
@@ -134,8 +134,8 @@ BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
       Tmat(tt,ii)=-SXvector(indx)
     end do
   end do
-  do i=1,n_core_orb
+  do i=1,n_core_inact_orb
-    ii=list_core(i)
+    ii=list_core_inact(i)
     do a=1,n_virt_orb
       aa=list_virt(a)
       indx+=1

src/casscf/tot_en.irp.f

@@ -10,19 +10,19 @@
    real*8                         :: e_one_all,e_two_all
    e_one_all=0.D0
    e_two_all=0.D0
-   do i=1,n_core_orb
+   do i=1,n_core_inact_orb
-     ii=list_core(i)
+     ii=list_core_inact(i)
      e_one_all+=2.D0*mo_one_e_integrals(ii,ii)
-     do j=1,n_core_orb
+     do j=1,n_core_inact_orb
-       jj=list_core(j)
+       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_orb
+       t3=t+n_core_inact_orb
        do u=1,n_act_orb
          uu=list_act(u)
-         u3=u+n_core_orb
+         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
@@ -34,9 +34,9 @@
        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_orb
+         v3=v+n_core_inact_orb
          do x=1,n_act_orb
-           x3=x+n_core_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
@@ -44,12 +44,12 @@
    end do
    ecore    =nuclear_repulsion
    ecore_bis=nuclear_repulsion
-   do i=1,n_core_orb
+   do i=1,n_core_inact_orb
-     ii=list_core(i)
+     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_orb
+     do j=1,n_core_inact_orb
-       jj=list_core(j)
+       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
@@ -61,14 +61,14 @@
    etwo_ter=0.D0
    do t=1,n_act_orb
      tt=list_act(t)
-     t3=t+n_core_orb
+     t3=t+n_core_inact_orb
      do u=1,n_act_orb
        uu=list_act(u)
-       u3=u+n_core_orb
+       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_orb
+       do i=1,n_core_inact_orb
-         ii=list_core(i)
+         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)           &
@@ -76,10 +76,10 @@
        end do
        do v=1,n_act_orb
          vv=list_act(v)
-         v3=v+n_core_orb
+         v3=v+n_core_inact_orb
          do x=1,n_act_orb
            xx=list_act(x)
-           x3=x+n_core_orb
+           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)

src/selectors_full/selectors.irp.f

@@ -38,22 +38,6 @@ END_PROVIDER
   END_DOC
   integer                        :: i,k
 
-!  if (threshold_selectors == 1.d0) then
-!
-!    do i=1,N_det_selectors
-!      do k=1,N_int
-!        psi_selectors(k,1,i) = psi_det(k,1,i)
-!        psi_selectors(k,2,i) = psi_det(k,2,i)
-!      enddo
-!    enddo
-!    do k=1,N_states
-!      do i=1,N_det_selectors
-!        psi_selectors_coef(i,k) = psi_coef(i,k)
-!      enddo
-!    enddo
-!
-!  else
-
   do i=1,N_det_selectors
     do k=1,N_int
       psi_selectors(k,1,i) = psi_det_sorted(k,1,i)
@@ -66,7 +50,6 @@ END_PROVIDER
     enddo
   enddo
 
-!  endif
 END_PROVIDER