Add possibility of single precision cholesky vectors

src/ao_two_e_ints/cholesky.irp.f

@@ -469,7 +469,7 @@ END_PROVIDER
      !$OMP PARALLEL DO PRIVATE(k,j)
      do k=1,rank
        do j=1,ao_num
-           cholesky_ao(1:ao_num,j,rank-k+1) = L((j-1_8)*ao_num+1_8:1_8*j*ao_num,rank-k+1)
+           cholesky_ao(1:ao_num,j,k) = L((j-1_8)*ao_num+1_8:1_8*j*ao_num,rank-k+1)
        enddo
      enddo
      !$OMP END PARALLEL DO

src/mo_two_e_ints/EZFIO.cfg

@@ -29,4 +29,10 @@ doc: Read/Write MO integrals with the long range interaction from/to disk [    W
 interface: ezfio,provider,ocaml
 default: None
 
+[mo_cholesky_double]
+type: logical
+doc: Use double precision to build integrals from Cholesky vectors
+interface: ezfio,provider,ocaml
+default: True
+
 

src/mo_two_e_ints/cholesky.irp.f

@@ -175,9 +175,9 @@ BEGIN_PROVIDER [ real, cholesky_mo_transp_sp, (cholesky_mo_num, mo_num, mo_num)
 
  integer :: i,j,k
  !$OMP PARALLEL DO PRIVATE(k)
- do k=1,cholesky_mo_num
+ do j=1,mo_num
-  do j=1,mo_num
+  do i=1,mo_num
-    do i=1,mo_num
+   do k=1,cholesky_mo_num
       cholesky_mo_transp_sp(k,i,j) = cholesky_mo_transp(k,i,j)
     enddo
   enddo

src/mo_two_e_ints/map_integrals.irp.f

@@ -9,15 +9,6 @@ BEGIN_PROVIDER [ logical, all_mo_integrals ]
   PROVIDE mo_two_e_integrals_in_map mo_integrals_cache mo_two_e_integrals_jj_exchange mo_two_e_integrals_jj_anti mo_two_e_integrals_jj big_array_exchange_integrals big_array_coulomb_integrals mo_one_e_integrals
 END_PROVIDER
 
-BEGIN_PROVIDER [ logical, mo_cholesky_double ]
-  implicit none
-  BEGIN_DOC
-! If true, use double precision to compute integrals from cholesky vectors
-  END_DOC
-  mo_cholesky_double = .True.
-END_PROVIDER
-
-
 !! MO Map
 !! ======
 
@@ -193,13 +184,14 @@ double precision function get_two_e_integral(i,j,k,l,map)
       if (mo_cholesky_double) then
         get_two_e_integral = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, cholesky_mo_transp(1,j,l), 1)
       else
-        get_two_e_integral = 0.d0
+        get_two_e_integral = sdot(cholesky_mo_num, cholesky_mo_transp_sp(1,i,k), 1, cholesky_mo_transp_sp(1,j,l), 1)
-        do isplit=1,4
+!        get_two_e_integral = 0.d0
-          get_two_e_integral = get_two_e_integral + &
+!        do isplit=1,4
-                               sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+!          get_two_e_integral = get_two_e_integral + &
-                                    cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, &
+!                               sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
-                                    cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1)
+!                                    cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, &
-        enddo
+!                                    cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1)
+!        enddo
       endif
 
     else
@@ -244,7 +236,7 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
   ii = ior(ii, j-mo_integrals_cache_min)
 
   if (do_mo_cholesky.and. .not.mo_cholesky_double) then
-    allocate(out_val_sp(sze))
+    allocate(out_val_sp(mo_num))
   endif
 
   if (iand(ii, -mo_integrals_cache_size) == 0) then
@@ -261,22 +253,27 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
             cholesky_mo_transp(1,j,l), 1, 0.d0, &
             out_val, 1)
         else
-          integer :: isplit
+          call sgemv('T', cholesky_mo_num, mo_integrals_cache_min-1, 1., &
-          call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+            cholesky_mo_transp_sp(1,1,k), cholesky_mo_num, &
-              mo_integrals_cache_min-1, 1., &
+            cholesky_mo_transp_sp(1,j,l), 1, 0., &
-              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
+            out_val_sp, 1)
-              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
-              out_val_sp, 1)
-            out_val(1:mo_integrals_cache_min-1) += out_val_sp(1:mo_integrals_cache_min-1)
           out_val(1:mo_integrals_cache_min-1) = out_val_sp(1:mo_integrals_cache_min-1)
-          do isplit=2,4
+!          integer :: isplit
-            call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+!          isplit=1
-              mo_integrals_cache_min-1, 1., &
+!          call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
-              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
+!              mo_integrals_cache_min-1, 1., &
-              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
+!              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
-              out_val_sp, 1)
+!              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
-            out_val(1:mo_integrals_cache_min-1) += out_val_sp(1:mo_integrals_cache_min-1)
+!              out_val_sp, 1)
-          enddo
+!          out_val(1:mo_integrals_cache_min-1) = out_val_sp(1:mo_integrals_cache_min-1)
+!          do isplit=2,4
+!              call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+!                mo_integrals_cache_min-1, 1., &
+!                cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
+!                cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
+!                out_val_sp, 1)
+!              out_val(1:mo_integrals_cache_min-1) += out_val_sp(1:mo_integrals_cache_min-1)
+!          enddo
         endif
 
       else
@@ -318,15 +315,26 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
              cholesky_mo_transp(1,j,l), 1, 0.d0, &
              out_val(mo_integrals_cache_max+1), 1)
         else
-          out_val = 0.d0
+          call sgemv('T', cholesky_mo_num, mo_num-mo_integrals_cache_max, 1., &
-          do isplit=1,4
+             cholesky_mo_transp_sp(1,mo_integrals_cache_max+1,k), cholesky_mo_num, &
-            call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+             cholesky_mo_transp_sp(1,j,l), 1, 0., &
-              mo_num-mo_integrals_cache_max, 1., &
+             out_val_sp(mo_integrals_cache_max+1), 1)
-              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),mo_integrals_cache_max+1,k), cholesky_mo_num, &
+          out_val(mo_integrals_cache_max+1:sze) = out_val_sp(mo_integrals_cache_max+1:sze)
-              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
+!          isplit=1
-              out_val_sp(mo_integrals_cache_max+1), 1)
+!          call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
-            out_val(mo_integrals_cache_max+1:sze) += out_val_sp(mo_integrals_cache_max+1:sze)
+!            mo_num-mo_integrals_cache_max, 1., &
-          enddo
+!            cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),mo_integrals_cache_max+1,k), cholesky_mo_num, &
+!            cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
+!            out_val_sp(mo_integrals_cache_max+1), 1)
+!          out_val(mo_integrals_cache_max+1:sze) = out_val_sp(mo_integrals_cache_max+1:sze)
+!          do isplit=2,4
+!            call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+!              mo_num-mo_integrals_cache_max, 1., &
+!              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),mo_integrals_cache_max+1,k), cholesky_mo_num, &
+!              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
+!              out_val_sp(mo_integrals_cache_max+1), 1)
+!            out_val(mo_integrals_cache_max+1:sze) += out_val_sp(mo_integrals_cache_max+1:sze)
+!          enddo
         endif
 
       else
@@ -366,15 +374,26 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
                cholesky_mo_transp(1,j,l), 1, 0.d0, &
                out_val, 1)
       else
-          out_val = 0.d0
+          call sgemv('T', cholesky_mo_num, sze, 1., &
-          do isplit=1,4
+               cholesky_mo_transp_sp(1,1,k), cholesky_mo_num, &
-            call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+               cholesky_mo_transp_sp(1,j,l), 1, 0., &
-              sze, 1., &
+               out_val_sp, 1)
-              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
+          out_val(1:sze) = out_val_sp(1:sze)
-              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
+!          isplit=1
-              out_val_sp, 1)
+!          call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
-            out_val(1:sze) += out_val_sp(1:sze)
+!            sze, 1., &
-          enddo
+!            cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
+!            cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
+!            out_val_sp, 1)
+!          out_val(1:sze) = out_val_sp(1:sze)
+!          do isplit=2,4
+!            call sgemv('T', cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+!              sze, 1., &
+!              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
+!              cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1, 0., &
+!              out_val_sp, 1)
+!            out_val(1:sze) += out_val_sp(1:sze)
+!          enddo
       endif
 
     else
@@ -459,23 +478,30 @@ subroutine get_mo_two_e_integrals_ij(k,l,sze,out_array,map)
              cholesky_mo_transp(1,1,l), cholesky_mo_num, 0.d0, &
              out_array, sze)
       else
-          integer :: isplit
+          real, allocatable :: out_array_sp(:,:)
-          double precision, allocatable :: out_array_sp(:,:)
           allocate(out_array_sp(sze,sze))
-          call sgemm('T', 'N', mo_num, mo_num, &
+          call sgemm('T', 'N', mo_num, mo_num, cholesky_mo_num, 1.0, &
-             cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), 1.d0, &
+             cholesky_mo_transp_sp(1,1,k), cholesky_mo_num, &
-             cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
+             cholesky_mo_transp_sp(1,1,l), cholesky_mo_num, 0.0, &
-             cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,l), cholesky_mo_num, 0.d0, &
              out_array_sp, sze)
           out_array(1:sze,1:sze) = out_array_sp(1:sze,1:sze)
-          do isplit=2,4
+!
-            call sgemm('T', 'N', mo_num, mo_num, &
+!          isplit=1
-             cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), 1.d0, &
+!          call sgemm('T', 'N', mo_num, mo_num, &
-             cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
+!             cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), 1., &
-             cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,l), cholesky_mo_num, 0.d0, &
+!             cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
-             out_array_sp, sze)
+!             cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,l), cholesky_mo_num, 0., &
-            out_array(1:sze,1:sze) = out_array(1:sze,1:sze) + out_array_sp(1:sze,1:sze)
+!             out_array_sp, sze)
-         enddo
+!          out_array(1:sze,1:sze) = out_array_sp(1:sze,1:sze)
+!          integer :: isplit
+!          do isplit=2,4
+!            call sgemm('T', 'N', mo_num, mo_num, &
+!             cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), 1., &
+!             cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,k), cholesky_mo_num, &
+!             cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),1,l), cholesky_mo_num, 0., &
+!             out_array_sp, sze)
+!            out_array(1:sze,1:sze) = out_array(1:sze,1:sze) + out_array_sp(1:sze,1:sze)
+!         enddo
          deallocate(out_array_sp)
       endif
 
@@ -631,12 +657,13 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
 
           integer :: isplit
           do i=1,mo_integrals_cache_min-1
-             out_val(i) = 0.d0
+             out_val(i) = sdot(cholesky_mo_num, cholesky_mo_transp_sp(1,i,k), 1, cholesky_mo_transp_sp(1,i,l), 1)
-             do isplit=1,4
+!             out_val(i) = 0.d0
-               out_val(i) += sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+!             do isplit=1,4
-                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, &
+!               out_val(i) += sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
-                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1)
+!                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, &
-             enddo
+!                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1)
+!             enddo
           enddo
 
           do i=mo_integrals_cache_min,mo_integrals_cache_max
@@ -644,12 +671,13 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
           enddo
 
           do i=mo_integrals_cache_max, sze
-             out_val(i) = 0.d0
+             out_val(i) = sdot(cholesky_mo_num, cholesky_mo_transp_sp(1,i,k), 1, cholesky_mo_transp_sp(1,i,l), 1)
-             do isplit=1,4
+!             out_val(i) = 0.d0
-               out_val(i) += sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+!             do isplit=1,4
-                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, &
+!               out_val(i) += sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
-                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1)
+!                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, &
-             enddo
+!                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1)
+!             enddo
           enddo
 
         endif
@@ -663,12 +691,13 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
           enddo
         else
           do i=1,sze
-             out_val(i) = 0.d0
+             out_val(i) = sdot(cholesky_mo_num, cholesky_mo_transp_sp(1,i,k), 1, cholesky_mo_transp_sp(1,i,l), 1)
-             do isplit=1,4
+!             out_val(i) = 0.d0
-               out_val(i) += sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
+!             do isplit=1,4
-                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, &
+!               out_val(i) += sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), &
-                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1)
+!                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, &
-             enddo
+!                                          cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1)
+!             enddo
           enddo
         endif