Handling of two different mo and ao integrals map

plugins/DFT_Utils/grid_density.irp.f

@@ -5,7 +5,7 @@
 
 BEGIN_PROVIDER [integer, n_points_radial_grid]
  implicit none
- n_points_radial_grid = 10000
+ n_points_radial_grid = 100
 END_PROVIDER 
 
 
@@ -188,6 +188,7 @@ END_PROVIDER
      call give_all_mos_at_r(r,mos_array)
      do m = 1, mo_tot_num
       do i = 1, mo_tot_num
+       if(dabs(one_body_dm_mo_alpha(i,m,i_state)).lt.1.d-10)cycle
        contrib = mos_array(i) * mos_array(m)
        one_body_dm_mo_alpha_at_grid_points(l,k,j,i_state) += one_body_dm_mo_alpha(i,m,i_state) * contrib
        one_body_dm_mo_beta_at_grid_points(l,k,j,i_state) += one_body_dm_mo_beta(i,m,i_state) * contrib

plugins/Integrals_erf/EZFIO.cfg

@@ -1,37 +1,3 @@
-[do_direct_integrals]
-type: logical
-doc: Compute integrals on the fly
-interface: ezfio,provider,ocaml
-default: False
-ezfio_name: direct
-
-[no_vvvv_integrals]
-type: logical
-doc: If True, computes all integrals except for the integrals having 4 virtual index
-interface: ezfio,provider,ocaml
-default: False
-ezfio_name: no_vvvv_integrals
-
-[no_ivvv_integrals]
-type: logical
-doc: Can be switched on only if  no_vvvv_integrals  is True, then do not computes the integrals having 3 virtual index and 1 belonging to the core inactive active orbitals
-interface: ezfio,provider,ocaml
-default: False
-ezfio_name: no_ivvv_integrals
-
-[no_vvv_integrals]
-type: logical
-doc: Can be switched on only if  no_vvvv_integrals  is True, then do not computes the integrals having 3 virtual orbitals
-interface: ezfio,provider,ocaml
-default: False
-ezfio_name: no_vvv_integrals
-
-[disk_access_mo_integrals]
-type: Disk_access
-doc: Read/Write MO integrals from/to disk [ Write | Read | None ] 
-interface: ezfio,provider,ocaml
-default: None
-
 [disk_access_ao_integrals_erf]
 type: Disk_access
 doc: Read/Write AO integrals with the long range interaction from/to disk [ Write | Read | None ] 
@@ -45,14 +11,6 @@ doc: Read/Write MO integrals with the long range interaction from/to disk [ Writ
 interface: ezfio,provider,ocaml
 default: None
 
-[disk_access_ao_integrals]
-type: Disk_access
-doc: Read/Write AO integrals from/to disk [ Write | Read | None ] 
-interface: ezfio,provider,ocaml
-default: None
-
-
-
 [ao_integrals_threshold]
 type: Threshold
 doc: If |<pq|rs>| < ao_integrals_threshold then <pq|rs> is zero

plugins/Integrals_erf/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Pseudo Bitmask ZMQ
+Pseudo Bitmask ZMQ Integrals_Bielec

plugins/Integrals_erf/ao_bi_integrals_erf.irp.f

@@ -294,7 +294,7 @@ subroutine compute_ao_bielec_integrals_erf(j,k,l,sze,buffer_value)
     buffer_value = 0._integral_kind
     return
   endif
-  if (ao_bielec_integral_schwartz(j,l) < thresh ) then
+  if (ao_bielec_integral_erf_schwartz(j,l) < thresh ) then
     buffer_value = 0._integral_kind
     return
   endif

plugins/Integrals_erf/ao_bielec_integrals_in_map_slave.irp.f

@@ -1,225 +0,0 @@
-subroutine ao_bielec_integrals_in_map_slave_tcp(i)
-  implicit none
-  integer, intent(in)            :: i
-  BEGIN_DOC
-! Computes a buffer of integrals. i is the ID of the current thread.
-  END_DOC
-  call ao_bielec_integrals_in_map_slave(0,i)
-end
-
-
-subroutine ao_bielec_integrals_in_map_slave_inproc(i)
-  implicit none
-  integer, intent(in)            :: i
-  BEGIN_DOC
-! Computes a buffer of integrals. i is the ID of the current thread.
-  END_DOC
-  call ao_bielec_integrals_in_map_slave(1,i)
-end
-
-
-subroutine push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id)
-  use f77_zmq
-  use map_module
-  implicit none
-  BEGIN_DOC
-! Push integrals in the push socket
-  END_DOC
-  integer(ZMQ_PTR), intent(in)   :: zmq_socket_push
-  integer, intent(in)            :: n_integrals
-  integer(key_kind), intent(in)  :: buffer_i(*)
-  real(integral_kind), intent(in) :: buffer_value(*)
-  integer, intent(in)            :: task_id
-  integer                        :: rc
-  
-  rc = f77_zmq_send( zmq_socket_push, n_integrals, 4, ZMQ_SNDMORE)
-  if (rc /= 4) then
-    print *, irp_here,  ': f77_zmq_send( zmq_socket_push, n_integrals, 4, ZMQ_SNDMORE)'
-    stop 'error'
-  endif
-  
-  rc = f77_zmq_send( zmq_socket_push, buffer_i, key_kind*n_integrals, ZMQ_SNDMORE)
-  if (rc /= key_kind*n_integrals) then
-    print *, irp_here,  ': f77_zmq_send( zmq_socket_push, buffer_i, key_kind*n_integrals, ZMQ_SNDMORE)'
-    stop 'error'
-  endif
-  
-  rc = f77_zmq_send( zmq_socket_push, buffer_value, integral_kind*n_integrals, ZMQ_SNDMORE)
-  if (rc /= integral_kind*n_integrals) then
-    print *, irp_here,  ': f77_zmq_send( zmq_socket_push, buffer_value, integral_kind*n_integrals, 0)'
-    stop 'error'
-  endif
-  
-  rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0)
-  if (rc /= 4) then
-    print *, irp_here,  ': f77_zmq_send( zmq_socket_push, task_id, 4, 0)'
-    stop 'error'
-  endif
-
-! Activate is zmq_socket_push is a REQ
-  integer :: idummy
-  rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
-  if (rc /= 4) then
-    print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
-    stop 'error'
-  endif
-end
-
-
-
-
-
-subroutine ao_bielec_integrals_in_map_slave(thread,iproc)
-  use map_module
-  use f77_zmq
-  implicit none
-  BEGIN_DOC
-! Computes a buffer of integrals
-  END_DOC
-
-  integer, intent(in)            :: thread, iproc
-
-  integer                        :: j,l,n_integrals
-  integer                        :: rc
-  real(integral_kind), allocatable :: buffer_value(:)
-  integer(key_kind), allocatable :: buffer_i(:)
-  
-  integer                        :: worker_id, task_id
-  character*(512)                :: task
-
-  integer(ZMQ_PTR),external      :: new_zmq_to_qp_run_socket
-  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket
-
-  integer(ZMQ_PTR), external     :: new_zmq_push_socket
-  integer(ZMQ_PTR)               :: zmq_socket_push
-
-  character*(64)                 :: state
-
-  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
-  zmq_socket_push      = new_zmq_push_socket(thread)
-
-  allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) )
-
-  call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
-
-  do 
-    call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task)
-    if (task_id == 0) exit
-    read(task,*) j, l
-    call compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) 
-    call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
-    call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id)
-  enddo
-
-
-  call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
-  deallocate( buffer_i, buffer_value )
-  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
-  call end_zmq_push_socket(zmq_socket_push,thread)
-
-end
-
-
-subroutine ao_bielec_integrals_in_map_collector
-  use map_module
-  use f77_zmq
-  implicit none
-  BEGIN_DOC
-! Collects results from the AO integral calculation
-  END_DOC
-
-  integer                        :: j,l,n_integrals
-  integer                        :: rc
-  
-  real(integral_kind), allocatable :: buffer_value(:)
-  integer(key_kind), allocatable :: buffer_i(:)
-  
-  integer(ZMQ_PTR),external      :: new_zmq_to_qp_run_socket
-  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket
-  
-  integer(ZMQ_PTR), external     :: new_zmq_pull_socket
-  integer(ZMQ_PTR)               :: zmq_socket_pull
-  
-  integer*8                      :: control, accu
-  integer                        :: task_id, more, sze
-
-  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
-  zmq_socket_pull = new_zmq_pull_socket()
-
-  sze = ao_num*ao_num
-  allocate ( buffer_i(sze), buffer_value(sze) )
-
-  accu = 0_8
-  more = 1
-  do while (more == 1)
-    
-    rc = f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0)
-    if (rc == -1) then
-      n_integrals = 0
-      return
-    endif
-    if (rc /= 4) then
-      print *, irp_here,  ': f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0)'
-      stop 'error'
-    endif
-    
-    if (n_integrals >= 0) then
-      
-      if (n_integrals > sze) then
-        deallocate (buffer_value, buffer_i)
-        sze = n_integrals
-        allocate (buffer_value(sze), buffer_i(sze))
-      endif
-
-      rc = f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0)
-      if (rc /= key_kind*n_integrals) then
-        print *,  rc, key_kind, n_integrals
-        print *, irp_here,  ': f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0)'
-        stop 'error'
-      endif
-      
-      rc = f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0)
-      if (rc /= integral_kind*n_integrals) then
-        print *, irp_here,  ': f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0)'
-        stop 'error'
-      endif
-
-      rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
-
-! Activate if zmq_socket_pull is a REP
-      rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
-      if (rc /= 4) then
-        print *,  irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
-        stop 'error'
-      endif
-
-      
-      call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
-      accu += n_integrals
-      if (task_id /= 0) then
-        call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more)
-      endif
-    endif
-
-  enddo
-
-  deallocate( buffer_i, buffer_value )
-
-  integer (map_size_kind) :: get_ao_map_size 
-  control = get_ao_map_size(ao_integrals_map)
-
-  if (control /= accu) then
-      print *, ''
-      print *, irp_here
-      print *, 'Control : ', control
-      print *, 'Accu    : ', accu
-      print *, 'Some integrals were lost during the parallel computation.'
-      print *, 'Try to reduce the number of threads.'
-      stop
-  endif
-
-  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
-  call end_zmq_pull_socket(zmq_socket_pull)
-
-end
-

plugins/Integrals_erf/gauss_legendre.irp.f

@@ -1,66 +0,0 @@
-BEGIN_PROVIDER [ integer, n_pt_max_integrals_16 ]
- implicit none
- BEGIN_DOC
- ! Aligned n_pt_max_integrals
- END_DOC
- integer, external :: align_double
- n_pt_max_integrals_16 = align_double(n_pt_max_integrals)
-END_PROVIDER
-
- BEGIN_PROVIDER [ double precision, gauleg_t2, (n_pt_max_integrals_16,n_pt_max_integrals/2) ]
-&BEGIN_PROVIDER [ double precision, gauleg_w, (n_pt_max_integrals_16,n_pt_max_integrals/2) ]
-   implicit none
-   BEGIN_DOC
-   ! t_w(i,1,k) = w(i)
-   ! t_w(i,2,k) = t(i)
-   END_DOC
-   integer                        :: i,j,l
-   l=0
-   do i = 2,n_pt_max_integrals,2
-     l = l+1
-     call gauleg(0.d0,1.d0,gauleg_t2(1,l),gauleg_w(1,l),i)
-     do j=1,i
-       gauleg_t2(j,l) *= gauleg_t2(j,l)
-     enddo
-   enddo
-   
-END_PROVIDER
- 
-subroutine gauleg(x1,x2,x,w,n)
-   implicit none
-   BEGIN_DOC
-   ! Gauss-Legendre
-   END_DOC
-   integer, intent(in)            :: n
-   double precision, intent(in)   :: x1, x2
-   double precision, intent (out) :: x(n),w(n)
-   double precision, parameter    :: eps=3.d-14
-   
-   integer                        :: m,i,j
-   double precision               :: xm, xl, z, z1, p1, p2, p3, pp, dn
-   m=(n+1)/2
-   xm=0.5d0*(x2+x1)
-   xl=0.5d0*(x2-x1)
-   dn = dble(n)
-   do i=1,m
-     z=dcos(3.141592654d0*(dble(i)-.25d0)/(dble(n)+.5d0))
-     z1 = z+1.d0
-     do while (dabs(z-z1) > eps)
-       p1=1.d0
-       p2=0.d0
-       do j=1,n
-         p3=p2
-         p2=p1
-         p1=(dble(j+j-1)*z*p2-dble(j-1)*p3)/j
-       enddo
-       pp=dn*(z*p1-p2)/(z*z-1.d0)
-       z1=z
-       z=z1-p1/pp
-     end do
-     x(i)=xm-xl*z
-     x(n+1-i)=xm+xl*z
-     w(i)=(xl+xl)/((1.d0-z*z)*pp*pp)
-     w(n+1-i)=w(i)
-   enddo
-end
- 

plugins/Integrals_erf/integrals_3_index.irp.f

@@ -1,22 +0,0 @@
- BEGIN_PROVIDER [double precision, big_array_coulomb_integrals, (mo_tot_num_align,mo_tot_num, mo_tot_num)]
-&BEGIN_PROVIDER [double precision, big_array_exchange_integrals,(mo_tot_num_align,mo_tot_num, mo_tot_num)]
- implicit none
- integer :: i,j,k,l
- double precision :: get_mo_bielec_integral
- double precision :: integral
-
- do k = 1, mo_tot_num
-  do i = 1, mo_tot_num
-   do j = 1, mo_tot_num
-     l = j
-     integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
-     big_array_coulomb_integrals(j,i,k) = integral
-     l = j
-     integral = get_mo_bielec_integral(i,j,l,k,mo_integrals_map)
-     big_array_exchange_integrals(j,i,k) = integral
-   enddo
-  enddo
- enddo
-
-
-END_PROVIDER 

plugins/Integrals_erf/integrals_3_index_erf.irp.f

@@ -0,0 +1,22 @@
+ BEGIN_PROVIDER [double precision, big_array_coulomb_integrals_erf, (mo_tot_num_align,mo_tot_num, mo_tot_num)]
+&BEGIN_PROVIDER [double precision, big_array_exchange_integrals_erf,(mo_tot_num_align,mo_tot_num, mo_tot_num)]
+ implicit none
+ integer :: i,j,k,l
+ double precision :: get_mo_bielec_integral_erf
+ double precision :: integral
+
+ do k = 1, mo_tot_num
+  do i = 1, mo_tot_num
+   do j = 1, mo_tot_num
+     l = j
+     integral = get_mo_bielec_integral_erf(i,j,k,l,mo_integrals_erf_map)
+     big_array_coulomb_integrals_erf(j,i,k) = integral
+     l = j
+     integral = get_mo_bielec_integral_erf(i,j,l,k,mo_integrals_erf_map)
+     big_array_exchange_integrals_erf(j,i,k) = integral
+   enddo
+  enddo
+ enddo
+
+
+END_PROVIDER 

plugins/Integrals_erf/map_integrals.irp.f

@@ -1,717 +0,0 @@
-use map_module
-
-!! AO Map
-!! ======
-
-BEGIN_PROVIDER [ type(map_type), ao_integrals_map ]
-  implicit none
-  BEGIN_DOC
-  ! AO integrals
-  END_DOC
-  integer(key_kind)              :: key_max
-  integer(map_size_kind)         :: sze
-  call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max)
-  sze = key_max
-  call map_init(ao_integrals_map,sze)
-  print*,  'AO map initialized : ', sze
-END_PROVIDER
-
-subroutine bielec_integrals_index(i,j,k,l,i1)
-  use map_module
-  implicit none
-  integer, intent(in)            :: i,j,k,l
-  integer(key_kind), intent(out) :: i1
-  integer(key_kind)              :: p,q,r,s,i2
-  p = min(i,k)
-  r = max(i,k)
-  p = p+ishft(r*r-r,-1)
-  q = min(j,l)
-  s = max(j,l)
-  q = q+ishft(s*s-s,-1)
-  i1 = min(p,q)
-  i2 = max(p,q)
-  i1 = i1+ishft(i2*i2-i2,-1)
-end
-
-subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
-  use map_module
-  implicit none
-  integer, intent(out)           :: i(8),j(8),k(8),l(8)
-  integer(key_kind), intent(in)  :: i1
-  integer(key_kind)              :: i2,i3
-  i = 0
-  i2   = ceiling(0.5d0*(dsqrt(8.d0*dble(i1)+1.d0)-1.d0))
-  l(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i2)+1.d0)-1.d0))
-  i3   = i1 - ishft(i2*i2-i2,-1)
-  k(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i3)+1.d0)-1.d0))
-  j(1) = int(i2 - ishft(l(1)*l(1)-l(1),-1),4)
-  i(1) = int(i3 - ishft(k(1)*k(1)-k(1),-1),4)
-
-              !ijkl
-  i(2) = i(1) !ilkj
-  j(2) = l(1)
-  k(2) = k(1)
-  l(2) = j(1)
-
-  i(3) = k(1) !kjil
-  j(3) = j(1)
-  k(3) = i(1)
-  l(3) = l(1)
-
-  i(4) = k(1) !klij
-  j(4) = l(1)
-  k(4) = i(1)
-  l(4) = j(1)
-
-  i(5) = j(1) !jilk
-  j(5) = i(1)
-  k(5) = l(1)
-  l(5) = k(1)
-
-  i(6) = j(1) !jkli
-  j(6) = k(1)
-  k(6) = l(1)
-  l(6) = i(1)
-
-  i(7) = l(1) !lijk
-  j(7) = i(1)
-  k(7) = j(1)
-  l(7) = k(1)
-
-  i(8) = l(1) !lkji
-  j(8) = k(1)
-  k(8) = j(1)
-  l(8) = i(1)
-
-  integer :: ii, jj
-  do ii=2,8
-    do jj=1,ii-1
-      if ( (i(ii) == i(jj)).and. &
-           (j(ii) == j(jj)).and. &
-           (k(ii) == k(jj)).and. &
-           (l(ii) == l(jj)) ) then
-         i(ii) = 0
-         exit
-      endif
-    enddo
-  enddo
-  do ii=1,8
-    if (i(ii) /= 0) then
-      call bielec_integrals_index(i(ii),j(ii),k(ii),l(ii),i2)
-      if (i1 /= i2) then
-        print *,  i1, i2
-        print *,  i(ii), j(jj), k(jj), l(jj)
-        stop 'bielec_integrals_index_reverse failed'
-      endif
-    endif
-  enddo
-
-
-end
-
- BEGIN_PROVIDER [ integer, ao_integrals_cache_min ]
-&BEGIN_PROVIDER [ integer, ao_integrals_cache_max ]
- implicit none
- BEGIN_DOC
- ! Min and max values of the AOs for which the integrals are in the cache
- END_DOC
- ao_integrals_cache_min = max(1,ao_num - 63)
- ao_integrals_cache_max = ao_num
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ]
- implicit none
- BEGIN_DOC
- ! Cache of AO integrals for fast access
- END_DOC
- PROVIDE ao_bielec_integrals_in_map
- integer                        :: i,j,k,l,ii
- integer(key_kind)              :: idx
- real(integral_kind)            :: integral
- !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
- do l=ao_integrals_cache_min,ao_integrals_cache_max
-   do k=ao_integrals_cache_min,ao_integrals_cache_max
-     do j=ao_integrals_cache_min,ao_integrals_cache_max
-       do i=ao_integrals_cache_min,ao_integrals_cache_max
-         !DIR$ FORCEINLINE
-         call bielec_integrals_index(i,j,k,l,idx)
-         !DIR$ FORCEINLINE
-         call map_get(ao_integrals_map,idx,integral)
-         ii = l-ao_integrals_cache_min
-         ii = ior( ishft(ii,6), k-ao_integrals_cache_min)
-         ii = ior( ishft(ii,6), j-ao_integrals_cache_min)
-         ii = ior( ishft(ii,6), i-ao_integrals_cache_min)
-         ao_integrals_cache(ii) = integral
-       enddo
-     enddo
-   enddo
- enddo
- !$OMP END PARALLEL DO
-
-END_PROVIDER
-
-
-double precision function get_ao_bielec_integral(i,j,k,l,map) result(result)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Gets one AO bi-electronic integral from the AO map
-  END_DOC
-  integer, intent(in)            :: i,j,k,l
-  integer(key_kind)              :: idx
-  type(map_type), intent(inout)  :: map
-  integer                        :: ii
-  real(integral_kind)            :: tmp
-  PROVIDE ao_bielec_integrals_in_map ao_integrals_cache ao_integrals_cache_min
-  !DIR$ FORCEINLINE
-  if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
-    tmp = 0.d0
-  else if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < ao_integrals_threshold) then
-    tmp = 0.d0
-  else
-    ii = l-ao_integrals_cache_min
-    ii = ior(ii, k-ao_integrals_cache_min)
-    ii = ior(ii, j-ao_integrals_cache_min)
-    ii = ior(ii, i-ao_integrals_cache_min)
-    if (iand(ii, -64) /= 0) then
-      !DIR$ FORCEINLINE
-      call bielec_integrals_index(i,j,k,l,idx)
-      !DIR$ FORCEINLINE
-      call map_get(map,idx,tmp)
-      tmp = tmp
-    else
-      ii = l-ao_integrals_cache_min
-      ii = ior( ishft(ii,6), k-ao_integrals_cache_min)
-      ii = ior( ishft(ii,6), j-ao_integrals_cache_min)
-      ii = ior( ishft(ii,6), i-ao_integrals_cache_min)
-      tmp = ao_integrals_cache(ii)
-    endif
-  endif
-  result = tmp
-end
-
-
-subroutine get_ao_bielec_integrals(j,k,l,sze,out_val)
-  use map_module
-  BEGIN_DOC
-  ! Gets multiple AO bi-electronic integral from the AO map .
-  ! All i are retrieved for j,k,l fixed.
-  END_DOC
-  implicit none
-  integer, intent(in)            :: j,k,l, sze
-  real(integral_kind), intent(out) :: out_val(sze)
-  
-  integer                        :: i
-  integer(key_kind)              :: hash
-  double precision               :: thresh
-  PROVIDE ao_bielec_integrals_in_map ao_integrals_map
-  thresh = ao_integrals_threshold
-  
-  if (ao_overlap_abs(j,l) < thresh) then
-    out_val = 0.d0
-    return
-  endif
-  
-  double precision :: get_ao_bielec_integral
-  do i=1,sze
-    out_val(i) = get_ao_bielec_integral(i,j,k,l,ao_integrals_map)
-  enddo
-  
-end
-
-subroutine get_ao_bielec_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Gets multiple AO bi-electronic integral from the AO map .
-  ! All non-zero i are retrieved for j,k,l fixed.
-  END_DOC
-  integer, intent(in)            :: j,k,l, sze
-  real(integral_kind), intent(out) :: out_val(sze)
-  integer, intent(out)           :: out_val_index(sze),non_zero_int
-  
-  integer                        :: i
-  integer(key_kind)              :: hash
-  double precision               :: thresh,tmp
-  PROVIDE ao_bielec_integrals_in_map
-  thresh = ao_integrals_threshold
-  
-  non_zero_int = 0
-  if (ao_overlap_abs(j,l) < thresh) then
-    out_val = 0.d0
-    return
-  endif
- 
-  non_zero_int = 0
-  do i=1,sze
-    integer, external :: ao_l4
-    double precision, external :: ao_bielec_integral
-    !DIR$ FORCEINLINE
-    if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh) then
-      cycle
-    endif
-    call bielec_integrals_index(i,j,k,l,hash)
-    call map_get(ao_integrals_map, hash,tmp)
-    if (dabs(tmp) < thresh ) cycle
-    non_zero_int = non_zero_int+1
-    out_val_index(non_zero_int) = i
-    out_val(non_zero_int) = tmp
-  enddo
-  
-end
-
-
-function get_ao_map_size()
-  implicit none
-  integer (map_size_kind) :: get_ao_map_size
-  BEGIN_DOC
-  ! Returns the number of elements in the AO map
-  END_DOC
-  get_ao_map_size = ao_integrals_map % n_elements
-end
-
-subroutine clear_ao_map
-  implicit none
-  BEGIN_DOC
-  ! Frees the memory of the AO map
-  END_DOC
-  call map_deinit(ao_integrals_map)
-  FREE ao_integrals_map
-end
-
-
-!! MO Map
-!! ======
-
-BEGIN_PROVIDER [ type(map_type), mo_integrals_map ]
-  implicit none
-  BEGIN_DOC
-  ! MO integrals
-  END_DOC
-  integer(key_kind)              :: key_max
-  integer(map_size_kind)         :: sze
-  call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max)
-  sze = key_max
-  call map_init(mo_integrals_map,sze)
-  print*, 'MO map initialized'
-END_PROVIDER
-
-subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Create new entry into AO map
-  END_DOC
-  
-  integer, intent(in)                :: n_integrals
-  integer(key_kind), intent(inout)   :: buffer_i(n_integrals)
-  real(integral_kind), intent(inout) :: buffer_values(n_integrals)
-  
-  call map_append(ao_integrals_map, buffer_i, buffer_values, n_integrals)
-end
-
-subroutine insert_into_mo_integrals_map(n_integrals,                 &
-      buffer_i, buffer_values, thr)
-  use map_module
-  implicit none
-  
-  BEGIN_DOC
-  ! Create new entry into MO map, or accumulate in an existing entry
-  END_DOC
-  
-  integer, intent(in)                :: n_integrals
-  integer(key_kind), intent(inout)   :: buffer_i(n_integrals)
-  real(integral_kind), intent(inout) :: buffer_values(n_integrals)
-  real(integral_kind), intent(in)    :: thr
-  call map_update(mo_integrals_map, buffer_i, buffer_values, n_integrals, thr)
-end
-
- BEGIN_PROVIDER [ integer, mo_integrals_cache_min ]
-&BEGIN_PROVIDER [ integer, mo_integrals_cache_max ]
- implicit none
- BEGIN_DOC
- ! Min and max values of the MOs for which the integrals are in the cache
- END_DOC
- mo_integrals_cache_min = max(1,elec_alpha_num - 31)
- mo_integrals_cache_max = min(mo_tot_num,mo_integrals_cache_min+63)
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0:64*64*64*64) ]
- implicit none
- BEGIN_DOC
- ! Cache of MO integrals for fast access
- END_DOC
- PROVIDE mo_bielec_integrals_in_map
- integer                        :: i,j,k,l
- integer                        :: ii
- integer(key_kind)              :: idx
- real(integral_kind)            :: integral
- FREE ao_integrals_cache
- !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
- do l=mo_integrals_cache_min,mo_integrals_cache_max
-   do k=mo_integrals_cache_min,mo_integrals_cache_max
-     do j=mo_integrals_cache_min,mo_integrals_cache_max
-       do i=mo_integrals_cache_min,mo_integrals_cache_max
-         !DIR$ FORCEINLINE
-         call bielec_integrals_index(i,j,k,l,idx)
-         !DIR$ FORCEINLINE
-         call map_get(mo_integrals_map,idx,integral)
-         ii = l-mo_integrals_cache_min
-         ii = ior( ishft(ii,6), k-mo_integrals_cache_min)
-         ii = ior( ishft(ii,6), j-mo_integrals_cache_min)
-         ii = ior( ishft(ii,6), i-mo_integrals_cache_min)
-         mo_integrals_cache(ii) = integral
-       enddo
-     enddo
-   enddo
- enddo
- !$OMP END PARALLEL DO
-
-END_PROVIDER
-
-
-double precision function get_mo_bielec_integral(i,j,k,l,map)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Returns one integral <ij|kl> in the MO basis
-  END_DOC
-  integer, intent(in)            :: i,j,k,l
-  integer(key_kind)              :: idx
-  integer                        :: ii
-  type(map_type), intent(inout)  :: map
-  real(integral_kind)            :: tmp
-  PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
-  ii = l-mo_integrals_cache_min
-  ii = ior(ii, k-mo_integrals_cache_min)
-  ii = ior(ii, j-mo_integrals_cache_min)
-  ii = ior(ii, i-mo_integrals_cache_min)
-  if (iand(ii, -64) /= 0) then
-    !DIR$ FORCEINLINE
-    call bielec_integrals_index(i,j,k,l,idx)
-    !DIR$ FORCEINLINE
-    call map_get(map,idx,tmp)
-    get_mo_bielec_integral = dble(tmp)
-  else
-    ii = l-mo_integrals_cache_min
-    ii = ior( ishft(ii,6), k-mo_integrals_cache_min)
-    ii = ior( ishft(ii,6), j-mo_integrals_cache_min)
-    ii = ior( ishft(ii,6), i-mo_integrals_cache_min)
-    get_mo_bielec_integral = mo_integrals_cache(ii)
-  endif
-end
-
-
-double precision function mo_bielec_integral(i,j,k,l)
-  implicit none
-  BEGIN_DOC
-  ! Returns one integral <ij|kl> in the MO basis
-  END_DOC
-  integer, intent(in)            :: i,j,k,l
-  double precision               :: get_mo_bielec_integral
-  PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
-  !DIR$ FORCEINLINE
-  PROVIDE mo_bielec_integrals_in_map
-  mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
-  return
-end
-
-subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Returns multiple integrals <ij|kl> in the MO basis, all
-  ! i for j,k,l fixed.
-  END_DOC
-  integer, intent(in)            :: j,k,l, sze
-  double precision, intent(out)  :: out_val(sze)
-  type(map_type), intent(inout)  :: map
-  integer                        :: i
-  integer(key_kind)              :: hash(sze)
-  real(integral_kind)            :: tmp_val(sze)
-  PROVIDE mo_bielec_integrals_in_map
-  
-  do i=1,sze
-    !DIR$ FORCEINLINE
-    call bielec_integrals_index(i,j,k,l,hash(i))
-  enddo
-  
-  if (key_kind == 8) then
-    call map_get_many(map, hash, out_val, sze)
-  else
-    call map_get_many(map, hash, tmp_val, sze)
-    ! Conversion to double precision 
-    do i=1,sze
-      out_val(i) = dble(tmp_val(i))
-    enddo
-  endif
-end
-
-subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Returns multiple integrals <ij|kl> in the MO basis, all
-  ! i(1)j(2) 1/r12 k(1)l(2)
-  ! i, j for k,l fixed.
-  END_DOC
-  integer, intent(in)            :: k,l, sze
-  double precision, intent(out)  :: out_array(sze,sze)
-  type(map_type), intent(inout)  :: map
-  integer                        :: i,j,kk,ll,m
-  integer(key_kind),allocatable  :: hash(:)
-  integer  ,allocatable          :: pairs(:,:), iorder(:)
-  real(integral_kind), allocatable :: tmp_val(:)
-
-  PROVIDE mo_bielec_integrals_in_map
-  allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), &
-  tmp_val(sze*sze))
-  
-  kk=0
-  out_array = 0.d0
-  do j=1,sze
-   do i=1,sze
-    kk += 1
-    !DIR$ FORCEINLINE
-    call bielec_integrals_index(i,j,k,l,hash(kk))
-    pairs(1,kk) = i
-    pairs(2,kk) = j
-    iorder(kk) = kk
-   enddo
-  enddo
-
-  logical :: integral_is_in_map
-  if (key_kind == 8) then
-    call i8radix_sort(hash,iorder,kk,-1)
-  else if (key_kind == 4) then
-    call iradix_sort(hash,iorder,kk,-1)
-  else if (key_kind == 2) then
-    call i2radix_sort(hash,iorder,kk,-1)
-  endif
-
-  call map_get_many(mo_integrals_map, hash, tmp_val, kk)
-
-  do ll=1,kk
-    m = iorder(ll)
-    i=pairs(1,m)
-    j=pairs(2,m)
-    out_array(i,j) = tmp_val(ll)
-  enddo  
-
-  deallocate(pairs,hash,iorder,tmp_val)
-end
-
-subroutine get_mo_bielec_integrals_coulomb_ii(k,l,sze,out_val,map)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Returns multiple integrals <ki|li> 
-  ! k(1)i(2) 1/r12 l(1)i(2) :: out_val(i1)
-  ! for k,l fixed.
-  END_DOC
-  integer, intent(in)            :: k,l, sze
-  double precision, intent(out)  :: out_val(sze)
-  type(map_type), intent(inout)  :: map
-  integer                        :: i
-  integer(key_kind)              :: hash(sze)
-  real(integral_kind)            :: tmp_val(sze)
-  PROVIDE mo_bielec_integrals_in_map
-  
-  integer :: kk
-  do i=1,sze
-    !DIR$ FORCEINLINE
-    call bielec_integrals_index(k,i,l,i,hash(i))
-  enddo
-  
-  if (key_kind == 8) then
-    call map_get_many(map, hash, out_val, sze)
-  else
-    call map_get_many(map, hash, tmp_val, sze)
-    ! Conversion to double precision 
-    do i=1,sze
-      out_val(i) = dble(tmp_val(i))
-    enddo
-  endif
-end
-
-subroutine get_mo_bielec_integrals_exch_ii(k,l,sze,out_val,map)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Returns multiple integrals <ki|il> 
-  ! k(1)i(2) 1/r12 i(1)l(2) :: out_val(i1)
-  ! for k,l fixed.
-  END_DOC
-  integer, intent(in)            :: k,l, sze
-  double precision, intent(out)  :: out_val(sze)
-  type(map_type), intent(inout)  :: map
-  integer                        :: i
-  integer(key_kind)              :: hash(sze)
-  real(integral_kind)            :: tmp_val(sze)
-  PROVIDE mo_bielec_integrals_in_map
-  
-  integer :: kk
-  do i=1,sze
-    !DIR$ FORCEINLINE
-    call bielec_integrals_index(k,i,i,l,hash(i))
-  enddo
-  
-  if (key_kind == 8) then
-    call map_get_many(map, hash, out_val, sze)
-  else
-    call map_get_many(map, hash, tmp_val, sze)
-    ! Conversion to double precision 
-    do i=1,sze
-      out_val(i) = dble(tmp_val(i))
-    enddo
-  endif
-end
-
-
-integer*8 function get_mo_map_size()
-  implicit none
-  BEGIN_DOC
-  ! Return the number of elements in the MO map
-  END_DOC
-  get_mo_map_size = mo_integrals_map % n_elements
-end
-
-BEGIN_TEMPLATE
-
-subroutine dump_$ao_integrals(filename)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Save to disk the $ao integrals
-  END_DOC
-  character*(*), intent(in)      :: filename
-  integer(cache_key_kind), pointer :: key(:)
-  real(integral_kind), pointer   :: val(:)
-  integer*8                      :: i,j, n
-  call ezfio_set_work_empty(.False.)
-  open(unit=66,file=filename,FORM='unformatted')
-  write(66) integral_kind, key_kind
-  write(66) $ao_integrals_map%sorted, $ao_integrals_map%map_size,    &
-      $ao_integrals_map%n_elements
-  do i=0_8,$ao_integrals_map%map_size
-    write(66) $ao_integrals_map%map(i)%sorted, $ao_integrals_map%map(i)%map_size,&
-        $ao_integrals_map%map(i)%n_elements
-  enddo
-  do i=0_8,$ao_integrals_map%map_size
-    key => $ao_integrals_map%map(i)%key
-    val => $ao_integrals_map%map(i)%value
-    n = $ao_integrals_map%map(i)%n_elements
-    write(66) (key(j), j=1,n), (val(j), j=1,n)
-  enddo
-  close(66)
-  
-end
-
-IRP_IF COARRAY
-subroutine communicate_$ao_integrals()
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Communicate the $ao integrals with co-array
-  END_DOC
-  integer(cache_key_kind), pointer :: key(:)
-  real(integral_kind), pointer   :: val(:)
-  integer*8                      :: i,j, k, nmax
-  integer*8, save                :: n[*]
-  integer                        :: copy_n
-
-  real(integral_kind), allocatable            :: buffer_val(:)[:]
-  integer(cache_key_kind), allocatable        :: buffer_key(:)[:]
-  real(integral_kind), allocatable            :: copy_val(:)
-  integer(key_kind), allocatable              :: copy_key(:)
-
-  n = 0_8
-  do i=0_8,$ao_integrals_map%map_size
-    n = max(n,$ao_integrals_map%map(i)%n_elements)
-  enddo
-  sync all
-  nmax = 0_8
-  do j=1,num_images()
-    nmax = max(nmax,n[j])
-  enddo
-  allocate( buffer_key(nmax)[*], buffer_val(nmax)[*])
-  allocate( copy_key(nmax), copy_val(nmax))
-  do i=0_8,$ao_integrals_map%map_size
-    key => $ao_integrals_map%map(i)%key
-    val => $ao_integrals_map%map(i)%value
-    n = $ao_integrals_map%map(i)%n_elements
-    do j=1,n
-      buffer_key(j) = key(j)
-      buffer_val(j) = val(j)
-    enddo
-    sync all
-    do j=1,num_images()
-      if (j /= this_image()) then
-        copy_n = n[j]
-        do k=1,copy_n
-          copy_val(k) = buffer_val(k)[j]
-          copy_key(k) = buffer_key(k)[j]
-          copy_key(k) = copy_key(k)+ishft(i,-map_shift)
-        enddo
-        call map_append($ao_integrals_map, copy_key, copy_val, copy_n )
-      endif
-    enddo
-    sync all
-  enddo
-  deallocate( buffer_key, buffer_val, copy_val, copy_key)
-  
-end
-IRP_ENDIF 
-
-
-integer function load_$ao_integrals(filename)
-  implicit none
-  BEGIN_DOC
-  ! Read from disk the $ao integrals
-  END_DOC
-  character*(*), intent(in)      :: filename
-  integer*8                      :: i
-  integer(cache_key_kind), pointer :: key(:)
-  real(integral_kind), pointer   :: val(:)
-  integer                        :: iknd, kknd
-  integer*8                      :: n, j
-  load_$ao_integrals = 1
-  open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
-  read(66,err=98,end=98) iknd, kknd
-  if (iknd /= integral_kind) then
-    print *,  'Wrong integrals kind in file :', iknd
-    stop 1
-  endif
-  if (kknd /= key_kind) then
-    print *,  'Wrong key kind in file :', kknd
-    stop 1
-  endif
-  read(66,err=98,end=98) $ao_integrals_map%sorted, $ao_integrals_map%map_size,&
-      $ao_integrals_map%n_elements
-  do i=0_8, $ao_integrals_map%map_size
-    read(66,err=99,end=99) $ao_integrals_map%map(i)%sorted,          &
-        $ao_integrals_map%map(i)%map_size, $ao_integrals_map%map(i)%n_elements
-    call cache_map_reallocate($ao_integrals_map%map(i),$ao_integrals_map%map(i)%map_size)
-  enddo
-  do i=0_8, $ao_integrals_map%map_size
-    key => $ao_integrals_map%map(i)%key
-    val => $ao_integrals_map%map(i)%value
-    n = $ao_integrals_map%map(i)%n_elements
-    read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n)
-  enddo
-  call map_sort($ao_integrals_map)
-  load_$ao_integrals = 0
-  return
-  99 continue
-  call map_deinit($ao_integrals_map)
-  98 continue
-  stop 'Problem reading $ao_integrals_map file in work/'
-  
-end
-
-SUBST [ ao_integrals_map, ao_integrals, ao_num ]
-ao_integrals_map ; ao_integrals ; ao_num ;;
-mo_integrals_map ; mo_integrals ; mo_tot_num ;;
-END_TEMPLATE

plugins/Integrals_erf/map_integrals_erf.irp.f

@@ -156,7 +156,7 @@ subroutine get_ao_bielec_integrals_erf_non_zero(j,k,l,sze,out_val,out_val_index,
     integer, external :: ao_l4
     double precision, external :: ao_bielec_integral_erf
     !DIR$ FORCEINLINE
-    if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh) then
+    if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thresh) then
       cycle
     endif
     call bielec_integrals_index(i,j,k,l,hash)
@@ -395,7 +395,7 @@ BEGIN_PROVIDER [ double precision, mo_integrals_erf_cache, (0:64*64*64*64) ]
  integer                        :: ii
  integer(key_kind)              :: idx
  real(integral_kind)            :: integral
- FREE ao_integrals_cache
+ FREE ao_integrals_erf_cache
  !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
  do l=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
    do k=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
@@ -478,7 +478,7 @@ subroutine get_mo_bielec_integrals_erf(j,k,l,sze,out_val,map)
   integer                        :: i
   integer(key_kind)              :: hash(sze)
   real(integral_kind)            :: tmp_val(sze)
-  PROVIDE mo_bielec_integrals_in_map
+  PROVIDE mo_bielec_integrals_erf_in_map
   
   do i=1,sze
     !DIR$ FORCEINLINE
@@ -564,7 +564,7 @@ subroutine get_mo_bielec_integrals_erf_coulomb_ii(k,l,sze,out_val,map)
   integer                        :: i
   integer(key_kind)              :: hash(sze)
   real(integral_kind)            :: tmp_val(sze)
-  PROVIDE mo_bielec_integrals_in_map
+  PROVIDE mo_bielec_integrals_erf_in_map
   
   integer :: kk
   do i=1,sze

plugins/Integrals_erf/mo_bi_integrals_erf.irp.f

@@ -587,8 +587,8 @@ END_PROVIDER
   
   do j=1,mo_tot_num
     do i=1,mo_tot_num
-      mo_bielec_integral_erf_jj(i,j) = get_mo_bielec_integral_erf(i,j,i,j,mo_integrals_map)
+      mo_bielec_integral_erf_jj(i,j) = get_mo_bielec_integral_erf(i,j,i,j,mo_integrals_erf_map)
-      mo_bielec_integral_erf_jj_exchange(i,j) = get_mo_bielec_integral_erf(i,j,j,i,mo_integrals_map)
+      mo_bielec_integral_erf_jj_exchange(i,j) = get_mo_bielec_integral_erf(i,j,j,i,mo_integrals_erf_map)
       mo_bielec_integral_erf_jj_anti(i,j) = mo_bielec_integral_erf_jj(i,j) - mo_bielec_integral_erf_jj_exchange(i,j)
     enddo
   enddo

plugins/Integrals_erf/providers_ao_standard.irp.f

@@ -1,116 +0,0 @@
-
-BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
-  implicit none
-  use f77_zmq
-  use map_module
-  BEGIN_DOC
-  !  Map of Atomic integrals
-  !     i(r1) j(r2) 1/r12 k(r1) l(r2)
-  END_DOC
-  
-  integer                        :: i,j,k,l
-  double precision               :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2
-  double precision               :: integral, wall_0
-  include 'Utils/constants.include.F'
-  
-  ! For integrals file
-  integer(key_kind),allocatable  :: buffer_i(:)
-  integer,parameter              :: size_buffer = 1024*64
-  real(integral_kind),allocatable :: buffer_value(:)
-  
-  integer                        :: n_integrals, rc
-  integer                        :: kk, m, j1, i1, lmax
-  character*(64)                 :: fmt
-  
-  integral = ao_bielec_integral(1,1,1,1)
-  
-  real                           :: map_mb
-  PROVIDE read_ao_integrals disk_access_ao_integrals
-  if (read_ao_integrals) then
-    print*,'Reading the AO integrals'
-      call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
-      print*, 'AO integrals provided'
-      ao_bielec_integrals_in_map = .True.
-      return
-  endif
-  
-  print*, 'Providing the AO integrals'
-  call wall_time(wall_0)
-  call wall_time(wall_1)
-  call cpu_time(cpu_1)
-
-  integer(ZMQ_PTR) :: zmq_to_qp_run_socket
-  call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals')
-
-  character(len=:), allocatable :: task
-  allocate(character(len=ao_num*12) :: task)
-  write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))'
-  do l=1,ao_num
-    write(task,fmt) (i,l, i=1,l)
-    call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task))
-  enddo
-  deallocate(task)
-  
-  call zmq_set_running(zmq_to_qp_run_socket)
-
-  PROVIDE nproc
-  !$OMP PARALLEL DEFAULT(private) num_threads(nproc+1)
-      i = omp_get_thread_num()
-      if (i==0) then
-        call ao_bielec_integrals_in_map_collector(i)
-      else
-        call ao_bielec_integrals_in_map_slave_inproc(i)
-      endif
-  !$OMP END PARALLEL
-
-  call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals')
-
-
-  print*, 'Sorting the map'
-  call map_sort(ao_integrals_map)
-  call cpu_time(cpu_2)
-  call wall_time(wall_2)
-  integer(map_size_kind)         :: get_ao_map_size, ao_map_size
-  ao_map_size = get_ao_map_size()
-  
-  print*, 'AO integrals provided:'
-  print*, ' Size of AO map :         ', map_mb(ao_integrals_map) ,'MB'
-  print*, ' Number of AO integrals :', ao_map_size
-  print*, ' cpu  time :',cpu_2 - cpu_1, 's'
-  print*, ' wall time :',wall_2 - wall_1, 's  ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )'
-  
-  ao_bielec_integrals_in_map = .True.
-
-  if (write_ao_integrals) then
-    call ezfio_set_work_empty(.False.)
-    call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
-    call ezfio_set_integrals_erf_disk_access_ao_integrals("Read")
-  endif
-  
-END_PROVIDER
- 
-BEGIN_PROVIDER [ double precision, ao_bielec_integral_schwartz,(ao_num,ao_num)  ]
-  implicit none
-  BEGIN_DOC
-  !  Needed to compute Schwartz inequalities
-  END_DOC
-  
-  integer                        :: i,k
-  double precision               :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2
-  
-  ao_bielec_integral_schwartz(1,1) = ao_bielec_integral(1,1,1,1)
-  !$OMP PARALLEL DO PRIVATE(i,k)                                     &
-      !$OMP DEFAULT(NONE)                                            &
-      !$OMP SHARED (ao_num,ao_bielec_integral_schwartz)              &
-      !$OMP SCHEDULE(dynamic)
-  do i=1,ao_num
-    do k=1,i
-      ao_bielec_integral_schwartz(i,k) = dsqrt(ao_bielec_integral(i,k,i,k))
-      ao_bielec_integral_schwartz(k,i) = ao_bielec_integral_schwartz(i,k)
-    enddo
-  enddo
-  !$OMP END PARALLEL DO
-  
-END_PROVIDER
-
-

plugins/Integrals_erf/qp_ao_erf_ints.irp.f

@@ -2,7 +2,7 @@ program qp_ao_ints
   use omp_lib
   implicit none
   BEGIN_DOC
-! Increments a running calculation to compute AO integrals
+! Increments a running calculation to compute AO integral_erfs
   END_DOC
   integer :: i
 
@@ -11,18 +11,18 @@ program qp_ao_ints
   zmq_context = f77_zmq_ctx_new ()
 
   ! Set the state of the ZMQ
-  zmq_state = 'ao_integrals'
+  zmq_state = 'ao_integral_erfs'
 
   ! Provide everything needed
-  double precision :: integral, ao_bielec_integral
+  double precision :: integral_erf, ao_bielec_integral_erf
-  integral = ao_bielec_integral(1,1,1,1)
+  integral_erf = ao_bielec_integral_erf(1,1,1,1)
 
   character*(64) :: state
   call wait_for_state(zmq_state,state)
   do while (state /= 'Stopped')
     !$OMP PARALLEL DEFAULT(PRIVATE) PRIVATE(i)
     i = omp_get_thread_num()
-    call ao_bielec_integrals_in_map_slave_tcp(i)
+    call ao_bielec_integrals_erf_in_map_slave_tcp(i)
     !$OMP END PARALLEL
     call wait_for_state(zmq_state,state)
   enddo

plugins/Integrals_erf/read_write.irp.f

@@ -1,51 +1,3 @@
-BEGIN_PROVIDER [ logical, read_ao_integrals ]
-&BEGIN_PROVIDER [ logical, read_mo_integrals ]
-&BEGIN_PROVIDER [ logical, write_ao_integrals ]
-&BEGIN_PROVIDER [ logical, write_mo_integrals ]
-
- BEGIN_DOC
-! One level of abstraction for disk_access_ao_integrals and disk_access_mo_integrals
- END_DOC
-implicit none
-
-    if (disk_access_ao_integrals.EQ.'Read') then
-        read_ao_integrals =  .True.
-        write_ao_integrals = .False.
-
-    else if  (disk_access_ao_integrals.EQ.'Write') then
-        read_ao_integrals = .False.
-        write_ao_integrals =  .True.
-    
-    else if (disk_access_ao_integrals.EQ.'None') then
-        read_ao_integrals = .False.
-        write_ao_integrals = .False.
-
-    else
-        print *, 'bielec_integrals/disk_access_ao_integrals has a wrong type'
-        stop 1
-
-    endif
-
-    if (disk_access_mo_integrals.EQ.'Read') then
-        read_mo_integrals =  .True.
-        write_mo_integrals = .False.
-
-    else if  (disk_access_mo_integrals.EQ.'Write') then
-        read_mo_integrals = .False.
-        write_mo_integrals =  .True.
-
-    else if (disk_access_mo_integrals.EQ.'None') then
-        read_mo_integrals = .False.
-        write_mo_integrals = .False.
-
-    else
-        print *, 'bielec_integrals/disk_access_mo_integrals has a wrong type'
-        stop 1
-
-    endif
-
-END_PROVIDER
-
 BEGIN_PROVIDER [ logical, read_ao_integrals_erf ]
 &BEGIN_PROVIDER [ logical, read_mo_integrals_erf ]
 &BEGIN_PROVIDER [ logical, write_ao_integrals_erf ]

plugins/Integrals_erf/test_bitmasks_integrals.irp.f

@@ -1,33 +0,0 @@
-program pouet
- implicit none
-
- call routine
-
-end
-
-subroutine routine 
- implicit none
-  integer(bit_kind)  :: mask_ijkl(N_int,4)
-  integer, allocatable           :: list_ijkl(:,:)
-  integer :: i,j
-  integer :: n_i,n_j,n_k,n_l
-  do i = 1,N_int 
-   mask_ijkl(i,1) = inact_bitmask(i,1)
-   mask_ijkl(i,2) = inact_bitmask(i,1)
-   mask_ijkl(i,3) = inact_bitmask(i,1)
-   mask_ijkl(i,4) = inact_bitmask(i,1)
-  enddo
-  allocate(list_ijkl(mo_tot_num,4))
-  call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int )
-  call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int )
-  call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int )
-  call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int )
-  print*,'n_i,n_j = ',n_i,n_j
-  print*,'n_k,n_l = ',n_k,n_l
-  do i =1, n_i
-   print*,list_ijkl(i,1), list_ijkl(i,2)
-  enddo
-  deallocate(list_ijkl)
- 
-
-end

plugins/Integrals_restart_DFT/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Integrals_Monoelec Integrals_erf Determinants DFT_Utils

plugins/Integrals_restart_DFT/README.rst

@@ -0,0 +1,12 @@
+==============
+core_integrals
+==============
+
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.

plugins/Integrals_restart_DFT/short_range_coulomb.irp.f

@@ -0,0 +1,79 @@
+BEGIN_PROVIDER [double precision, density_matrix_read, (mo_tot_num, mo_tot_num)]
+ implicit none
+ integer :: i,j,k,l
+ logical                        :: exists
+ call ezfio_has_determinants_density_matrix_mo_disk(exists)
+ if(exists)then
+  print*, 'reading the density matrix from input'
+  call ezfio_get_determinants_density_matrix_mo_disk(exists)
+  print*, 'reading done'
+ else 
+  print*, 'no density matrix found in EZFIO file ...'
+  print*, 'stopping ..'
+  stop
+ endif
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [double precision, effective_short_range_operator, (mo_tot_num,mo_tot_num)]
+ implicit none
+ integer :: i,j,k,l,m,n
+ double precision :: get_mo_bielec_integral,get_mo_bielec_integral_erf
+ double precision :: integral, integral_erf
+ effective_short_range_operator = 0.d0
+ do i = 1, mo_tot_num
+  do j = 1, mo_tot_num
+   if(dabs(one_body_dm_mo(i,j)).le.1.d-10)cycle
+   do k = 1, mo_tot_num
+    do l = 1, mo_tot_num
+     integral = get_mo_bielec_integral(i,k,j,l,mo_integrals_map)
+!    integral_erf = get_mo_bielec_integral_erf(i,k,j,l,mo_integrals_erf_map)
+     effective_short_range_operator(l,k) += one_body_dm_mo(i,j) * integral
+    enddo
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [double precision, effective_one_e_potential, (mo_tot_num_align, mo_tot_num,N_states)]
+ implicit none
+ integer :: i,j,i_state
+ effective_one_e_potential = 0.d0
+ do i_state = 1, N_states
+  do i = 1, mo_tot_num
+   do j = 1, mo_tot_num
+    effective_one_e_potential(i,j,i_state) = effective_short_range_operator(i,j) + mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) & 
+                                   + 0.5d0 * (lda_ex_potential_alpha_ao(i,j,i_state) + lda_ex_potential_beta_ao(i,j,i_state))
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
+subroutine save_one_e_effective_potential  
+ implicit none
+ double precision, allocatable :: tmp(:,:)
+ allocate(tmp(size(effective_one_e_potential,1),size(effective_one_e_potential,2)))
+ integer :: i,j
+ do i = 1, mo_tot_num
+  do j = 1, mo_tot_num
+   tmp(i,j) = effective_one_e_potential(i,j,1)
+  enddo
+ enddo
+ call write_one_e_integrals('mo_one_integral', tmp,      &
+      size(tmp,1), size(tmp,2))
+ call ezfio_set_integrals_monoelec_disk_access_only_mo_one_integrals("Read")
+ deallocate(tmp)
+
+end
+
+subroutine save_erf_bi_elec_integrals
+ implicit none
+ integer :: i,j,k,l
+ PROVIDE mo_bielec_integrals_erf_in_map
+ call ezfio_set_work_empty(.False.)
+ call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_erf_map)
+ call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
+end

plugins/Integrals_restart_DFT/write_integrals_restart_dft.irp.f

@@ -0,0 +1,18 @@
+program write_integrals
+ implicit none
+ read_wf = .true.
+ touch read_wf
+ disk_access_only_mo_one_integrals = "None"
+ touch disk_access_only_mo_one_integrals
+ disk_access_mo_integrals = "None"
+ touch disk_access_mo_integrals
+ call routine
+
+end
+
+subroutine routine
+ implicit none
+ call save_one_e_effective_potential
+ call save_erf_bi_elec_integrals
+
+end

plugins/Slater_rules_DFT/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Determinants Integrals_restart_DFT Davidson 

plugins/Slater_rules_DFT/README.rst

@@ -0,0 +1,12 @@
+================
+Slater_rules_DFT
+================
+
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.

plugins/Slater_rules_DFT/Slater_rules_DFT.main.irp.f

@@ -0,0 +1,38 @@
+program Slater_rules_DFT
+  implicit none
+  BEGIN_DOC
+! TODO
+  END_DOC
+  print *, '  _/                              '
+  print *, ' -:\_?,     _Jm####La             '
+  print *, 'J"(:" >  _]#AZ#Z#UUZ##,           '
+  print *, '_,::./   %(|i%12XmX1*1XL      _?, '
+  print *, '  \..\ _\(vmWQwodY+ia%lnL  _",/ ( '
+  print *, '   .:< ]J=mQD?WXn<uQWmmvd, -.-:=!'
+  print *, '   "{Z jC]QW|=3Zv)Bi3BmXv3  =   _7'
+  print *, '    ]h[Z6)WQ;)jZs]C;|$BZv+, : ./  '
+  print *, '    -#sJX%$Wmm#ev]hinW#Xi:` c ;   '
+  print *, '     #X#X23###1}vI$WWmX1>|,)nr"   '
+  print *, '     4XZ#Xov1v}=)vnXAX1nnv;1n"    '
+  print *, '     ]XX#ZXoovvvivnnnlvvo2*i7     '
+  print *, '     "23Z#1S2oo2XXSnnnoSo2>v"     '
+  print *, '      miX#L -~`""!!1}oSoe|i7      '
+  print *, '      4cn#m,        v221=|v[      '
+  print *, '      ]hI3Zma,;..__wXSe=+vo       '
+  print *, '      ]Zov*XSUXXZXZXSe||vo2       '
+  print *, '      ]Z#><iiii|i||||==vn2(       '
+  print *, '      ]Z#i<ii||+|=||=:{no2[       '
+  print *, '      ]ZUsiiiiivi|=||=vo22[       '
+  print *, '      ]XZvlliiIi|i=|+|vooo        '
+  print *, '      =v1llli||||=|||||lii(       '
+  print *, '      ]iillii||||||||=>=|<        '
+  print *, '      -ziiiii||||||+||==+>        '
+  print *, '       -%|+++||=|=+|=|==/         '
+  print *, '        -a>====+|====-:-          '
+  print *, '          "~,- --   /-            '
+  print *, '            -.     )>             '
+  print *, '           .~      +-             '
+  print *, '           . .... : .             '
+  print *, '            -------~              '
+  print *, ''
+end

plugins/Slater_rules_DFT/energy.irp.f

@@ -0,0 +1,7 @@
+! BEGIN_PROVIDER [double precision, energy_total]
+!&BEGIN_PROVIDER [double precision, energy_one_electron]
+!&BEGIN_PROVIDER [double precision, energy_hartree]
+!&BEGIN_PROVIDER [double precision, energy]
+! implicit none
+!
+!END_PROVIDER 

plugins/Slater_rules_DFT/slater_rules_erf.irp.f

@@ -0,0 +1,445 @@
+
+subroutine i_H_j_erf(key_i,key_j,Nint,hij)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns <i|H|j> where i and j are determinants
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: key_i(Nint,2), key_j(Nint,2)
+  double precision, intent(out)  :: hij
+  
+  integer                        :: exc(0:2,2,2)
+  integer                        :: degree
+  double precision               :: get_mo_bielec_integral_erf
+  integer                        :: m,n,p,q
+  integer                        :: i,j,k
+  integer                        :: occ(Nint*bit_kind_size,2)
+  double precision               :: diag_H_mat_elem_erf, phase,phase_2
+  integer                        :: n_occ_ab(2)
+  PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_map big_array_exchange_integrals_erf
+  
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+  ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
+  ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
+  ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
+  ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
+  
+  hij = 0.d0
+  !DIR$ FORCEINLINE
+  call get_excitation_degree(key_i,key_j,degree,Nint)
+  integer :: spin
+  select case (degree)
+    case (2)
+      call get_double_excitation(key_i,key_j,exc,phase,Nint)
+      if (exc(0,1,1) == 1) then
+        ! Mono alpha, mono beta
+        if(exc(1,1,1) == exc(1,2,2) )then
+         hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
+        else if (exc(1,2,1) ==exc(1,1,2))then
+         hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
+        else
+         hij = phase*get_mo_bielec_integral_erf(                          &
+             exc(1,1,1),                                              &
+             exc(1,1,2),                                              &
+             exc(1,2,1),                                              &
+             exc(1,2,2) ,mo_integrals_erf_map)
+        endif
+      else if (exc(0,1,1) == 2) then
+        ! Double alpha
+        hij = phase*(get_mo_bielec_integral_erf(                         &
+            exc(1,1,1),                                              &
+            exc(2,1,1),                                              &
+            exc(1,2,1),                                              &
+            exc(2,2,1) ,mo_integrals_erf_map) -                          &
+            get_mo_bielec_integral_erf(                                  &
+            exc(1,1,1),                                              &
+            exc(2,1,1),                                              &
+            exc(2,2,1),                                              &
+            exc(1,2,1) ,mo_integrals_erf_map) )
+      else if (exc(0,1,2) == 2) then
+        ! Double beta
+        hij = phase*(get_mo_bielec_integral_erf(                         &
+            exc(1,1,2),                                              &
+            exc(2,1,2),                                              &
+            exc(1,2,2),                                              &
+            exc(2,2,2) ,mo_integrals_erf_map) -                          &
+            get_mo_bielec_integral_erf(                                  &
+            exc(1,1,2),                                              &
+            exc(2,1,2),                                              &
+            exc(2,2,2),                                              &
+            exc(1,2,2) ,mo_integrals_erf_map) )
+      endif
+    case (1)
+      call get_mono_excitation(key_i,key_j,exc,phase,Nint)
+      !DIR$ FORCEINLINE
+      call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
+      if (exc(0,1,1) == 1) then
+        ! Mono alpha
+        m = exc(1,1,1)
+        p = exc(1,2,1)
+        spin = 1
+        do i = 1, n_occ_ab(1)
+         hij += -big_array_exchange_integrals_erf(occ(i,1),m,p) + big_array_coulomb_integrals_erf(occ(i,1),m,p)
+        enddo
+        do i = 1, n_occ_ab(2)
+         hij += big_array_coulomb_integrals_erf(occ(i,2),m,p)
+        enddo
+      else
+        ! Mono beta
+        m = exc(1,1,2)
+        p = exc(1,2,2)
+        spin = 2
+        do i = 1, n_occ_ab(2)
+         hij += -big_array_exchange_integrals_erf(occ(i,2),m,p) + big_array_coulomb_integrals_erf(occ(i,2),m,p)
+        enddo
+        do i = 1, n_occ_ab(1)
+         hij += big_array_coulomb_integrals_erf(occ(i,1),m,p)
+        enddo
+      endif
+      hij = hij + mo_nucl_elec_integral(m,p) + mo_kinetic_integral(m,p)
+      hij = hij * phase
+    case (0)
+      hij = diag_H_mat_elem_erf(key_i,Nint)
+  end select
+end
+
+double precision function diag_H_mat_elem_erf(key_i,Nint)
+ implicit none
+ integer(bit_kind), intent(in) :: key_i(N_int,2)
+ integer, intent(in)  :: Nint
+ integer :: i,j
+ integer                        :: occ(Nint*bit_kind_size,2)
+ integer                        :: n_occ_ab(2)
+ call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
+ diag_H_mat_elem_erf = 0.d0
+ ! alpha - alpha
+ do i = 1, n_occ_ab(1)
+  diag_H_mat_elem_erf += mo_nucl_elec_integral(occ(i,1),mo_nucl_elec_integral(i,1))
+  do j = i+1, n_occ_ab(1)
+   diag_H_mat_elem_erf += mo_bielec_integral_erf_jj_anti(occ(i,1),occ(j,1))
+  enddo
+ enddo
+
+ ! beta - beta 
+ do i = 1, n_occ_ab(2)
+  diag_H_mat_elem_erf += mo_nucl_elec_integral(occ(i,2),mo_nucl_elec_integral(i,2))
+  do j = i+1, n_occ_ab(2)
+   diag_H_mat_elem_erf += mo_bielec_integral_erf_jj_anti(occ(i,2),occ(j,2))
+  enddo
+ enddo
+
+ ! alpha - beta 
+ do i = 1, n_occ_ab(1)
+  do j = 1, n_occ_ab(2)
+   diag_H_mat_elem_erf += mo_bielec_integral_erf_jj(occ(i,1),occ(j,2))
+  enddo
+ enddo
+
+end
+
+
+
+subroutine i_H_j_erf_and_short_coulomb(key_i,key_j,Nint,hij)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns <i|H|j> where i and j are determinants
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: key_i(Nint,2), key_j(Nint,2)
+  double precision, intent(out)  :: hij
+  
+  integer                        :: exc(0:2,2,2)
+  integer                        :: degree
+  double precision               :: get_mo_bielec_integral_erf
+  integer                        :: m,n,p,q
+  integer                        :: i,j,k
+  integer                        :: occ(Nint*bit_kind_size,2)
+  double precision               :: diag_H_mat_elem_erf, phase,phase_2
+  integer                        :: n_occ_ab(2)
+  PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_map big_array_exchange_integrals_erf
+  
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+  ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
+  ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
+  ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
+  ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
+  
+  hij = 0.d0
+  !DIR$ FORCEINLINE
+  call get_excitation_degree(key_i,key_j,degree,Nint)
+  integer :: spin
+  select case (degree)
+    case (2)
+      call get_double_excitation(key_i,key_j,exc,phase,Nint)
+      if (exc(0,1,1) == 1) then
+        ! Mono alpha, mono beta
+        if(exc(1,1,1) == exc(1,2,2) )then
+         hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
+        else if (exc(1,2,1) ==exc(1,1,2))then
+         hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
+        else
+         hij = phase*get_mo_bielec_integral_erf(                          &
+             exc(1,1,1),                                              &
+             exc(1,1,2),                                              &
+             exc(1,2,1),                                              &
+             exc(1,2,2) ,mo_integrals_erf_map)
+        endif
+      else if (exc(0,1,1) == 2) then
+        ! Double alpha
+        hij = phase*(get_mo_bielec_integral_erf(                         &
+            exc(1,1,1),                                              &
+            exc(2,1,1),                                              &
+            exc(1,2,1),                                              &
+            exc(2,2,1) ,mo_integrals_erf_map) -                          &
+            get_mo_bielec_integral_erf(                                  &
+            exc(1,1,1),                                              &
+            exc(2,1,1),                                              &
+            exc(2,2,1),                                              &
+            exc(1,2,1) ,mo_integrals_erf_map) )
+      else if (exc(0,1,2) == 2) then
+        ! Double beta
+        hij = phase*(get_mo_bielec_integral_erf(                         &
+            exc(1,1,2),                                              &
+            exc(2,1,2),                                              &
+            exc(1,2,2),                                              &
+            exc(2,2,2) ,mo_integrals_erf_map) -                          &
+            get_mo_bielec_integral_erf(                                  &
+            exc(1,1,2),                                              &
+            exc(2,1,2),                                              &
+            exc(2,2,2),                                              &
+            exc(1,2,2) ,mo_integrals_erf_map) )
+      endif
+    case (1)
+      call get_mono_excitation(key_i,key_j,exc,phase,Nint)
+      !DIR$ FORCEINLINE
+      call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
+      if (exc(0,1,1) == 1) then
+        ! Mono alpha
+        m = exc(1,1,1)
+        p = exc(1,2,1)
+        spin = 1
+        do i = 1, n_occ_ab(1)
+         hij += -big_array_exchange_integrals_erf(occ(i,1),m,p) + big_array_coulomb_integrals_erf(occ(i,1),m,p)
+        enddo
+        do i = 1, n_occ_ab(2)
+         hij += big_array_coulomb_integrals_erf(occ(i,2),m,p)
+        enddo
+      else
+        ! Mono beta
+        m = exc(1,1,2)
+        p = exc(1,2,2)
+        spin = 2
+        do i = 1, n_occ_ab(2)
+         hij += -big_array_exchange_integrals_erf(occ(i,2),m,p) + big_array_coulomb_integrals_erf(occ(i,2),m,p)
+        enddo
+        do i = 1, n_occ_ab(1)
+         hij += big_array_coulomb_integrals_erf(occ(i,1),m,p)
+        enddo
+      endif
+      hij = hij + mo_nucl_elec_integral(m,p) + mo_kinetic_integral(m,p) + effective_short_range_operator(m,p)
+      hij = hij * phase
+    case (0)
+      hij = diag_H_mat_elem_erf(key_i,Nint)
+  end select
+end
+
+double precision function diag_H_mat_elem_erf_and_short_coulomb(key_i,Nint)
+ implicit none
+ integer(bit_kind), intent(in) :: key_i(N_int,2)
+ integer, intent(in)  :: Nint
+ integer :: i,j
+ integer                        :: occ(Nint*bit_kind_size,2)
+ integer                        :: n_occ_ab(2)
+
+ call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
+ diag_H_mat_elem_erf_and_short_coulomb = 0.d0
+ ! alpha - alpha
+ do i = 1, n_occ_ab(1)
+  diag_H_mat_elem_erf_and_short_coulomb += mo_nucl_elec_integral(occ(i,1),mo_nucl_elec_integral(i,1)) + mo_kinetic_integral(occ(i,1),mo_nucl_elec_integral(i,1)) &
+                                         + effective_short_range_operator(occ(i,1),occ(i,1))
+  do j = i+1, n_occ_ab(1)
+   diag_H_mat_elem_erf_and_short_coulomb += mo_bielec_integral_erf_jj_anti(occ(i,1),occ(j,1))
+  enddo
+ enddo
+
+ ! beta - beta 
+ do i = 1, n_occ_ab(2)
+  diag_H_mat_elem_erf_and_short_coulomb += mo_nucl_elec_integral(occ(i,2),mo_nucl_elec_integral(i,2)) + mo_kinetic_integral(occ(i,2),mo_nucl_elec_integral(i,2)) & 
+                                         + effective_short_range_operator(occ(i,2),occ(i,2))
+  do j = i+1, n_occ_ab(2)
+   diag_H_mat_elem_erf_and_short_coulomb += mo_bielec_integral_erf_jj_anti(occ(i,2),occ(j,2))
+  enddo
+ enddo
+
+ ! alpha - beta 
+ do i = 1, n_occ_ab(1)
+  do j = 1, n_occ_ab(2)
+   diag_H_mat_elem_erf_and_short_coulomb += mo_bielec_integral_erf_jj(occ(i,1),occ(j,2))
+  enddo
+ enddo
+
+end
+
+
+subroutine i_H_j_erf_component(key_i,key_j,Nint,hij_core,hij_hartree,hij_erf,hij_total)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns <i|H|j> where i and j are determinants
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: key_i(Nint,2), key_j(Nint,2)
+  double precision, intent(out)  :: hij_core
+  double precision, intent(out)  :: hij_hartree
+  double precision, intent(out)  :: hij_erf
+  double precision, intent(out)  :: hij_total
+  
+  integer                        :: exc(0:2,2,2)
+  integer                        :: degree
+  double precision               :: get_mo_bielec_integral_erf
+  integer                        :: m,n,p,q
+  integer                        :: i,j,k
+  integer                        :: occ(Nint*bit_kind_size,2)
+  double precision               :: diag_H_mat_elem_erf, phase,phase_2
+  integer                        :: n_occ_ab(2)
+  PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_map big_array_exchange_integrals_erf
+  
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+  ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
+  ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
+  ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
+  ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
+  
+  hij_core = 0.d0
+  hij_hartree = 0.d0
+  hij_erf = 0.d0
+ 
+  !DIR$ FORCEINLINE
+  call get_excitation_degree(key_i,key_j,degree,Nint)
+  integer :: spin
+  select case (degree)
+    case (2)
+      call get_double_excitation(key_i,key_j,exc,phase,Nint)
+      if (exc(0,1,1) == 1) then
+        ! Mono alpha, mono beta
+        if(exc(1,1,1) == exc(1,2,2) )then
+         hij_erf = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
+        else if (exc(1,2,1) ==exc(1,1,2))then
+         hij_erf = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
+        else
+         hij_erf = phase*get_mo_bielec_integral_erf(                          &
+             exc(1,1,1),                                              &
+             exc(1,1,2),                                              &
+             exc(1,2,1),                                              &
+             exc(1,2,2) ,mo_integrals_erf_map)
+        endif
+      else if (exc(0,1,1) == 2) then
+        ! Double alpha
+        hij_erf = phase*(get_mo_bielec_integral_erf(                         &
+            exc(1,1,1),                                              &
+            exc(2,1,1),                                              &
+            exc(1,2,1),                                              &
+            exc(2,2,1) ,mo_integrals_erf_map) -                          &
+            get_mo_bielec_integral_erf(                                  &
+            exc(1,1,1),                                              &
+            exc(2,1,1),                                              &
+            exc(2,2,1),                                              &
+            exc(1,2,1) ,mo_integrals_erf_map) )
+      else if (exc(0,1,2) == 2) then
+        ! Double beta
+        hij_erf = phase*(get_mo_bielec_integral_erf(                         &
+            exc(1,1,2),                                              &
+            exc(2,1,2),                                              &
+            exc(1,2,2),                                              &
+            exc(2,2,2) ,mo_integrals_erf_map) -                          &
+            get_mo_bielec_integral_erf(                                  &
+            exc(1,1,2),                                              &
+            exc(2,1,2),                                              &
+            exc(2,2,2),                                              &
+            exc(1,2,2) ,mo_integrals_erf_map) )
+      endif
+    case (1)
+      call get_mono_excitation(key_i,key_j,exc,phase,Nint)
+      !DIR$ FORCEINLINE
+      call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
+      if (exc(0,1,1) == 1) then
+        ! Mono alpha
+        m = exc(1,1,1)
+        p = exc(1,2,1)
+        spin = 1
+        do i = 1, n_occ_ab(1)
+         hij_erf += -big_array_exchange_integrals_erf(occ(i,1),m,p) + big_array_coulomb_integrals_erf(occ(i,1),m,p)
+        enddo
+        do i = 1, n_occ_ab(2)
+         hij_erf += big_array_coulomb_integrals_erf(occ(i,2),m,p)
+        enddo
+      else
+        ! Mono beta
+        m = exc(1,1,2)
+        p = exc(1,2,2)
+        spin = 2
+        do i = 1, n_occ_ab(2)
+         hij_erf += -big_array_exchange_integrals_erf(occ(i,2),m,p) + big_array_coulomb_integrals_erf(occ(i,2),m,p)
+        enddo
+        do i = 1, n_occ_ab(1)
+         hij_erf += big_array_coulomb_integrals_erf(occ(i,1),m,p)
+        enddo
+      endif
+      hij_core = mo_nucl_elec_integral(m,p) + mo_kinetic_integral(m,p) 
+      hij_hartree =  effective_short_range_operator(m,p)
+      hij_total = (hij_erf + hij_core + hij_hartree) * phase
+    case (0)
+      call diag_H_mat_elem_erf_component(key_i,hij_core,hij_hartree,hij_erf,hij_total,Nint)
+  end select
+end
+
+subroutine diag_H_mat_elem_erf_component(key_i,hij_core,hij_hartree,hij_erf,hij_total,Nint)
+ implicit none
+ integer(bit_kind), intent(in) :: key_i(N_int,2)
+ integer, intent(in)  :: Nint
+ double precision, intent(out)  :: hij_core
+ double precision, intent(out)  :: hij_hartree
+ double precision, intent(out)  :: hij_erf
+ double precision, intent(out)  :: hij_total
+ integer :: i,j
+ integer                        :: occ(Nint*bit_kind_size,2)
+ integer                        :: n_occ_ab(2)
+
+ call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
+ hij_core = 0.d0
+ hij_hartree = 0.d0
+ hij_erf = 0.d0
+ ! alpha - alpha
+ do i = 1, n_occ_ab(1)
+  hij_core += mo_nucl_elec_integral(occ(i,1),mo_nucl_elec_integral(i,1)) + mo_kinetic_integral(occ(i,1),mo_nucl_elec_integral(i,1))  
+  hij_hartree += effective_short_range_operator(occ(i,1),occ(i,1))
+  do j = i+1, n_occ_ab(1)
+   hij_erf += mo_bielec_integral_erf_jj_anti(occ(i,1),occ(j,1))
+  enddo
+ enddo
+
+ ! beta - beta 
+ do i = 1, n_occ_ab(2)
+  hij_core += mo_nucl_elec_integral(occ(i,2),mo_nucl_elec_integral(i,2)) + mo_kinetic_integral(occ(i,2),mo_nucl_elec_integral(i,2))   
+  hij_hartree += effective_short_range_operator(occ(i,2),occ(i,2))
+  do j = i+1, n_occ_ab(2)
+   hij_erf += mo_bielec_integral_erf_jj_anti(occ(i,2),occ(j,2))
+  enddo
+ enddo
+
+ ! alpha - beta 
+ do i = 1, n_occ_ab(1)
+  do j = 1, n_occ_ab(2)
+   hij_erf += mo_bielec_integral_erf_jj(occ(i,1),occ(j,2))
+  enddo
+ enddo
+ hij_total = hij_erf + hij_hartree + hij_core
+
+end
+
+

src/Integrals_Monoelec/EZFIO.cfg

@@ -4,6 +4,14 @@ doc: Read/Write MO one-electron integrals from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
 default: None
 
+
+[disk_access_only_mo_one_integrals]
+type: Disk_access
+doc: Read/Write MO for only the total one-electron integrals which can be anything [ Write | Read | None ] 
+interface: ezfio,provider,ocaml
+default: None
+
+
 [disk_access_ao_one_integrals]
 type: Disk_access
 doc: Read/Write AO one-electron integrals from/to disk [ Write | Read | None ] 

src/Integrals_Monoelec/mo_mono_ints.irp.f

@@ -6,7 +6,8 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to
   ! sum of the kinetic and nuclear electronic potential
   END_DOC
   print*,'Providing the mono electronic integrals'
-  if (read_mo_one_integrals) then
+  if (read_only_mo_one_integrals) then
+    print*, 'Reading the mono electronic integrals from disk'
     call read_one_e_integrals('mo_one_integral', mo_mono_elec_integral,      &
         size(mo_mono_elec_integral,1), size(mo_mono_elec_integral,2))
     print *,  'MO N-e integrals read from disk'
@@ -14,7 +15,7 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to
    do j = 1, mo_tot_num
      do i = 1, mo_tot_num
        mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + &
-          mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j)
+                                    mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j)
      enddo
    enddo
   endif

src/Integrals_Monoelec/read_write.irp.f

@@ -1,5 +1,6 @@
  BEGIN_PROVIDER [ logical, read_ao_one_integrals ]
 &BEGIN_PROVIDER [ logical, read_mo_one_integrals ]
+&BEGIN_PROVIDER [ logical, read_only_mo_one_integrals ]
 &BEGIN_PROVIDER [ logical, write_ao_one_integrals ]
 &BEGIN_PROVIDER [ logical, write_mo_one_integrals ]
    
@@ -21,10 +22,14 @@
      write_ao_one_integrals = .False.
      
    else
-     print *, 'bielec_integrals/disk_access_ao_integrals has a wrong type'
+     print *, 'monoelec_integrals/disk_access_ao_integrals has a wrong type'
      stop 1
      
    endif
+ 
+   if (disk_access_only_mo_one_integrals.EQ.'Read')then
+     read_only_mo_one_integrals = .True.
+   endif
    
    if (disk_access_mo_one_integrals.EQ.'Read') then
      read_mo_one_integrals =  .True.
@@ -39,7 +44,7 @@
      write_mo_one_integrals = .False.
      
    else
-     print *, 'bielec_integrals/disk_access_mo_integrals has a wrong type'
+     print *, 'monoelec_integrals/disk_access_mo_integrals has a wrong type'
      stop 1
      
    endif