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

etc/qp.rc

@@ -190,21 +190,8 @@ _qp_Complete()
             ;;
         esac;;
       set_file)
-        # Caching the search results to reduce repeated find calls
+        compopt -o nospace
-        if [[ -z "$QP_FILE_CACHE" || "$CACHE_DIR" != "$PWD" ]]; then
+        COMPREPLY=( $(compgen -d -- "$cur") )
-          CACHE_DIR="$PWD"
-          QP_FILE_CACHE=$(find . -type f -name .version -exec dirname {} \; | sed 's/\/\.version$//')
-        fi
-
-        # Support for relative paths
-        prefix=$(dirname "${cur}")
-        if [[ "$prefix" != "." ]]; then
-          dirs=$(echo "$QP_FILE_CACHE" | grep "^$prefix")
-        else
-          dirs="$QP_FILE_CACHE"
-        fi
-
-        COMPREPLY=( $(compgen -W "$dirs" -- "$cur") )
         return 0
         ;;
       plugins)

ocaml/Input_ao_basis.ml

@@ -13,8 +13,6 @@ module Ao_basis : sig
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
-      ao_normalized   : bool;
-      primitives_normalized : bool;
     } [@@deriving sexp]
   ;;
   val read : unit -> t option
@@ -36,8 +34,6 @@ end = struct
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
-      ao_normalized   : bool;
-      primitives_normalized : bool;
     } [@@deriving sexp]
   ;;
 
@@ -133,24 +129,6 @@ end = struct
     Ezfio.get_ao_basis_ao_cartesian ()
   ;;
 
-  let read_ao_normalized () =
-    if not (Ezfio.has_ao_basis_ao_normalized()) then
-       get_default "ao_normalized"
-       |> bool_of_string
-       |> Ezfio.set_ao_basis_ao_normalized
-    ;
-    Ezfio.get_ao_basis_ao_normalized ()
-  ;;
-
-  let read_primitives_normalized () =
-    if not (Ezfio.has_ao_basis_primitives_normalized()) then
-       get_default "primitives_normalized"
-       |> bool_of_string
-       |> Ezfio.set_ao_basis_primitives_normalized
-    ;
-    Ezfio.get_ao_basis_primitives_normalized ()
-  ;;
-
   let to_long_basis b =
     let ao_num = AO_number.to_int b.ao_num in
     let gto_array = Array.init (AO_number.to_int b.ao_num)
@@ -213,8 +191,6 @@ end = struct
          ao_coef         ;
          ao_expo         ;
          ao_cartesian    ;
-         ao_normalized   ;
-         primitives_normalized ;
        } = b
      in
      write_md5 b ;
@@ -247,8 +223,6 @@ end = struct
      ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ;
 
      Ezfio.set_ao_basis_ao_cartesian(ao_cartesian);
-     Ezfio.set_ao_basis_ao_normalized(ao_normalized);
-     Ezfio.set_ao_basis_primitives_normalized(primitives_normalized);
 
      let ao_coef =
       Array.to_list ao_coef
@@ -280,8 +254,6 @@ end = struct
             ao_coef         = read_ao_coef () ;
             ao_expo         = read_ao_expo () ;
             ao_cartesian    = read_ao_cartesian () ;
-            ao_normalized   = read_ao_normalized () ;
-            primitives_normalized   = read_primitives_normalized () ;
           }
         in
         to_md5 result
@@ -392,8 +364,6 @@ end = struct
       { ao_basis = name ;
         ao_num ; ao_prim_num ; ao_prim_num_max ; ao_nucl ;
         ao_power ; ao_coef ; ao_expo ; ao_cartesian ;
-        ao_normalized = bool_of_string @@ get_default "ao_normalized";
-        primitives_normalized = bool_of_string @@ get_default "primitives_normalized";
         }
   ;;
 
@@ -448,14 +418,6 @@ Cartesian coordinates (6d,10f,...) ::
 
   ao_cartesian = %s
 
-Use normalized primitive functions ::
-
-  primitives_normalized = %s
-
-Use normalized basis functions ::
-
-  ao_normalized = %s
-
 Basis set (read-only) ::
 
 %s
@@ -469,8 +431,6 @@ Basis set (read-only) ::
 
 "   (AO_basis_name.to_string b.ao_basis)
     (string_of_bool b.ao_cartesian)
-    (string_of_bool b.primitives_normalized)
-    (string_of_bool b.ao_normalized)
     (Basis.to_string short_basis
        |> String_ext.split ~on:'\n'
        |> list_map (fun x-> "  "^x)
@@ -507,8 +467,6 @@ ao_power                = %s
 ao_coef                 = %s
 ao_expo                 = %s
 ao_cartesian            = %s
-ao_normalized           = %s
-primitives_normalized   = %s
 md5                     = %s
 "
     (AO_basis_name.to_string b.ao_basis)
@@ -525,8 +483,6 @@ md5                     = %s
     (b.ao_expo  |> Array.to_list |> list_map AO_expo.to_string
       |> String.concat ", ")
     (b.ao_cartesian |> string_of_bool)
-    (b.ao_normalized |> string_of_bool)
-    (b.primitives_normalized |> string_of_bool)
     (to_md5 b |> MD5.to_string )
 
   ;;

scripts/qp_import_trexio.py

@@ -191,11 +191,14 @@ def write_ezfio(trexio_filename, filename):
             ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
 
 
-            shell_factor = trexio.read_basis_shell_factor(trexio_file)
             prim_factor  = trexio.read_basis_prim_factor(trexio_file)
-            for i,p in enumerate(prim_factor):
+            ezfio.set_basis_prim_normalization_factor(prim_factor)
-                coefficient[i] *= prim_factor[i]
+            ezfio.set_basis_primitives_normalized(True)
-            ezfio.set_ao_basis_primitives_normalized(True)
+            ezfio.set_basis_ao_normalized(False)
+
+            shell_factor = trexio.read_basis_shell_factor(trexio_file)
+            for i, shell_idx in enumerate(shell_index):
+               coefficient[i] *= shell_factor[shell_idx]
             ezfio.set_basis_prim_coef(coefficient)
 
         elif basis_type.lower() == "numerical":
@@ -247,7 +250,6 @@ def write_ezfio(trexio_filename, filename):
             ezfio.set_basis_shell_index([x+1 for x in shell_index])
             ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
 
-            shell_factor = trexio.read_basis_shell_factor(trexio_file)
         else:
            raise TypeError
 
@@ -317,13 +319,18 @@ def write_ezfio(trexio_filename, filename):
                 exponent.append(expo[i])
                 num_prim.append(num_prim0[i])
 
-        print (len(coefficient), ao_num)
         assert (len(coefficient) == ao_num)
+
         ezfio.set_ao_basis_ao_power(power_x + power_y + power_z)
         ezfio.set_ao_basis_ao_prim_num(num_prim)
 
         prim_num_max = max( [ len(x) for x in coefficient ] )
 
+        ao_normalization = trexio.read_ao_normalization(trexio_file_cart)
+        for i, coef in enumerate(coefficient):
+           for j in range(len(coef)):
+              coef[j] *= ao_normalization[i]
+
         for i in range(ao_num):
             coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)]
             exponent   [i] += [0. for j in range(len(exponent[i]), prim_num_max)]
@@ -338,7 +345,6 @@ def write_ezfio(trexio_filename, filename):
                 coef.append(coefficient[j])
                 expo.append(exponent[j])
 
-#        ezfio.set_ao_basis_ao_prim_num_max(prim_num_max)
         ezfio.set_ao_basis_ao_coef(coef)
         ezfio.set_ao_basis_ao_expo(expo)
 
@@ -388,14 +394,6 @@ def write_ezfio(trexio_filename, filename):
 
       # Read coefs from temporary cartesian file created in the AO section
       MoMatrix = trexio.read_mo_coefficient(trexio_file_cart)
-
-      # Renormalize MO coefs if needed
-      if trexio.has_ao_normalization(trexio_file_cart):
-        ezfio.set_ao_basis_ao_normalized(False)
-        norm = trexio.read_ao_normalization(trexio_file_cart)
-#        for j in range(mo_num):
-#           for i,f in enumerate(norm):
-#              MoMatrix[i,j] *= f
       ezfio.set_mo_basis_mo_coef(MoMatrix)
 
       mo_occ = [ 0. for i in range(mo_num) ]

src/ao_basis/EZFIO.cfg

@@ -55,18 +55,6 @@ doc: If |true|, use |AOs| in Cartesian coordinates (6d,10f,...)
 interface: ezfio, provider
 default: false
 
-[ao_normalized]
-type: logical
-doc: Use normalized basis functions
-interface: ezfio, provider
-default: true
-
-[primitives_normalized]
-type: logical
-doc: Use normalized primitive functions
-interface: ezfio, provider
-default: true
-
 [use_cgtos]
 type: logical
 doc: If true, use cgtos for AO integrals

src/ao_basis/aos.irp.f

@@ -82,10 +82,11 @@ END_PROVIDER
     enddo
     ao_coef_normalization_factor(i) = 1.d0/dsqrt(norm)
 
-    if (.not.ao_normalized) then
+    if (ao_normalized) then
       do j=1,ao_prim_num(i)
         ao_coef_normalized(i,j) = ao_coef_normalized(i,j) * ao_coef_normalization_factor(i)
       enddo
+    else
       ao_coef_normalization_factor(i) = 1.d0
     endif
   enddo

src/basis/EZFIO.cfg

@@ -73,3 +73,15 @@ size: (basis.prim_num)
 interface: ezfio, provider
 
 
+[primitives_normalized]
+type: logical
+doc: If true, assume primitive basis functions are normalized
+interface: ezfio, provider, ocaml
+default: true
+
+[ao_normalized]
+type: logical
+doc: If true, normalize the basis functions
+interface: ezfio, provider, ocaml
+default: false
+

src/basis/basis.irp.f

@@ -6,6 +6,11 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
 
   logical                        :: has
   PROVIDE ezfio_filename
+  if (.not.ao_normalized) then
+    shell_normalization_factor = 1.d0
+    return
+  endif
+
   if (mpi_master) then
     if (size(shell_normalization_factor) == 0) return
 
@@ -70,6 +75,12 @@ BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ]
 
   logical                        :: has
   PROVIDE ezfio_filename
+
+  if (.not.primitives_normalized) then
+    prim_normalization_factor(:) = 1.d0
+    return
+  endif
+
   if (mpi_master) then
     if (size(prim_normalization_factor) == 0) return
 
@@ -95,9 +106,9 @@ BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ]
 
         do k=1, prim_num
           if (shell_index(k) /= i) cycle
-            call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), &
+          call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), &
-              powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+            powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
-            prim_normalization_factor(k) = 1.d0/dsqrt(norm)
+          prim_normalization_factor(k) = 1.d0/dsqrt(norm)
         enddo
       enddo
 

src/casscf_cipsi/bielec.irp.f

@@ -1,7 +1,7 @@
 BEGIN_PROVIDER [real*8, bielec_PQxx_array, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
   BEGIN_DOC
-  ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! 
+  ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!!
-  ! 
+  !
   ! Replaced by the Cholesky-based function bielec_PQxx
   !
   ! bielec_PQxx_array : integral (pq|xx) with p,q arbitrary, x core or active
@@ -13,10 +13,10 @@ BEGIN_PROVIDER [real*8, bielec_PQxx_array, (mo_num, mo_num,n_core_inact_act_orb,
   print*,''
   print*,'Providing bielec_PQxx_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!'
   print*,''
-  
+
   bielec_PQxx_array(:,:,:,:) = 0.d0
-  PROVIDE mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
-  
+
   !$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_array, &
@@ -61,8 +61,8 @@ END_PROVIDER
 
 BEGIN_PROVIDER [real*8, bielec_PxxQ_array, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
   BEGIN_DOC
-  ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! 
+  ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!!
-  ! 
+  !
   ! Replaced by the Cholesky-based function bielec_PxxQ
   !
   ! bielec_PxxQ_array : integral (px|xq) with p,q arbitrary, x core or active
@@ -72,11 +72,11 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_array, (mo_num,n_core_inact_act_orb,n_core_i
   integer                        :: i,j,ii,jj,p,q,i3,j3,t3,v3
   double precision, allocatable  :: integrals_array(:,:)
   real*8                         :: mo_two_e_integral
-  
+
   print*,''
   print*,'Providing bielec_PxxQ_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!'
   print*,''
-  PROVIDE mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
   bielec_PxxQ_array = 0.d0
 
   !$OMP PARALLEL DEFAULT(NONE) &
@@ -85,7 +85,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_array, (mo_num,n_core_inact_act_orb,n_core_i
   !$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)
@@ -143,17 +143,17 @@ 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
-  ! 
+  !
   ! This array can be stored anyway. Ex: 50 active orbitals, 1500 MOs ==> 8x50^3x1500 = 1.5 Gb
   END_DOC
   implicit none
   integer                        :: i,j,k,p,t,u,v
   double precision, external     :: mo_two_e_integral
-  double precision :: wall0, wall1 
+  double precision :: wall0, wall1
   call wall_time(wall0)
   print*,'Providing bielecCI'
-  PROVIDE mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
-  
+
   !$OMP PARALLEL DO DEFAULT(NONE) &
   !$OMP PRIVATE(i,j,k,p,t,u,v) &
   !$OMP SHARED(mo_num,n_act_orb,list_act,bielecCI)

src/casscf_cipsi/hessian.irp.f

@@ -380,7 +380,7 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
  !  c-v |     X  X 
  !  a-v |        X
 
-  provide mo_two_e_integrals_in_map
+  provide all_mo_integrals
 
  hessmat = 0.d0
 

src/casscf_cipsi/hessian_old.irp.f

@@ -13,18 +13,18 @@ BEGIN_PROVIDER [real*8, hessmat_old, (nMonoEx,nMonoEx)]
   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)
@@ -38,7 +38,7 @@ BEGIN_PROVIDER [real*8, hessmat_old, (nMonoEx,nMonoEx)]
       hessmat_old(jndx,indx)=res
     end do
   end do
-  
+
 END_PROVIDER
 
 subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res)
@@ -75,9 +75,9 @@ subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res)
   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
@@ -174,10 +174,10 @@ subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res)
       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)]
@@ -190,26 +190,26 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)]
   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
+  provide all_mo_integrals
-  
+
   !$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 ) 
+  ! (DOUBLY OCCUPIED ---> ACT )
   do i=1,n_core_inact_orb
     do t=1,n_act_orb
       indx = t + (i-1)*n_act_orb
@@ -226,14 +226,14 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)]
           jndx+=1
         end do
       end do
-      ! (DOUBLY OCCUPIED ---> VIRTUAL) 
+      ! (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) 
+      ! (ACTIVE ---> VIRTUAL)
       do u=1,n_act_orb
         do a=1,n_virt_orb
           hessmat_peter(jndx,indx)=hessmat_itua(i,t,u,a)
@@ -243,15 +243,15 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)]
     end do
   end do
   !$OMP END DO NOWAIT
-  
+
   indx_shift = n_core_inact_orb*n_act_orb
   !$OMP DO
-  ! (DOUBLY OCCUPIED ---> VIRTUAL) 
+  ! (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) 
+      ! (DOUBLY OCCUPIED ---> VIRTUAL)
       do j=i,n_core_inact_orb
         if (i.eq.j) then
           bstart=a
@@ -263,7 +263,7 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)]
           jndx+=1
         end do
       end do
-      ! (ACT ---> VIRTUAL) 
+      ! (ACT ---> VIRTUAL)
       do t=1,n_act_orb
         do b=1,n_virt_orb
           hessmat_peter(jndx,indx)=hessmat_iatb(i,a,t,b)
@@ -273,15 +273,15 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)]
     end do
   end do
   !$OMP END DO NOWAIT
-  
+
   indx_shift += n_core_inact_orb*n_virt_orb
-  !$OMP DO 
+  !$OMP DO
-  ! (ACT ---> VIRTUAL) 
+  ! (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) 
+      ! (ACT ---> VIRTUAL)
       do u=t,n_act_orb
         if (t.eq.u) then
           bstart=a
@@ -295,16 +295,16 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)]
       end do
     end do
   end do
-  !$OMP 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/ccsd/ccsd_space_orb_sub.irp.f

@@ -35,7 +35,7 @@ subroutine run_ccsd_space_orb
     PROVIDE cholesky_mo_transp
     FREE cholesky_ao
   else
-    PROVIDE mo_two_e_integrals_in_map
+    PROVIDE all_mo_integrals
   endif
 
   det = psi_det(:,:,cc_ref)

src/cipsi/stochastic_cipsi.irp.f

@@ -1,11 +1,11 @@
-subroutine run_stochastic_cipsi(Ev,PT2) 
+subroutine run_stochastic_cipsi(Ev,PT2)
   use selection_types
   implicit none
   BEGIN_DOC
 ! Selected Full Configuration Interaction with Stochastic selection and PT2.
   END_DOC
   integer                        :: i,j,k
-  double precision, intent(out)  :: Ev(N_states), PT2(N_states) 
+  double precision, intent(out)  :: Ev(N_states), PT2(N_states)
   double precision, allocatable  :: zeros(:)
   integer                        :: to_select
   type(pt2_type)                 :: pt2_data, pt2_data_err
@@ -14,7 +14,7 @@ subroutine run_stochastic_cipsi(Ev,PT2)
 
   double precision :: rss
   double precision, external :: memory_of_double
-  PROVIDE H_apply_buffer_allocated distributed_davidson mo_two_e_integrals_in_map
+  PROVIDE H_apply_buffer_allocated distributed_davidson all_mo_integrals
 
   threshold_generators = 1.d0
   SOFT_TOUCH threshold_generators

src/cipsi_utils/pt2_stoch_routines.irp.f

@@ -165,7 +165,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       state_average_weight(pt2_stoch_istate) = 1.d0
       TOUCH state_average_weight pt2_stoch_istate selection_weight
 
-      PROVIDE nproc pt2_F mo_two_e_integrals_in_map mo_one_e_integrals pt2_w
+      PROVIDE nproc pt2_F all_mo_integrals mo_one_e_integrals pt2_w
       PROVIDE psi_selectors pt2_u pt2_J pt2_R
       call new_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2')
 

src/cipsi_utils/slave_cipsi.irp.f

@@ -14,7 +14,7 @@ subroutine run_slave_cipsi
 end
 
 subroutine provide_everything
-  PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag
+  PROVIDE H_apply_buffer_allocated all_mo_integrals psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag
   PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context
   PROVIDE psi_det psi_coef threshold_generators state_average_weight
   PROVIDE N_det_selectors pt2_stoch_istate N_det selection_weight pseudo_sym

src/davidson/diagonalization_h_dressed.irp.f

@@ -31,7 +31,7 @@ subroutine davidson_diag_h(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nint,
   ASSERT (sze > 0)
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
-  PROVIDE mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
   allocate(H_jj(sze))
 
   H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)

src/davidson/diagonalization_hcsf_dressed.irp.f

@@ -30,7 +30,7 @@ subroutine davidson_diag_h_csf(dets_in,u_in,dim_in,energies,sze,sze_csf,N_st,N_s
   ASSERT (sze > 0)
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
-  PROVIDE mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
   allocate(H_jj(sze))
 
   H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -62,7 +62,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
   ASSERT (sze > 0)
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
-  PROVIDE mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
   allocate(H_jj(sze))
 
   H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)

src/davidson/diagonalization_nonsym_h_dressed.irp.f

@@ -42,7 +42,7 @@ subroutine davidson_diag_nonsym_h(dets_in, u_in, dim_in, energies, sze, N_st, N_
   ASSERT (sze  > 0)
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
-  PROVIDE mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
 
   allocate(H_jj(sze))
 

src/determinants/h_apply_nozmq.template.f

@@ -17,7 +17,7 @@ subroutine $subroutine($params_main)
   double precision, allocatable  :: fock_diag_tmp(:,:)
 
   $initialization
-  PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators
+  PROVIDE H_apply_buffer_allocated all_mo_integrals psi_det_generators psi_coef_generators
 
 
   call wall_time(wall_0)

src/determinants/slater_rules.irp.f

@@ -521,7 +521,7 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
   integer                        :: occ(Nint*bit_kind_size,2)
   double precision               :: diag_H_mat_elem, phase
   integer                        :: n_occ_ab(2)
-  PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
+  PROVIDE all_mo_integrals
 
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
@@ -623,7 +623,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
   integer                        :: occ(Nint*bit_kind_size,2)
   double precision               :: diag_H_mat_elem, phase
   integer                        :: n_occ_ab(2)
-  PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
+  PROVIDE all_mo_integrals
 
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
@@ -724,7 +724,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble,phase)
   integer                        :: n_occ_ab(2)
   logical                        :: has_mipi(Nint*bit_kind_size)
   double precision               :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
-  PROVIDE mo_two_e_integrals_in_map mo_integrals_map
+  PROVIDE all_mo_integrals
 
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
@@ -2227,7 +2227,7 @@ subroutine i_H_j_single_spin(key_i,key_j,Nint,spin,hij)
   integer                        :: exc(0:2,2)
   double precision               :: phase
 
-  PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
 
   call get_single_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint)
   call get_single_excitation_from_fock(key_i,key_j,exc(1,1),exc(1,2),spin,phase,hij)
@@ -2248,7 +2248,7 @@ subroutine i_H_j_double_spin(key_i,key_j,Nint,hij)
   double precision               :: phase
   double precision, external     :: get_two_e_integral
 
-  PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
   call get_double_excitation_spin(key_i,key_j,exc,phase,Nint)
   hij = phase*(get_two_e_integral(                             &
       exc(1,1),                                                    &
@@ -2277,7 +2277,7 @@ subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij)
   double precision               :: phase, phase2
   double precision, external     :: get_two_e_integral
 
-  PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
 
   call get_single_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint)
   call get_single_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint)

src/determinants/slater_rules_wee_mono.irp.f

@@ -13,7 +13,7 @@ subroutine i_Wee_j_single(key_i,key_j,Nint,spin,hij)
   integer                        :: exc(0:2,2)
   double precision               :: phase
 
-  PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
 
   call get_single_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint)
   call single_excitation_wee(key_i,key_j,exc(1,1),exc(1,2),spin,phase,hij)
@@ -285,7 +285,7 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij)
   integer                        :: occ(Nint*bit_kind_size,2)
   double precision               :: diag_H_mat_elem, phase,phase_2
   integer                        :: n_occ_ab(2)
-  PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals ref_bitmask_two_e_energy
+  PROVIDE all_mo_integrals ref_bitmask_two_e_energy
 
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)

src/determinants/utils.irp.f

@@ -4,12 +4,11 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
  BEGIN_DOC
  ! |H| matrix on the basis of the Slater determinants defined by psi_det
  END_DOC
- integer :: i,j,k
+ integer :: i,j
  double precision :: hij
- integer :: degree(N_det),idx(0:N_det)
+ PROVIDE all_mo_integrals
- call  i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij)
  print*,'Providing the H_matrix_all_dets ...'
- !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hij,degree,idx,k) &
+ !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hij) &
  !$OMP SHARED (N_det, psi_det, N_int,H_matrix_all_dets)
  do i =1,N_det
    do j = i, N_det
@@ -31,7 +30,7 @@ BEGIN_PROVIDER [ double precision, H_matrix_diag_all_dets,(N_det) ]
  integer :: i
  double precision :: hij
  integer :: degree(N_det)
- call  i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij)
+ PROVIDE all_mo_integrals
  !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,hij,degree) &
  !$OMP SHARED (N_det, psi_det, N_int,H_matrix_diag_all_dets)
  do i =1,N_det

src/dressing/dress_slave.irp.f

@@ -15,7 +15,7 @@ subroutine dress_slave
 end
 
 subroutine provide_everything
-  PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
+  PROVIDE H_apply_buffer_allocated all_mo_integrals psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
 end
 
 subroutine run_wf

src/dressing/dress_stoch_routines.irp.f

@@ -258,7 +258,7 @@ subroutine ZMQ_dress(E, dress, delta_out, delta_s2_out, relative_error)
     state_average_weight(dress_stoch_istate) = 1.d0
     TOUCH state_average_weight dress_stoch_istate
 
-    provide nproc mo_two_e_integrals_in_map mo_one_e_integrals psi_selectors pt2_F pt2_N_teeth dress_M_m
+    provide nproc all_mo_integrals mo_one_e_integrals psi_selectors pt2_F pt2_N_teeth dress_M_m
 
     print *, '========== ================= ================= ================='
     print *, ' Samples        Energy         Stat. Error         Seconds      '

src/fci/fci.irp.f

@@ -36,8 +36,9 @@ program fci
   !
   END_DOC
 
+  PROVIDE all_mo_integrals
   if (.not.is_zmq_slave) then
-    PROVIDE psi_det psi_coef mo_two_e_integrals_in_map
+    PROVIDE psi_det psi_coef
 
     write(json_unit,json_array_open_fmt) 'fci'
 
@@ -55,7 +56,7 @@ program fci
     call json_close
 
   else
-    PROVIDE mo_two_e_integrals_in_map pt2_min_parallel_tasks
+    PROVIDE pt2_min_parallel_tasks
 
     call run_slave_cipsi
 

src/fci/pt2.irp.f

@@ -11,7 +11,7 @@ program pt2
   !
   ! The main option for the |PT2| correction is the
   ! :option:`perturbation pt2_relative_error` which is the relative
-  ! stochastic error on the |PT2| to reach before stopping the 
+  ! stochastic error on the |PT2| to reach before stopping the
   ! sampling.
   !
   END_DOC
@@ -19,7 +19,7 @@ program pt2
      read_wf = .True.
      threshold_generators = 1.d0
      SOFT_TOUCH read_wf threshold_generators
-     PROVIDE mo_two_e_integrals_in_map
+     PROVIDE all_mo_integrals
      PROVIDE psi_energy
      call run
   else
@@ -32,22 +32,22 @@ subroutine run
   use selection_types
   integer                        :: i,j,k
   logical, external              :: detEq
-  
+
   type(pt2_type)                 :: pt2_data, pt2_data_err
   integer                        :: degree
   integer                        :: n_det_before, to_select
   double precision               :: threshold_davidson_in
-  
+
   double precision               :: relative_error
   double precision, allocatable  :: E_CI_before(:)
-  
+
   allocate ( E_CI_before(N_states))
   call pt2_alloc(pt2_data, N_states)
   call pt2_alloc(pt2_data_err, N_states)
-  
+
   E_CI_before(:) = psi_energy(:) + nuclear_repulsion
   relative_error=PT2_relative_error
-  
+
   if (do_pt2) then
      call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, relative_error, 0) ! Stochastic PT2
   else
@@ -56,7 +56,7 @@ subroutine run
 
   call print_summary(psi_energy_with_nucl_rep(1:N_states), &
        pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2)
-  
+
   call save_energy(E_CI_before, pt2_data % pt2)
   call pt2_dealloc(pt2_data)
   call pt2_dealloc(pt2_data_err)

src/mo_optimization/debug_gradient_list_opt.irp.f

@@ -19,7 +19,7 @@
 
 
 program debug_gradient_list
-  
+
   implicit none
 
   ! Variables
@@ -30,12 +30,12 @@ program debug_gradient_list
   double precision              :: threshold
   double precision              :: max_error, max_elem, norm
   integer                       :: nb_error
-  
+
   m = dim_list_act_orb
-  ! Definition of n  
+  ! Definition of n
   n = m*(m-1)/2
 
-  PROVIDE mo_two_e_integrals_in_map ! Verifier pour suppression
+  PROVIDE all_mo_integrals ! Verifier pour suppression
 
   ! Allocation
   allocate(v_grad(n), v_grad2(n))
@@ -44,15 +44,15 @@ program debug_gradient_list
 
   call diagonalize_ci ! Verifier pour suppression
 
-  ! Gradient  
+  ! Gradient
   call gradient_list_opt(n,m,list_act,v_grad,max_elem,norm)
   call first_gradient_list_opt(n,m,list_act,v_grad2)
-  
+
-  
+
   v_grad = v_grad - v_grad2
   nb_error = 0
-  max_error = 0d0 
+  max_error = 0d0
-  threshold = 1d-12 
+  threshold = 1d-12
 
   do i = 1, n
     if (ABS(v_grad(i)) > threshold) then
@@ -65,9 +65,9 @@ program debug_gradient_list
 
     endif
   enddo
- 
+
   print*,''
-  print*,'Check the gradient' 
+  print*,'Check the gradient'
   print*,'Threshold:', threshold
   print*,'Nb error:', nb_error
   print*,'Max error:', max_error

src/mo_optimization/debug_gradient_opt.irp.f

@@ -19,7 +19,7 @@
 
 
 program debug_gradient
-  
+
   implicit none
 
   ! Variables
@@ -30,27 +30,27 @@ program debug_gradient
   double precision              :: threshold
   double precision              :: max_error, max_elem
   integer                       :: nb_error
-  
+
-  ! Definition of n  
+  ! Definition of n
   n = mo_num*(mo_num-1)/2
 
-  PROVIDE mo_two_e_integrals_in_map ! Check for suppression
+  PROVIDE all_mo_integrals ! Check for suppression
 
   ! Allocation
   allocate(v_grad(n), v_grad2(n))
 
   ! Calculation
 
-  call diagonalize_ci 
+  call diagonalize_ci
 
-  ! Gradient  
+  ! Gradient
   call first_gradient_opt(n,v_grad)
   call gradient_opt(n,v_grad2,max_elem)
-  
+
   v_grad = v_grad - v_grad2
   nb_error = 0
-  max_error = 0d0 
+  max_error = 0d0
-  threshold = 1d-12 
+  threshold = 1d-12
 
   do i = 1, n
     if (ABS(v_grad(i)) > threshold) then
@@ -63,9 +63,9 @@ program debug_gradient
 
     endif
   enddo
- 
+
   print*,''
-  print*,'Check the gradient' 
+  print*,'Check the gradient'
   print*,'Threshold :', threshold
   print*,'Nb error :', nb_error
   print*,'Max error :', max_error

src/mo_optimization/debug_hessian_list_opt.irp.f

@@ -43,18 +43,18 @@ program debug_hessian_list_opt
   double precision              :: max_error, max_error_H
   integer                       :: nb_error, nb_error_H
   double precision              :: threshold
-  
+
   m = dim_list_act_orb !mo_num
 
-  ! Definition of n  
+  ! Definition of n
   n = m*(m-1)/2
 
-  PROVIDE mo_two_e_integrals_in_map 
+  PROVIDE all_mo_integrals
 
-  ! Hessian 
+  ! Hessian
   if (optimization_method == 'full') then
     print*,'Use the full hessian matrix'
-    allocate(H(n,n),H2(n,n))  
+    allocate(H(n,n),H2(n,n))
     allocate(h_f(m,m,m,m),h_f2(m,m,m,m))
 
     call hessian_list_opt(n,m,list_act,H,h_f)
@@ -65,7 +65,7 @@ program debug_hessian_list_opt
     h_f = h_f - h_f2
     H = H - H2
     max_error = 0d0
-    nb_error = 0    
+    nb_error = 0
     threshold = 1d-12
 
     do l = 1, m
@@ -99,7 +99,7 @@ program debug_hessian_list_opt
 
         endif
       enddo
-    enddo 
+    enddo
 
     ! Deallocation
     deallocate(H, H2, h_f, h_f2)
@@ -110,26 +110,26 @@ program debug_hessian_list_opt
     allocate(H(n,1),H2(n,1))
     call diag_hessian_list_opt(n,m,list_act,H)
     call first_diag_hessian_list_opt(n,m,list_act,H2)
-    
+
     H = H - H2
-  
+
     max_error_H = 0d0
     nb_error_H = 0
- 
+
     do i = 1, n
       if (ABS(H(i,1)) > threshold) then
         print*, H(i,1)
         nb_error_H = nb_error_H + 1
- 
+
         if (ABS(H(i,1)) > ABS(max_error_H)) then
           max_error_H = H(i,1)
         endif
- 
+
       endif
     enddo
-   
+
  endif
-  
+
   print*,''
   if (optimization_method == 'full') then
     print*,'Check of the full hessian'
@@ -140,8 +140,8 @@ program debug_hessian_list_opt
   else
     print*,'Check of the diagonal hessian'
   endif
-   
+
   print*,'Nb error_H:', nb_error_H
   print*,'Max error_H:', max_error_H
- 
+
 end program

src/mo_optimization/debug_hessian_opt.irp.f

@@ -36,20 +36,20 @@ program debug_hessian
   double precision              :: max_error, max_error_H
   integer                       :: nb_error, nb_error_H
   double precision              :: threshold
-  
+
-  ! Definition of n  
+  ! Definition of n
   n = mo_num*(mo_num-1)/2
 
-  PROVIDE mo_two_e_integrals_in_map 
+  PROVIDE all_mo_integrals
 
   ! Allocation
-  allocate(H(n,n),H2(n,n))  
+  allocate(H(n,n),H2(n,n))
   allocate(h_f(mo_num,mo_num,mo_num,mo_num),h_f2(mo_num,mo_num,mo_num,mo_num))
 
   ! Calculation
-  
+
-  ! Hessian 
+  ! Hessian
-  if (optimization_method == 'full') then 
+  if (optimization_method == 'full') then
 
     print*,'Use the full hessian matrix'
     call hessian_opt(n,H,h_f)
@@ -59,7 +59,7 @@ program debug_hessian
     h_f = h_f - h_f2
     H = H - H2
     max_error = 0d0
-    nb_error = 0    
+    nb_error = 0
     threshold = 1d-12
 
     do l = 1, mo_num
@@ -93,7 +93,7 @@ program debug_hessian
 
        endif
      enddo
-   enddo 
+   enddo
 
   elseif (optimization_method == 'diag') then
 
@@ -128,43 +128,43 @@ program debug_hessian
     enddo
 
     h=H-H2
-  
+
     max_error_H = 0d0
     nb_error_H = 0
- 
+
     do j = 1, n
       do i = 1, n
         if (ABS(H(i,j)) > threshold) then
           print*, H(i,j)
           nb_error_H = nb_error_H + 1
- 
+
           if (ABS(H(i,j)) > ABS(max_error_H)) then
             max_error_H = H(i,j)
           endif
- 
+
         endif
       enddo
     enddo
-  
+
   else
     print*,'Unknown optimization_method, please select full, diag'
     call abort
   endif
-  
+
   print*,''
   if (optimization_method == 'full') then
     print*,'Check the full hessian'
   else
     print*,'Check the diagonal hessian'
   endif
-   
+
   print*,'Threshold :', threshold
   print*,'Nb error :', nb_error
   print*,'Max error :', max_error
   print*,''
   print*,'Nb error_H :', nb_error_H
   print*,'Max error_H :', max_error_H
- 
+
   ! Deallocation
   deallocate(H,H2,h_f,h_f2)
 

src/mo_optimization_utils/routine_opt_mos.irp.f

@@ -9,7 +9,7 @@ subroutine run_optimization_mos_CIPSI
   logical :: not_converged
   character (len=100) :: filename
 
-  PROVIDE psi_det psi_coef mo_two_e_integrals_in_map ao_pseudo_integrals
+  PROVIDE psi_det psi_coef all_mo_integrals ao_pseudo_integrals
   allocate(Ev(N_states),PT2(N_states))
 
   not_converged = .True.
@@ -30,7 +30,7 @@ subroutine run_optimization_mos_CIPSI
       print*,'======================'
       print*,' Cipsi step:', nb_iter
       print*,'======================'
-      print*,'' 
+      print*,''
       print*,'********** cipsi step **********'
       ! cispi calculation
       call run_stochastic_cipsi(Ev,PT2)
@@ -70,7 +70,7 @@ subroutine run_optimization_mos_CIPSI
           print*, 'The program will exit'
           exit
       endif
-      
+
       ! To double the number of determinants in the wf
       N_det_max = int(dble(n_det * 2)*0.9)
       TOUCH N_det_max

src/mo_optimization_utils/run_orb_opt_trust_v2.irp.f

@@ -93,7 +93,7 @@ subroutine run_orb_opt_trust_v2
   integer,allocatable           :: tmp_list(:), key(:)
   double precision, allocatable :: tmp_m_x(:,:),tmp_R(:,:), tmp_x(:), W(:,:), e_val(:)
 
-  PROVIDE mo_two_e_integrals_in_map ci_energy psi_det psi_coef
+  PROVIDE all_mo_integrals ci_energy psi_det psi_coef
 
 ! Allocation
 
@@ -110,17 +110,17 @@ allocate(tmp_R(m,m), tmp_m_x(m,m), tmp_x(tmp_n))
 allocate(e_val(tmp_n),key(tmp_n),v_grad(tmp_n))
 
 ! Method
-!    There are three different methods : 
+!    There are three different methods :
 !    - the "full" hessian, which uses all the elements of the hessian
 !      matrix"
 !    - the "diagonal" hessian, which uses only the diagonal elements of the
 !      hessian
-!    - without the hessian (hessian = identity matrix) 
+!    - without the hessian (hessian = identity matrix)
 
 
 !Display the method
  print*, 'Method :', optimization_method
-if (optimization_method == 'full') then 
+if (optimization_method == 'full') then
   print*, 'Full hessian'
   allocate(H(tmp_n,tmp_n), h_f(m,m,m,m),W(tmp_n,tmp_n))
   tmp_n2 = tmp_n
@@ -147,13 +147,13 @@ print*, 'Absolute value of the hessian:', absolute_eig
 ! - We diagonalize the hessian
 ! - We compute a step and loop to reduce the radius of the
 !   trust region (and the size of the step by the way) until the step is
-!   accepted 
+!   accepted
-! - We repeat the process until the convergence 
+! - We repeat the process until the convergence
 !   NB: the convergence criterion can be changed
 
 
 ! Loop until the convergence of the optimization
-! call diagonalize_ci 
+! call diagonalize_ci
 
 !### Initialization ###
 nb_iter = 0
@@ -183,14 +183,14 @@ do while (not_converged)
     ! Full hessian
     call hessian_list_opt(tmp_n, m, tmp_list, H, h_f)
 
-    ! Diagonalization of the hessian 
+    ! Diagonalization of the hessian
     call diagonalization_hessian(tmp_n, H, e_val, w)
 
   elseif (optimization_method == 'diag') then
-    ! Diagonal hessian 
+    ! Diagonal hessian
     call diag_hessian_list_opt(tmp_n, m, tmp_list, H)
   else
-    ! Identity matrix 
+    ! Identity matrix
     do tmp_i = 1, tmp_n
       H(tmp_i,1) = 1d0
     enddo
@@ -212,7 +212,7 @@ do while (not_converged)
   endif
 
   ! Init before the internal loop
-  cancel_step = .True. ! To enter in the loop just after 
+  cancel_step = .True. ! To enter in the loop just after
   nb_cancel = 0
   nb_sub_iter = 0
 
@@ -225,15 +225,15 @@ do while (not_converged)
     print*,'Max elem grad:', max_elem_grad
     print*,'-----------------------------'
 
-    ! Hessian,gradient,Criterion -> x 
+    ! Hessian,gradient,Criterion -> x
-    call trust_region_step_w_expected_e(tmp_n,tmp_n2,H,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,tmp_x,must_exit) 
+    call trust_region_step_w_expected_e(tmp_n,tmp_n2,H,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,tmp_x,must_exit)
 
     if (must_exit) then
       print*,'step_in_trust_region sends: Exit'
       exit
     endif
 
-    ! 1D tmp -> 2D tmp 
+    ! 1D tmp -> 2D tmp
     call vec_to_mat_v2(tmp_n, m, tmp_x, tmp_m_x)
 
     ! Rotation matrix for the active MOs
@@ -250,14 +250,14 @@ do while (not_converged)
     call sub_to_full_rotation_matrix(m, tmp_list, tmp_R, R)
 
     ! MO rotations
-    call apply_mo_rotation(R, prev_mos)   
+    call apply_mo_rotation(R, prev_mos)
 
     ! Update of the energy before the diagonalization of the hamiltonian
     call clear_mo_map
-    TOUCH mo_coef psi_det psi_coef ci_energy two_e_dm_mo 
+    TOUCH mo_coef psi_det psi_coef ci_energy two_e_dm_mo
     call state_average_energy(criterion)
 
-    ! Criterion -> step accepted or rejected 
+    ! Criterion -> step accepted or rejected
     call trust_region_is_step_cancelled(nb_iter, prev_criterion, criterion, criterion_model, rho, cancel_step)
 
     ! Cancellation of the step if necessary
@@ -283,7 +283,7 @@ do while (not_converged)
   ! To exit the external loop if must_exit = .True.
   if (must_exit) then
     exit
-  endif 
+  endif
 
   ! Step accepted, nb iteration + 1
   nb_iter = nb_iter + 1
@@ -295,7 +295,7 @@ do while (not_converged)
   endif
   if (nb_iter >= optimization_max_nb_iter) then
     print*,'Not converged: nb_iter >= optimization_max_nb_iter'
-    not_converged = .False. 
+    not_converged = .False.
   endif
 
   if (.not. not_converged) then

src/mo_two_e_ints/map_integrals.irp.f

@@ -1,5 +1,15 @@
 use map_module
 
+BEGIN_PROVIDER [ logical, all_mo_integrals ]
+  implicit none
+  BEGIN_DOC
+! Used to provide everything needed before using MO integrals
+! PROVIDE all_mo_integrals
+  END_DOC
+  PROVIDE mo_two_e_integrals_in_map mo_integrals_cache mo_two_e_integrals_jj_exchange mo_two_e_integrals_jj_anti mo_two_e_integrals_jj big_array_exchange_integrals big_array_coulomb_integrals mo_one_e_integrals
+END_PROVIDER
+
+
 !! MO Map
 !! ======
 
@@ -35,20 +45,24 @@ end
  BEGIN_PROVIDER [ integer, mo_integrals_cache_min ]
 &BEGIN_PROVIDER [ integer, mo_integrals_cache_max ]
 &BEGIN_PROVIDER [ integer, mo_integrals_cache_size ]
+&BEGIN_PROVIDER [ integer*8, mo_integrals_cache_size_8 ]
  implicit none
  BEGIN_DOC
  ! Min and max values of the MOs for which the integrals are in the cache
  END_DOC
 
- mo_integrals_cache_size  = 2**mo_integrals_cache_shift
+ mo_integrals_cache_size  = shiftl(1,mo_integrals_cache_shift)
+ mo_integrals_cache_size_8  = shiftl(1_8, mo_integrals_cache_shift*4)
+
 
  mo_integrals_cache_min = max(1,elec_alpha_num - (mo_integrals_cache_size/2 - 1) )
  mo_integrals_cache_max = min(mo_num, mo_integrals_cache_min + mo_integrals_cache_size - 1)
- print *, 'MO integrals cache: (', mo_integrals_cache_min, ', ', mo_integrals_cache_max, ')'
+ print *, 'MO integrals cache: (', mo_integrals_cache_min, ', ', mo_integrals_cache_max, '), ', &
+          shiftr(mo_integrals_cache_size_8, 17), 'MiB'
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:(1_8*mo_integrals_cache_size)**4) ]
+BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:mo_integrals_cache_size_8) ]
  implicit none
  BEGIN_DOC
  ! Cache of MO integrals for fast access
@@ -67,8 +81,7 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:(1_8*mo_integrals_ca
      do k=mo_integrals_cache_min,mo_integrals_cache_max
          ii = int(l-mo_integrals_cache_min,8)
          ii = ior( shiftl(ii,mo_integrals_cache_shift), int(k-mo_integrals_cache_min,8))
-         ii = shiftl(ii,mo_integrals_cache_shift)
+         ii = shiftl(ii,2*mo_integrals_cache_shift)
-         ii = shiftl(ii,mo_integrals_cache_shift)
          call dgemm('T','N', mo_integrals_cache_max-mo_integrals_cache_min+1, &
                              mo_integrals_cache_max-mo_integrals_cache_min+1, &
            cholesky_mo_num, 1.d0, &
@@ -130,7 +143,7 @@ double precision function get_two_e_integral(i,j,k,l,map)
   END_DOC
   integer, intent(in)            :: i,j,k,l
   integer(key_kind)              :: idx
-  integer                        :: ii
+  integer                        :: ii, kk
   type(map_type), intent(inout)  :: map
   real(integral_kind)            :: tmp
 
@@ -166,8 +179,11 @@ double precision function get_two_e_integral(i,j,k,l,map)
 
       double precision, external :: ddot
       get_two_e_integral = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, cholesky_mo_transp(1,j,l), 1)
-!       double precision, external :: get_from_mo_cholesky_cache
+
-!       get_two_e_integral = get_from_mo_cholesky_cache(i,j,k,l,.False.)
+!      get_two_e_integral = 0.d0
+!      do kk=1,cholesky_mo_num
+!        get_two_e_integral = get_two_e_integral + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l)
+!      enddo
 
     else
 
@@ -328,7 +344,7 @@ double precision function mo_two_e_integral(i,j,k,l)
   END_DOC
   integer, intent(in)            :: i,j,k,l
   double precision               :: get_two_e_integral
-  PROVIDE mo_two_e_integrals_in_map mo_integrals_cache
+  PROVIDE all_mo_integrals
   !DIR$ FORCEINLINE
   mo_two_e_integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
   return
@@ -503,11 +519,46 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
 
     if (do_mo_cholesky) then
 
-      double precision, external :: ddot
+      if ( (k>=mo_integrals_cache_min).and.(k<=mo_integrals_cache_max).and. &
-      do i=1,sze
+           (l>=mo_integrals_cache_min).and.(l<=mo_integrals_cache_max) ) then
-        out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, &
+
-                                           cholesky_mo_transp(1,i,l), 1)
+        double precision, external :: ddot
-      enddo
+        integer :: kk
+
+        do i=1,mo_integrals_cache_min-1
+          out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, &
+                                             cholesky_mo_transp(1,i,l), 1)
+!          out_val(i) = 0.d0
+!          do kk=1,cholesky_mo_num
+!            out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l)
+!          enddo
+        enddo
+
+        do i=mo_integrals_cache_min,mo_integrals_cache_max
+          out_val(i) = get_two_e_integral_cache(i,i,k,l)
+        enddo
+
+        do i=mo_integrals_cache_max, sze
+          out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, &
+                                             cholesky_mo_transp(1,i,l), 1)
+!          out_val(i) = 0.d0
+!          do kk=1,cholesky_mo_num
+!            out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l)
+!          enddo
+        enddo
+
+      else
+
+        do i=1,sze
+          out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, &
+                                             cholesky_mo_transp(1,i,l), 1)
+!          out_val(i) = 0.d0
+!          do kk=1,cholesky_mo_num
+!            out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l)
+!          enddo
+        enddo
+
+      endif
 
     else
 

src/mu_of_r/f_hf_utils.irp.f

@@ -1,41 +1,41 @@
 BEGIN_PROVIDER [double precision, two_e_int_hf_f, (n_basis_orb,n_basis_orb,n_max_occ_val_orb_for_hf,n_max_occ_val_orb_for_hf)]
  implicit none
  BEGIN_DOC
-! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) 
+! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space)
-! 
+!
-! needed to compute the function f_{HF}(r_1,r_2) 
+! needed to compute the function f_{HF}(r_1,r_2)
 !
 ! two_e_int_hf_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis"
  END_DOC
  integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
  double precision :: get_two_e_integral
- PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
+ PROVIDE all_mo_integrals big_array_exchange_integrals
- do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1 
+ do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1
   m = list_valence_orb_for_hf(orb_m,1)
-  do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2 
+  do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2
    n = list_valence_orb_for_hf(orb_n,1)
-   do orb_i = 1, n_basis_orb ! electron 1 
+   do orb_i = 1, n_basis_orb ! electron 1
     i = list_basis(orb_i)
-    do orb_j = 1, n_basis_orb ! electron 2 
+    do orb_j = 1, n_basis_orb ! electron 2
      j = list_basis(orb_j)
      !              2       1     2     1
-     two_e_int_hf_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) 
+     two_e_int_hf_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map)
     enddo
    enddo
   enddo
  enddo
-END_PROVIDER 
+END_PROVIDER
 
 subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens)
  implicit none
  BEGIN_DOC
-! f_HF_val_ab(r1,r2) = function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) 
+! f_HF_val_ab(r1,r2) = function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018)
-! 
-! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931   2937)
-! 
-! two_bod_dens = on-top pair density of the HF wave function 
 !
-! < HF | wee_{\alpha\beta} | HF > =  \int (r1,r2) f_HF_ab(r1,r2) excluding all contributions from "core" "electrons" 
+! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931   2937)
+!
+! two_bod_dens = on-top pair density of the HF wave function
+!
+! < HF | wee_{\alpha\beta} | HF > =  \int (r1,r2) f_HF_ab(r1,r2) excluding all contributions from "core" "electrons"
  END_DOC
  double precision, intent(in) :: r1(3), r2(3)
  double precision, intent(out):: f_HF_val_ab,two_bod_dens
@@ -50,8 +50,8 @@ subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens)
  allocate(mos_array_valence_r1(n_basis_orb) , mos_array_valence_r2(n_basis_orb), mos_array_r1(mo_num), mos_array_r2(mo_num))
  allocate(mos_array_valence_hf_r1(n_occ_val_orb_for_hf(1)) , mos_array_valence_hf_r2(n_occ_val_orb_for_hf(2)) )
  ! You get all orbitals in r_1 and r_2
- call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r1,mos_array_r1)
- call give_all_mos_at_r(r2,mos_array_r2) 
+ call give_all_mos_at_r(r2,mos_array_r2)
  ! You extract the occupied ALPHA/BETA orbitals belonging to the space \mathcal{A}
  do i_m = 1, n_occ_val_orb_for_hf(1)
   mos_array_valence_hf_r1(i_m) = mos_array_r1(list_valence_orb_for_hf(i_m,1))
@@ -61,7 +61,7 @@ subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens)
  enddo
 
  ! You extract the orbitals belonging to the space \mathcal{B}
- do i_m = 1, n_basis_orb 
+ do i_m = 1, n_basis_orb
   mos_array_valence_r1(i_m) = mos_array_r1(list_basis(i_m))
   mos_array_valence_r2(i_m) = mos_array_r2(list_basis(i_m))
  enddo
@@ -70,24 +70,24 @@ subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens)
  f_HF_val_ab = 0.d0
  two_bod_dens = 0.d0
  ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space
- do m = 1, n_occ_val_orb_for_hf(1)! electron 1 
+ do m = 1, n_occ_val_orb_for_hf(1)! electron 1
  ! You browse all OCCUPIED BETA  electrons in the \mathcal{A} space
-  do n = 1, n_occ_val_orb_for_hf(2)! electron 2 
+  do n = 1, n_occ_val_orb_for_hf(2)! electron 2
    ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2)
-   two_bod_dens += mos_array_valence_hf_r1(m) * mos_array_valence_hf_r1(m) * mos_array_valence_hf_r2(n) * mos_array_valence_hf_r2(n) 
+   two_bod_dens += mos_array_valence_hf_r1(m) * mos_array_valence_hf_r1(m) * mos_array_valence_hf_r2(n) * mos_array_valence_hf_r2(n)
-   ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space 
+   ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space
    do i = 1, n_basis_orb
     do j = 1, n_basis_orb
      !                                             2 1 2 1
-     f_HF_val_ab +=  two_e_int_hf_f(j,i,n,m)  & 
+     f_HF_val_ab +=  two_e_int_hf_f(j,i,n,m)  &
-     * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m)  & 
+     * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m)  &
-     * mos_array_valence_r2(j) * mos_array_valence_hf_r2(n)    
+     * mos_array_valence_r2(j) * mos_array_valence_hf_r2(n)
     enddo
    enddo
   enddo
  enddo
  ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
- f_HF_val_ab *= 2.d0 
+ f_HF_val_ab *= 2.d0
  two_bod_dens *= 2.d0
 end
 
@@ -95,14 +95,14 @@ end
 subroutine integral_f_HF_valence_ab(r1,int_f_HF_val_ab)
  implicit none
  BEGIN_DOC
-! in_f_HF_val_ab(r_1) = \int dr_2 f_{\Psi^B}(r_1,r_2) 
+! in_f_HF_val_ab(r_1) = \int dr_2 f_{\Psi^B}(r_1,r_2)
-! 
+!
-! where f_{\Psi^B}(r_1,r_2) is defined by Eq. (22) of J. Chem. Phys. 149, 194301 (2018) 
+! where f_{\Psi^B}(r_1,r_2) is defined by Eq. (22) of J. Chem. Phys. 149, 194301 (2018)
 !
 ! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931   2937)
 !
-! Such function can be used to test if the f_HF_val_ab(r_1,r_2) is correctly built. 
+! Such function can be used to test if the f_HF_val_ab(r_1,r_2) is correctly built.
-! 
+!
 ! < HF | wee_{\alpha\beta} | HF > =  \int (r1) int_f_HF_val_ab(r_1)
  END_DOC
  double precision, intent(in) :: r1(3)
@@ -114,28 +114,28 @@ subroutine integral_f_HF_valence_ab(r1,int_f_HF_val_ab)
  double precision, allocatable  :: mos_array_valence_r1(:)
  double precision, allocatable  :: mos_array_valence_hf_r1(:)
  double precision :: get_two_e_integral
- call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r1,mos_array_r1)
  allocate(mos_array_valence_r1( n_basis_orb ))
  allocate(mos_array_valence_hf_r1( n_occ_val_orb_for_hf(1) ) )
  do i_m = 1, n_occ_val_orb_for_hf(1)
   mos_array_valence_hf_r1(i_m) = mos_array_r1(list_valence_orb_for_hf(i_m,1))
  enddo
- do i_m = 1, n_basis_orb 
+ do i_m = 1, n_basis_orb
   mos_array_valence_r1(i_m) = mos_array_r1(list_basis(i_m))
  enddo
 
  int_f_HF_val_ab = 0.d0
  ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space
- do m = 1, n_occ_val_orb_for_hf(1)! electron 1 
+ do m = 1, n_occ_val_orb_for_hf(1)! electron 1
  ! You browse all OCCUPIED BETA  electrons in the \mathcal{A} space
-  do n = 1, n_occ_val_orb_for_hf(2)! electron 2 
+  do n = 1, n_occ_val_orb_for_hf(2)! electron 2
-   ! You browse all ORBITALS in the \mathacal{B} space 
+   ! You browse all ORBITALS in the \mathacal{B} space
    do i = 1, n_basis_orb
-   ! due to integration in real-space and the use of orthonormal MOs, a Kronecker delta_jn shoes up 
+   ! due to integration in real-space and the use of orthonormal MOs, a Kronecker delta_jn shoes up
     j = n
      !                                             2 1 2 1
-     int_f_HF_val_ab +=  two_e_int_hf_f(j,i,n,m)  & 
+     int_f_HF_val_ab +=  two_e_int_hf_f(j,i,n,m)  &
-     * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m) 
+     * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m)
    enddo
   enddo
  enddo

src/mu_of_r/f_psi_i_a_v_utils.irp.f

@@ -1,7 +1,7 @@
 subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens)
  implicit none
  BEGIN_DOC
-! contribution from purely inactive orbitals to f_{\Psi^B}(r_1,r_2) for a CAS wave function 
+! contribution from purely inactive orbitals to f_{\Psi^B}(r_1,r_2) for a CAS wave function
  END_DOC
  double precision, intent(in) :: r1(3),r2(3)
  double precision, intent(out):: f_ii_val_ab,two_bod_dens
@@ -9,13 +9,13 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens)
  integer :: i_i,i_j
  double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:)
  double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:)
- double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) 
+ double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:)
  double precision :: get_two_e_integral
  ! You get all orbitals in r_1 and r_2
  allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
- call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r1,mos_array_r1)
- call give_all_mos_at_r(r2,mos_array_r2) 
+ call give_all_mos_at_r(r2,mos_array_r2)
- ! You extract the inactive orbitals 
+ ! You extract the inactive orbitals
  allocate(mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) )
  do i_m = 1, n_inact_orb
   mos_array_inact_r1(i_m) = mos_array_r1(list_inact(i_m))
@@ -26,7 +26,7 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens)
 
  ! You extract the orbitals belonging to the space \mathcal{B}
  allocate(mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) )
- do i_m = 1, n_basis_orb 
+ do i_m = 1, n_basis_orb
   mos_array_basis_r1(i_m) = mos_array_r1(list_basis(i_m))
   mos_array_basis_r2(i_m) = mos_array_r2(list_basis(i_m))
  enddo
@@ -34,30 +34,30 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens)
  f_ii_val_ab = 0.d0
  two_bod_dens = 0.d0
  ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space
- do m = 1, n_inact_orb ! electron 1 
+ do m = 1, n_inact_orb ! electron 1
  ! You browse all OCCUPIED BETA  electrons in the \mathcal{A} space
-  do n = 1, n_inact_orb ! electron 2 
+  do n = 1, n_inact_orb ! electron 2
    ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2)
-   two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_r2(n) 
+   two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_r2(n)
-   ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space 
+   ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space
    do i = 1, n_basis_orb
     do j = 1, n_basis_orb
      !                              2 1 2 1
-     f_ii_val_ab +=  two_e_int_ii_f(j,i,n,m)   * mos_array_inact_r1(m) * mos_array_basis_r1(i)  & 
+     f_ii_val_ab +=  two_e_int_ii_f(j,i,n,m)   * mos_array_inact_r1(m) * mos_array_basis_r1(i)  &
-                                               * mos_array_inact_r2(n) * mos_array_basis_r2(j)    
+                                               * mos_array_inact_r2(n) * mos_array_basis_r2(j)
     enddo
    enddo
   enddo
  enddo
  ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
- f_ii_val_ab  *= 2.d0 
+ f_ii_val_ab  *= 2.d0
  two_bod_dens *= 2.d0
 end
 
 
 subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
  BEGIN_DOC
-! contribution from inactive and active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function 
+! contribution from inactive and active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function
  END_DOC
  implicit none
  integer, intent(in) :: istate
@@ -65,7 +65,7 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
  double precision, intent(out):: f_ia_val_ab,two_bod_dens
  integer :: i,orb_i,a,orb_a,n,m,b
  double precision :: rho
- double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) 
+ double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:)
  double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:)
  double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:)
  double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:)
@@ -75,10 +75,10 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
  two_bod_dens = 0.d0
  ! You get all orbitals in r_1 and r_2
  allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
- call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r1,mos_array_r1)
- call give_all_mos_at_r(r2,mos_array_r2) 
+ call give_all_mos_at_r(r2,mos_array_r2)
 
- ! You extract the inactive orbitals 
+ ! You extract the inactive orbitals
  allocate( mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) )
  do i = 1, n_inact_orb
   mos_array_inact_r1(i) = mos_array_r1(list_inact(i))
@@ -87,7 +87,7 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
   mos_array_inact_r2(i) = mos_array_r2(list_inact(i))
  enddo
 
- ! You extract the active orbitals 
+ ! You extract the active orbitals
  allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) )
  do i= 1, n_act_orb
   mos_array_act_r1(i) = mos_array_r1(list_act(i))
@@ -109,11 +109,11 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
  ! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2)
  allocate(rho_tilde(n_inact_orb,n_act_orb))
  two_bod_dens = 0.d0
- do a = 1, n_act_orb 
+ do a = 1, n_act_orb
   do i = 1, n_inact_orb
    rho_tilde(i,a) = 0.d0
-   do b = 1, n_act_orb 
+   do b = 1, n_act_orb
-    rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate) 
+    rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate)
     two_bod_dens += mos_array_inact_r1(i) * mos_array_inact_r1(i) * mos_array_act_r2(a) * mos_array_act_r2(b) * rho
     rho_tilde(i,a) += rho * mos_array_inact_r1(i) * mos_array_act_r2(b)
    enddo
@@ -125,12 +125,12 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
  allocate( v_tilde(n_inact_orb,n_act_orb)   )
  allocate( integrals_array(mo_num,mo_num) )
  v_tilde = 0.d0
- do a = 1, n_act_orb 
+ do a = 1, n_act_orb
   orb_a = list_act(a)
   do i = 1, n_inact_orb
    v_tilde(i,a) = 0.d0
    orb_i = list_inact(i)
-!   call get_mo_two_e_integrals_ij(orb_i,orb_a,mo_num,integrals_array,mo_integrals_map) 
+!   call get_mo_two_e_integrals_ij(orb_i,orb_a,mo_num,integrals_array,mo_integrals_map)
    do m = 1, n_basis_orb
     do n = 1, n_basis_orb
 !     v_tilde(i,a) += integrals_array(n,m) * mos_array_basis_r2(n) * mos_array_basis_r1(m)
@@ -146,14 +146,14 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
   enddo
  enddo
  ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
- f_ia_val_ab *= 2.d0 
+ f_ia_val_ab *= 2.d0
  two_bod_dens *= 2.d0
 end
 
 
 subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
  BEGIN_DOC
-! contribution from purely active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function 
+! contribution from purely active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function
  END_DOC
  implicit none
  integer, intent(in) :: istate
@@ -161,7 +161,7 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
  double precision, intent(out):: f_aa_val_ab,two_bod_dens
  integer :: i,orb_i,a,orb_a,n,m,b,c,d
  double precision :: rho
- double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) 
+ double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:)
  double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:)
  double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:)
  double precision, allocatable :: integrals_array(:,:),rho_tilde(:,:),v_tilde(:,:)
@@ -170,10 +170,10 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
  two_bod_dens = 0.d0
  ! You get all orbitals in r_1 and r_2
  allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
- call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r1,mos_array_r1)
- call give_all_mos_at_r(r2,mos_array_r2) 
+ call give_all_mos_at_r(r2,mos_array_r2)
 
- ! You extract the active orbitals 
+ ! You extract the active orbitals
  allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) )
  do i= 1, n_act_orb
   mos_array_act_r1(i) = mos_array_r1(list_act(i))
@@ -201,8 +201,8 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
    do c = 1, n_act_orb  ! 1
     do d = 1, n_act_orb  ! 2
      rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate)
-     rho_tilde(b,a) += rho 
+     rho_tilde(b,a) += rho
-     two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) 
+     two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b)
     enddo
    enddo
   enddo
@@ -212,7 +212,7 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
  ! v_tilde(i,a) = \sum_{m,n} phi_m(1) * phi_n(2) < i a | m n >
  allocate( v_tilde(n_act_orb,n_act_orb)   )
  v_tilde = 0.d0
- do a = 1, n_act_orb 
+ do a = 1, n_act_orb
   do b = 1, n_act_orb
    v_tilde(b,a) = 0.d0
    do m = 1, n_basis_orb
@@ -235,92 +235,92 @@ end
 BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act_orb,n_act_orb)]
  implicit none
  BEGIN_DOC
-! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) 
+! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space)
-! 
+!
-! needed to compute the function f_{ii}(r_1,r_2) 
+! needed to compute the function f_{ii}(r_1,r_2)
 !
 ! two_e_int_aa_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis"
  END_DOC
  integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
  double precision :: integrals_array(mo_num,mo_num),get_two_e_integral
- PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
+ PROVIDE all_mo_integrals
- do orb_m = 1, n_act_orb ! electron 1 
+ do orb_m = 1, n_act_orb ! electron 1
   m = list_act(orb_m)
-  do orb_n = 1, n_act_orb ! electron 2 
+  do orb_n = 1, n_act_orb ! electron 2
    n = list_act(orb_n)
-   call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) 
+   call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map)
-   do orb_i = 1, n_basis_orb ! electron 1 
+   do orb_i = 1, n_basis_orb ! electron 1
     i = list_basis(orb_i)
-    do orb_j = 1, n_basis_orb ! electron 2 
+    do orb_j = 1, n_basis_orb ! electron 2
      j = list_basis(orb_j)
      !              2       1     2     1
-     two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) 
+     two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map)
-!     two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) 
+!     two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i)
     enddo
    enddo
   enddo
  enddo
-END_PROVIDER 
+END_PROVIDER
 
 BEGIN_PROVIDER [double precision, two_e_int_ia_f, (n_basis_orb,n_basis_orb,n_inact_orb,n_act_orb)]
  implicit none
  BEGIN_DOC
-! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) 
+! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space)
-! 
+!
-! needed to compute the function f_{ia}(r_1,r_2) 
+! needed to compute the function f_{ia}(r_1,r_2)
 !
 ! two_e_int_aa_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis"
  END_DOC
  integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
  double precision :: integrals_array(mo_num,mo_num),get_two_e_integral
- PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
+ PROVIDE all_mo_integrals
- do orb_m = 1, n_act_orb ! electron 1 
+ do orb_m = 1, n_act_orb ! electron 1
   m = list_act(orb_m)
-  do orb_n = 1, n_inact_orb ! electron 2 
+  do orb_n = 1, n_inact_orb ! electron 2
    n = list_inact(orb_n)
-   call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) 
+   call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map)
-   do orb_i = 1, n_basis_orb ! electron 1 
+   do orb_i = 1, n_basis_orb ! electron 1
     i = list_basis(orb_i)
-    do orb_j = 1, n_basis_orb ! electron 2 
+    do orb_j = 1, n_basis_orb ! electron 2
      j = list_basis(orb_j)
      !              2       1     2     1
-!     two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) 
+!     two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map)
-     two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) 
+     two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i)
     enddo
    enddo
   enddo
  enddo
-END_PROVIDER 
+END_PROVIDER
 
 BEGIN_PROVIDER [double precision, two_e_int_ii_f, (n_basis_orb,n_basis_orb,n_inact_orb,n_inact_orb)]
  implicit none
  BEGIN_DOC
-! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) 
+! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space)
-! 
+!
-! needed to compute the function f_{ii}(r_1,r_2) 
+! needed to compute the function f_{ii}(r_1,r_2)
 !
 ! two_e_int_ii_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis"
  END_DOC
  integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
  double precision :: get_two_e_integral,integrals_array(mo_num,mo_num)
- PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
+ PROVIDE all_mo_integrals
- do orb_m = 1, n_inact_orb ! electron 1 
+ do orb_m = 1, n_inact_orb ! electron 1
   m = list_inact(orb_m)
-  do orb_n = 1, n_inact_orb ! electron 2 
+  do orb_n = 1, n_inact_orb ! electron 2
    n = list_inact(orb_n)
-   call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) 
+   call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map)
-   do orb_i = 1, n_basis_orb ! electron 1 
+   do orb_i = 1, n_basis_orb ! electron 1
     i = list_basis(orb_i)
-    do orb_j = 1, n_basis_orb ! electron 2 
+    do orb_j = 1, n_basis_orb ! electron 2
      j = list_basis(orb_j)
      !              2       1     2     1
-!     two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) 
+!     two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map)
-     two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) 
+     two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i)
     enddo
    enddo
   enddo
  enddo
-END_PROVIDER 
+END_PROVIDER
 
 
 subroutine give_mu_of_r_cas(r,istate,mu_of_r,f_psi,n2_psi)
@@ -343,13 +343,13 @@ subroutine give_mu_of_r_cas(r,istate,mu_of_r,f_psi,n2_psi)
  ! active-active part of f_psi(r1,r2)
  call give_f_aa_val_ab(r,r,f_aa_val_ab,two_bod_dens_aa,istate)
 
- f_psi  = f_ii_val_ab     + f_ia_val_ab     + f_aa_val_ab 
+ f_psi  = f_ii_val_ab     + f_ia_val_ab     + f_aa_val_ab
  n2_psi = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa
  if(n2_psi.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * n2_psi.lt.0.d0)then
    w_psi   = 1.d+10
- else 
+ else
    w_psi  = f_psi /  n2_psi
  endif
  mu_of_r  = w_psi * sqpi * 0.5d0
- 
+
 end

src/mu_of_r/mu_of_r_conditions.irp.f

@@ -1,12 +1,12 @@
 
  BEGIN_PROVIDER [double precision, mu_of_r_prov, (n_points_final_grid,N_states) ]
- implicit none 
+ implicit none
  BEGIN_DOC
- ! general variable for mu(r) 
+ ! general variable for mu(r)
  !
- ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) 
+ ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
  !
- ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals 
+ ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
  !
  ! in the two-body density matrix are excluded
  END_DOC
@@ -31,7 +31,7 @@
     mu_of_r_prov(ipoint,istate) =  mu_of_r_hf_sparse(ipoint)
    else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
     mu_of_r_prov(ipoint,istate) =  mu_of_r_psi_cas(ipoint,istate)
-   else 
+   else
     print*,'you requested the following mu_of_r_potential'
     print*,mu_of_r_potential
     print*,'which does not correspond to any of the options for such keyword'
@@ -42,23 +42,23 @@
 
  if (write_mu_of_r) then
    print*,'Writing mu(r) on disk ...'
-   call ezfio_set_mu_of_r_io_mu_of_r('Read')                                                                               
+   call ezfio_set_mu_of_r_io_mu_of_r('Read')
-   call ezfio_set_mu_of_r_mu_of_r_disk(mu_of_r_prov)                                                                               
+   call ezfio_set_mu_of_r_mu_of_r_disk(mu_of_r_prov)
  endif
 
  call wall_time(wall1)
  print*,'Time to provide mu_of_r = ',wall1-wall0
- END_PROVIDER 
+ END_PROVIDER
 
 
  BEGIN_PROVIDER [double precision, mu_of_r_hf, (n_points_final_grid) ]
- implicit none 
+ implicit none
  BEGIN_DOC
  ! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO)
  !
  ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B
  !
- ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals 
+ ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
  !
  ! in the two-body density matrix are excluded
  END_DOC
@@ -70,14 +70,14 @@
  sqpi = dsqrt(dacos(-1.d0))
  !$OMP PARALLEL DO &
  !$OMP DEFAULT (NONE)  &
- !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & 
+ !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) &
- !$OMP ShARED (n_points_final_grid,mu_of_r_hf,f_hf_cholesky,on_top_hf_grid,sqpi) 
+ !$OMP ShARED (n_points_final_grid,mu_of_r_hf,f_hf_cholesky,on_top_hf_grid,sqpi)
  do ipoint = 1, n_points_final_grid
   f_hf   = f_hf_cholesky(ipoint)
   on_top = on_top_hf_grid(ipoint)
   if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then
     w_hf   = 1.d+10
-  else 
+  else
     w_hf  = f_hf /  on_top
   endif
   mu_of_r_hf(ipoint) =  w_hf * sqpi * 0.5d0
@@ -85,16 +85,16 @@
  !$OMP END PARALLEL DO
  call wall_time(wall1)
  print*,'Time to provide mu_of_r_hf = ',wall1-wall0
- END_PROVIDER 
+ END_PROVIDER
 
  BEGIN_PROVIDER [double precision, mu_of_r_hf_sparse, (n_points_final_grid) ]
- implicit none 
+ implicit none
  BEGIN_DOC
  ! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO)
  !
  ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B
  !
- ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals 
+ ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
  !
  ! in the two-body density matrix are excluded
  END_DOC
@@ -106,14 +106,14 @@
  PROVIDE f_hf_cholesky_sparse on_top_hf_grid
  !$OMP PARALLEL DO &
  !$OMP DEFAULT (NONE)  &
- !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & 
+ !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) &
- !$OMP ShARED (n_points_final_grid,mu_of_r_hf_sparse,f_hf_cholesky_sparse,on_top_hf_grid,sqpi) 
+ !$OMP ShARED (n_points_final_grid,mu_of_r_hf_sparse,f_hf_cholesky_sparse,on_top_hf_grid,sqpi)
  do ipoint = 1, n_points_final_grid
   f_hf   = f_hf_cholesky_sparse(ipoint)
   on_top = on_top_hf_grid(ipoint)
   if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then
     w_hf   = 1.d+10
-  else 
+  else
     w_hf  = f_hf /  on_top
   endif
   mu_of_r_hf_sparse(ipoint) =  w_hf * sqpi * 0.5d0
@@ -121,36 +121,36 @@
  !$OMP END PARALLEL DO
  call wall_time(wall1)
  print*,'Time to provide mu_of_r_hf_sparse = ',wall1-wall0
- END_PROVIDER 
+ END_PROVIDER
 
  BEGIN_PROVIDER [double precision, mu_of_r_hf_old, (n_points_final_grid) ]
- implicit none 
+ implicit none
  BEGIN_DOC
  ! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO)
  !
  ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B
  !
- ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals 
+ ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
  !
  ! in the two-body density matrix are excluded
  END_DOC
  integer :: ipoint
  double precision :: wall0,wall1,f_hf,on_top,w_hf,sqpi
- PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
+ PROVIDE all_mo_integrals
  print*,'providing mu_of_r_hf_old ...'
  call wall_time(wall0)
  sqpi = dsqrt(dacos(-1.d0))
- provide f_psi_hf_ab 
+ provide f_psi_hf_ab
  !$OMP PARALLEL DO &
  !$OMP DEFAULT (NONE)  &
- !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & 
+ !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) &
- !$OMP ShARED (n_points_final_grid,mu_of_r_hf_old,f_psi_hf_ab,on_top_hf_mu_r,sqpi) 
+ !$OMP ShARED (n_points_final_grid,mu_of_r_hf_old,f_psi_hf_ab,on_top_hf_mu_r,sqpi)
  do ipoint = 1, n_points_final_grid
   f_hf   = f_psi_hf_ab(ipoint)
   on_top = on_top_hf_mu_r(ipoint)
   if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then
     w_hf   = 1.d+10
-  else 
+  else
     w_hf  = f_hf /  on_top
   endif
   mu_of_r_hf_old(ipoint) =  w_hf * sqpi * 0.5d0
@@ -158,17 +158,17 @@
  !$OMP END PARALLEL DO
  call wall_time(wall1)
  print*,'Time to provide mu_of_r_hf_old = ',wall1-wall0
- END_PROVIDER 
+ END_PROVIDER
 
 
  BEGIN_PROVIDER [double precision, mu_of_r_psi_cas, (n_points_final_grid,N_states) ]
- implicit none 
+ implicit none
  BEGIN_DOC
  ! mu(r) computed with a wave function developped in an active space
  !
- ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) 
+ ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
  !
- ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals 
+ ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
  !
  ! in the one- and two-body density matrix are excluded
  END_DOC
@@ -181,15 +181,15 @@
  provide f_psi_cas_ab
  !$OMP PARALLEL DO &
  !$OMP DEFAULT (NONE)  &
- !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) & 
+ !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) &
- !$OMP SHARED (n_points_final_grid,mu_of_r_psi_cas,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states) 
+ !$OMP SHARED (n_points_final_grid,mu_of_r_psi_cas,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states)
  do istate = 1, N_states
   do ipoint = 1, n_points_final_grid
-   f_psi  = f_psi_cas_ab(ipoint,istate) 
+   f_psi  = f_psi_cas_ab(ipoint,istate)
-   on_top = on_top_cas_mu_r(ipoint,istate) 
+   on_top = on_top_cas_mu_r(ipoint,istate)
    if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then
      w_psi   = 1.d+10
-   else 
+   else
      w_psi  = f_psi /  on_top
    endif
    mu_of_r_psi_cas(ipoint,istate) = w_psi * sqpi * 0.5d0
@@ -198,15 +198,15 @@
  !$OMP END PARALLEL DO
  call wall_time(wall1)
  print*,'Time to provide mu_of_r_psi_cas = ',wall1-wall0
- END_PROVIDER 
+ END_PROVIDER
 
 
  BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)]
  implicit none
  BEGIN_DOC
- ! average value of mu(r) weighted with the total one-e density and divided by the number of electrons 
+ ! average value of mu(r) weighted with the total one-e density and divided by the number of electrons
  !
- ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals 
+ ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
  !
  ! in the one- and two-body density matrix are excluded
  END_DOC
@@ -216,12 +216,12 @@
  do istate = 1, N_states
   do ipoint = 1, n_points_final_grid
    weight =final_weight_at_r_vector(ipoint)
-   density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) & 
+   density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) &
            + one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
    if(mu_of_r_prov(ipoint,istate).gt.1.d+09)cycle
    mu_average_prov(istate) += mu_of_r_prov(ipoint,istate) * weight * density
   enddo
   mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate)
  enddo
- END_PROVIDER 
+ END_PROVIDER
 

src/scf_utils/roothaan_hall_scf.irp.f

@@ -222,13 +222,13 @@ END_DOC
   endif
 
 
-  ! Identify degenerate MOs and force them on the axes
+  ! Identify degenerate MOs and force them to be on the axes
   allocate(S(ao_num,ao_num))
   i=1
   do while (i<mo_num)
-    j=i
+    j=i+1
     m=1
-    do while ( (j+1<mo_num).and.(fock_matrix_diag_mo(j+1)-fock_matrix_diag_mo(i) < 1.d-8) )
+    do while ( (j<=mo_num).and.(fock_matrix_diag_mo(j)-fock_matrix_diag_mo(i) < 1.d-5) )
       j += 1
       m += 1
     enddo
@@ -236,7 +236,7 @@ END_DOC
       call dgemm('N','T',ao_num,ao_num,m,1.d0,mo_coef(1,i),size(mo_coef,1),mo_coef(1,i),size(mo_coef,1),0.d0,S,size(S,1))
       call pivoted_cholesky( S, m, -1.d0, ao_num, mo_coef(1,i))
     endif
-    i = j+1
+    i = j
   enddo
 
   if(do_mom)then

src/tools/fcidump.irp.f

@@ -41,7 +41,7 @@ program fcidump
   integer(key_kind), allocatable :: keys(:)
   double precision, allocatable  :: values(:)
   integer(cache_map_size_kind)   :: n_elements, n_elements_max
-  PROVIDE mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
 
   double precision :: get_two_e_integral, integral
 

src/tools/fcidump_pyscf.irp.f

@@ -41,7 +41,7 @@ program fcidump_pyscf
   integer(key_kind), allocatable :: keys(:)
   double precision, allocatable  :: values(:)
   integer(cache_map_size_kind)   :: n_elements, n_elements_max
-  PROVIDE mo_two_e_integrals_in_map
+  PROVIDE all_mo_integrals
 
   double precision :: get_two_e_integral, integral
 

src/tools/four_idx_transform.irp.f

@@ -18,6 +18,6 @@ program four_idx_transform
   io_mo_two_e_integrals = 'Write'
   SOFT_TOUCH io_mo_two_e_integrals
   if (.true.) then
-    PROVIDE mo_two_e_integrals_in_map
+    PROVIDE all_mo_integrals
   endif
 end

src/trexio/export_trexio_routines.irp.f

@@ -271,11 +271,7 @@ subroutine export_trexio(update,full_path)
      call trexio_assert(rc, TREXIO_SUCCESS)
 
      allocate(factor(shell_num))
-!     if (ao_normalized) then
+     factor(1:shell_num) = shell_normalization_factor(1:shell_num)
-       factor(1:shell_num) = shell_normalization_factor(1:shell_num)
-!     else
-!       factor(1:shell_num) = 1.d0
-!     endif
      rc = trexio_write_basis_shell_factor(f(1), factor)
      call trexio_assert(rc, TREXIO_SUCCESS)
 
@@ -291,11 +287,12 @@ subroutine export_trexio(update,full_path)
     call trexio_assert(rc, TREXIO_SUCCESS)
 
     allocate(factor(prim_num))
-!    if (primitives_normalized) then
+    if (primitives_normalized) then
       factor(1:prim_num) = prim_normalization_factor(1:prim_num)
-!    else
+    else
-!      factor(1:prim_num) = 1.d0
+      factor(1:prim_num) = 1.d0
-!    endif
+    endif
+
     rc = trexio_write_basis_prim_factor(f(1), factor)
     call trexio_assert(rc, TREXIO_SUCCESS)
     deallocate(factor)
@@ -324,14 +321,10 @@ subroutine export_trexio(update,full_path)
     C_A(3) = 0.d0
 
     allocate(factor(ao_num))
-    if (ao_normalized) then
+    do i=1,ao_num
-      do i=1,ao_num
+      l = ao_first_of_shell(ao_shell(i))
-        l = ao_first_of_shell(ao_shell(i))
+      factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0))
-        factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0))
+    enddo
-      enddo
-    else
-      factor(:) = 1.d0
-    endif
     rc = trexio_write_ao_normalization(f(1), factor)
     call trexio_assert(rc, TREXIO_SUCCESS)
     deallocate(factor)
@@ -505,7 +498,7 @@ subroutine export_trexio(update,full_path)
 
   if (export_mo_two_e_ints) then
     print *, 'MO two-e integrals'
-    PROVIDE mo_two_e_integrals_in_map
+    PROVIDE all_mo_integrals
 
     double precision, external :: mo_two_e_integral
 

src/two_body_rdm/act_2_rdm.irp.f

@@ -4,34 +4,34 @@
  BEGIN_DOC
 !                             12 12
 !                 1 2 1 2 == <ij|kl>
-! act_2_rdm_ab_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of alpha/beta+beta/alpha electrons 
+! act_2_rdm_ab_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of alpha/beta+beta/alpha electrons
-! 
+!
 !   <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>
 !
 ! + <Psi_{istate}| a^{\dagger}_{i \beta} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \beta} |Psi_{istate}>
 !
-! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
+! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
 !
 ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2
 !
-! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
+! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
 !
- END_DOC 
+ END_DOC
  integer :: ispin
  double precision :: wall_1, wall_2
- character*(128) :: name_file 
+ character*(128) :: name_file
  name_file = 'act_2_rdm_ab_mo'
  ! condition for alpha/beta spin
  print*,''
  print*,'Providing act_2_rdm_ab_mo '
- ispin = 3 
+ ispin = 3
  act_2_rdm_ab_mo = 0.d0
- provide mo_two_e_integrals_in_map
+ provide all_mo_integrals
  call wall_time(wall_1)
  if(read_two_body_rdm_ab)then
   print*,'Reading act_2_rdm_ab_mo from disk ...'
   call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_ab_mo,name_file)
- else 
+ else
   call orb_range_2_rdm_openmp(act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  endif
  if(write_two_body_rdm_ab)then
@@ -42,36 +42,36 @@
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_ab_mo',wall_2 - wall_1
  act_2_rdm_ab_mo *= 2.d0
- END_PROVIDER 
+ END_PROVIDER
 
 
  BEGIN_PROVIDER [double precision, act_2_rdm_aa_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
  implicit none
  BEGIN_DOC
-! act_2_rdm_aa_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of ALPHA/ALPHA electrons 
+! act_2_rdm_aa_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of ALPHA/ALPHA electrons
-! 
+!
 ! <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \alpha} |Psi_{istate}>
 !
-! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
+! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
 !
 ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)
 !
-! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
+! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
- END_DOC 
+ END_DOC
  integer :: ispin
  double precision :: wall_1, wall_2
  ! condition for alpha/beta spin
  print*,''
  print*,'Providing act_2_rdm_aa_mo '
- character*(128) :: name_file 
+ character*(128) :: name_file
  name_file = 'act_2_rdm_aa_mo'
- ispin = 1 
+ ispin = 1
  act_2_rdm_aa_mo = 0.d0
  call wall_time(wall_1)
  if(read_two_body_rdm_aa)then
   print*,'Reading act_2_rdm_aa_mo from disk ...'
   call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_aa_mo,name_file)
- else 
+ else
   call orb_range_2_rdm_openmp(act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  endif
  if(write_two_body_rdm_aa)then
@@ -83,36 +83,36 @@
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_aa_mo',wall_2 - wall_1
  act_2_rdm_aa_mo *= 2.d0
- END_PROVIDER 
+ END_PROVIDER
- 
+
 
  BEGIN_PROVIDER [double precision, act_2_rdm_bb_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
  implicit none
  BEGIN_DOC
-! act_2_rdm_bb_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of BETA/BETA electrons 
+! act_2_rdm_bb_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of BETA/BETA electrons
-! 
+!
 ! <Psi_{istate}| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi_{istate}>
 !
-! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
+! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
 !
 ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)
 !
-! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
+! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
- END_DOC 
+ END_DOC
  integer :: ispin
  double precision :: wall_1, wall_2
  ! condition for beta/beta spin
  print*,''
  print*,'Providing act_2_rdm_bb_mo '
- character*(128) :: name_file 
+ character*(128) :: name_file
  name_file = 'act_2_rdm_bb_mo'
- ispin = 2 
+ ispin = 2
  act_2_rdm_bb_mo = 0.d0
  call wall_time(wall_1)
  if(read_two_body_rdm_bb)then
   print*,'Reading act_2_rdm_bb_mo from disk ...'
   call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_bb_mo,name_file)
- else 
+ else
   call orb_range_2_rdm_openmp(act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  endif
  if(write_two_body_rdm_bb)then
@@ -124,35 +124,35 @@
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_bb_mo',wall_2 - wall_1
  act_2_rdm_bb_mo *= 2.d0
- END_PROVIDER 
+ END_PROVIDER
 
  BEGIN_PROVIDER [double precision, act_2_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
  implicit none
  BEGIN_DOC
-! act_2_rdm_spin_trace_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM 
+! act_2_rdm_spin_trace_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM
-! 
+!
 ! \sum_{\sigma,\sigma'}<Psi_{istate}| a^{\dagger}_{i \sigma} a^{\dagger}_{j \sigma'} a_{l \sigma'} a_{k \sigma} |Psi_{istate}>
 !
-! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
+! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
 !
 ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec}^{act} * (N_{elec}^{act} - 1)
 !
-! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
+! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
- END_DOC 
+ END_DOC
  integer :: ispin
  double precision :: wall_1, wall_2
  ! condition for beta/beta spin
  print*,''
  print*,'Providing act_2_rdm_spin_trace_mo '
- character*(128) :: name_file 
+ character*(128) :: name_file
  name_file = 'act_2_rdm_spin_trace_mo'
- ispin = 4 
+ ispin = 4
  act_2_rdm_spin_trace_mo = 0.d0
  call wall_time(wall_1)
  if(read_two_body_rdm_spin_trace)then
   print*,'Reading act_2_rdm_spin_trace_mo from disk ...'
   call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_spin_trace_mo,name_file)
- else 
+ else
   call orb_range_2_rdm_openmp(act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  endif
  if(write_two_body_rdm_spin_trace)then
@@ -164,4 +164,4 @@
  act_2_rdm_spin_trace_mo *= 2.d0
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_spin_trace_mo',wall_2 - wall_1
- END_PROVIDER 
+ END_PROVIDER

src/two_rdm_routines/davidson_like_2rdm.irp.f

@@ -2,9 +2,9 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
    use bitmasks
    implicit none
    BEGIN_DOC
-   ! if ispin == 1 :: alpha/alpha 2rdm 
+   ! if ispin == 1 :: alpha/alpha 2rdm
-   !          == 2 :: beta /beta  2rdm 
+   !          == 2 :: beta /beta  2rdm
-   !          == 3 :: alpha/beta + beta/alpha  2rdm 
+   !          == 3 :: alpha/beta + beta/alpha  2rdm
    !          == 4 :: spin traced 2rdm :: aa + bb + ab + ba
    !
    ! Assumes that the determinants are in psi_det
@@ -19,7 +19,7 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
    integer                        :: k
    double precision, allocatable  :: u_t(:,:)
    !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
-   PROVIDE mo_two_e_integrals_in_map
+   PROVIDE all_mo_integrals
    allocate(u_t(N_st,N_det))
    do k=1,N_st
      call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
@@ -30,14 +30,14 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
        u_t,                                                          &
        size(u_t, 1),                                                 &
        N_det, N_st)
-   
+
    call orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,1,N_det,0,1)
    deallocate(u_t)
-   
+
    do k=1,N_st
      call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
    enddo
-   
+
 end
 
 subroutine orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
@@ -52,11 +52,11 @@ subroutine orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_
    integer, intent(in)             :: dim1,norb,list_orb(norb),ispin
    double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
    double precision, intent(in)   :: u_t(N_st,N_det)
-   
+
    integer :: k
-   
+
    PROVIDE N_int
-   
+
    select case (N_int)
      case (1)
        call orb_range_2_rdm_openmp_work_1(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
@@ -70,9 +70,9 @@ subroutine orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_
        call orb_range_2_rdm_openmp_work_N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
    end select
 end
- 
+
- 
+
- 
+
 
  BEGIN_TEMPLATE
 subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
@@ -81,9 +81,9 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
    implicit none
    BEGIN_DOC
    ! Computes the two rdm for the N_st vectors |u_t>
-   ! if ispin == 1 :: alpha/alpha 2rdm 
+   ! if ispin == 1 :: alpha/alpha 2rdm
-   !          == 2 :: beta /beta  2rdm 
+   !          == 2 :: beta /beta  2rdm
-   !          == 3 :: alpha/beta  2rdm 
+   !          == 3 :: alpha/beta  2rdm
    !          == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba))
    ! The 2rdm will be computed only on the list of orbitals list_orb, which contains norb
    ! Default should be 1,N_det,0,1 for istart,iend,ishift,istep
@@ -92,7 +92,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
    double precision, intent(in)   :: u_t(N_st,N_det)
    integer, intent(in)            :: dim1,norb,list_orb(norb),ispin
    double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
-   
+
    integer(omp_lock_kind)          :: lock_2rdm
    integer                        :: i,j,k,l
    integer                        :: k_a, k_b, l_a, l_b
@@ -116,7 +116,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
    integer, allocatable           :: keys(:,:)
    double precision, allocatable  :: values(:,:)
    integer                        :: nkeys,sze_buff
-   alpha_alpha = .False.   
+   alpha_alpha = .False.
    beta_beta   = .False.
    alpha_beta  = .False.
    spin_trace  = .False.
@@ -134,27 +134,27 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
     stop
    endif
 
-   
+
    PROVIDE N_int
 
    call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)
-   sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 
+   sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60
-   list_orb_reverse = -1000 
+   list_orb_reverse = -1000
    do i = 1, norb
-    list_orb_reverse(list_orb(i)) = i 
+    list_orb_reverse(list_orb(i)) = i
    enddo
    maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
    allocate(idx0(maxab))
-   
+
    do i=1,maxab
      idx0(i) = i
    enddo
    call omp_init_lock(lock_2rdm)
-   
+
    ! Prepare the array of all alpha single excitations
    ! -------------------------------------------------
-   
+
-   PROVIDE N_int nthreads_davidson elec_alpha_num 
+   PROVIDE N_int nthreads_davidson elec_alpha_num
     !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson)      &
             !$OMP   SHARED(psi_bilinear_matrix_rows, N_det,lock_2rdm,&
             !$OMP          psi_bilinear_matrix_columns,              &
@@ -166,7 +166,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
             !$OMP          psi_bilinear_matrix_order_transp_reverse, &
             !$OMP          psi_bilinear_matrix_columns_loc,          &
             !$OMP          psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, &
-            !$OMP          istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, & 
+            !$OMP          istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, &
             !$OMP          ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask)     &
             !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1,  &
             !$OMP          lcol, lrow, l_a, l_b,                     &
@@ -174,35 +174,35 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
             !$OMP          tmp_det2, idx, l, kcol_prev,              &
             !$OMP          singles_a, n_singles_a, singles_b,        &
             !$OMP          n_singles_b, nkeys, keys, values)
-   
+
    ! Alpha/Beta double excitations
    ! =============================
    nkeys = 0
-   allocate( keys(4,sze_buff), values(n_st,sze_buff)) 
+   allocate( keys(4,sze_buff), values(n_st,sze_buff))
    allocate( buffer($N_int,maxab),                                   &
        singles_a(maxab),                                             &
        singles_b(maxab),                                             &
        doubles(maxab),                                               &
        idx(maxab))
-   
+
    kcol_prev=-1
-   
+
    ASSERT (iend <= N_det)
    ASSERT (istart > 0)
    ASSERT (istep  > 0)
-   
+
     !$OMP DO SCHEDULE(dynamic,64)
    do k_a=istart+ishift,iend,istep
-     
+
      krow = psi_bilinear_matrix_rows(k_a)
      ASSERT (krow <= N_det_alpha_unique)
-     
+
      kcol = psi_bilinear_matrix_columns(k_a)
      ASSERT (kcol <= N_det_beta_unique)
-     
+
      tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-     
+
      if (kcol /= kcol_prev) then
        call get_all_spin_singles_$N_int(                             &
            psi_det_beta_unique, idx0,                                &
@@ -210,58 +210,58 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
            singles_b, n_singles_b)
      endif
      kcol_prev = kcol
-     
+
      ! Loop over singly excited beta columns
      ! -------------------------------------
-     
+
      do i=1,n_singles_b
        lcol = singles_b(i)
-       
+
        tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
-       
+
        l_a = psi_bilinear_matrix_columns_loc(lcol)
        ASSERT (l_a <= N_det)
-       
+
        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
          lrow = psi_bilinear_matrix_rows(l_a)
          ASSERT (lrow <= N_det_alpha_unique)
-         
+
          buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
-         
+
          ASSERT (l_a <= N_det)
          idx(j) = l_a
          l_a = l_a+1
        enddo
        j = j-1
-       
+
        call get_all_spin_singles_$N_int(                             &
            buffer, idx, tmp_det(1,1), j,                             &
            singles_a, n_singles_a )
-       
+
        ! Loop over alpha singles
        ! -----------------------
-       
+
       if(alpha_beta.or.spin_trace)then
        do k = 1,n_singles_a
          l_a = singles_a(k)
          ASSERT (l_a <= N_det)
-         
+
          lrow = psi_bilinear_matrix_rows(l_a)
          ASSERT (lrow <= N_det_alpha_unique)
-         
+
          tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
 !         print*,'nkeys before = ',nkeys
          do l= 1, N_states
            c_1(l) = u_t(l,l_a) * u_t(l,k_a)
          enddo
          if(alpha_beta)then
-          ! only ONE contribution 
+          ! only ONE contribution
           if (nkeys+1 .ge. sze_buff) then
             call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
             nkeys = 0
           endif
          else if (spin_trace)then
-          ! TWO contributions 
+          ! TWO contributions
           if (nkeys+2 .ge. sze_buff) then
             call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
             nkeys = 0
@@ -271,42 +271,42 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
 
        enddo
       endif
-       
+
       call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
       nkeys = 0
      enddo
-     
+
    enddo
      !$OMP END DO
-   
+
      !$OMP DO SCHEDULE(dynamic,64)
    do k_a=istart+ishift,iend,istep
-     
+
-     
+
      ! Single and double alpha exitations
      ! ===================================
-     
+
-     
+
      ! Initial determinant is at k_a in alpha-major representation
      ! -----------------------------------------------------------------------
-     
+
      krow = psi_bilinear_matrix_rows(k_a)
      ASSERT (krow <= N_det_alpha_unique)
-     
+
      kcol = psi_bilinear_matrix_columns(k_a)
      ASSERT (kcol <= N_det_beta_unique)
-     
+
      tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-     
+
      ! Initial determinant is at k_b in beta-major representation
      ! ----------------------------------------------------------------------
-     
+
      k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
      ASSERT (k_b <= N_det)
-     
+
      spindet(1:$N_int) = tmp_det(1:$N_int,1)
-     
+
      ! Loop inside the beta column to gather all the connected alphas
      lcol = psi_bilinear_matrix_columns(k_a)
      l_a = psi_bilinear_matrix_columns_loc(lcol)
@@ -316,28 +316,28 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
        if (lcol /= kcol) exit
        lrow = psi_bilinear_matrix_rows(l_a)
        ASSERT (lrow <= N_det_alpha_unique)
-       
+
        buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
        idx(i) = l_a
        l_a = l_a+1
      enddo
      i = i-1
-     
+
      call get_all_spin_singles_and_doubles_$N_int(                   &
          buffer, idx, spindet, i,                                    &
          singles_a, doubles, n_singles_a, n_doubles )
-     
+
      ! Compute Hij for all alpha singles
      ! ----------------------------------
-     
+
      tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
      do i=1,n_singles_a
        l_a = singles_a(i)
        ASSERT (l_a <= N_det)
-       
+
        lrow = psi_bilinear_matrix_rows(l_a)
        ASSERT (lrow <= N_det_alpha_unique)
-       
+
        tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
        do l= 1, N_states
          c_1(l) = u_t(l,l_a) * u_t(l,k_a)
@@ -356,23 +356,23 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
        endif
        call orb_range_off_diag_single_to_all_states_aa_dm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
        endif
-       
+
      enddo
-     
+
      call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
      nkeys = 0
-     
+
      ! Compute Hij for all alpha doubles
      ! ----------------------------------
-     
+
      if(alpha_alpha.or.spin_trace)then
       do i=1,n_doubles
         l_a = doubles(i)
         ASSERT (l_a <= N_det)
-        
+
         lrow = psi_bilinear_matrix_rows(l_a)
         ASSERT (lrow <= N_det_alpha_unique)
-        
+
         do l= 1, N_states
           c_1(l) = u_t(l,l_a) * u_t(l,k_a)
         enddo
@@ -385,29 +385,29 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
      endif
      call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
      nkeys = 0
-     
+
-     
+
      ! Single and double beta excitations
      ! ==================================
-     
+
-     
+
      ! Initial determinant is at k_a in alpha-major representation
      ! -----------------------------------------------------------------------
-     
+
      krow = psi_bilinear_matrix_rows(k_a)
      kcol = psi_bilinear_matrix_columns(k_a)
-     
+
      tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-     
+
      spindet(1:$N_int) = tmp_det(1:$N_int,2)
-     
+
      ! Initial determinant is at k_b in beta-major representation
      ! -----------------------------------------------------------------------
-     
+
      k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
      ASSERT (k_b <= N_det)
-     
+
      ! Loop inside the alpha row to gather all the connected betas
      lrow = psi_bilinear_matrix_transp_rows(k_b)
      l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
@@ -417,28 +417,28 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
        if (lrow /= krow) exit
        lcol = psi_bilinear_matrix_transp_columns(l_b)
        ASSERT (lcol <= N_det_beta_unique)
-       
+
        buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
        idx(i) = l_b
        l_b = l_b+1
      enddo
      i = i-1
-     
+
      call get_all_spin_singles_and_doubles_$N_int(                   &
          buffer, idx, spindet, i,                                    &
          singles_b, doubles, n_singles_b, n_doubles )
-     
+
      ! Compute Hij for all beta singles
      ! ----------------------------------
-     
+
      tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      do i=1,n_singles_b
        l_b = singles_b(i)
        ASSERT (l_b <= N_det)
-       
+
        lcol = psi_bilinear_matrix_transp_columns(l_b)
        ASSERT (lcol <= N_det_beta_unique)
-       
+
        tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
        l_a = psi_bilinear_matrix_transp_order(l_b)
        do l= 1, N_states
@@ -461,18 +461,18 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
      enddo
      call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
      nkeys = 0
-     
+
      ! Compute Hij for all beta doubles
      ! ----------------------------------
-     
+
      if(beta_beta.or.spin_trace)then
       do i=1,n_doubles
         l_b = doubles(i)
         ASSERT (l_b <= N_det)
-        
+
         lcol = psi_bilinear_matrix_transp_columns(l_b)
         ASSERT (lcol <= N_det_beta_unique)
-        
+
         l_a = psi_bilinear_matrix_transp_order(l_b)
         do l= 1, N_states
           c_1(l) = u_t(l,l_a) * u_t(l,k_a)
@@ -484,59 +484,59 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
         call orb_range_off_diag_double_to_all_states_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
 !        print*,'to do orb_range_off_diag_double_to_2_rdm_bb_dm_buffer'
         ASSERT (l_a <= N_det)
-        
+
       enddo
      endif
      call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
      nkeys = 0
-     
+
-     
+
      ! Diagonal contribution
      ! =====================
-     
+
-     
+
      ! Initial determinant is at k_a in alpha-major representation
      ! -----------------------------------------------------------------------
-     
+
      krow = psi_bilinear_matrix_rows(k_a)
      ASSERT (krow <= N_det_alpha_unique)
-     
+
      kcol = psi_bilinear_matrix_columns(k_a)
      ASSERT (kcol <= N_det_beta_unique)
-     
+
      tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-     
+
      double precision, external     :: diag_wee_mat_elem, diag_S_mat_elem
-     
+
      double precision               :: c_1(N_states)
      do l = 1, N_states
        c_1(l) = u_t(l,k_a) * u_t(l,k_a)
      enddo
-     
+
      call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
      nkeys = 0
      call orb_range_diag_to_all_states_2_rdm_dm_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
      call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
      nkeys = 0
-     
+
    end do
     !$OMP END DO
    deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values)
     !$OMP END PARALLEL
-   
+
 end
- 
+
  SUBST [ N_int ]
- 
+
  1;;
  2;;
  3;;
  4;;
  N_int;;
- 
+
  END_TEMPLATE
- 
+
 
 subroutine update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
  use omp_lib
@@ -545,7 +545,7 @@ subroutine update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,loc
  integer, intent(in) :: keys(4,nkeys)
  double precision, intent(in) :: values(n_st,nkeys)
  double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,n_st)
- 
+
  integer(omp_lock_kind),intent(inout):: lock_2rdm
 
  integer :: istate

src/two_rdm_routines/davidson_like_trans_2rdm.irp.f

@@ -2,12 +2,12 @@ subroutine orb_range_2_trans_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N
    use bitmasks
    implicit none
    BEGIN_DOC
-   ! if ispin == 1 :: alpha/alpha 2_rdm 
+   ! if ispin == 1 :: alpha/alpha 2_rdm
-   !          == 2 :: beta /beta  2_rdm 
+   !          == 2 :: beta /beta  2_rdm
-   !          == 3 :: alpha/beta + beta/alpha  2trans_rdm 
+   !          == 3 :: alpha/beta + beta/alpha  2trans_rdm
    !          == 4 :: spin traced 2_rdm :: aa + bb + ab + ba
    !
-   ! notice that here it is the TRANSITION RDM THAT IS COMPUTED 
+   ! notice that here it is the TRANSITION RDM THAT IS COMPUTED
    !
    ! THE DIAGONAL PART IS THE USUAL ONE FOR A GIVEN STATE
    ! Assumes that the determinants are in psi_det
@@ -22,7 +22,7 @@ subroutine orb_range_2_trans_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N
    integer                        :: k
    double precision, allocatable  :: u_t(:,:)
    !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
-   PROVIDE mo_two_e_integrals_in_map
+   PROVIDE all_mo_integrals
    allocate(u_t(N_st,N_det))
    do k=1,N_st
      call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
@@ -33,14 +33,14 @@ subroutine orb_range_2_trans_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N
        u_t,                                                          &
        size(u_t, 1),                                                 &
        N_det, N_st)
-   
+
    call orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,1,N_det,0,1)
    deallocate(u_t)
-   
+
    do k=1,N_st
      call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
    enddo
-   
+
 end
 
 subroutine orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
@@ -55,11 +55,11 @@ subroutine orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,
    integer, intent(in)             :: dim1,norb,list_orb(norb),ispin
    double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st)
    double precision, intent(in)    :: u_t(N_st,N_det)
-   
+
    integer :: k
-   
+
    PROVIDE N_int
-   
+
    select case (N_int)
      case (1)
        call orb_range_2_trans_rdm_openmp_work_1(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
@@ -73,8 +73,8 @@ subroutine orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,
        call orb_range_2_trans_rdm_openmp_work_N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
    end select
 end
- 
+
- 
+
  BEGIN_TEMPLATE
 subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
    use bitmasks
@@ -82,9 +82,9 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
    implicit none
    BEGIN_DOC
    ! Computes the two trans_rdm for the N_st vectors |u_t>
-   ! if ispin == 1 :: alpha/alpha 2trans_rdm 
+   ! if ispin == 1 :: alpha/alpha 2trans_rdm
-   !          == 2 :: beta /beta  2trans_rdm 
+   !          == 2 :: beta /beta  2trans_rdm
-   !          == 3 :: alpha/beta  2trans_rdm 
+   !          == 3 :: alpha/beta  2trans_rdm
    !          == 4 :: spin traced 2trans_rdm :: aa + bb + 0.5 (ab + ba))
    ! The 2trans_rdm will be computed only on the list of orbitals list_orb, which contains norb
    ! Default should be 1,N_det,0,1 for istart,iend,ishift,istep
@@ -93,7 +93,7 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
    double precision, intent(in)   :: u_t(N_st,N_det)
    integer, intent(in)            :: dim1,norb,list_orb(norb),ispin
    double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st)
-   
+
    integer(omp_lock_kind)          :: lock_2trans_rdm
    integer                        :: i,j,k,l
    integer                        :: k_a, k_b, l_a, l_b
@@ -118,7 +118,7 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
    double precision, allocatable  :: values(:,:,:)
    integer                        :: nkeys,sze_buff
    integer :: ll
-   alpha_alpha = .False.   
+   alpha_alpha = .False.
    beta_beta   = .False.
    alpha_beta  = .False.
    spin_trace  = .False.
@@ -136,27 +136,27 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
     stop
    endif
 
-   
+
    PROVIDE N_int
 
    call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)
-   sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 
+   sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60
-   list_orb_reverse = -1000 
+   list_orb_reverse = -1000
    do i = 1, norb
-    list_orb_reverse(list_orb(i)) = i 
+    list_orb_reverse(list_orb(i)) = i
    enddo
    maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
    allocate(idx0(maxab))
-   
+
    do i=1,maxab
      idx0(i) = i
    enddo
    call omp_init_lock(lock_2trans_rdm)
-   
+
    ! Prepare the array of all alpha single excitations
    ! -------------------------------------------------
-   
+
-   PROVIDE N_int nthreads_davidson elec_alpha_num 
+   PROVIDE N_int nthreads_davidson elec_alpha_num
     !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson)      &
             !$OMP   SHARED(psi_bilinear_matrix_rows, N_det,lock_2trans_rdm,&
             !$OMP          psi_bilinear_matrix_columns,              &
@@ -168,7 +168,7 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
             !$OMP          psi_bilinear_matrix_order_transp_reverse, &
             !$OMP          psi_bilinear_matrix_columns_loc,          &
             !$OMP          psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, &
-            !$OMP          istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, & 
+            !$OMP          istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, &
             !$OMP          ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask)     &
             !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1,  &
             !$OMP          lcol, lrow, l_a, l_b,                     &
@@ -176,35 +176,35 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
             !$OMP          tmp_det2, idx, l, kcol_prev,              &
             !$OMP          singles_a, n_singles_a, singles_b,        &
             !$OMP          n_singles_b, nkeys, keys, values)
-   
+
    ! Alpha/Beta double excitations
    ! =============================
    nkeys = 0
-   allocate( keys(4,sze_buff), values(n_st,n_st,sze_buff)) 
+   allocate( keys(4,sze_buff), values(n_st,n_st,sze_buff))
    allocate( buffer($N_int,maxab),                                   &
        singles_a(maxab),                                             &
        singles_b(maxab),                                             &
        doubles(maxab),                                               &
        idx(maxab))
-   
+
    kcol_prev=-1
-   
+
    ASSERT (iend <= N_det)
    ASSERT (istart > 0)
    ASSERT (istep  > 0)
-   
+
     !$OMP DO SCHEDULE(dynamic,64)
    do k_a=istart+ishift,iend,istep
-     
+
      krow = psi_bilinear_matrix_rows(k_a)
      ASSERT (krow <= N_det_alpha_unique)
-     
+
      kcol = psi_bilinear_matrix_columns(k_a)
      ASSERT (kcol <= N_det_beta_unique)
-     
+
      tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-     
+
      if (kcol /= kcol_prev) then
        call get_all_spin_singles_$N_int(                             &
            psi_det_beta_unique, idx0,                                &
@@ -212,45 +212,45 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
            singles_b, n_singles_b)
      endif
      kcol_prev = kcol
-     
+
      ! Loop over singly excited beta columns
      ! -------------------------------------
-     
+
      do i=1,n_singles_b
        lcol = singles_b(i)
-       
+
        tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
-       
+
        l_a = psi_bilinear_matrix_columns_loc(lcol)
        ASSERT (l_a <= N_det)
-       
+
        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
          lrow = psi_bilinear_matrix_rows(l_a)
          ASSERT (lrow <= N_det_alpha_unique)
-         
+
          buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
-         
+
          ASSERT (l_a <= N_det)
          idx(j) = l_a
          l_a = l_a+1
        enddo
        j = j-1
-       
+
        call get_all_spin_singles_$N_int(                             &
            buffer, idx, tmp_det(1,1), j,                             &
            singles_a, n_singles_a )
-       
+
        ! Loop over alpha singles
        ! -----------------------
-       
+
       if(alpha_beta.or.spin_trace)then
        do k = 1,n_singles_a
          l_a = singles_a(k)
          ASSERT (l_a <= N_det)
-         
+
          lrow = psi_bilinear_matrix_rows(l_a)
          ASSERT (lrow <= N_det_alpha_unique)
-         
+
          tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
 !         print*,'nkeys before = ',nkeys
          do ll = 1, N_states
@@ -259,13 +259,13 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
           enddo
          enddo
          if(alpha_beta)then
-          ! only ONE contribution 
+          ! only ONE contribution
           if (nkeys+1 .ge. sze_buff) then
             call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
             nkeys = 0
           endif
          else if (spin_trace)then
-          ! TWO contributions 
+          ! TWO contributions
           if (nkeys+2 .ge. sze_buff) then
             call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
             nkeys = 0
@@ -275,42 +275,42 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
 
        enddo
       endif
-       
+
       call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
       nkeys = 0
      enddo
-     
+
    enddo
      !$OMP END DO
-   
+
      !$OMP DO SCHEDULE(dynamic,64)
    do k_a=istart+ishift,iend,istep
-     
+
-     
+
      ! Single and double alpha exitations
      ! ===================================
-     
+
-     
+
      ! Initial determinant is at k_a in alpha-major representation
      ! -----------------------------------------------------------------------
-     
+
      krow = psi_bilinear_matrix_rows(k_a)
      ASSERT (krow <= N_det_alpha_unique)
-     
+
      kcol = psi_bilinear_matrix_columns(k_a)
      ASSERT (kcol <= N_det_beta_unique)
-     
+
      tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-     
+
      ! Initial determinant is at k_b in beta-major representation
      ! ----------------------------------------------------------------------
-     
+
      k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
      ASSERT (k_b <= N_det)
-     
+
      spindet(1:$N_int) = tmp_det(1:$N_int,1)
-     
+
      ! Loop inside the beta column to gather all the connected alphas
      lcol = psi_bilinear_matrix_columns(k_a)
      l_a = psi_bilinear_matrix_columns_loc(lcol)
@@ -320,28 +320,28 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
        if (lcol /= kcol) exit
        lrow = psi_bilinear_matrix_rows(l_a)
        ASSERT (lrow <= N_det_alpha_unique)
-       
+
        buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
        idx(i) = l_a
        l_a = l_a+1
      enddo
      i = i-1
-     
+
      call get_all_spin_singles_and_doubles_$N_int(                   &
          buffer, idx, spindet, i,                                    &
          singles_a, doubles, n_singles_a, n_doubles )
-     
+
      ! Compute Hij for all alpha singles
      ! ----------------------------------
-     
+
      tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
      do i=1,n_singles_a
        l_a = singles_a(i)
        ASSERT (l_a <= N_det)
-       
+
        lrow = psi_bilinear_matrix_rows(l_a)
        ASSERT (lrow <= N_det_alpha_unique)
-       
+
        tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
        do ll= 1, N_states
         do l= 1, N_states
@@ -362,23 +362,23 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
        endif
        call orb_range_off_diag_single_to_all_states_aa_trans_rdm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
        endif
-       
+
      enddo
-     
+
      call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
      nkeys = 0
-     
+
      ! Compute Hij for all alpha doubles
      ! ----------------------------------
-     
+
      if(alpha_alpha.or.spin_trace)then
       do i=1,n_doubles
         l_a = doubles(i)
         ASSERT (l_a <= N_det)
-        
+
         lrow = psi_bilinear_matrix_rows(l_a)
         ASSERT (lrow <= N_det_alpha_unique)
-        
+
         do ll= 1, N_states
          do l= 1, N_states
            c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a)
@@ -393,29 +393,29 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
      endif
      call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
      nkeys = 0
-     
+
-     
+
      ! Single and double beta excitations
      ! ==================================
-     
+
-     
+
      ! Initial determinant is at k_a in alpha-major representation
      ! -----------------------------------------------------------------------
-     
+
      krow = psi_bilinear_matrix_rows(k_a)
      kcol = psi_bilinear_matrix_columns(k_a)
-     
+
      tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-     
+
      spindet(1:$N_int) = tmp_det(1:$N_int,2)
-     
+
      ! Initial determinant is at k_b in beta-major representation
      ! -----------------------------------------------------------------------
-     
+
      k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
      ASSERT (k_b <= N_det)
-     
+
      ! Loop inside the alpha row to gather all the connected betas
      lrow = psi_bilinear_matrix_transp_rows(k_b)
      l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
@@ -425,28 +425,28 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
        if (lrow /= krow) exit
        lcol = psi_bilinear_matrix_transp_columns(l_b)
        ASSERT (lcol <= N_det_beta_unique)
-       
+
        buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
        idx(i) = l_b
        l_b = l_b+1
      enddo
      i = i-1
-     
+
      call get_all_spin_singles_and_doubles_$N_int(                   &
          buffer, idx, spindet, i,                                    &
          singles_b, doubles, n_singles_b, n_doubles )
-     
+
      ! Compute Hij for all beta singles
      ! ----------------------------------
-     
+
      tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      do i=1,n_singles_b
        l_b = singles_b(i)
        ASSERT (l_b <= N_det)
-       
+
        lcol = psi_bilinear_matrix_transp_columns(l_b)
        ASSERT (lcol <= N_det_beta_unique)
-       
+
        tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
        l_a = psi_bilinear_matrix_transp_order(l_b)
        do ll= 1, N_states
@@ -471,18 +471,18 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
      enddo
      call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
      nkeys = 0
-     
+
      ! Compute Hij for all beta doubles
      ! ----------------------------------
-     
+
      if(beta_beta.or.spin_trace)then
       do i=1,n_doubles
         l_b = doubles(i)
         ASSERT (l_b <= N_det)
-        
+
         lcol = psi_bilinear_matrix_transp_columns(l_b)
         ASSERT (lcol <= N_det_beta_unique)
-        
+
         l_a = psi_bilinear_matrix_transp_order(l_b)
         do ll= 1, N_states
          do l= 1, N_states
@@ -496,61 +496,61 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb
         call orb_range_off_diag_double_to_all_states_trans_rdm_bb_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
 !        print*,'to do orb_range_off_diag_double_to_2_trans_rdm_bb_dm_buffer'
         ASSERT (l_a <= N_det)
-        
+
       enddo
      endif
      call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
      nkeys = 0
-     
+
-     
+
      ! Diagonal contribution
      ! =====================
-     
+
-     
+
      ! Initial determinant is at k_a in alpha-major representation
      ! -----------------------------------------------------------------------
-     
+
      krow = psi_bilinear_matrix_rows(k_a)
      ASSERT (krow <= N_det_alpha_unique)
-     
+
      kcol = psi_bilinear_matrix_columns(k_a)
      ASSERT (kcol <= N_det_beta_unique)
-     
+
      tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
      tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-     
+
      double precision, external     :: diag_wee_mat_elem, diag_S_mat_elem
-     
+
      double precision               :: c_1(N_states,N_states)
      do ll = 1, N_states
       do l = 1, N_states
         c_1(l,ll) = u_t(ll,k_a) * u_t(l,k_a)
       enddo
      enddo
-     
+
      call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
      nkeys = 0
      call orb_range_diag_to_all_states_2_rdm_trans_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
      call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
      nkeys = 0
-     
+
    end do
     !$OMP END DO
    deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values)
     !$OMP END PARALLEL
-   
+
 end
- 
+
  SUBST [ N_int ]
- 
+
  1;;
  2;;
  3;;
  4;;
  N_int;;
- 
+
  END_TEMPLATE
- 
+
 subroutine update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
  use omp_lib
  implicit none
@@ -558,7 +558,7 @@ subroutine update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_arr
  integer, intent(in) :: keys(4,nkeys)
  double precision, intent(in) :: values(n_st,n_st,nkeys)
  double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,n_st,n_st)
- 
+
  integer(omp_lock_kind),intent(inout):: lock_2rdm
 
  integer :: i,h1,h2,p1,p2,istate,jstate

src/utils_cc/mo_integrals_cc.irp.f

@@ -77,7 +77,7 @@ subroutine gen_v_space(n1,n2,n3,n4,list1,list2,list3,list4,v)
   else
     double precision              :: get_two_e_integral
 
-    PROVIDE mo_two_e_integrals_in_map
+    PROVIDE all_mo_integrals
 
     !$OMP PARALLEL &
     !$OMP SHARED(n1,n2,n3,n4,list1,list2,list3,list4,v,mo_integrals_map) &
@@ -161,7 +161,7 @@ BEGIN_PROVIDER [double precision, cc_space_v, (mo_num,mo_num,mo_num,mo_num)]
     integer          :: i,j,k,l
     double precision :: get_two_e_integral
 
-    PROVIDE mo_two_e_integrals_in_map
+    PROVIDE all_mo_integrals
 
     !$OMP PARALLEL &
     !$OMP SHARED(cc_space_v,mo_num,mo_integrals_map) &
@@ -194,7 +194,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_oooo, (cc_nOa, cc_nOa, cc_nOa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_oooo,1)
     n2 = size(cc_space_v_oooo,2)
     n3 = size(cc_space_v_oooo,3)
@@ -237,7 +237,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_vooo, (cc_nVa, cc_nOa, cc_nOa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_vooo,1)
     n2 = size(cc_space_v_vooo,2)
     n3 = size(cc_space_v_vooo,3)
@@ -281,7 +281,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_ovoo, (cc_nOa, cc_nVa, cc_nOa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_ovoo,1)
     n2 = size(cc_space_v_ovoo,2)
     n3 = size(cc_space_v_ovoo,3)
@@ -315,7 +315,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_oovo, (cc_nOa, cc_nOa, cc_nVa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_oovo,1)
     n2 = size(cc_space_v_oovo,2)
     n3 = size(cc_space_v_oovo,3)
@@ -349,7 +349,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_ooov, (cc_nOa, cc_nOa, cc_nOa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_oovo,1)
     n2 = size(cc_space_v_oovo,2)
     n3 = size(cc_space_v_oovo,3)
@@ -383,7 +383,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_vvoo, (cc_nVa, cc_nVa, cc_nOa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_vvoo,1)
     n2 = size(cc_space_v_vvoo,2)
     n3 = size(cc_space_v_vvoo,3)
@@ -426,7 +426,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_vovo, (cc_nVa, cc_nOa, cc_nVa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_vovo,1)
     n2 = size(cc_space_v_vovo,2)
     n3 = size(cc_space_v_vovo,3)
@@ -469,7 +469,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_voov, (cc_nVa, cc_nOa, cc_nOa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_voov,1)
     n2 = size(cc_space_v_voov,2)
     n3 = size(cc_space_v_voov,3)
@@ -503,7 +503,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_ovvo, (cc_nOa, cc_nVa, cc_nVa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_ovvo,1)
     n2 = size(cc_space_v_ovvo,2)
     n3 = size(cc_space_v_ovvo,3)
@@ -537,7 +537,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_ovov, (cc_nOa, cc_nVa, cc_nOa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_ovov,1)
     n2 = size(cc_space_v_ovov,2)
     n3 = size(cc_space_v_ovov,3)
@@ -571,7 +571,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_oovv, (cc_nOa, cc_nOa, cc_nVa, cc_n
 
     integer :: i1, i2, i3, i4
     integer :: n1, n2, n3, n4
-    
+
     n1 = size(cc_space_v_oovv,1)
     n2 = size(cc_space_v_oovv,2)
     n3 = size(cc_space_v_oovv,3)