From f42ffa77849704e300e626cb8212d573a68bf206 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 14 Sep 2017 12:43:18 +0200 Subject: [PATCH] Fixed integral8 --- plugins/Full_CI_ZMQ/selection.irp.f | 75 ++++++++---------------- src/Determinants/EZFIO.cfg | 2 +- src/Integrals_Bielec/map_integrals.irp.f | 15 +++-- 3 files changed, 34 insertions(+), 58 deletions(-) diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index acf19392..f404d069 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/plugins/Full_CI_ZMQ/selection.irp.f @@ -9,29 +9,6 @@ BEGIN_PROVIDER [ integer, fragment_count ] END_PROVIDER -double precision function integral8(i,j,k,l) - implicit none - - integer, intent(in) :: i,j,k,l - double precision, external :: get_mo_bielec_integral - integer :: ii - ii = l-mo_integrals_cache_min - ii = ior(ii, k-mo_integrals_cache_min) - ii = ior(ii, j-mo_integrals_cache_min) - ii = ior(ii, i-mo_integrals_cache_min) - if (iand(ii, -64) /= 0) then - integral8 = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) - else - ii = l-mo_integrals_cache_min - ii = ior( ishft(ii,6), k-mo_integrals_cache_min) - ii = ior( ishft(ii,6), j-mo_integrals_cache_min) - ii = ior( ishft(ii,6), i-mo_integrals_cache_min) - integral8 = mo_integrals_cache(ii) - endif -end function - - - subroutine assert(cond, msg) character(*), intent(in) :: msg logical, intent(in) :: cond @@ -135,7 +112,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2) integer :: i, j, h1, h2, p1, p2, sfix, hfix, pfix, hmob, pmob, puti double precision :: hij - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) integer, parameter :: turn2(2) = (/2,1/) @@ -148,7 +125,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) if(bannedOrb(puti)) cycle p1 = p(turn3_2(1,i), sp) p2 = p(turn3_2(2,i), sp) - hij = integral8(p1, p2, h1, h2) - integral8(p2, p1, h1, h2) + hij = mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2, p1, h1, h2) hij *= get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2) vect(:, puti) += hij * coefs end do @@ -161,7 +138,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) puti = p(j, sp) if(bannedOrb(puti)) cycle pmob = p(turn2(j), sp) - hij = integral8(pfix, pmob, hfix, hmob) + hij = mo_bielec_integral(pfix, pmob, hfix, hmob) hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix) vect(:, puti) += hij * coefs end do @@ -173,7 +150,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) p2 = p(2,sfix) h1 = h(1,sfix) h2 = h(2,sfix) - hij = (integral8(p1,p2,h1,h2) - integral8(p2,p1,h1,h2)) + hij = (mo_bielec_integral(p1,p2,h1,h2) - mo_bielec_integral(p2,p1,h1,h2)) hij *= get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2) vect(:, puti) += hij * coefs end if @@ -198,7 +175,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) logical, allocatable :: lbanned(:) integer(bit_kind) :: det(N_int, 2) double precision :: hij - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral allocate (lbanned(mo_tot_num)) lbanned = bannedOrb @@ -217,13 +194,13 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) do i=1,hole-1 if(lbanned(i)) cycle - hij = (integral8(p1, p2, i, hole) - integral8(p2, p1, i, hole)) + hij = (mo_bielec_integral(p1, p2, i, hole) - mo_bielec_integral(p2, p1, i, hole)) hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2) vect(1:N_states,i) += hij * coefs(1:N_states) end do do i=hole+1,mo_tot_num if(lbanned(i)) cycle - hij = (integral8(p1, p2, hole, i) - integral8(p2, p1, hole, i)) + hij = (mo_bielec_integral(p1, p2, hole, i) - mo_bielec_integral(p2, p1, hole, i)) hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2) vect(1:N_states,i) += hij * coefs(1:N_states) end do @@ -235,7 +212,7 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) p2 = p(1, sh) do i=1,mo_tot_num if(lbanned(i)) cycle - hij = integral8(p1, p2, i, hole) + hij = mo_bielec_integral(p1, p2, i, hole) hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2) vect(1:N_states,i) += hij * coefs(1:N_states) end do @@ -722,7 +699,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral integer :: i, j, tip, ma, mi, puti, putj integer :: h1, h2, p1, p2, i1, i2 @@ -757,7 +734,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) h1 = h(1, ma) h2 = h(2, ma) - hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) + hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) if(ma == 1) then mat(:, putj, puti) += coefs * hij else @@ -776,7 +753,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(banned(puti,putj,bant)) cycle p1 = p(turn2(i), 1) - hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) + hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) mat(:, puti, putj) += coefs * hij end do end do @@ -796,7 +773,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) i2 = turn2d(2, i, j) p1 = p(i1, ma) p2 = p(i2, ma) - hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) + hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) mat(:, puti, putj) += coefs * hij end do end do @@ -810,7 +787,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(banned(puti,putj,1)) cycle p2 = p(i, ma) - hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2) + hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2) mat(:, min(puti, putj), max(puti, putj)) += coefs * hij end do else ! tip == 4 @@ -821,7 +798,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) p2 = p(2, mi) h1 = h(1, mi) h2 = h(2, mi) - hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2) + hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2) mat(:, puti, putj) += coefs * hij end if end if @@ -841,7 +818,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp double precision :: hij, tmp_row(N_states, mo_tot_num), tmp_row2(N_states, mo_tot_num) - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral logical :: ok logical, allocatable :: lbanned(:,:) @@ -881,12 +858,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) tmp_row = 0d0 do putj=1, hfix-1 if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle - hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) + hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) tmp_row(1:N_states,putj) += hij * coefs(1:N_states) end do do putj=hfix+1, mo_tot_num if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle - hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) + hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) tmp_row(1:N_states,putj) += hij * coefs(1:N_states) end do @@ -906,13 +883,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) !p1 fixed putj = p1 if(.not. banned(putj,puti,bant)) then - hij = integral8(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix) + hij = mo_bielec_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix) tmp_row(:,puti) += hij * coefs end if putj = p2 if(.not. banned(putj,puti,bant)) then - hij = integral8(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix) + hij = mo_bielec_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix) tmp_row2(:,puti) += hij * coefs end if end do @@ -934,12 +911,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) tmp_row = 0d0 do putj=1,hfix-1 if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle - hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) + hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) tmp_row(:,putj) += hij * coefs end do do putj=hfix+1,mo_tot_num if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle - hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) + hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) tmp_row(:,putj) += hij * coefs end do @@ -957,13 +934,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(lbanned(puti,ma)) cycle putj = p2 if(.not. banned(puti,putj,1)) then - hij = integral8(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1) + hij = mo_bielec_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1) tmp_row(:,puti) += hij * coefs end if putj = p1 if(.not. banned(puti,putj,1)) then - hij = integral8(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2) + hij = mo_bielec_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2) tmp_row2(:,puti) += hij * coefs end if end do @@ -1015,7 +992,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer :: i, j, s, h1, h2, p1, p2, puti, putj double precision :: hij, phase - double precision, external :: get_phase_bi, integral8 + double precision, external :: get_phase_bi, mo_bielec_integral logical :: ok integer :: bant @@ -1035,7 +1012,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) - hij = integral8(p1, p2, h1, h2) * phase + hij = mo_bielec_integral(p1, p2, h1, h2) * phase end if mat(:, p1, p2) += coefs(:) * hij end do @@ -1052,7 +1029,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int) call i_h_j(gen, det, N_int, hij) else - hij = (integral8(p1, p2, puti, putj) - integral8(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2) + hij = (mo_bielec_integral(p1, p2, puti, putj) - mo_bielec_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2) end if mat(:, puti, putj) += coefs(:) * hij end do diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index f4c5d866..6bf6faff 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -50,7 +50,7 @@ default: 0.99 type: Threshold doc: Thresholds on selectors (fraction of the norm) interface: ezfio,provider,ocaml -default: 0.999 +default: 0.995 [n_int] interface: ezfio diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 00f43785..996f8464 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -356,13 +356,13 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:128_8*128_8*128_8*12 FREE ao_integrals_cache !$OMP PARALLEL DO PRIVATE (i,j,k,l,i4,j4,k4,l4,idx,ii,integral) do l=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + l4 = int(l,4) do k=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + k4 = int(k,4) do j=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + j4 = int(j,4) do i=mo_integrals_cache_min_8,mo_integrals_cache_max_8 i4 = int(i,4) - j4 = int(j,4) - k4 = int(k,4) - l4 = int(l,4) !DIR$ FORCEINLINE call bielec_integrals_index(i4,j4,k4,l4,idx) !DIR$ FORCEINLINE @@ -399,17 +399,16 @@ double precision function get_mo_bielec_integral(i,j,k,l,map) ii = ior(ii, j-mo_integrals_cache_min) ii = ior(ii, i-mo_integrals_cache_min) if (iand(ii, -128) /= 0) then -! if (.True.) then !DIR$ FORCEINLINE call bielec_integrals_index(i,j,k,l,idx) !DIR$ FORCEINLINE call map_get(map,idx,tmp) get_mo_bielec_integral = dble(tmp) else - ii_8 = l-mo_integrals_cache_min_8 - ii_8 = ior( ishft(ii_8,7), k-mo_integrals_cache_min_8) - ii_8 = ior( ishft(ii_8,7), j-mo_integrals_cache_min_8) - ii_8 = ior( ishft(ii_8,7), i-mo_integrals_cache_min_8) + ii_8 = int(l,8)-mo_integrals_cache_min_8 + ii_8 = ior( ishft(ii_8,7), int(k,8)-mo_integrals_cache_min_8) + ii_8 = ior( ishft(ii_8,7), int(j,8)-mo_integrals_cache_min_8) + ii_8 = ior( ishft(ii_8,7), int(i,8)-mo_integrals_cache_min_8) get_mo_bielec_integral = mo_integrals_cache(ii_8) endif end