Accelerated cache

src/mo_two_e_ints/cholesky.irp.f

@@ -4,7 +4,7 @@ BEGIN_PROVIDER [ logical, do_mo_cholesky ]
  ! If True, use Cholesky vectors for MO integrals
  END_DOC
  do_mo_cholesky = do_ao_cholesky
- do_mo_cholesky = .False.
+! do_mo_cholesky = .False.
 END_PROVIDER
 
 BEGIN_PROVIDER [ integer, cholesky_mo_num ]

src/mo_two_e_ints/four_idx_novvvv.irp.f

@@ -1,180 +0,0 @@
-BEGIN_PROVIDER [ double precision, mo_coef_novirt, (ao_num,n_core_inact_act_orb) ]
- implicit none
- BEGIN_DOC
- ! MO coefficients without virtual MOs
- END_DOC
- integer :: j,jj
-
- do j=1,n_core_inact_act_orb
-   jj = list_core_inact_act(j)
-   mo_coef_novirt(:,j) = mo_coef(:,jj)
- enddo
-
-END_PROVIDER
-
-subroutine ao_to_mo_novirt(A_ao,LDA_ao,A_mo,LDA_mo)
-  implicit none
-  BEGIN_DOC
-  ! Transform A from the |AO| basis to the |MO| basis excluding virtuals
-  !
-  ! $C^\dagger.A_{ao}.C$
-  END_DOC
-  integer, intent(in)            :: LDA_ao,LDA_mo
-  double precision, intent(in)   :: A_ao(LDA_ao,ao_num)
-  double precision, intent(out)  :: A_mo(LDA_mo,n_core_inact_act_orb)
-  double precision, allocatable  :: T(:,:)
-
-  allocate ( T(ao_num,n_core_inact_act_orb) )
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
-
-  call dgemm('N','N', ao_num, n_core_inact_act_orb, ao_num,          &
-      1.d0, A_ao,LDA_ao,                                             &
-      mo_coef_novirt, size(mo_coef_novirt,1),                        &
-      0.d0, T, size(T,1))
-  
-  call dgemm('T','N', n_core_inact_act_orb, n_core_inact_act_orb, ao_num,&
-      1.d0, mo_coef_novirt,size(mo_coef_novirt,1),                   &
-      T, ao_num,                                                     &
-      0.d0, A_mo, size(A_mo,1))
-  
-  deallocate(T)
-end
-
-
-subroutine four_idx_novvvv
-  print*,'********'
-  print*,'********'
-  print*,'********'
-  print*,'WARNING :: Using four_idx_novvvv, and we are not sure that this routine is not bugged ...'
-  print*,'********'
-  print*,'********'
-  print*,'********'
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Retransform MO integrals for next CAS-SCF step
-  END_DOC
-  print*,'Using partial transformation'
-  print*,'It will not transform all integrals with at least 3 indices within the virtuals'
-  integer                        :: i,j,k,l,n_integrals
-  double precision, allocatable  :: f(:,:,:), f2(:,:,:), d(:,:), T(:,:,:,:), T2(:,:,:,:)
-  double precision, external     :: get_ao_two_e_integral
-  integer(key_kind), allocatable :: idx(:)
-  real(integral_kind), allocatable :: values(:)
-
-  integer                        :: p,q,r,s
-  double precision               :: c
-  allocate( T(n_core_inact_act_orb,n_core_inact_act_orb,ao_num,ao_num) , &
-            T2(n_core_inact_act_orb,n_core_inact_act_orb,ao_num,ao_num) )
-  
-  !$OMP PARALLEL DEFAULT(NONE)                                       &
-  !$OMP SHARED(mo_num,ao_num,T,n_core_inact_act_orb, mo_coef_transp, &
-  !$OMP   mo_integrals_threshold,mo_coef,mo_integrals_map,           &
-  !$OMP   list_core_inact_act,T2,ao_integrals_map)            &
-  !$OMP PRIVATE(i,j,k,l,p,q,r,s,idx,values,n_integrals,              &
-  !$OMP   f,f2,d,c)
-  allocate(f(ao_num,ao_num,ao_num), f2(ao_num,ao_num,ao_num), d(mo_num,mo_num), &
-           idx(mo_num*mo_num), values(mo_num*mo_num) )
-  
-  ! <aa|vv>
-  !$OMP DO
-  do s=1,ao_num
-    do r=1,ao_num
-      do q=1,ao_num
-        do p=1,r
-          f (p,q,r) = get_ao_two_e_integral(p,q,r,s,ao_integrals_map)
-          f (r,q,p) = f(p,q,r)
-        enddo
-      enddo
-    enddo
-    do r=1,ao_num
-      do q=1,ao_num
-        do p=1,ao_num
-          f2(p,q,r) = f(p,r,q)
-        enddo
-      enddo
-    enddo
-    ! f (p,q,r) = <pq|rs>
-    ! f2(p,q,r) = <pr|qs>
-
-    do r=1,ao_num
-      call ao_to_mo_novirt(f (1,1,r),size(f ,1),T (1,1,r,s),size(T,1))
-      call ao_to_mo_novirt(f2(1,1,r),size(f2,1),T2(1,1,r,s),size(T,1))
-    enddo
-    ! T (i,j,p,q) = <ij|rs>
-    ! T2(i,j,p,q) = <ir|js>
-
-  enddo
-  !$OMP END DO 
-
-  !$OMP DO
-  do j=1,n_core_inact_act_orb
-    do i=1,n_core_inact_act_orb
-      do s=1,ao_num
-        do r=1,ao_num
-          f (r,s,1) = T (i,j,r,s)
-          f2(r,s,1) = T2(i,j,r,s)
-        enddo
-      enddo
-      call ao_to_mo(f ,size(f ,1),d,size(d,1))
-      n_integrals = 0
-      do l=1,mo_num
-        do k=1,mo_num
-          n_integrals+=1
-          call two_e_integrals_index(list_core_inact_act(i),list_core_inact_act(j),k,l,idx(n_integrals))
-          values(n_integrals) = d(k,l)
-        enddo
-      enddo
-      call map_append(mo_integrals_map, idx, values, n_integrals)
-
-      call ao_to_mo(f2,size(f2,1),d,size(d,1))
-      n_integrals = 0
-      do l=1,mo_num
-        do k=1,mo_num
-          n_integrals+=1
-          call two_e_integrals_index(list_core_inact_act(i),k,list_core_inact_act(j),l,idx(n_integrals))
-          values(n_integrals) = d(k,l)
-        enddo
-      enddo
-      call map_append(mo_integrals_map, idx, values, n_integrals)
-    enddo
-  enddo
-  !$OMP END DO
-  deallocate(f,f2,d,idx,values)
-  
-  !$OMP END PARALLEL
-
-  deallocate(T,T2)
-  
-  
-  call map_sort(mo_integrals_map)
-  call map_unique(mo_integrals_map)
-  call map_shrink(mo_integrals_map,real(mo_integrals_threshold,integral_kind))
-  
-end
-
-subroutine four_idx_novvvv2
-  use bitmasks
-  implicit none
-  integer                        :: i
-  integer(bit_kind)              :: mask_ijkl(N_int,4)
-
-    print*, '<ix|ix>'
-    do i = 1,N_int
-      mask_ijkl(i,1) =  core_inact_act_bitmask_4(i,1)
-      mask_ijkl(i,2) =  full_ijkl_bitmask_4(i,1)
-      mask_ijkl(i,3) =  core_inact_act_bitmask_4(i,1)
-      mask_ijkl(i,4) =  full_ijkl_bitmask_4(i,1)
-    enddo
-    call add_integrals_to_map(mask_ijkl)
-
-    print*, '<ii|vv>'
-    do i = 1,N_int
-      mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1)
-      mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1)
-      mask_ijkl(i,3) = virt_bitmask(i,1)
-      mask_ijkl(i,4) = virt_bitmask(i,1)
-    enddo
-    call add_integrals_to_map(mask_ijkl)
-
-end

src/mo_two_e_ints/map_integrals.irp.f

@@ -48,6 +48,8 @@ end
    mo_integrals_cache_shift = 7  ! 7 = log(128). Max 7
  endif
 
+!mo_integrals_cache_shift = 2  ! 5 = log(32).
+
  mo_integrals_cache_size  = 2**mo_integrals_cache_shift
 
  mo_integrals_cache_min = max(1,elec_alpha_num - (mo_integrals_cache_size/2 - 1) )
@@ -112,6 +114,24 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0:mo_integrals_cache_siz
 END_PROVIDER
 
 
+double precision function get_two_e_integral_cache(i,j,k,l)
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! Returns one integral <ij|kl> in the MO basis taken from the cache
+  END_DOC
+  integer, intent(in)            :: i,j,k,l
+  integer                        :: ii
+
+  ii = l-mo_integrals_cache_min
+  ii = ior( shiftl(ii,mo_integrals_cache_shift), k-mo_integrals_cache_min)
+  ii = ior( shiftl(ii,mo_integrals_cache_shift), j-mo_integrals_cache_min)
+  ii = ior( shiftl(ii,mo_integrals_cache_shift), i-mo_integrals_cache_min)
+  get_two_e_integral_cache = mo_integrals_cache(ii)
+
+end
+
+
 double precision function get_two_e_integral(i,j,k,l,map)
   use map_module
   implicit none
@@ -123,7 +143,6 @@ double precision function get_two_e_integral(i,j,k,l,map)
   integer                        :: ii
   type(map_type), intent(inout)  :: map
   real(integral_kind)            :: tmp
-  integer                        :: kk
 
   PROVIDE mo_two_e_integrals_in_map mo_integrals_cache do_mo_cholesky
 
@@ -145,13 +164,9 @@ double precision function get_two_e_integral(i,j,k,l,map)
   ii = ior(ii, i-mo_integrals_cache_min)
 
   if (iand(ii, -mo_integrals_cache_size) == 0) then
-    ! Integrals is in the cache
 
-    ii = l-mo_integrals_cache_min
+    double precision, external :: get_two_e_integral_cache
-    ii = ior( shiftl(ii,mo_integrals_cache_shift), k-mo_integrals_cache_min)
+    get_two_e_integral = get_two_e_integral_cache(i,j,k,l)
-    ii = ior( shiftl(ii,mo_integrals_cache_shift), j-mo_integrals_cache_min)
-    ii = ior( shiftl(ii,mo_integrals_cache_shift), i-mo_integrals_cache_min)
-    get_two_e_integral = mo_integrals_cache(ii)
 
   else
 
@@ -199,7 +214,6 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
   double precision, intent(out)  :: out_val(sze)
   type(map_type), intent(inout)  :: map
   integer                        :: i
-  double precision, external :: get_two_e_integral
 
   integer                        :: ii
   real(integral_kind)            :: tmp
@@ -256,13 +270,7 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
 
     endif
 
-
+    call get_mo_two_e_integrals_cache(j,k,l,sze,out_val)
-    ii = l-mo_integrals_cache_min
-    ii = ior( shiftl(ii, mo_integrals_cache_shift), k-mo_integrals_cache_min)
-    ii = ior( shiftl(ii, mo_integrals_cache_shift), j-mo_integrals_cache_min)
-    ii = shiftl(ii, mo_integrals_cache_shift)
-    out_val(mo_integrals_cache_min:mo_integrals_cache_max) = &
-       mo_integrals_cache(ii:ii+mo_integrals_cache_max-mo_integrals_cache_min)
 
     if (mo_integrals_cache_max < mo_num) then
 
@@ -333,6 +341,26 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
 
 end
 
+subroutine get_mo_two_e_integrals_cache(j,k,l,sze,out_val)
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! Returns multiple integrals <ij|kl> in the MO basis, all
+  ! i for j,k,l fixed, all integrals from the cache
+  END_DOC
+  integer, intent(in)            :: j,k,l, sze
+  double precision, intent(out)  :: out_val(sze)
+  integer                        :: ii
+
+  ii = l-mo_integrals_cache_min
+  ii = ior( shiftl(ii, mo_integrals_cache_shift), k-mo_integrals_cache_min)
+  ii = ior( shiftl(ii, mo_integrals_cache_shift), j-mo_integrals_cache_min)
+  ii = shiftl(ii, mo_integrals_cache_shift)
+  out_val(mo_integrals_cache_min:mo_integrals_cache_max) = &
+       mo_integrals_cache(ii:ii+mo_integrals_cache_max-mo_integrals_cache_min)
+
+end
+
 subroutine get_mo_two_e_integrals_ij(k,l,sze,out_array,map)
   use map_module
   implicit none
@@ -347,16 +375,32 @@ subroutine get_mo_two_e_integrals_ij(k,l,sze,out_array,map)
   integer                        :: j
   real(integral_kind), allocatable :: tmp_val(:)
 
-  if (do_mo_cholesky) then
+  if ( (mo_integrals_cache_min>1).or.(mo_integrals_cache_max<mo_num) ) then
-    call dgemm('T', 'N', mo_num, mo_num, cholesky_mo_num, 1.d0, &
+
-       cholesky_mo_transp(1,1,k), cholesky_mo_num, &
+    if (do_mo_cholesky) then
-       cholesky_mo_transp(1,1,l), cholesky_mo_num, 0.d0, &
+
-       out_array, sze)
+      call dgemm('T', 'N', mo_num, mo_num, cholesky_mo_num, 1.d0, &
+         cholesky_mo_transp(1,1,k), cholesky_mo_num, &
+         cholesky_mo_transp(1,1,l), cholesky_mo_num, 0.d0, &
+         out_array, sze)
+
+    else
+
+      do j=1,sze
+        call get_mo_two_e_integrals(j,k,l,sze,out_array(1,j),map)
+      enddo
+
+    endif
+
   else
+
+    double precision, external :: get_two_e_integral_cache
     do j=1,sze
-      call get_mo_two_e_integrals(j,k,l,sze,out_array(1,j),map)
+      call get_mo_two_e_integrals_cache(j,k,l,sze,out_array(1,j))
     enddo
+
   endif
+
 end
 
 subroutine get_mo_two_e_integrals_i1j1(k,l,sze,out_array,map)
@@ -373,17 +417,28 @@ subroutine get_mo_two_e_integrals_i1j1(k,l,sze,out_array,map)
   integer                        :: j
   PROVIDE mo_two_e_integrals_in_map
 
-  if (do_mo_cholesky) then
+  if ( (mo_integrals_cache_min>1).or.(mo_integrals_cache_max<mo_num) ) then
 
-    call dgemv('T', cholesky_mo_num, mo_num*mo_num, 1.d0, &
+    if (do_mo_cholesky) then
-       cholesky_mo_transp(1,1,1), cholesky_mo_num, &
+
-       cholesky_mo_transp(1,k,l), 1, 0.d0, &
+      call dgemv('T', cholesky_mo_num, mo_num*mo_num, 1.d0, &
-       out_array, 1)
+         cholesky_mo_transp(1,1,1), cholesky_mo_num, &
+         cholesky_mo_transp(1,k,l), 1, 0.d0, &
+         out_array, 1)
+
+    else
+
+      do j=1,sze
+        call get_mo_two_e_integrals(k,j,l,sze,out_array(1,j),map)
+      enddo
+
+    endif
 
   else
 
+    double precision, external :: get_two_e_integral_cache
     do j=1,sze
-      call get_mo_two_e_integrals(k,j,l,sze,out_array(1,j),map)
+      call get_mo_two_e_integrals_cache(k,j,l,sze,out_array(1,j))
     enddo
 
   endif
@@ -406,21 +461,33 @@ subroutine get_mo_two_e_integrals_coulomb_ii(k,l,sze,out_val,map)
   double precision, external     :: get_two_e_integral
   PROVIDE mo_two_e_integrals_in_map
 
-  if (do_mo_cholesky) then
+  if ( (mo_integrals_cache_min>1).or.(mo_integrals_cache_max<mo_num) ) then
 
-    call dgemv('T', cholesky_mo_num, mo_num, 1.d0, &
+    if (do_mo_cholesky) then
-       cholesky_mo_transp(1,1,1), cholesky_mo_num*(mo_num+1), &
+
-       cholesky_mo_transp(1,k,l), 1, 0.d0, &
+      call dgemv('T', cholesky_mo_num, mo_num, 1.d0, &
-       out_val, 1)
+         cholesky_mo_transp(1,1,1), cholesky_mo_num*(mo_num+1), &
+         cholesky_mo_transp(1,k,l), 1, 0.d0, &
+         out_val, 1)
+
+    else
+
+      do i=1,sze
+        out_val(i) = get_two_e_integral(i,k,i,l,map)
+      enddo
+
+    endif
 
   else
 
+    double precision, external :: get_two_e_integral_cache
     do i=1,sze
-      out_val(i) = get_two_e_integral(k,i,l,i,map)
+      out_val(i) = get_two_e_integral_cache(i,k,i,l)
     enddo
 
   endif
 
+
 end
 
 subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
@@ -438,9 +505,33 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
   double precision, external     :: get_two_e_integral
   PROVIDE mo_two_e_integrals_in_map
 
-  do i=1,sze
+  if ( (mo_integrals_cache_min>1).or.(mo_integrals_cache_max<mo_num) ) then
-    out_val(i) = get_two_e_integral(k,i,i,l,map)
+
-  enddo
+    if (do_mo_cholesky) then
+
+      double precision, external :: ddot
+      do i=1,sze
+        out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, &
+                                           cholesky_mo_transp(1,i,l), 1)
+      enddo
+
+    else
+
+      do i=1,sze
+        out_val(i) = get_two_e_integral(i,i,k,l,map)
+      enddo
+
+    endif
+
+  else
+
+    double precision, external :: get_two_e_integral_cache
+    do i=1,sze
+      out_val(i) = get_two_e_integral_cache(i,i,k,l)
+    enddo
+
+  endif
+
 
 end
 

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -15,10 +15,6 @@
 ! - All integrals are stored in the map or cache map
 ! - 1,2,3-index arrays are built from the map
 !
-! TODO:
-! - get_mo_integrals using cholesky
-! - get_mo_integralss using cholesky
-! - get_mo_integralss in PT2
 
 subroutine mo_two_e_integrals_index(i,j,k,l,i1)
   use map_module