Added files

TODO

@@ -41,7 +41,6 @@
   * Examples : subroutine example_module
 
 * Config file for Cray
-* Noms des modules et fichiers en minuscule
 * Implicit none obligatoire
 
 * Doc des executables

src/cis/cis.irp.f

@@ -1,5 +1,8 @@
 program cis
   implicit none
+  BEGIN_DOC
+! Configuration Interaction with Single excitations.
+  END_DOC
   read_wf = .False.
   SOFT_TOUCH read_wf
   call run

src/cisd/cisd.irp.f

@@ -1,5 +1,8 @@
-program cis
+program cisd
   implicit none
+  BEGIN_DOC
+! Configuration Interaction with Single and Double excitations.
+  END_DOC
   read_wf = .False.
   SOFT_TOUCH read_wf
   call run

src/determinants/density_matrix.irp.f

@@ -489,73 +489,3 @@ END_PROVIDER
 END_PROVIDER
 
 
- BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_old, (mo_tot_num,mo_tot_num,N_states) ]
-&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_old, (mo_tot_num,mo_tot_num,N_states) ]
-   implicit none
-   BEGIN_DOC
-   ! Alpha and beta one-body density matrix for each state
-   END_DOC
-
-   integer                        :: j,k,l,m
-   integer                        :: occ(N_int*bit_kind_size,2)
-   double precision               :: ck, cl, ckl
-   double precision               :: phase
-   integer                        :: h1,h2,p1,p2,s1,s2, degree
-   integer                        :: exc(0:2,2,2),n_occ(2)
-   double precision, allocatable  :: tmp_a(:,:,:), tmp_b(:,:,:)
-
-     one_body_dm_mo_alpha_old = 0.d0
-     one_body_dm_mo_beta_old  = 0.d0
-     !$OMP PARALLEL DEFAULT(NONE)                                         &
-        !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
-        !$OMP  tmp_a, tmp_b, n_occ)&
-        !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
-        !$OMP  elec_beta_num,one_body_dm_mo_alpha_old,one_body_dm_mo_beta_old,N_det,&
-        !$OMP  mo_tot_num)
-     allocate(tmp_a(mo_tot_num,mo_tot_num,N_states), tmp_b(mo_tot_num,mo_tot_num,N_states) )
-     tmp_a = 0.d0
-     tmp_b = 0.d0
-     !$OMP DO SCHEDULE(dynamic)
-     do k=1,N_det
-       call bitstring_to_list_ab(psi_det(1,1,k), occ, n_occ, N_int)
-       do m=1,N_states
-         ck = psi_coef(k,m)*psi_coef(k,m) 
-         do l=1,elec_alpha_num
-           j = occ(l,1)
-           tmp_a(j,j,m) += ck
-         enddo
-         do l=1,elec_beta_num
-           j = occ(l,2)
-           tmp_b(j,j,m) += ck
-         enddo
-       enddo
-       do l=1,k-1
-         call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int)
-         if (degree /= 1) then
-           cycle
-         endif
-         call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
-         call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-         do m=1,N_states
-           ckl = psi_coef(k,m) * psi_coef(l,m) * phase 
-           if (s1==1) then
-             tmp_a(h1,p1,m) += ckl
-             tmp_a(p1,h1,m) += ckl
-           else
-             tmp_b(h1,p1,m) += ckl
-             tmp_b(p1,h1,m) += ckl
-           endif
-         enddo
-       enddo
-     enddo
-     !$OMP END DO NOWAIT
-     !$OMP CRITICAL
-     one_body_dm_mo_alpha_old(:,:,:) = one_body_dm_mo_alpha_old(:,:,:) + tmp_a(:,:,:)
-     !$OMP END CRITICAL
-     !$OMP CRITICAL
-     one_body_dm_mo_beta_old(:,:,:)  = one_body_dm_mo_beta_old(:,:,:)  + tmp_b(:,:,:)
-     !$OMP END CRITICAL
-     deallocate(tmp_a,tmp_b)
-     !$OMP END PARALLEL
-
-END_PROVIDER

src/ezfio_files/NEED

@@ -0,0 +1 @@
+mpi

src/ezfio_files/README.rst

@@ -0,0 +1,8 @@
+===========
+Ezfio_files
+===========
+
+This modules essentially contains the name of the |EZFIO| directory in the
+:c:data:`ezfio_filename` variable. This is read as the first argument of the
+command-line, or as the :envvar:`QP_INPUT` environment variable.
+

src/ezfio_files/ezfio.irp.f

@@ -0,0 +1,53 @@
+BEGIN_PROVIDER [ character*(128), ezfio_filename ]
+  implicit none
+  BEGIN_DOC
+  ! Name of EZFIO file. It is obtained from the QPACKAGE_INPUT environment
+  ! variable if it is set, or as the 1st argument of the command line.
+  END_DOC
+
+  PROVIDE mpi_initialized
+  
+  ! Get the QPACKAGE_INPUT environment variable
+  call getenv('QPACKAGE_INPUT',ezfio_filename)
+  if (ezfio_filename == '') then
+    ! Get from the command line
+    integer                        :: iargc
+    call getarg(0,ezfio_filename)
+    if (iargc() /= 1) then
+      print *, 'Missing EZFIO file name in the command line:'
+      print *, trim(ezfio_filename)//' <ezfio_file>'
+      stop 1
+    endif
+    call getarg(1,ezfio_filename)
+  endif
+
+  ! Check that file exists
+  logical :: exists
+  inquire(file=trim(ezfio_filename)//'/ezfio/creation',exist=exists)
+  if (.not.exists) then
+    print *, 'Error: file '//trim(ezfio_filename)//' does not exist'
+    stop 1
+  endif
+
+  call ezfio_set_file(ezfio_filename)
+
+  ! Adjust out-of-memory killer flag such that the current process will be
+  ! killed first by the OOM killer, allowing compute nodes to survive 
+  integer :: getpid
+  character*(64) :: command, pidc
+  write(pidc,*) getpid()
+  write(command,*) 'echo 15 > /proc//'//trim(adjustl(pidc))//'/oom_adj'
+  call system(command)
+
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ character*(128), ezfio_work_dir ]
+ implicit none
+ BEGIN_DOC
+ ! EZFIO/work/
+ END_DOC
+ call ezfio_set_work_empty(.False.)
+ ezfio_work_dir = trim(ezfio_filename)//'/work/'
+END_PROVIDER
+

src/ezfio_files/get_unit_and_open.irp.f

@@ -0,0 +1,58 @@
+integer function getUnitAndOpen(f,mode)
+  implicit none
+
+  BEGIN_DOC
+!  :f:
+!     file name
+!
+!  :mode:
+!     'R' : READ, UNFORMATTED
+!     'W' : WRITE, UNFORMATTED
+!     'r' : READ, FORMATTED
+!     'w' : WRITE, FORMATTED
+!     'a' : APPEND, FORMATTED
+!     'x' : READ/WRITE, FORMATTED
+!
+  END_DOC
+
+  character*(*)      :: f
+  character*(128)    :: new_f
+  integer            :: iunit
+  logical            :: is_open, exists
+  character          :: mode
+
+  is_open = .True.
+  iunit = 10
+  new_f = f
+  do while (is_open)
+    inquire(unit=iunit,opened=is_open)
+    if (.not.is_open) then
+      getUnitAndOpen = iunit
+    endif
+    iunit = iunit+1
+  enddo
+  if (mode.eq.'r') then
+    inquire(file=f,exist=exists)
+    if (.not.exists) then
+      open(unit=getUnitAndOpen,file=f,status='NEW',action='WRITE',form='FORMATTED')
+      close(unit=getUnitAndOpen)
+    endif
+    open(unit=getUnitAndOpen,file=f,status='OLD',action='READ',form='FORMATTED')
+  else if (mode.eq.'R') then
+    inquire(file=f,exist=exists)
+    if (.not.exists) then
+      open(unit=getUnitAndOpen,file=f,status='NEW',action='WRITE',form='UNFORMATTED')
+      close(unit=getUnitAndOpen)
+    endif
+    open(unit=getUnitAndOpen,file=f,status='OLD',action='READ',form='UNFORMATTED')
+  else if (mode.eq.'W') then
+    open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='WRITE',form='UNFORMATTED')
+  else if (mode.eq.'w') then
+    open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='WRITE',form='FORMATTED')
+  else if (mode.eq.'a') then
+    open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='WRITE',position='APPEND',form='FORMATTED')
+  else if (mode.eq.'x') then
+    open(unit=getUnitAndOpen,file=new_f,form='FORMATTED')
+  endif
+end function getUnitAndOpen
+

src/ezfio_files/output.irp.f

@@ -0,0 +1,85 @@
+ BEGIN_PROVIDER [ double precision, output_wall_time_0 ]
+&BEGIN_PROVIDER [ double precision, output_cpu_time_0 ]
+  implicit none
+  BEGIN_DOC
+  ! Initial CPU and wall times when printing in the output files
+  END_DOC
+  call cpu_time(output_wall_time_0)
+  call wall_time(output_wall_time_0)
+END_PROVIDER
+ 
+ 
+subroutine write_time(iunit)
+  implicit none
+  BEGIN_DOC
+  ! Write a time stamp in the output for chronological reconstruction
+  END_DOC
+  integer, intent(in)            :: iunit
+  double precision               :: wt, ct
+  if (.not.mpi_master) then
+    return
+  endif
+  write(6,*)
+  call print_memory_usage()
+  call cpu_time(ct)
+  call wall_time(wt)
+  write(6,'(A,F15.6,A,F15.6,A)') &
+    '.. >>>>> [ WALL TIME: ', wt-output_wall_time_0, &
+        ' s ] [ CPU  TIME: ', ct-output_cpu_time_0, ' s ] <<<<< ..'
+  write(6,*)
+end
+ 
+subroutine write_double(iunit,value,label)
+  implicit none
+  BEGIN_DOC
+  ! Write a double precision value in output
+  END_DOC
+  if (.not.mpi_master) then
+    return
+  endif
+  integer, intent(in)            :: iunit
+  double precision               :: value
+  character*(*)                  :: label
+  character*(64), parameter      :: f = '(A50,G24.16)'
+  character*(50)                 :: newlabel
+  write(newlabel,'(A,A)') '* ',trim(label)
+  write(6,f) newlabel, value
+end
+
+
+subroutine write_int(iunit,value,label)
+  implicit none
+  BEGIN_DOC
+  ! Write an integer value in output
+  END_DOC
+  if (.not.mpi_master) then
+    return
+  endif
+  integer, intent(in)            :: iunit
+  integer                        :: value
+  character*(*)                  :: label
+  character*(64), parameter      :: f = '(A50,I16)'
+  character*(50)                 :: newlabel
+  write(newlabel,'(A,A)') '* ',trim(label)
+  write(6,f) newlabel, value
+end
+
+
+subroutine write_bool(iunit,value,label)
+  implicit none
+  BEGIN_DOC
+  ! Write an logical value in output
+  END_DOC
+  if (.not.mpi_master) then
+    return
+  endif
+  integer, intent(in)            :: iunit
+  logical                        :: value
+  character*(*)                  :: label
+  character*(64), parameter      :: f = '(A50,L1)'
+  character*(50)                 :: newlabel
+  write(newlabel,'(A,A)') '* ',trim(label)
+  write(6,f) newlabel, value
+end
+
+

src/ezfio_files/q_package.ezfio_config

@@ -0,0 +1,6 @@
+work
+  empty     logical
+
+save
+  empty     logical
+

src/fci/fci.irp.f

@@ -1,5 +1,8 @@
-program fci_zmq
+program fci
   implicit none
+  BEGIN_DOC
+! Selected Full Configuration Interaction.
+  END_DOC
   integer                        :: i,j,k
   double precision, allocatable  :: pt2(:), variance(:), norm(:), rpt2(:)
   integer                        :: n_det_before, to_select

src/fci/pt2.irp.f

@@ -1,5 +1,9 @@
-program pt2_stoch
+program pt2
   implicit none
+  BEGIN_DOC
+  ! Second order perturbative correction to the wave function contained in the
+  ! EZFIO directory.
+  END_DOC
   read_wf = .True.
   SOFT_TOUCH read_wf
   PROVIDE mo_bielec_integrals_in_map

src/four_idx/four_idx_transform.irp.f

@@ -1,7 +1,7 @@
-program four_idx
+program four_idx_transform
   implicit none
   BEGIN_DOC
-! 4-index transformation from AO to MO integrals
+! 4-index transformation of two-electron integrals from AO to MO integrals
   END_DOC
   
   disk_access_mo_integrals = 'Write'

src/generators_cas/NEED

@@ -0,0 +1 @@
+determinants

src/generators_cas/README.rst

@@ -0,0 +1,12 @@
+==============
+Generators_CAS
+==============
+
+Module defining the generator determinants as those belonging to a |CAS|.
+The |MOs| belonging to the |CAS| are those which were set as active with
+the :ref:`qp_set_mo_class` command.
+
+This module is intended to be included in the :file:`NEED` file to define
+generators on a |CAS|.
+
+

src/generators_cas/generators.irp.f

@@ -0,0 +1,83 @@
+use bitmasks
+
+BEGIN_PROVIDER [ integer, N_det_generators ]
+  implicit none
+  BEGIN_DOC
+  ! Number of generator detetrminants
+  END_DOC
+  integer                        :: i,k,l
+  logical                        :: good
+  integer, external :: number_of_holes,number_of_particles
+  call write_time(6)
+  N_det_generators = 0
+  do i=1,N_det
+    good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 )
+    if (good) then
+      N_det_generators += 1
+    endif
+  enddo
+  N_det_generators = max(N_det_generators,1)
+  call write_int(6,N_det_generators,'Number of generators')
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_order, (psi_det_size) ]
+  implicit none
+  BEGIN_DOC
+  ! For Single reference wave functions, the generator is the
+  ! Hartree-Fock determinant
+  END_DOC
+  integer                        :: i, k, l, m
+  logical                        :: good
+  integer, external :: number_of_holes,number_of_particles
+  integer, allocatable :: nongen(:)
+  integer :: inongen
+
+  allocate(nongen(N_det))
+
+  inongen = 0
+  m=0
+  do i=1,N_det
+    good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 )
+    if (good) then
+      m = m+1
+      psi_det_sorted_gen_order(i) = m
+      do k=1,N_int
+        psi_det_generators(k,1,m) = psi_det_sorted(k,1,i)
+        psi_det_generators(k,2,m) = psi_det_sorted(k,2,i)
+      enddo
+      psi_coef_generators(m,:) = psi_coef_sorted(i,:)
+    else
+      inongen += 1
+      nongen(inongen) = i
+    endif
+  enddo
+  
+  psi_det_sorted_gen(:,:,:N_det_generators) = psi_det_generators(:,:,:N_det_generators)
+  psi_coef_sorted_gen(:N_det_generators, :) = psi_coef_generators(:N_det_generators, :)
+  do i=1,inongen
+    psi_det_sorted_gen_order(nongen(i)) = N_det_generators+i
+    psi_det_sorted_gen(:,:,N_det_generators+i) = psi_det_sorted(:,:,nongen(i))
+    psi_coef_sorted_gen(N_det_generators+i, :) = psi_coef_sorted(nongen(i),:)
+  end do
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, size_select_max]
+  implicit none
+  BEGIN_DOC
+  ! Size of the select_max array
+  END_DOC
+  size_select_max = 10000
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ]
+  implicit none
+  BEGIN_DOC
+  ! Memo to skip useless selectors
+  END_DOC
+  select_max = huge(1.d0)
+END_PROVIDER
+

src/generators_full/NEED

@@ -0,0 +1,2 @@
+determinants
+hartree_fock

src/generators_full/README.rst

@@ -0,0 +1,9 @@
+===============
+Generators_full
+===============
+
+Module defining the generator determinants as all the determinants of the
+variational space.
+
+This module is intended to be included in the :file:`NEED` file to define
+a full set of generators.

src/generators_full/generators.irp.f

@@ -0,0 +1,82 @@
+use bitmasks
+
+BEGIN_PROVIDER [ integer, N_det_generators ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the number of generators is 1 : the
+ ! Hartree-Fock determinant
+ END_DOC
+ integer :: i
+ double precision :: norm
+ call write_time(6)
+ norm = 1.d0
+ N_det_generators = N_det
+ do i=1,N_det
+   norm = norm - psi_average_norm_contrib_sorted(i)
+   if (norm - 1.d-10 < 1.d0 - threshold_generators) then
+     N_det_generators = i
+     exit
+   endif
+ enddo
+ N_det_generators = max(N_det_generators,1)
+ call write_int(6,N_det_generators,'Number of generators')
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det)
+ psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted(1:N_det,1:N_states)
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_order,     (psi_det_size)  ]
+
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_det_sorted_gen = psi_det_sorted
+ psi_coef_sorted_gen = psi_coef_sorted
+ psi_det_sorted_gen_order = psi_det_sorted_order
+END_PROVIDER
+
+
+BEGIN_PROVIDER [integer, degree_max_generators]
+ implicit none
+ BEGIN_DOC
+! Max degree of excitation (respect to HF) of the generators
+ END_DOC
+ integer :: i,degree
+  degree_max_generators = 0
+  do i = 1, N_det_generators
+   call get_excitation_degree(HF_bitmask,psi_det_generators(1,1,i),degree,N_int)
+   if(degree .gt. degree_max_generators)then
+    degree_max_generators = degree
+   endif
+  enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ integer, size_select_max]
+ implicit none
+ BEGIN_DOC
+ ! Size of the select_max array
+ END_DOC
+ size_select_max = 10000
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ]
+ implicit none
+ BEGIN_DOC
+ ! Memo to skip useless selectors
+ END_DOC
+ select_max = huge(1.d0)
+END_PROVIDER
+

src/integrals_bielec/EZFIO.cfg

@@ -0,0 +1,55 @@
+[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]
+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
+interface: ezfio,provider,ocaml
+default: 1.e-15
+ezfio_name: threshold_ao
+
+[mo_integrals_threshold]
+type: Threshold
+doc: If | <ij|kl> | < `mo_integrals_threshold` then <ij|kl> is zero
+interface: ezfio,provider,ocaml
+default: 1.e-15
+ezfio_name: threshold_mo
+
+
+[no_vvvv_integrals]
+type: logical
+doc: If `True`, computes all integrals except for the integrals having 4 virtual indices
+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 does not compute the integrals with 3 virtual indices 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 does not compute the integrals with 3 virtual orbitals
+interface: ezfio,provider,ocaml
+default: False
+ezfio_name: no_vvv_integrals
+
+[do_direct_integrals]
+type: logical
+doc: Compute integrals on the fly (very slow, only for debugging)
+interface: ezfio,provider,ocaml
+default: False
+ezfio_name: direct
+

src/integrals_bielec/NEED

@@ -0,0 +1,7 @@
+ao_one_e_integrals
+mo_one_e_integrals
+pseudo
+bitmask
+zmq
+ao_basis
+mo_basis

src/integrals_bielec/README.rst

@@ -0,0 +1,21 @@
+================
+Integrals_Bielec
+================
+
+Here, all two-electron integrals (:math:`1/r_{12}`) are computed.
+As they have 4 indices and many are zero, they are stored in a map, as defined
+in :file:`Utils/map_module.f90`.
+
+To fetch an |AO| integral, use the
+`get_ao_bielec_integral(i,j,k,l,ao_integrals_map)` function, and
+to fetch an |MO| integral, use
+`get_mo_bielec_integral(i,j,k,l,mo_integrals_map)` or
+`mo_bielec_integral(i,j,k,l)`.
+
+The conventions are:
+
+* For |AO| integrals : (ik|jl) = (11|22)
+* For |MO| integrals : <ij|kl> = <12|12>
+
+
+

src/integrals_bielec/ao_bielec_integrals_in_map_slave.irp.f

@@ -0,0 +1,244 @@
+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
+
+IRP_IF ZMQ_PUSH
+IRP_ELSE
+  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
+IRP_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()
+
+  integer, external :: connect_to_taskserver
+  if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then
+    call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
+    return
+  endif
+
+  zmq_socket_push      = new_zmq_push_socket(thread)
+
+  allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) )
+
+
+  do 
+    integer, external :: get_task_from_taskserver
+    if (get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task) == -1) then
+      exit
+    endif
+    if (task_id == 0) exit
+    read(task,*) j, l
+    integer, external :: task_done_to_taskserver
+    call compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) 
+    if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) == -1) then
+        stop 'Unable to send task_done'
+    endif
+    call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id)
+  enddo
+
+  integer, external :: disconnect_from_taskserver
+  if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then
+    continue
+  endif
+  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(zmq_socket_pull)
+  use map_module
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Collects results from the AO integral calculation
+  END_DOC
+
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
+  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*8                      :: control, accu, sze
+  integer                        :: task_id, more
+
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_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)
+
+IRP_IF ZMQ_PUSH
+IRP_ELSE
+      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
+IRP_ENDIF
+
+      
+      call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
+      accu += n_integrals
+      if (task_id /= 0) then
+        integer, external :: zmq_delete_task
+        if (zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) == -1) then
+          stop 'Unable to delete task'
+        endif
+      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)
+
+end
+

src/integrals_bielec/core_quantities.irp.f

@@ -0,0 +1,35 @@
+BEGIN_PROVIDER [double precision, core_energy]
+ implicit none
+ integer :: i,j,k,l
+ core_energy = 0.d0
+ do i = 1, n_core_orb
+  j = list_core(i)
+  core_energy += 2.d0 * mo_mono_elec_integral(j,j) + mo_bielec_integral_jj(j,j)
+  do k = i+1, n_core_orb
+   l = list_core(k)
+   core_energy += 2.d0 * (2.d0 * mo_bielec_integral_jj(j,l) - mo_bielec_integral_jj_exchange(j,l))
+  enddo
+ enddo
+ core_energy += nuclear_repulsion
+
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, core_fock_operator, (mo_tot_num,mo_tot_num)]
+ implicit none
+ integer :: i,j,k,l,m,n
+ double precision :: get_mo_bielec_integral
+ BEGIN_DOC 
+! this is the contribution to the Fock operator from the core electrons
+ END_DOC
+ core_fock_operator = 0.d0
+ do i = 1, n_act_orb
+  j = list_act(i)
+  do k = 1, n_act_orb
+   l = list_act(k)
+   do m = 1, n_core_orb
+    n = list_core(m)
+    core_fock_operator(j,l) += 2.d0 * get_mo_bielec_integral(j,n,l,n,mo_integrals_map) - get_mo_bielec_integral(j,n,n,l,mo_integrals_map)
+   enddo
+  enddo
+ enddo
+END_PROVIDER

src/integrals_bielec/gauss_legendre.irp.f

@@ -0,0 +1,57 @@
+ BEGIN_PROVIDER [ double precision, gauleg_t2, (n_pt_max_integrals,n_pt_max_integrals/2) ]
+&BEGIN_PROVIDER [ double precision, gauleg_w, (n_pt_max_integrals,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
+ 

src/integrals_bielec/integrals_3_index.irp.f

@@ -0,0 +1,22 @@
+ BEGIN_PROVIDER [double precision, big_array_coulomb_integrals, (mo_tot_num,mo_tot_num, mo_tot_num)]
+&BEGIN_PROVIDER [double precision, big_array_exchange_integrals,(mo_tot_num,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 

src/integrals_bielec/map_integrals.irp.f

@@ -0,0 +1,694 @@
+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+shiftr(r*r-r,1)
+  q = min(j,l)
+  s = max(j,l)
+  q = q+shiftr(s*s-s,1)
+  i1 = min(p,q)
+  i2 = max(p,q)
+  i1 = i1+shiftr(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 - shiftr(i2*i2-i2,1)
+  k(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i3)+1.d0)-1.d0))
+  j(1) = int(i2 - shiftr(l(1)*l(1)-l(1),1),4)
+  i(1) = int(i3 - shiftr(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(ii), k(ii), l(ii)
+        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( shiftl(ii,6), k-ao_integrals_cache_min)
+         ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
+         ii = ior( shiftl(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)
+    else
+      ii = l-ao_integrals_cache_min
+      ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
+      ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
+      ii = ior( shiftl(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: ', sze
+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*4, mo_integrals_cache_min ]
+&BEGIN_PROVIDER [ integer*4, mo_integrals_cache_max ]
+&BEGIN_PROVIDER [ integer*8, mo_integrals_cache_min_8 ]
+&BEGIN_PROVIDER [ integer*8, mo_integrals_cache_max_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_min_8 = max(1_8,elec_alpha_num - 63_8)
+ mo_integrals_cache_max_8 = min(int(mo_tot_num,8),mo_integrals_cache_min_8+127_8)
+ mo_integrals_cache_min   = max(1,elec_alpha_num - 63)
+ mo_integrals_cache_max   = min(mo_tot_num,mo_integrals_cache_min+127)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:128_8*128_8*128_8*128_8) ]
+ implicit none
+ BEGIN_DOC
+ ! Cache of MO integrals for fast access
+ END_DOC
+ PROVIDE mo_bielec_integrals_in_map
+ integer*8                      :: i,j,k,l
+ integer*4                      :: i4,j4,k4,l4
+ integer*8                      :: ii
+ integer(key_kind)              :: idx
+ real(integral_kind)            :: integral
+ FREE ao_integrals_cache
+ !$OMP PARALLEL DO PRIVATE (i,j,k,l,i4,j4,k4,l4,idx,ii,integral)
+ do l=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+   l4 = int(l,4)
+   do k=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+     k4 = int(k,4)
+     do j=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+       j4 = int(j,4)
+       do i=mo_integrals_cache_min_8,mo_integrals_cache_max_8
+         i4 = int(i,4)
+         !DIR$ FORCEINLINE
+         call bielec_integrals_index(i4,j4,k4,l4,idx)
+         !DIR$ FORCEINLINE
+         call map_get(mo_integrals_map,idx,integral)
+         ii = l-mo_integrals_cache_min_8
+         ii = ior( shiftl(ii,7), k-mo_integrals_cache_min_8)
+         ii = ior( shiftl(ii,7), j-mo_integrals_cache_min_8)
+         ii = ior( shiftl(ii,7), i-mo_integrals_cache_min_8)
+         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
+  integer*8                      :: ii_8
+  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, -128) /= 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_8 = int(l,8)-mo_integrals_cache_min_8
+    ii_8 = ior( shiftl(ii_8,7), int(k,8)-mo_integrals_cache_min_8)
+    ii_8 = ior( shiftl(ii_8,7), int(j,8)-mo_integrals_cache_min_8)
+    ii_8 = ior( shiftl(ii_8,7), int(i,8)-mo_integrals_cache_min_8)
+    get_mo_bielec_integral = mo_integrals_cache(ii_8)
+  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
+  PROVIDE mo_bielec_integrals_in_map
+  !DIR$ FORCEINLINE
+  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
+  double precision, external :: get_mo_bielec_integral
+  PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
+  
+  integer                        :: ii, ii0
+  integer*8                      :: ii_8, ii0_8
+  real(integral_kind)            :: tmp
+  integer(key_kind)              :: i1, idx
+  integer(key_kind)              :: p,q,r,s,i2
+  PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
+
+  ii0 = l-mo_integrals_cache_min
+  ii0 = ior(ii0, k-mo_integrals_cache_min)
+  ii0 = ior(ii0, j-mo_integrals_cache_min)
+
+  ii0_8 = int(l,8)-mo_integrals_cache_min_8
+  ii0_8 = ior( shiftl(ii0_8,7), int(k,8)-mo_integrals_cache_min_8)
+  ii0_8 = ior( shiftl(ii0_8,7), int(j,8)-mo_integrals_cache_min_8)
+
+  q = min(j,l)
+  s = max(j,l)
+  q = q+shiftr(s*s-s,1)
+
+  do i=1,sze
+    ii = ior(ii0, i-mo_integrals_cache_min)
+    if (iand(ii, -128) == 0) then
+      ii_8 = ior( shiftl(ii0_8,7), int(i,8)-mo_integrals_cache_min_8)
+      out_val(i) = mo_integrals_cache(ii_8)
+    else
+      p = min(i,k)
+      r = max(i,k)
+      p = p+shiftr(r*r-r,1)
+      i1 = min(p,q)
+      i2 = max(p,q)
+      idx = i1+shiftr(i2*i2-i2,1)
+      !DIR$ FORCEINLINE
+      call map_get(map,idx,tmp)
+      out_val(i) = dble(tmp)
+    endif
+  enddo
+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 (integral_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 (integral_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
+  if (.not.mpi_master) then
+    return
+  endif
+  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
+
+
+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

src/integrals_bielec/read_write.irp.f

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