Accelerated Cholesky

src/ao_two_e_ints/EZFIO.cfg

@@ -11,6 +11,12 @@ interface: ezfio,provider,ocaml
 default: 1.e-15
 ezfio_name: threshold_ao
 
+[ao_cholesky_threshold]
+type: Threshold
+doc: If | (ii|jj) | < `ao_cholesky_threshold` then (ii|jj) is zero
+interface: ezfio,provider,ocaml
+default: 1.e-12
+
 [do_direct_integrals]
 type: logical
 doc: Compute integrals on the fly (very slow, only for debugging)

src/ao_two_e_ints/cholesky.irp.f

@@ -4,29 +4,7 @@ BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
  ! Number of Cholesky vectors in AO basis
  END_DOC
 
- integer :: i,j,k,l
+ cholesky_ao_num_guess = ao_num*ao_num / 2
- double precision :: xnorm0, x, integral
- double precision, external :: ao_two_e_integral
-
- cholesky_ao_num_guess = 0
- xnorm0 = 0.d0
- x = 0.d0
- do j=1,ao_num
-   do i=1,ao_num
-     integral = ao_two_e_integral(i,i,j,j)
-     if (integral > ao_integrals_threshold) then
-       cholesky_ao_num_guess += 1
-     else
-       x += integral
-     endif
-   enddo
- enddo
- print *, 'Cholesky decomposition of AO integrals'
- print *, '--------------------------------------'
- print *, ''
- print *, 'Estimated Error: ', x
- print *, 'Guess size: ', cholesky_ao_num_guess, '(', 100.d0*dble(cholesky_ao_num_guess)/dble(ao_num*ao_num), ' %)'
-
 END_PROVIDER
 
  BEGIN_PROVIDER [ integer, cholesky_ao_num ]
@@ -39,7 +17,7 @@ END_PROVIDER
  END_DOC
 
  type(c_ptr) :: ptr
- integer :: fd, i,j,k,l, rank
+ integer :: fd, i,j,k,l,m,rank
  double precision, pointer :: ao_integrals(:,:,:,:)
  double precision, external :: ao_two_e_integral
 
@@ -49,24 +27,49 @@ END_PROVIDER
    8, fd, .False., ptr)
  call c_f_pointer(ptr, ao_integrals, (/ao_num, ao_num, ao_num, ao_num/))
 
- double precision :: integral
+ print*, 'Providing the AO integrals (Cholesky)'
+ call wall_time(wall_1)
+ call cpu_time(cpu_1)
+
+ ao_integrals = 0.d0
+
+ double precision :: integral, cpu_1, cpu_2, wall_1, wall_2
  logical, external :: ao_two_e_integral_zero
- !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k,l, integral) SCHEDULE(dynamic)
+
- do l=1,ao_num
+ !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l, integral, wall_2)
-  do j=1,l
+ do m=0,9
-   do k=1,ao_num
+   do l=1+m,ao_num,10
-    do i=1,k
+     !$OMP DO SCHEDULE(dynamic)
-     if (ao_two_e_integral_zero(i,j,k,l)) cycle
+     do j=1,l
-     integral = ao_two_e_integral(i,k,j,l)
+       do k=1,ao_num
-     ao_integrals(i,k,j,l) = integral
+         do i=1,min(k,j)
-     ao_integrals(k,i,j,l) = integral
+           if (ao_two_e_integral_zero(i,j,k,l)) cycle
-     ao_integrals(i,k,l,j) = integral
+           integral = ao_two_e_integral(i,k,j,l)
-     ao_integrals(k,i,l,j) = integral
+           ao_integrals(i,k,j,l) = integral
-    enddo
+           ao_integrals(k,i,j,l) = integral
+           ao_integrals(i,k,l,j) = integral
+           ao_integrals(k,i,l,j) = integral
+           ao_integrals(j,l,i,k) = integral
+           ao_integrals(j,l,k,i) = integral
+           ao_integrals(l,j,i,k) = integral
+           ao_integrals(l,j,k,i) = integral
+         enddo
+       enddo
+     enddo
+     !$OMP END DO NOWAIT
    enddo
-  enddo
+   !$OMP MASTER
+   call wall_time(wall_2)
+   print '(F10.2,'' %  in'', 4X, I10, '' s.'')', (m+1) * 10, wall_2-wall_1
+   !$OMP END MASTER
  enddo
- !$OMP END PARALLEL DO
+ !$OMP END PARALLEL
+
+ call wall_time(wall_2)
+ call cpu_time(cpu_2)
+ print*, 'AO integrals provided:'
+ 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)), ' )'
 
  ! Call Lapack
  cholesky_ao_num = cholesky_ao_num_guess

src/ccsd/ccsd_space_orb_sub.irp.f

@@ -16,7 +16,7 @@ subroutine run_ccsd_space_orb
   double precision, allocatable :: all_err(:,:), all_t(:,:)
   integer, allocatable          :: list_occ(:), list_vir(:)
   integer(bit_kind)             :: det(N_int,2)
-  integer                       :: nO, nV, nOa, nOb, nVa, nVb, n_spin(4)
+  integer                       :: nO, nV, nOa, nVa
 
   PROVIDE mo_two_e_integrals_in_map
 
@@ -24,12 +24,8 @@ subroutine run_ccsd_space_orb
   print*,'Reference determinant:'
   call print_det(det,N_int)
 
-  ! Extract number of occ/vir alpha/beta spin orbitals
+  nOa = cc_nOa
-  !call extract_n_spin(det,n_spin)
+  nVa = cc_nVa
-  nOa = cc_nOa !n_spin(1)
-  nOb = cc_nOb !n_spin(2)
-  nVa = cc_nVa !n_spin(3)
-  nVb = cc_nVb !n_spin(4)
 
   ! Check that the reference is a closed shell determinant
   if (cc_ref_is_open_shell) then