From f58df5e81669226e728f744308593db5a5e4cad0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 7 Jun 2024 16:09:53 +0200 Subject: [PATCH] Added do_mo_cholesky --- src/ao_two_e_ints/EZFIO.cfg | 2 +- src/ao_two_e_ints/cholesky.irp.f | 6 +- src/mo_two_e_ints/EZFIO.cfg | 6 ++ src/mo_two_e_ints/cholesky.irp.f | 117 +++++++++++++++------- src/mo_two_e_ints/four_idx_novvvv.irp.f | 9 -- src/mo_two_e_ints/integrals_3_index.irp.f | 2 +- src/mo_two_e_ints/mo_bi_integrals.irp.f | 2 +- src/trexio/import_trexio_integrals.irp.f | 5 + 8 files changed, 97 insertions(+), 52 deletions(-) diff --git a/src/ao_two_e_ints/EZFIO.cfg b/src/ao_two_e_ints/EZFIO.cfg index c2e083a3..a985149e 100644 --- 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 diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index 6778d5c7..a1cd8e5b 100644 --- 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 diff --git a/src/mo_two_e_ints/EZFIO.cfg b/src/mo_two_e_ints/EZFIO.cfg index 088a2416..49a2952c 100644 --- 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 ] diff --git a/src/mo_two_e_ints/cholesky.irp.f b/src/mo_two_e_ints/cholesky.irp.f index 0d0989d7..d3affd68 100644 --- 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 diff --git a/src/mo_two_e_ints/four_idx_novvvv.irp.f b/src/mo_two_e_ints/four_idx_novvvv.irp.f index 2be09689..80af35dc 100644 --- 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 diff --git a/src/mo_two_e_ints/integrals_3_index.irp.f b/src/mo_two_e_ints/integrals_3_index.irp.f index eb05da84..c0dab506 100644 --- 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)) diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index 0e77b6a2..cb3f4bc6 100644 --- 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 diff --git a/src/trexio/import_trexio_integrals.irp.f b/src/trexio/import_trexio_integrals.irp.f index 8c6b79d7..5a6b3c03 100644 --- 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