Added do_mo_cholesky

2024-06-27 23:52:24 +02:00 · 2024-06-07 16:09:53 +02:00 · 2024-06-07 16:09:53 +02:00 · f58df5e816
commit f58df5e816
parent af8973770e
8 changed files with 97 additions and 52 deletions
--- a/src/ao_two_e_ints/EZFIO.cfg
+++ b/src/ao_two_e_ints/EZFIO.cfg
@ -6,7 +6,7 @@ default: None

 [io_ao_cholesky]
 type: Disk_access
-doc: Read/Write |AO| integrals from/to disk [ Write | Read | None ] 
+doc: Read/Write |AO| Cholesky integrals from/to disk [ Write | Read | None ] 
 interface: ezfio,provider,ocaml
 default: None

--- a/src/ao_two_e_ints/cholesky.irp.f
+++ b/src/ao_two_e_ints/cholesky.irp.f
@ -76,7 +76,7 @@ END_PROVIDER
   deallocate(cholesky_ao)

   if (read_ao_cholesky) then
-     print *,  'Reading Cholesky vectors from disk...'
+     print *,  'Reading Cholesky AO vectors from disk...'
     iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R')
     read(iunit) rank
     allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr)
@ -486,7 +486,7 @@ END_PROVIDER
     cholesky_ao_num = rank

     if (write_ao_cholesky) then
-       print *,  'Writing Cholesky vectors to disk...'
+       print *,  'Writing Cholesky AO vectors to disk...'
       iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'W')
       write(iunit) rank
       write(iunit) cholesky_ao
@ -499,7 +499,7 @@ END_PROVIDER
   print *, 'Rank  : ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
   print *,  ''
   call wall_time(wall1)
-   print*,'Time to provide AO cholesky vectors = ',wall1-wall0
+   print*,'Time to provide AO cholesky vectors = ',(wall1-wall0)/60.d0, ' min'

 END_PROVIDER

--- a/src/mo_two_e_ints/EZFIO.cfg
+++ b/src/mo_two_e_ints/EZFIO.cfg
@ -1,3 +1,9 @@
+[io_mo_cholesky]
+type: Disk_access
+doc: Read/Write |MO| Cholesky integrals from/to disk [ Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
+
 [io_mo_two_e_integrals]
 type: Disk_access
 doc: Read/Write |MO| integrals from/to disk [ Write | Read | None ] 
--- a/src/mo_two_e_ints/cholesky.irp.f
+++ b/src/mo_two_e_ints/cholesky.irp.f
@ -1,55 +1,98 @@
+BEGIN_PROVIDER [ logical, do_mo_cholesky ]
+ implicit none
+ BEGIN_DOC
+ ! If True, use Cholesky vectors for MO integrals
+ END_DOC
+ do_mo_cholesky = do_ao_cholesky
+END_PROVIDER
+
 BEGIN_PROVIDER [ integer, cholesky_mo_num ]
 implicit none
 BEGIN_DOC
 ! Number of Cholesky vectors in MO basis
 END_DOC
- cholesky_mo_num = cholesky_ao_num
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, cholesky_mo, (mo_num, mo_num, cholesky_mo_num) ]
- implicit none
- BEGIN_DOC
- ! Cholesky vectors in MO basis
- END_DOC
-
- integer :: k, i, j
-
- call set_multiple_levels_omp(.False.)
- !$OMP PARALLEL DO PRIVATE(k)
- do k=1,cholesky_mo_num
-  do j=1,mo_num
-    do i=1,mo_num
-      cholesky_mo(i,j,k) = cholesky_mo_transp(k,i,j)
-    enddo
-  enddo
- enddo
- !$OMP END PARALLEL DO
-
+ integer, external              :: getUnitAndOpen
+ integer                        :: iunit
+ if (read_mo_cholesky) then
+   iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_mo_transp', 'R')
+   read(iunit) cholesky_mo_num
+   close(iunit)
+ else
+   cholesky_mo_num = cholesky_ao_num
+ endif
 END_PROVIDER

+!BEGIN_PROVIDER [ double precision, cholesky_mo, (mo_num, mo_num, cholesky_mo_num) ]
+! implicit none
+! BEGIN_DOC
+! ! Cholesky vectors in MO basis
+! END_DOC
+!
+! integer :: k, i, j
+!
+! call set_multiple_levels_omp(.False.)
+! !$OMP PARALLEL DO PRIVATE(k)
+! do k=1,cholesky_mo_num
+!  do j=1,mo_num
+!    do i=1,mo_num
+!      cholesky_mo(i,j,k) = cholesky_mo_transp(k,i,j)
+!    enddo
+!  enddo
+! enddo
+! !$OMP END PARALLEL DO
+!
+!END_PROVIDER
+!
 BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_mo_num, mo_num, mo_num) ]
 implicit none
 BEGIN_DOC
- ! Cholesky vectors in MO basis
+ ! Cholesky vectors in MO basis. Warning: it is transposed wrt cholesky_ao:
+ !
+ ! -  cholesky_ao        is (ao_num^2 x cholesky_ao_num)
+ !
+ ! - cholesky_mo_transp is (cholesky_mo_num x mo_num^2)
 END_DOC

 double precision, allocatable :: X(:,:,:)
 double precision :: wall0, wall1
- integer :: ierr
- print *, 'AO->MO Transformation of Cholesky vectors'
-  call wall_time(wall0)
+ integer, external              :: getUnitAndOpen
+ integer                        :: iunit, ierr, rank

- allocate(X(mo_num,cholesky_mo_num,ao_num), stat=ierr)
- if (ierr /= 0) then
-   print *, irp_here, ': Allocation failed'
- endif
- call dgemm('T','N', ao_num*cholesky_mo_num, mo_num, ao_num, 1.d0, &
-     cholesky_ao, ao_num, mo_coef, ao_num, 0.d0, X, ao_num*cholesky_mo_num)
- call dgemm('T','N', cholesky_mo_num*mo_num, mo_num, ao_num, 1.d0, &
-     X, ao_num, mo_coef, ao_num, 0.d0, cholesky_mo_transp, cholesky_mo_num*mo_num)
- deallocate(X)
-  call wall_time(wall1)
- print*,'Time for AO->MO Cholesky vectors = ',wall1-wall0
+  if (read_mo_cholesky) then
+      print *,  'Reading Cholesky MO vectors from disk...'
+     iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_mo_transp', 'R')
+     read(iunit) rank
+     if (cholesky_mo_num /= rank) then
+        stop 'inconsistent rank'
+     endif
+     read(iunit) cholesky_mo_transp
+     close(iunit)
+  else
+     print *, 'AO->MO Transformation of Cholesky vectors'
+     call wall_time(wall0)
+
+     allocate(X(mo_num,cholesky_mo_num,ao_num), stat=ierr)
+     if (ierr /= 0) then
+       print *, irp_here, ': Allocation failed'
+     endif
+     call dgemm('T','N', ao_num*cholesky_mo_num, mo_num, ao_num, 1.d0, &
+         cholesky_ao, ao_num, mo_coef, ao_num, 0.d0, X, ao_num*cholesky_mo_num)
+     call dgemm('T','N', cholesky_mo_num*mo_num, mo_num, ao_num, 1.d0, &
+         X, ao_num, mo_coef, ao_num, 0.d0, cholesky_mo_transp, cholesky_mo_num*mo_num)
+     deallocate(X)
+     call wall_time(wall1)
+     print*,'Time to provide MO cholesky vectors = ',(wall1-wall0)/60.d0, ' min'
+
+
+     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_transp
+       close(iunit)
+       call ezfio_set_mo_two_e_ints_io_mo_cholesky('Read')
+     endif
+  endif

 END_PROVIDER

--- a/src/mo_two_e_ints/four_idx_novvvv.irp.f
+++ b/src/mo_two_e_ints/four_idx_novvvv.irp.f
@ -1,12 +1,3 @@
-!BEGIN_PROVIDER [ logical, no_vvvv_integrals  ]
-!  implicit none
-!  BEGIN_DOC
-! If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices
-!  END_DOC
-!
-!  no_vvvv_integrals = .False.
-!END_PROVIDER
-
 BEGIN_PROVIDER [ double precision, mo_coef_novirt, (ao_num,n_core_inact_act_orb) ]
 implicit none
 BEGIN_DOC
--- a/src/mo_two_e_ints/integrals_3_index.irp.f
+++ b/src/mo_two_e_ints/integrals_3_index.irp.f
@ -10,7 +10,7 @@
 double precision :: get_two_e_integral
 double precision :: integral

- if (do_ao_cholesky) then
+ if (do_mo_cholesky) then

    double precision, allocatable :: buffer_jj(:,:), buffer(:,:,:)
    allocate(buffer_jj(cholesky_mo_num,mo_num), buffer(mo_num,mo_num,mo_num))
--- a/src/mo_two_e_ints/mo_bi_integrals.irp.f
+++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f
@ -1362,7 +1362,7 @@ END_PROVIDER
  double precision               :: get_two_e_integral


-  if (do_ao_cholesky) then
+  if (do_mo_cholesky) then
    double precision, allocatable :: buffer(:,:)
    allocate (buffer(cholesky_mo_num,mo_num))
    do k=1,cholesky_mo_num
--- a/src/trexio/import_trexio_integrals.irp.f
+++ b/src/trexio/import_trexio_integrals.irp.f
@ -41,6 +41,11 @@ subroutine run(f)
  integer         , allocatable :: Vi(:,:)
  double precision              :: s

+! TODO:
+! - If Cholesky AO in trexio file, read cholesky ao vectors
+! - If Cholesky MO in trexio file, read cholesky mo vectors
+! - If Cholesky MO not in trexio file, force do_cholesky_mo to False
+
  if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
    rc = trexio_read_nucleus_repulsion(f, s)
    if (rc /= TREXIO_SUCCESS) then