From 5150497318b800ca332b9e407790112369362432 Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Tue, 20 Feb 2018 15:51:53 +0100 Subject: [PATCH] working but super slow mrcc_sto --- plugins/Full_CI_ZMQ/selection.irp.f | 4 +- plugins/dress_zmq/alpha_factory.irp.f | 44 ++++++++----- plugins/dress_zmq/dressing_vector.irp.f | 13 ++-- plugins/mrcc_sto/mrcc_sto.irp.f | 85 +++++++++++++++++++------ plugins/mrcepa0/dressing.irp.f | 22 ++++++- plugins/mrcepa0/dressing_vector.irp.f | 9 +++ 6 files changed, 136 insertions(+), 41 deletions(-) diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index 378e51c4..acda9fa6 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/plugins/Full_CI_ZMQ/selection.irp.f @@ -71,7 +71,6 @@ subroutine select_connected(i_generator,E0,pt2,b,subset) hole_mask(k,2) = iand(generators_bitmask(k,2,s_hole,l), psi_det_generators(k,2,i_generator)) particle_mask(k,1) = iand(generators_bitmask(k,1,s_part,l), not(psi_det_generators(k,1,i_generator)) ) particle_mask(k,2) = iand(generators_bitmask(k,2,s_part,l), not(psi_det_generators(k,2,i_generator)) ) - enddo call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b,subset) enddo @@ -299,7 +298,8 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), particle_mask(k,1)) particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), particle_mask(k,2)) enddo - + + integer :: N_holes(2), N_particles(2) integer :: hole_list(N_int*bit_kind_size,2) integer :: particle_list(N_int*bit_kind_size,2) diff --git a/plugins/dress_zmq/alpha_factory.irp.f b/plugins/dress_zmq/alpha_factory.irp.f index a70c22ae..7f2952b7 100644 --- a/plugins/dress_zmq/alpha_factory.irp.f +++ b/plugins/dress_zmq/alpha_factory.irp.f @@ -31,6 +31,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index integer :: h1,h2,s1,s2,s3,i1,i2,ib,sp,k,i,j,nt,ii,n integer(bit_kind) :: hole(N_int,2), particle(N_int,2), mask(N_int, 2), pmask(N_int, 2) + integer(bit_kind) :: mmask(N_int, 2) logical :: fullMatch, ok integer(bit_kind) :: mobMask(N_int, 2), negMask(N_int, 2) @@ -58,11 +59,24 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index monoBdo = .true. do k=1,N_int - hole (k,1) = iand(psi_det_generators(k,1,i_generator), generators_bitmask(k,1,s_hole,bitmask_index)) - hole (k,2) = iand(psi_det_generators(k,2,i_generator), generators_bitmask(k,2,s_hole,bitmask_index)) - particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), generators_bitmask(k,1,s_part,bitmask_index)) - particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), generators_bitmask(k,2,s_part,bitmask_index)) + !hole (k,1) = iand(psi_det_generators(k,1,i_generator), generators_bitmask(k,1,s_hole,bitmask_index)) + !hole (k,2) = iand(psi_det_generators(k,2,i_generator), generators_bitmask(k,2,s_hole,bitmask_index)) + !particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), generators_bitmask(k,1,s_part,bitmask_index)) + !particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), generators_bitmask(k,2,s_part,bitmask_index)) + hole (k,1) = iand(psi_det_generators(k,1,i_generator), full_ijkl_bitmask(k)) + hole (k,2) = iand(psi_det_generators(k,2,i_generator), full_ijkl_bitmask(k)) + particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), full_ijkl_bitmask(k)) + particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), full_ijkl_bitmask(k)) + enddo + + !if(i_generator == 34) then + ! call debug_det(psi_det_generators(1,1,34), N_int) + ! call debug_det(generators_bitmask(1,1,s_part,bitmask_index), N_int) + ! call debug_det(particle, N_int) + ! print *, "dddddddddddd" + ! stop + !end if integer :: N_holes(2), N_particles(2) integer :: hole_list(N_int*bit_kind_size,2) @@ -335,10 +349,10 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index if(siz > size(abuf)) stop "buffer too small in alpha_factory" abuf = 0 call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, indexes_end, abuf, interesting) - indexes_end(:,:) -= 1 - do i=1,siz - if(abuf(i) < 1 .or. abuf(i) > N_det) stop "foireous abuf" - end do + !indexes_end(:,:) -= 1 + !do i=1,siz + ! if(abuf(i) < 1 .or. abuf(i) > N_det) stop "foireous abuf" + !end do !print *, "IND1", indexes(1,:) !print *, "IND2", indexes_end(1,:) !stop @@ -374,7 +388,7 @@ subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexe allocate(labuf(N_det), putten(N_det)) putten = .false. - st1 = indexes_end(0,0) + st1 = indexes_end(0,0)-1 !! if(st1 > 0) labuf(:st1) = abuf(:st1) st1 += 1 @@ -382,19 +396,19 @@ subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexe s1 = 1 s2 = 2 lindex(:, 1) = indexes(1:,0) - lindex_end(:,1) = indexes_end(1:,0) + lindex_end(:,1) = indexes_end(1:,0)-1 lindex(:, 2) = indexes(0, 1:) - lindex_end(:, 2) = indexes_end(0, 1:) + lindex_end(:, 2) = indexes_end(0, 1:)-1 else if(sp == 2) then s1 = 2 s2 = 2 lindex(:, 2) = indexes(0, 1:) - lindex_end(:, 2) = indexes_end(0, 1:) + lindex_end(:, 2) = indexes_end(0, 1:)-1 else if(sp == 1) then s1 = 1 s2 = 1 lindex(:, 1) = indexes(1:, 0) - lindex_end(:,1) = indexes_end(1:, 0) + lindex_end(:,1) = indexes_end(1:, 0)-1 end if do i=1,mo_tot_num @@ -432,8 +446,8 @@ subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexe end if if(indexes(i,j) /= 0) then - st4 = st3 + 1 + indexes_end(i,j)-indexes(i,j) - labuf(st3:st4-1) = abuf(indexes(i,j):indexes_end(i,j)) + st4 = st3 + 1 + indexes_end(i,j)-indexes(i,j) -1!! + labuf(st3:st4-1) = abuf(indexes(i,j):indexes_end(i,j)-1) !! else st4 = st3 end if diff --git a/plugins/dress_zmq/dressing_vector.irp.f b/plugins/dress_zmq/dressing_vector.irp.f index db1f428b..0d24ed94 100644 --- a/plugins/dress_zmq/dressing_vector.irp.f +++ b/plugins/dress_zmq/dressing_vector.irp.f @@ -1,3 +1,9 @@ + BEGIN_PROVIDER [ integer, nalp ] +&BEGIN_PROVIDER [ integer, ninc ] + nalp = 0 + ninc = 0 + END_PROVIDER + BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ] &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ] implicit none @@ -11,9 +17,6 @@ double precision :: f, tmp double precision, external :: u_dot_v - print *, "DELTA_IJ", delta_ij(1,:10,1) - print *, "DELTA_IJ div", delta_ij(1,:10,1) / psi_coef(dressed_column_idx(1),1) - !stop do k=1,N_states l = dressed_column_idx(k) f = 1.d0/psi_coef(l,k) @@ -27,6 +30,8 @@ tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) dressing_column_s(l,k) -= tmp * f enddo - !stop + print *, "NALP", nalp + print *, "NINC", ninc + print *, "DELTA_IJ", dressing_column_h(:10,1) END_PROVIDER diff --git a/plugins/mrcc_sto/mrcc_sto.irp.f b/plugins/mrcc_sto/mrcc_sto.irp.f index 37b40a0e..cde31adb 100644 --- a/plugins/mrcc_sto/mrcc_sto.irp.f +++ b/plugins/mrcc_sto/mrcc_sto.irp.f @@ -4,18 +4,11 @@ program mrcc_sto BEGIN_DOC ! TODO END_DOC - print *, "========================" - print *, "========================" - print *, "========================" - print *, "MRCC_STO not implemented - acts as a unittest for dress_zmq" - print *, "========================" - print *, "========================" - print *, "========================" call dress_zmq() end - BEGIN_PROVIDER [ integer, idx_non_ref_from_sorted, (N_det) ] -&BEGIN_PROVIDER [ integer, psi_from_sorted, (N_det) ] + BEGIN_PROVIDER [ integer, psi_from_sorted, (N_det) ] +&BEGIN_PROVIDER [ integer, idx_non_ref_from_sorted, (N_det) ] implicit none integer :: i,inpsisor @@ -30,10 +23,17 @@ end end do END_PROVIDER + subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha) use bitmasks implicit none - + BEGIN_DOC + !delta_ij_loc(:,:,1) : dressing column for H + !delta_ij_loc(:,:,2) : dressing column for S2 + !minilist : indices of determinants connected to alpha ( in psi_det_sorted ) + !n_minilist : size of minilist + !alpha : alpha determinant + END_DOC integer(bit_kind), intent(in) :: alpha(N_int,2) integer,intent(in) :: minilist(n_minilist), n_minilist double precision, intent(inout) :: delta_ij_loc(N_states,N_det,2) @@ -49,22 +49,36 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha) double precision :: ci_inv(N_states) integer :: exc(0:2,2,2) integer :: h1,h2,p1,p2,s1,s2 - integer(bit_kind) :: tmp_det(N_int,2) + integer(bit_kind) :: tmp_det(N_int,2), ctrl integer :: i_state, k_sd, l_sd, m_sd, ll_sd, i_I double precision, allocatable :: hij_cache(:), sij_cache(:) double precision :: Delta_E_inv(N_states) double precision :: sdress, hdress double precision :: c0(N_states) - logical :: ok + logical :: ok, ok2 + integer :: old_ninc + double precision :: shdress + if(n_minilist == 1) return + + shdress = 0d0 + old_ninc = ninc + if (perturbative_triples) then PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat endif + + do i_I=1,N_det_ref + call get_excitation_degree(alpha,psi_ref(1,1,i_I),degree1,N_int) + if(degree1 <= 2) return + end do + allocate(hij_cache(N_det), sij_cache(N_det)) allocate (dIa_hla(N_states,N_det), dIa_sla(N_states,N_det)) allocate (idx_alpha(0:n_minilist)) + do i_state=1,N_states c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state) enddo @@ -81,18 +95,31 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha) exit end if end do + + !if(ok) then + ! call get_excitation(psi_det_sorted(1,1,k_sd),alpha,exc,degree1,phase,N_int) + ! if(degree1 == 0 .or. degree1 > 2) stop "minilist error" + ! call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2) + ! + ! if(h1 > 10 .or. p1 < 7 .or. p1 == 8 .or. p1 == 9) ok = .false. + ! if(ok .and. degree1 == 2) then + ! if(h2 > 10 .or. p2 < 7 .or. p2 == 8 .or. p2 == 9) ok = .false. + ! end if + ! !if(degree1 == 0 .or. degree1 > 2) stop "minilist error" + ! !iand(xor(psi_det_sorted(i,2,k_sd), alpha(i,2)), alpha(i,2)) + !end if - if( xor(ok, idx_non_ref_from_sorted(k_sd) > 0)) stop "BUGUE" + !if( xor(ok, idx_non_ref_from_sorted(k_sd) > 0)) stop "BUGUE" if(ok) then ll_sd += 1 idx_alpha(ll_sd) = k_sd ! call i_h_j(alpha,psi_non_ref(1,1,idx_alpha(l_sd)),N_int,hij_cache(k_sd)) ! call get_s2(alpha,psi_non_ref(1,1,idx_alpha(l_sd)),N_int,sij_cache(k_sd)) - call i_h_j(alpha,psi_det_sorted(1,1,k_sd),N_int,hij_cache(k_sd)) - call get_s2(alpha,psi_det_sorted(1,1,k_sd),N_int,sij_cache(k_sd)) + call i_h_j(alpha,psi_det_sorted(1,1,k_sd),N_int,hij_cache(k_sd)) + call get_s2(alpha,psi_det_sorted(1,1,k_sd),N_int,sij_cache(k_sd)) end if enddo - + if(ll_sd <= 1) return idx_alpha(0) = ll_sd @@ -126,9 +153,12 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha) enddo call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, N_int) + ok2 = .false. do i_state=1,N_states dIK(i_state) = dij(i_I, idx_non_ref_from_sorted(idx_alpha(k_sd)), i_state) + if(dIK(i_state) /= 0d0) ok2 = .true. enddo + if(.not. ok2) cycle ! do i_state=1,N_states @@ -174,6 +204,12 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha) enddo enddo + ok2 = .false. + do i_state=1,N_states + if(dIa(i_state) /= 0d0) ok2 = .true. + enddo + if(.not. ok2) cycle + do i_state=1,N_states ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state) enddo @@ -192,9 +228,15 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha) !print *, i_state, idx_alpha(l_sd) k_sd = idx_alpha(l_sd) m_sd = psi_from_sorted(k_sd) - if(psi_det(1,1,m_sd) /= psi_det_sorted(1,1,k_sd)) stop "psi_from_sorted foireous" - hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state) - sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state) + hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state)! * c0(i_state) + sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state)! * c0(i_state) + !!$OMP ATOMIC + !shdress += 1d0 + nalp += 1 + if(hdress /= 0d0) then + ninc = ninc + 1 + !print *, "grepme2", hdress, shdress + end if !$OMP ATOMIC delta_ij_loc(i_state,m_sd,1) += hdress !$OMP ATOMIC @@ -203,6 +245,11 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha) enddo enddo enddo + + !if(ninc /= old_ninc) then + ! nalp = nalp + 1 + ! !print "(A8,I20,I20,E15.5)", "grepme", alpha(1,1), alpha(1,2), shdress + !end if end subroutine diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index a7f40b5b..1fc024b5 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -343,7 +343,10 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b !double precision, external :: get_dij, get_dij_index double precision :: Delta_E_inv(N_states) double precision, intent(inout) :: contrib(N_states) - double precision :: sdress, hdress + double precision :: sdress, hdress, shdress + integer :: old_ninc + + old_ninc = ninc if (perturbative_triples) then PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat @@ -408,6 +411,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b do i_alpha=1,N_tq + old_ninc = ninc + shdress = 0d0 + if(key_mask(1,1) /= 0) then call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) @@ -545,6 +551,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b k_sd = idx_alpha(l_sd) hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) + !!$OMP ATOMIC + if(hdress /= 0d0) ninc = ninc + 1 + shdress += hdress !$OMP ATOMIC contrib(i_state) += hdress * psi_coef(dressed_column_idx(i_state), i_state) * psi_non_ref_coef(k_sd, i_state) !$OMP ATOMIC @@ -554,7 +563,18 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b enddo enddo enddo + if(ninc /= old_ninc) then + nalp = nalp + 1 + !print "(A8,I20,I20,E15.5)", "grepme", tq(1,1,i_alpha), tq(1,2,i_alpha), shdress + !if(tq(1,1,i_alpha) == 1007 .and. tq(1,2,i_alpha) == 17301943) then + ! print *, "foinder", i_generator + ! call debug_det(psi_det_generators(1,1, i_generator), N_int) + ! call debug_det(tq(1,1,i_alpha), N_int) + ! stop + ! end if + end if enddo + deallocate (dIa_hla,dIa_sla,hij_cache,sij_cache) deallocate(miniList, idx_miniList) end diff --git a/plugins/mrcepa0/dressing_vector.irp.f b/plugins/mrcepa0/dressing_vector.irp.f index 6e6602ea..22394a92 100644 --- a/plugins/mrcepa0/dressing_vector.irp.f +++ b/plugins/mrcepa0/dressing_vector.irp.f @@ -1,3 +1,10 @@ + + BEGIN_PROVIDER [ integer, nalp ] +&BEGIN_PROVIDER [ integer, ninc ] + nalp = 0 + ninc = 0 +END_PROVIDER + BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ] &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ] implicit none @@ -24,6 +31,8 @@ tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) dressing_column_s(l,k) -= tmp * f enddo + print *, "NALP", nalp + print *, "NINC", ninc print *, "DRESS", dressing_column_h(:10,1) ! stop END_PROVIDER