From eb6e1d4339fe6342da40c56184119c0545b97507 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 13 Feb 2025 10:23:52 +0100 Subject: [PATCH 1/4] Introduced qp_bug.irp.f --- src/determinants/spindeterminants.irp.f | 22 ++++++++++++++++------ src/utils/bug.irp.f | 23 +++++++++++++++++++++++ 2 files changed, 39 insertions(+), 6 deletions(-) create mode 100644 src/utils/bug.irp.f diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index 01979c02..2f497bd7 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -197,7 +197,9 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint) enddo i += 1 - ASSERT (i <= N_det_alpha_unique) + if (i> N_det_alpha_unique) then + call qp_bug(irp_here, i, 'i> N_det_alpha_unique') + endif !DIR$ FORCEINLINE do while (spin_det_search_key(psi_det_alpha_unique(1,i),Nint) == det_ref) @@ -219,12 +221,15 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint) endif i += 1 if (i > N_det_alpha_unique) then - ASSERT (get_index_in_psi_det_alpha_unique > 0) - return + exit endif enddo + if (get_index_in_psi_det_alpha_unique <= 0) then + call qp_bug(irp_here, get_index_in_psi_det_alpha_unique, 'get_index_in_psi_det_alpha_unique <= 0') + endif + end integer function get_index_in_psi_det_beta_unique(key,Nint) @@ -277,7 +282,9 @@ integer function get_index_in_psi_det_beta_unique(key,Nint) enddo i += 1 - ASSERT (i <= N_det_beta_unique) + if (i > N_det_beta_unique) then + call qp_bug(irp_here, i, 'i> N_det_beta_unique') + endif !DIR$ FORCEINLINE do while (spin_det_search_key(psi_det_beta_unique(1,i),Nint) == det_ref) @@ -299,12 +306,15 @@ integer function get_index_in_psi_det_beta_unique(key,Nint) endif i += 1 if (i > N_det_beta_unique) then - ASSERT (get_index_in_psi_det_beta_unique > 0) - return + exit endif enddo + if (get_index_in_psi_det_beta_unique <= 0) then + call qp_bug(irp_here, i, 'get_index_in_psi_det_beta_unique <= 0') + endif + end diff --git a/src/utils/bug.irp.f b/src/utils/bug.irp.f new file mode 100644 index 00000000..0e2ad551 --- /dev/null +++ b/src/utils/bug.irp.f @@ -0,0 +1,23 @@ +subroutine qp_bug(from, code, message) + implicit none + BEGIN_DOC +! This routine prints a bug report + END_DOC + character*(*) :: from + integer :: code + character*(*) :: message + + print *, '' + print *, '=======================' + print *, 'Bug in Quantum Package!' + print *, '=======================' + print *, '' + print *, ' from: ', trim(from) + print *, ' code: ', code + print *, ' info: ', trim(message) + print *, '' + print *, 'Please report this bug at https://github.com/QuantumPackage/qp2/issues' + print *, 'with your output file attached.' + print *, '' + stop -1 +end subroutine qp_bug From 70ad9f31b390b5931392a52ae375fab7c1e9f241 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 13 Feb 2025 11:58:46 +0100 Subject: [PATCH 2/4] Accelerated PT2-cholesky with dgemm --- external/ezfio | 2 +- src/cipsi/selection.irp.f | 28 +++++++++++++------ src/mo_two_e_ints/map_integrals.irp.f | 40 ++++++++++++++++++++++----- 3 files changed, 54 insertions(+), 16 deletions(-) diff --git a/external/ezfio b/external/ezfio index d02132ea..dba01c4f 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit d02132ea79217c16fd24242e8f8b8a6c3ff68091 +Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3 diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index af7420c8..8c22ec85 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -1478,17 +1478,21 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer, parameter :: bant=1 - double precision, allocatable :: hij_cache1(:), hij_cache2(:) - allocate (hij_cache1(mo_num),hij_cache2(mo_num)) +! double precision, allocatable :: hij_cache1(:), hij_cache2(:) + double precision, allocatable :: hij_cache1(:,:), hij_cache2(:,:) +! allocate (hij_cache1(mo_num),hij_cache2(mo_num)) PROVIDE mo_integrals_threshold if(sp == 3) then ! AB + + allocate(hij_cache1(mo_num,mo_num)) + h1 = p(1,1) h2 = p(1,2) + call get_mo_two_e_integrals_ij(h2,h1,mo_num,hij_cache1,mo_integrals_map) do p1=1, mo_num if(bannedOrb(p1, 1)) cycle - call get_mo_two_e_integrals(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map) do p2=1, mo_num if(bannedOrb(p2,2)) cycle if(banned(p1, p2, bant)) cycle ! rentable? @@ -1497,7 +1501,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call i_h_j(gen, det, N_int, hij) else phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) - hij = hij_cache1(p2) * phase + hij = hij_cache1(p2,p1) * phase end if if (dabs(hij) < mo_integrals_threshold) cycle !DIR$ LOOP COUNT AVG(4) @@ -1507,13 +1511,18 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) end do end do + deallocate(hij_cache1) + else ! AA BB + + allocate(hij_cache1(mo_num,mo_num),hij_cache2(mo_num,mo_num)) + p1 = p(1,sp) p2 = p(2,sp) + call get_mo_two_e_integrals_ij(p2,p1,mo_num,hij_cache1,mo_integrals_map) + call get_mo_two_e_integrals_ij(p1,p2,mo_num,hij_cache2,mo_integrals_map) do puti=1, mo_num if (bannedOrb(puti, sp)) cycle - call get_mo_two_e_integrals(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map) - call get_mo_two_e_integrals(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map) do putj=puti+1, mo_num if(bannedOrb(putj, sp)) cycle if(banned(puti, putj, bant)) cycle ! rentable? @@ -1522,7 +1531,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call i_h_j(gen, det, N_int, hij) if (dabs(hij) < mo_integrals_threshold) cycle else - hij = hij_cache1(putj) - hij_cache2(putj) + hij = hij_cache1(putj,puti) - hij_cache2(putj,puti) if (dabs(hij) < mo_integrals_threshold) cycle hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) end if @@ -1532,9 +1541,12 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) enddo end do end do + + deallocate(hij_cache1,hij_cache2) + end if - deallocate(hij_cache1,hij_cache2) +! deallocate(hij_cache1,hij_cache2) end diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 5a005d7b..8f3c72d0 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -262,8 +262,14 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map) out_val, 1) else integer :: isplit - out_val = 0.d0 - do isplit=1,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) + 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, & @@ -442,16 +448,36 @@ subroutine get_mo_two_e_integrals_ij(k,l,sze,out_array,map) double precision, intent(out) :: out_array(sze,sze) type(map_type), intent(inout) :: map integer :: j - real(integral_kind), allocatable :: tmp_val(:) if ( (mo_integrals_cache_min>1).or.(mo_integrals_cache_max Date: Thu, 13 Feb 2025 15:05:27 +0100 Subject: [PATCH 3/4] Add possibility of single precision cholesky vectors --- src/ao_two_e_ints/cholesky.irp.f | 2 +- src/mo_two_e_ints/EZFIO.cfg | 6 + src/mo_two_e_ints/cholesky.irp.f | 6 +- src/mo_two_e_ints/map_integrals.irp.f | 193 +++++++++++++++----------- 4 files changed, 121 insertions(+), 86 deletions(-) diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index eb2f092a..d15ebdc3 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/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 diff --git a/src/mo_two_e_ints/EZFIO.cfg b/src/mo_two_e_ints/EZFIO.cfg index da9d8fc9..f26bfb61 100644 --- a/src/mo_two_e_ints/EZFIO.cfg +++ b/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 + diff --git a/src/mo_two_e_ints/cholesky.irp.f b/src/mo_two_e_ints/cholesky.irp.f index 835110de..062e52e2 100644 --- a/src/mo_two_e_ints/cholesky.irp.f +++ b/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 i=1,mo_num + do j=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 diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 8f3c72d0..bc83fbdd 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/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 - do isplit=1,4 - get_two_e_integral = get_two_e_integral + & - sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), & - cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, & - cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),j,l), 1) - enddo + get_two_e_integral = sdot(cholesky_mo_num, cholesky_mo_transp_sp(1,i,k), 1, cholesky_mo_transp_sp(1,j,l), 1) +! get_two_e_integral = 0.d0 +! do isplit=1,4 +! get_two_e_integral = get_two_e_integral + & +! sdot(cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), & +! cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,k), 1, & +! 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_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) + call sgemv('T', cholesky_mo_num, mo_integrals_cache_min-1, 1., & + cholesky_mo_transp_sp(1,1,k), cholesky_mo_num, & + cholesky_mo_transp_sp(1,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) - 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 +! integer :: isplit +! isplit=1 +! 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) +! 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 - do isplit=1,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 + call sgemv('T', cholesky_mo_num, mo_num-mo_integrals_cache_max, 1., & + cholesky_mo_transp_sp(1,mo_integrals_cache_max+1,k), cholesky_mo_num, & + cholesky_mo_transp_sp(1,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) +! isplit=1 +! 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) +! 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 - do isplit=1,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 + call sgemv('T', cholesky_mo_num, sze, 1., & + cholesky_mo_transp_sp(1,1,k), cholesky_mo_num, & + cholesky_mo_transp_sp(1,j,l), 1, 0., & + out_val_sp, 1) + out_val(1:sze) = out_val_sp(1:sze) +! isplit=1 +! 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) +! 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 - double precision, allocatable :: out_array_sp(:,:) + real, allocatable :: out_array_sp(:,:) allocate(out_array_sp(sze,sze)) - call sgemm('T', 'N', mo_num, mo_num, & - cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), 1.d0, & - 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.d0, & + call sgemm('T', 'N', mo_num, mo_num, cholesky_mo_num, 1.0, & + cholesky_mo_transp_sp(1,1,k), cholesky_mo_num, & + cholesky_mo_transp_sp(1,1,l), cholesky_mo_num, 0.0, & 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, & - cholesky_mo_num_split(isplit+1) - cholesky_mo_num_split(isplit), 1.d0, & - 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.d0, & - out_array_sp, sze) - out_array(1:sze,1:sze) = out_array(1:sze,1:sze) + out_array_sp(1:sze,1:sze) - enddo +! +! isplit=1 +! 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_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 - do isplit=1,4 - 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,k), 1, & - cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1) - enddo + out_val(i) = sdot(cholesky_mo_num, cholesky_mo_transp_sp(1,i,k), 1, cholesky_mo_transp_sp(1,i,l), 1) +! out_val(i) = 0.d0 +! do isplit=1,4 +! 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,k), 1, & +! 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 - do isplit=1,4 - 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,k), 1, & - cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1) - enddo + out_val(i) = sdot(cholesky_mo_num, cholesky_mo_transp_sp(1,i,k), 1, cholesky_mo_transp_sp(1,i,l), 1) +! out_val(i) = 0.d0 +! do isplit=1,4 +! 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,k), 1, & +! 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 - do isplit=1,4 - 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,k), 1, & - cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1) - enddo + out_val(i) = sdot(cholesky_mo_num, cholesky_mo_transp_sp(1,i,k), 1, cholesky_mo_transp_sp(1,i,l), 1) +! out_val(i) = 0.d0 +! do isplit=1,4 +! 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,k), 1, & +! cholesky_mo_transp_sp(cholesky_mo_num_split(isplit),i,l), 1) +! enddo enddo endif From bce68e7461440405f2cdaf8462a97163ed89899a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 13 Feb 2025 16:06:05 +0100 Subject: [PATCH 4/4] Simplified hij_cache in pt2 --- src/cipsi/selection.irp.f | 26 +++++++++----------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 8c22ec85..99bc7013 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -1478,19 +1478,17 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer, parameter :: bant=1 -! double precision, allocatable :: hij_cache1(:), hij_cache2(:) - double precision, allocatable :: hij_cache1(:,:), hij_cache2(:,:) -! allocate (hij_cache1(mo_num),hij_cache2(mo_num)) + double precision, allocatable :: hij_cache(:,:) PROVIDE mo_integrals_threshold - if(sp == 3) then ! AB + allocate(hij_cache(mo_num,mo_num)) - allocate(hij_cache1(mo_num,mo_num)) + if(sp == 3) then ! AB h1 = p(1,1) h2 = p(1,2) - call get_mo_two_e_integrals_ij(h2,h1,mo_num,hij_cache1,mo_integrals_map) + call get_mo_two_e_integrals_ij(h2,h1,mo_num,hij_cache,mo_integrals_map) do p1=1, mo_num if(bannedOrb(p1, 1)) cycle do p2=1, mo_num @@ -1501,7 +1499,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call i_h_j(gen, det, N_int, hij) else phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) - hij = hij_cache1(p2,p1) * phase + hij = hij_cache(p2,p1) * phase end if if (dabs(hij) < mo_integrals_threshold) cycle !DIR$ LOOP COUNT AVG(4) @@ -1511,16 +1509,11 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) end do end do - deallocate(hij_cache1) - else ! AA BB - allocate(hij_cache1(mo_num,mo_num),hij_cache2(mo_num,mo_num)) - p1 = p(1,sp) p2 = p(2,sp) - call get_mo_two_e_integrals_ij(p2,p1,mo_num,hij_cache1,mo_integrals_map) - call get_mo_two_e_integrals_ij(p1,p2,mo_num,hij_cache2,mo_integrals_map) + call get_mo_two_e_integrals_ij(p2,p1,mo_num,hij_cache,mo_integrals_map) do puti=1, mo_num if (bannedOrb(puti, sp)) cycle do putj=puti+1, mo_num @@ -1531,7 +1524,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call i_h_j(gen, det, N_int, hij) if (dabs(hij) < mo_integrals_threshold) cycle else - hij = hij_cache1(putj,puti) - hij_cache2(putj,puti) + hij = hij_cache(putj,puti) - hij_cache(puti,putj) if (dabs(hij) < mo_integrals_threshold) cycle hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) end if @@ -1542,11 +1535,10 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) end do end do - deallocate(hij_cache1,hij_cache2) - end if -! deallocate(hij_cache1,hij_cache2) + deallocate(hij_cache) + end