beginning to clean erf integrals

plugins/Integrals_erf/EZFIO.cfg

@@ -32,16 +32,16 @@ doc: Read/Write MO integrals from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
 default: None
 
-[disk_access_ao_integrals_standard]
+[disk_access_ao_integrals_erf]
 type: Disk_access
-doc: Read/Write AO integrals_standard from/to disk [ Write | Read | None ] 
+doc: Read/Write AO integrals with the long range interaction from/to disk [ Write | Read | None ] 
 interface: ezfio,provider,ocaml
 default: None
 
 
-[disk_access_mo_integrals_standard]
+[disk_access_mo_integrals_erf]
 type: Disk_access
-doc: Read/Write MO integrals_standard from/to disk [ Write | Read | None ] 
+doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ] 
 interface: ezfio,provider,ocaml
 default: None
 
@@ -74,11 +74,3 @@ interface: ezfio,provider,ocaml
 default: 0.5
 ezfio_name: mu_erf
 
-
-[long_range]
-type: logical
-doc: if true, compute all the integrals using the long range interaction
-interface: ezfio,provider,ocaml
-default: False
-ezfio_name: long_range
-

plugins/Integrals_erf/map_integrals_long_range.irp.f

@@ -28,6 +28,7 @@ END_PROVIDER
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ]
+  use map_module
  implicit none
  BEGIN_DOC
  ! Cache of AO integrals for fast access
@@ -42,7 +43,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ]
      do j=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
        do i=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
          !DIR$ FORCEINLINE
-         call bielec_integrals_erf_index(i,j,k,l,idx)
+         call bielec_integrals_index(i,j,k,l,idx)
          !DIR$ FORCEINLINE
          call map_get(ao_integrals_erf_map,idx,integral)
          ii = l-ao_integrals_erf_cache_min
@@ -343,7 +344,21 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_erf_map ]
   print*, 'MO map initialized'
 END_PROVIDER
 
-subroutine insert_into_mo_integrals_map(n_integrals,                 &
+subroutine insert_into_ao_integrals_erf_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_erf_map, buffer_i, buffer_values, n_integrals)
+end
+
+subroutine insert_into_mo_integrals_erf_map(n_integrals,                 &
       buffer_i, buffer_values, thr)
   use map_module
   implicit none

plugins/Integrals_erf/mo_bi_integrals.irp.f

@@ -122,7 +122,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
   if (write_mo_integrals) then
     call ezfio_set_work_empty(.False.)
     call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
-    call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
+    call ezfio_set_integrals_erf_disk_access_mo_integrals("Read")
   endif
   
 END_PROVIDER
@@ -1065,8 +1065,6 @@ end
   ! mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij
   END_DOC
   implicit none
-  long_range = .false.
-  touch long_range
 ! call routine_mo_bielec_integral_jj_from_ao
 
 ! END_PROVIDER 
@@ -1199,8 +1197,6 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_exchange_from_ao, (mo_tot_num_align,mo_tot_num) ]
 &BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_anti_from_ao, (mo_tot_num_align,mo_tot_num) ]
   implicit none
-  long_range = .false.
-  touch long_range
   BEGIN_DOC
   ! mo_bielec_integral_vv_from_ao(i,j) = J_ij
   ! mo_bielec_integral_vv_exchange_from_ao(i,j) = J_ij

plugins/Integrals_erf/mo_bi_integrals_erf.irp.f

@@ -30,14 +30,14 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_erf_in_map ]
   END_DOC
   
   mo_bielec_integrals_erf_in_map = .True.
-! if (read_mo_integrals) then
+  if (read_mo_integrals_erf) then
-!   print*,'Reading the MO integrals'
+    print*,'Reading the MO integrals_erf'
-!   call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
+    call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
-!   print*, 'MO integrals provided'
+    print*, 'MO integrals_erf provided'
-!   return
+    return
-! else
+  else
-    PROVIDE ao_bielec_integrals_in_map
+    PROVIDE ao_bielec_integrals_erf_in_map
-! endif
+  endif
   
  !if(no_vvvv_integrals)then
  !  integer                        :: i,j,k,l
@@ -117,13 +117,13 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_erf_in_map ]
  !  endif
  !  
  !else
-    call add_integrals_to_map(full_ijkl_bitmask_4)
+    call add_integrals_erf_to_map(full_ijkl_bitmask_4)
- !endif
- !if (write_mo_integrals) then
- !  call ezfio_set_work_empty(.False.)
- !  call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
- !  call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
  !endif
+  if (write_mo_integrals_erf) then
+    call ezfio_set_work_empty(.False.)
+    call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
+    call ezfio_set_integrals_erf_disk_access_mo_integrals_erf("Read")
+  endif
   
 END_PROVIDER
 
@@ -160,7 +160,7 @@ subroutine add_integrals_erf_to_map(mask_ijkl)
   integer                        :: i2,i3,i4
   double precision,parameter     :: thr_coef = 1.d-10
   
-  PROVIDE ao_bielec_integrals_in_map  mo_coef
+  PROVIDE ao_bielec_integrals_erf_in_map  mo_coef
   
   !Get list of MOs for i,j,k and l
   !-------------------------------

plugins/Integrals_erf/providers_ao_erf.irp.f

@@ -25,22 +25,22 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ]
   integral = ao_bielec_integral_erf(1,1,1,1)
   
   real                           :: map_mb
-! PROVIDE read_ao_integrals disk_access_ao_integrals
+  PROVIDE read_ao_integrals_erf disk_access_ao_integrals_erf
-! if (read_ao_integrals) then
+  if (read_ao_integrals_erf) then
-!   print*,'Reading the AO integrals'
+    print*,'Reading the AO integrals_erf'
-!     call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
+      call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
-!     print*, 'AO integrals provided'
+      print*, 'AO integrals_erf provided'
-!     ao_bielec_integrals_in_map = .True.
+      ao_bielec_integrals_erf_in_map = .True.
-!     return
+      return
-! endif
+  endif
   
-  print*, 'Providing the AO integrals'
+  print*, 'Providing the AO integrals_erf'
   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')
+  call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals_erf')
 
   character(len=:), allocatable :: task
   allocate(character(len=ao_num*12) :: task)
@@ -63,7 +63,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ]
       endif
   !$OMP END PARALLEL
 
-  call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals')
+  call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals_erf')
 
 
   print*, 'Sorting the map'
@@ -81,11 +81,11 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ]
   
   ao_bielec_integrals_erf_in_map = .True.
 
-! if (write_ao_integrals) then
+  if (write_ao_integrals_erf) then
-!   call ezfio_set_work_empty(.False.)
+    call ezfio_set_work_empty(.False.)
-!   call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_erf_map)
+    call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
-!   call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
+    call ezfio_set_integrals_erf_disk_access_ao_integrals_erf("Read")
-! endif
+  endif
   
 END_PROVIDER
  

plugins/Integrals_erf/providers_ao_standard.irp.f

@@ -84,7 +84,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
   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_bielec_disk_access_ao_integrals("Read")
+    call ezfio_set_integrals_erf_disk_access_ao_integrals("Read")
   endif
   
 END_PROVIDER

plugins/Integrals_erf/read_write.irp.f

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

src/Determinants/EZFIO.cfg

@@ -119,3 +119,9 @@ doc: Maximum number of determinants for which the full H matrix is stored. Be ca
 interface: ezfio,provider,ocaml
 default: 90000
 
+[density_matrix_mo_disk]
+type: double precision
+doc: coefficient of the ith ao on the jth mo
+interface: ezfio
+size: (mo_basis.mo_tot_num,mo_basis.mo_tot_num)
+

src/Determinants/density_matrix.irp.f

@@ -15,6 +15,21 @@
    enddo
 END_PROVIDER
 
+
+subroutine save_density_matrix_mo
+ implicit none
+ double precision, allocatable :: dm(:,:)
+ allocate(dm(mo_tot_num,mo_tot_num))
+ integer :: i,j
+ do i = 1, mo_tot_num
+  do j = 1, mo_tot_num
+   dm(i,j) = one_body_dm_mo_alpha_average(i,j)
+  enddo
+ enddo
+ call ezfio_set_determinants_density_matrix_mo_disk(dm)
+
+end
+
  BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ]
  implicit none 
  integer :: i,j,ispin,istate

src/Integrals_Bielec/EZFIO.cfg

@@ -52,18 +52,3 @@ interface: ezfio,provider,ocaml
 default: 1.e-15
 ezfio_name: threshold_mo
 
-[mu_erf]
-type: double precision
-doc: cutting of the interaction in the range separated model
-interface: ezfio,provider,ocaml
-default: 0.5
-ezfio_name: mu_erf
-
-
-[long_range]
-type: logical
-doc: if true, compute all the integrals using the long range interaction
-interface: ezfio,provider,ocaml
-default: False
-ezfio_name: long_range
-

src/Integrals_Bielec/ao_bi_integrals.irp.f

@@ -158,7 +158,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
         q_inv = 1.d0/qq
         schwartz_kl(s,r) = general_primitive_integral(dim1,          &
             Q_new,Q_center,fact_q,qq,q_inv,iorder_q,                 &
-            Q_new,Q_center,fact_q,qq,q_inv,iorder_q)      &
+            Q_new,Q_center,fact_q,qq,q_inv,iorder_q)                 &
             * coef2 
         schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r))
       enddo
@@ -471,15 +471,9 @@ double precision function general_primitive_integral(dim,            &
   
   ! Gaussian Product
   ! ----------------
-  double precision :: p_plus_q
-  if(long_range)then
-   p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
-  else 
-   p_plus_q = (p+q)
-  endif
-  pq = p_inv*0.5d0*q_inv
   
-  pq_inv = 0.5d0/p_plus_q
+  pq = p_inv*0.5d0*q_inv
+  pq_inv = 0.5d0/(p+q)
   p10_1 = q*pq  ! 1/(2p)
   p01_1 = p*pq  ! 1/(2q)
   pq_inv_2 = pq_inv+pq_inv
@@ -554,7 +548,7 @@ double precision function general_primitive_integral(dim,            &
     return
   endif
   
-  rho = p*q *pq_inv_2  ! le rho qui va bien
+  rho = p*q *pq_inv_2
   dist =  (P_center(1) - Q_center(1))*(P_center(1) - Q_center(1)) +  &
       (P_center(2) - Q_center(2))*(P_center(2) - Q_center(2)) +      &
       (P_center(3) - Q_center(3))*(P_center(3) - Q_center(3))
@@ -580,8 +574,7 @@ double precision function general_primitive_integral(dim,            &
   double precision               :: rint_sum
   accu = accu + rint_sum(n_pt_out,const,d1)
   
-  ! change p+q in dsqrt
+  general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p+q)
-  general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p_plus_q)
 end
 
 
@@ -627,16 +620,10 @@ double precision function ERI(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y,
   ASSERT (gama >= 0.d0)
   p = alpha + beta
   q = delta + gama
-  double precision :: p_plus_q
-  if(long_range)then
-   p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
-  else 
-   p_plus_q = (p+q)
-  endif
   ASSERT (p+q >= 0.d0)
   n_pt =  ishft( nx+ny+nz,1 )
   
-  coeff = pi_5_2 / (p * q * dsqrt(p_plus_q))
+  coeff = pi_5_2 / (p * q * dsqrt(p+q))
   if (n_pt == 0) then
     ERI = coeff
     return
@@ -663,21 +650,15 @@ subroutine integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q
   integer                        :: i, n_iter, n_pt, j
   double precision               :: I_f, pq_inv, p10_1, p10_2, p01_1, p01_2,rho,pq_inv_2
   integer :: ix,iy,iz, jx,jy,jz, sx,sy,sz
-  double precision :: p_plus_q
   
   j = ishft(n_pt,-1)
   ASSERT (n_pt > 1)
-  if(long_range)then
+  pq_inv = 0.5d0/(p+q)
-   p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
-  else 
-   p_plus_q = (p+q)
-  endif
-  pq_inv = 0.5d0/(p_plus_q)
   pq_inv_2 = pq_inv + pq_inv
   p10_1 = 0.5d0/p
   p01_1 = 0.5d0/q
-  p10_2 = 0.5d0 *  q /(p * p_plus_q)
+  p10_2 = 0.5d0 *  q /(p * q + p * p)
-  p01_2 = 0.5d0 *  p /(q * p_plus_q)
+  p01_2 = 0.5d0 *  p /(q * q + q * p)
   double precision               :: B00(n_pt_max_integrals)
   double precision               :: B10(n_pt_max_integrals), B01(n_pt_max_integrals)
   double precision               :: t1(n_pt_max_integrals), t2(n_pt_max_integrals)

src/Integrals_Monoelec/mo_mono_ints.irp.f

@@ -6,10 +6,23 @@ 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'
-  do j = 1, mo_tot_num
+  if (read_mo_one_integrals) then
-    do i = 1, mo_tot_num
+    call read_one_e_integrals('mo_one_integral', mo_mono_elec_integral,      &
-      mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + &
+        size(mo_mono_elec_integral,1), size(mo_mono_elec_integral,2))
-         mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j)
+    print *,  'MO N-e integrals read from disk'
-    enddo
+  else
-  enddo
+   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)
+     enddo
+   enddo
+  endif
+
+! if (write_mo_one_integrals) then
+!   call write_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 written to disk'
+! endif
+
 END_PROVIDER