Clean Cholesky MO

src/mo_two_e_ints/EZFIO.cfg

@@ -17,12 +17,6 @@ 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 3 or 4 virtual indices
-interface: ezfio,provider,ocaml
-default: false
-
 [io_mo_two_e_integrals_erf]
 type: Disk_access
 doc: Read/Write MO integrals with the long range interaction from/to disk [    Write | Read | None ]

src/mo_two_e_ints/cholesky.irp.f

@@ -68,7 +68,11 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_mo_num, mo_num,
      read(iunit) cholesky_mo_transp
      close(iunit)
   else
+     print *, ''
      print *, 'AO->MO Transformation of Cholesky vectors'
+     print *, '-----------------------------------------'
+     print *, ''
+
      call wall_time(wall0)
 
      allocate(X(mo_num,cholesky_mo_num,ao_num), stat=ierr)
@@ -87,7 +91,7 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_mo_num, mo_num,
      if (write_mo_cholesky) then
        print *,  'Writing Cholesky MO vectors to disk...'
        iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_mo_transp', 'W')
-       write(iunit) rank
+       write(iunit) cholesky_mo_num
        write(iunit) cholesky_mo_transp
        close(iunit)
        call ezfio_set_mo_two_e_ints_io_mo_cholesky('Read')

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -39,29 +39,16 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
     return
   endif
 
-  if (.not. do_direct_integrals) then
-    PROVIDE ao_two_e_integrals_in_map
-  endif
-
-  print *,  ''
-  print *,  'AO -> MO integrals transformation'
-  print *,  '---------------------------------'
-  print *,  ''
-
   call wall_time(wall_1)
   call cpu_time(cpu_1)
 
-  if(no_vvvv_integrals)then
+  if (do_mo_cholesky) then
-    call four_idx_novvvv_old
+    call add_integrals_to_map_cholesky
   else
-    if (do_ao_cholesky) then
+    if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then
-      call add_integrals_to_map_cholesky
+      call four_idx_dgemm
     else
-      if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then
+      call add_integrals_to_map(full_ijkl_bitmask_4)
-        call four_idx_dgemm
-      else
-        call add_integrals_to_map(full_ijkl_bitmask_4)
-      endif
     endif
   endif
 
@@ -92,8 +79,15 @@ subroutine four_idx_dgemm
   double precision, allocatable :: a1(:,:,:,:)
   double precision, allocatable :: a2(:,:,:,:)
 
+  PROVIDE ao_two_e_integrals_in_map mo_coef
+
+  print *, ''
+  print *, 'DGEMM-based AO->MO Transformation'
+  print *, '---------------------------------'
+  print *, ''
+
   if (ao_num > 1289) then
-    print *,  irp_here, ': Integer overflow in ao_num**3'
+    print *,  irp_here, ': Integer overflow in ao_num**3. Set do_ao_cholesky=.True.'
   endif
 
   allocate (a1(ao_num,ao_num,ao_num,ao_num))
@@ -213,6 +207,12 @@ subroutine add_integrals_to_map(mask_ijkl)
 
   PROVIDE ao_two_e_integrals_in_map  mo_coef
 
+
+  print *, ''
+  print *, 'Sparse AO->MO Transformation'
+  print *, '----------------------------'
+  print *, ''
+
   !Get list of MOs for i,j,k and l
   !-------------------------------
 
@@ -469,6 +469,7 @@ subroutine add_integrals_to_map_cholesky
   integer(key_kind)  , allocatable :: buffer_i(:)
   real(integral_kind), allocatable :: buffer_value(:)
 
+  PROVIDE cholesky_mo_transp
   call set_multiple_levels_omp(.False.)
 
   !$OMP PARALLEL DEFAULT(SHARED) &

src/mo_two_e_ints/no_vvvv.irp.f

@@ -1,88 +0,0 @@
-
-subroutine four_idx_novvvv_old
-  use map_module
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-  ! Retransform MO integrals for next CAS-SCF step
-  END_DOC
-  integer(bit_kind)              :: mask_ijkl(N_int,4)
-  integer(bit_kind)              :: mask_ijk(N_int,3)
-
-  print*,'Using partial transformation'
-  print*,'It will not transform all integrals with at least 3 indices within the virtuals'
-    integer                        :: i,j,k,l
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I I I I !!!!!!!!!!!!!!!!!!!!
-    ! (core+inact+act) ^ 4
-    ! <ii|ii>
-    print*, ''
-    print*, '<ii|ii>'
-    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) =  core_inact_act_bitmask_4(i,1)
-      mask_ijkl(i,4) =  core_inact_act_bitmask_4(i,1)
-    enddo
-    call add_integrals_to_map(mask_ijkl)
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I I V V !!!!!!!!!!!!!!!!!!!!
-    ! (core+inact+act) ^ 2  (virt) ^2
-    ! <iv|iv>  = J_iv
-    print*, ''
-    print*, '<iv|iv>'
-    do i = 1,N_int
-      mask_ijkl(i,1) =  core_inact_act_bitmask_4(i,1)
-      mask_ijkl(i,2) =  virt_bitmask(i,1)
-      mask_ijkl(i,3) =  core_inact_act_bitmask_4(i,1)
-      mask_ijkl(i,4) =  virt_bitmask(i,1)
-    enddo
-    call add_integrals_to_map(mask_ijkl)
-
-    ! (core+inact+act) ^ 2  (virt) ^2
-    ! <ii|vv> = (iv|iv)
-    print*, ''
-    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)
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!!
-!    if(.not.no_vvv_integrals)then
-      print*, ''
-      print*, '<rv|sv> and <rv|vs>'
-      do i = 1,N_int
-        mask_ijk(i,1) =  virt_bitmask(i,1)
-        mask_ijk(i,2) =  virt_bitmask(i,1)
-        mask_ijk(i,3) =  virt_bitmask(i,1)
-      enddo
-      call add_integrals_to_map_three_indices(mask_ijk)
-!    endif
-
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I I I V !!!!!!!!!!!!!!!!!!!!
-    ! (core+inact+act) ^ 3  (virt) ^1
-    ! <iv|ii>
-    print*, ''
-    print*, '<iv|ii>'
-    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) =  core_inact_act_bitmask_4(i,1)
-      mask_ijkl(i,4) =  virt_bitmask(i,1)
-    enddo
-    call add_integrals_to_map(mask_ijkl)
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I V V V !!!!!!!!!!!!!!!!!!!!
-    ! (core+inact+act) ^ 1  (virt) ^3
-    ! <iv|vv>
-!    if(.not.no_ivvv_integrals)then
-      print*, ''
-      print*, '<iv|vv>'
-      do i = 1,N_int
-        mask_ijkl(i,1) =  core_inact_act_bitmask_4(i,1)
-        mask_ijkl(i,2) =  virt_bitmask(i,1)
-        mask_ijkl(i,3) =  virt_bitmask(i,1)
-        mask_ijkl(i,4) =  virt_bitmask(i,1)
-      enddo
-      call add_integrals_to_map_no_exit_34(mask_ijkl)
-end