From 7ef8ee171da211f4b6f973c1f638f7adc90534b2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 29 Jul 2015 15:35:01 +0200 Subject: [PATCH] Accelerated get_excitation_degree_vector --- plugins/MRCC_Utils_new/mrcc_amplitudes.irp.f | 170 ++++++++++--------- src/Determinants/slater_rules.irp.f | 46 +++-- 2 files changed, 118 insertions(+), 98 deletions(-) diff --git a/plugins/MRCC_Utils_new/mrcc_amplitudes.irp.f b/plugins/MRCC_Utils_new/mrcc_amplitudes.irp.f index 8d1bce33..6746bee1 100644 --- a/plugins/MRCC_Utils_new/mrcc_amplitudes.irp.f +++ b/plugins/MRCC_Utils_new/mrcc_amplitudes.irp.f @@ -1,85 +1,93 @@ subroutine get_excitation_operators_for_one_ref(det_ref,i_state,ndetnonref,N_connect_ref,excitation_operators,amplitudes_phase_less,index_connected) - use bitmasks - implicit none - integer(bit_kind), intent(in) :: det_ref(N_int,2) - integer, intent(in) :: i_state,ndetnonref - integer*2, intent(out) :: excitation_operators(5,ndetnonref) - integer, intent(out) :: index_connected(ndetnonref) - integer, intent(out) :: N_connect_ref - double precision, intent(out) :: amplitudes_phase_less(ndetnonref) - - integer :: i,j,k,l,degree,h1,p1,h2,p2,s1,s2 - integer :: exc(0:2,2,2) - double precision :: phase,hij - BEGIN_DOC -! This subroutine provides all the amplitudes and excitation operators -! that one needs to go from the reference to the non reference wave function -! you enter with det_ref that is a reference determinant -! -! N_connect_ref is the number of determinants belonging to psi_non_ref -! that are connected to det_ref. -! -! amplitudes_phase_less(i) = amplitude phase less t_{I->i} = * lambda_mrcc(i) * phase(I->i) -! -! excitation_operators(:,i) represents the holes and particles that -! link the ith connected determinant to det_ref -! if :: -! excitation_operators(5,i) = 2 :: double excitation alpha -! excitation_operators(5,i) = -2 :: double excitation beta -!!! excitation_operators(1,i) :: hole 1 -!!! excitation_operators(2,i) :: particle 1 -!!! excitation_operators(3,i) :: hole 2 -!!! excitation_operators(4,i) :: particle 2 -! else if :: -! excitation_operators(5,i) = 1 :: single excitation alpha -!!! excitation_operators(1,i) :: hole 1 -!!! excitation_operators(2,i) :: particle 1 -! else if :: -! excitation_operators(5,i) = -1 :: single excitation beta -!!! excitation_operators(3,i) :: hole 1 -!!! excitation_operators(4,i) :: particle 1 -! else if :: -!!! excitation_operators(5,i) = 0 :: double excitation alpha/beta -!!! excitation_operators(1,i) :: hole 1 alpha -!!! excitation_operators(2,i) :: particle 1 alpha -!!! excitation_operators(3,i) :: hole 2 beta -!!! excitation_operators(4,i) :: particle 2 beta - END_DOC - N_connect_ref = 0 - do i = 1, ndetnonref - call i_H_j_phase_out(det_ref,psi_non_ref(1,1,i),N_int,hij,phase,exc,degree) -! if(dabs(hij).le.mo_integrals_threshold)cycle - N_connect_ref +=1 - index_connected(N_connect_ref) = i - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - amplitudes_phase_less(N_connect_ref) = hij * lambda_mrcc(i_state,i) !*phase - - if(degree==2)then - excitation_operators(1,N_connect_ref) = h1 - excitation_operators(2,N_connect_ref) = p1 - excitation_operators(3,N_connect_ref) = h2 - excitation_operators(4,N_connect_ref) = p2 - if(s1==s2.and.s1==1)then ! double alpha - excitation_operators(5,N_connect_ref)= 2 - elseif(s1==s2.and.s1==2)then ! double beta - excitation_operators(5,N_connect_ref)=-2 - else - excitation_operators(5,N_connect_ref)= 0 ! double alpha/beta - endif - elseif(degree==1)then - if(s1==1)then ! mono alpha - excitation_operators(5,N_connect_ref)= 1 - excitation_operators(1,N_connect_ref) = h1 - excitation_operators(2,N_connect_ref) = p1 - else ! mono beta - excitation_operators(5,N_connect_ref)=-1 - excitation_operators(3,N_connect_ref) = h1 - excitation_operators(4,N_connect_ref) = p1 - endif - else - N_connect_ref-=1 - endif + use bitmasks + implicit none + integer(bit_kind), intent(in) :: det_ref(N_int,2) + integer, intent(in) :: i_state,ndetnonref + integer*2, intent(out) :: excitation_operators(5,ndetnonref) + integer, intent(out) :: index_connected(ndetnonref) + integer, intent(out) :: N_connect_ref + double precision, intent(out) :: amplitudes_phase_less(ndetnonref) - enddo + integer :: i,j,k,l,degree,h1,p1,h2,p2,s1,s2 + integer :: exc(0:2,2,2) + double precision :: phase,hij + BEGIN_DOC + ! This subroutine provides all the amplitudes and excitation operators + ! that one needs to go from the reference to the non reference wave function + ! you enter with det_ref that is a reference determinant + ! + ! N_connect_ref is the number of determinants belonging to psi_non_ref + ! that are connected to det_ref. + ! + ! amplitudes_phase_less(i) = amplitude phase less t_{I->i} = * lambda_mrcc(i) * phase(I->i) + ! + ! excitation_operators(:,i) represents the holes and particles that + ! link the ith connected determinant to det_ref + ! if :: + ! excitation_operators(5,i) = 2 :: double excitation alpha + ! excitation_operators(5,i) = -2 :: double excitation beta + !!! excitation_operators(1,i) :: hole 1 + !!! excitation_operators(2,i) :: particle 1 + !!! excitation_operators(3,i) :: hole 2 + !!! excitation_operators(4,i) :: particle 2 + ! else if :: + ! excitation_operators(5,i) = 1 :: single excitation alpha + !!! excitation_operators(1,i) :: hole 1 + !!! excitation_operators(2,i) :: particle 1 + ! else if :: + ! excitation_operators(5,i) = -1 :: single excitation beta + !!! excitation_operators(3,i) :: hole 1 + !!! excitation_operators(4,i) :: particle 1 + ! else if :: + !!! excitation_operators(5,i) = 0 :: double excitation alpha/beta + !!! excitation_operators(1,i) :: hole 1 alpha + !!! excitation_operators(2,i) :: particle 1 alpha + !!! excitation_operators(3,i) :: hole 2 beta + !!! excitation_operators(4,i) :: particle 2 beta + END_DOC + N_connect_ref = 0 + do i = 1, ndetnonref + call i_H_j_phase_out(det_ref,psi_non_ref(1,1,i),N_int,hij,phase,exc,degree) + if (dabs(hij) <= mo_integrals_threshold) then + cycle + endif + N_connect_ref +=1 + index_connected(N_connect_ref) = i + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + amplitudes_phase_less(N_connect_ref) = hij * lambda_mrcc(i_state,i) !*phase + + if (degree==2) then + excitation_operators(1,N_connect_ref) = h1 + excitation_operators(2,N_connect_ref) = p1 + excitation_operators(3,N_connect_ref) = h2 + excitation_operators(4,N_connect_ref) = p2 + if(s1==s2.and.s1==1)then ! double alpha + excitation_operators(5,N_connect_ref) = 2 + elseif(s1==s2.and.s1==2)then ! double beta + excitation_operators(5,N_connect_ref) = -2 + else ! double alpha/beta + excitation_operators(5,N_connect_ref) = 0 + endif + + else if(degree==1) then + + if(s1==1)then ! mono alpha + excitation_operators(5,N_connect_ref) = 1 + excitation_operators(1,N_connect_ref) = h1 + excitation_operators(2,N_connect_ref) = p1 + else ! mono beta + excitation_operators(5,N_connect_ref) = -1 + excitation_operators(3,N_connect_ref) = h1 + excitation_operators(4,N_connect_ref) = p1 + endif + + else + + N_connect_ref-=1 + + endif + + enddo + end diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index b836d7c9..d0157ecc 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -960,7 +960,7 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) integer, intent(out) :: degree(sze) integer, intent(out) :: idx(0:sze) - integer :: i,l + integer :: i,l,d ASSERT (Nint > 0) ASSERT (sze > 0) @@ -970,9 +970,12 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) !DIR$ LOOP COUNT (1000) do i=1,sze - degree(l) = ishft(popcnt(xor( key1(1,1,i), key2(1,1))) + & + d = ishft(popcnt(xor( key1(1,1,i), key2(1,1))) + & popcnt(xor( key1(1,2,i), key2(1,2))),-1) - if (degree(l) < 3) then + if (d > 2) then + cycle + else + degree(l) = d idx(l) = i l = l+1 endif @@ -982,13 +985,16 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) !DIR$ LOOP COUNT (1000) do i=1,sze - degree(l) = ishft(popcnt(xor( key1(1,1,i), key2(1,1))) + & + d = ishft(popcnt(xor( key1(1,1,i), key2(1,1))) + & popcnt(xor( key1(1,2,i), key2(1,2))) + & popcnt(xor( key1(2,1,i), key2(2,1))) + & popcnt(xor( key1(2,2,i), key2(2,2))),-1) - if (degree(l) < 3) then - idx(l) = i - l = l+1 + if (d > 2) then + cycle + else + degree(l) = d + idx(l) = i + l = l+1 endif enddo @@ -996,15 +1002,18 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) !DIR$ LOOP COUNT (1000) do i=1,sze - degree(l) = ishft( popcnt(xor( key1(1,1,i), key2(1,1))) + & + d = ishft( popcnt(xor( key1(1,1,i), key2(1,1))) + & popcnt(xor( key1(1,2,i), key2(1,2))) + & popcnt(xor( key1(2,1,i), key2(2,1))) + & popcnt(xor( key1(2,2,i), key2(2,2))) + & popcnt(xor( key1(3,1,i), key2(3,1))) + & popcnt(xor( key1(3,2,i), key2(3,2))),-1) - if (degree(l) < 3) then - idx(l) = i - l = l+1 + if (d > 2) then + cycle + else + degree(l) = d + idx(l) = i + l = l+1 endif enddo @@ -1012,16 +1021,19 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) !DIR$ LOOP COUNT (1000) do i=1,sze - degree(l) = 0 + d = 0 !DEC$ LOOP COUNT MIN(4) do l=1,Nint - degree(l) = degree(l)+ popcnt(xor( key1(l,1,i), key2(l,1))) +& + d = d + popcnt(xor( key1(l,1,i), key2(l,1))) +& popcnt(xor( key1(l,2,i), key2(l,2))) enddo - degree(l) = ishft(degree(l),-1) - if (degree(l) < 3) then - idx(l) = i - l = l+1 + d = ishft(d,-1) + if (d > 2) then + cycle + else + degree(l) = d + idx(l) = i + l = l+1 endif enddo