From b1c7c121b21da7ece830289801931eed1d7f36cf Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 5 Jul 2019 15:39:27 +0200 Subject: [PATCH] working on pert rdms --- src/cipsi/pert_rdm_providers.irp.f | 166 ++++++++++++++++++ src/cipsi/selection.irp.f | 148 +--------------- src/cipsi/update_2rdm.irp.f | 52 ++++++ .../orb_range_routines_openmp.irp.f | 14 +- 4 files changed, 226 insertions(+), 154 deletions(-) create mode 100644 src/cipsi/pert_rdm_providers.irp.f diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f new file mode 100644 index 00000000..ad5355a2 --- /dev/null +++ b/src/cipsi/pert_rdm_providers.irp.f @@ -0,0 +1,166 @@ + +use bitmasks + +BEGIN_PROVIDER [logical , pert_2rdm ] + implicit none + pert_2rdm = .False. +END_PROVIDER + +BEGIN_PROVIDER [integer, n_orb_pert_rdm] + implicit none + n_orb_pert_rdm = n_act_orb +END_PROVIDER + +BEGIN_PROVIDER [integer, list_orb_reverse_pert_rdm, (mo_num)] + implicit none + list_orb_reverse_pert_rdm = list_act_reverse + +END_PROVIDER + +BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)] + implicit none + list_orb_pert_rdm = list_act + +END_PROVIDER + +subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf, psi_det_connection, psi_coef_connection, n_det_connection) + use bitmasks + use selection_types + implicit none + + integer, intent(in) :: n_det_connection + double precision, intent(in) :: psi_coef_connection(n_det_connection,N_states) + integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) + integer, intent(in) :: i_generator, sp, h1, h2 + double precision, intent(in) :: mat(N_states, mo_num, mo_num) + logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num) + double precision, intent(in) :: fock_diag_tmp(mo_num) + double precision, intent(in) :: E0(N_states) + double precision, intent(inout) :: pt2(N_states) + double precision, intent(inout) :: variance(N_states) + double precision, intent(inout) :: norm(N_states) + type(selection_buffer), intent(inout) :: buf + logical :: ok + integer :: s1, s2, p1, p2, ib, j, istate + integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) + double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states) + double precision, external :: diag_H_mat_elem_fock + double precision :: E_shift + + logical, external :: detEq + double precision, allocatable :: values(:) + integer, allocatable :: keys(:,:) + integer :: nkeys + integer :: sze_buff + sze_buff = 5 * mo_num ** 2 + allocate(keys(4,sze_buff),values(sze_buff)) + nkeys = 0 + if(sp == 3) then + s1 = 1 + s2 = 2 + else + s1 = sp + s2 = sp + end if + call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) + E_shift = 0.d0 + + if (h0_type == 'SOP') then + j = det_to_occ_pattern(i_generator) + E_shift = psi_det_Hii(i_generator) - psi_occ_pattern_Hii(j) + endif + + do p1=1,mo_num + if(bannedOrb(p1, s1)) cycle + ib = 1 + if(sp /= 3) ib = p1+1 + + do p2=ib,mo_num + +! ----- +! /!\ Generating only single excited determinants doesn't work because a +! determinant generated by a single excitation may be doubly excited wrt +! to a determinant of the future. In that case, the determinant will be +! detected as already generated when generating in the future with a +! double excitation. +! +! if (.not.do_singles) then +! if ((h1 == p1) .or. (h2 == p2)) then +! cycle +! endif +! endif +! +! if (.not.do_doubles) then +! if ((h1 /= p1).and.(h2 /= p2)) then +! cycle +! endif +! endif +! ----- + + if(bannedOrb(p2, s2)) cycle + if(banned(p1,p2)) cycle + + + if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle + call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) + + if (do_only_cas) then + integer, external :: number_of_holes, number_of_particles + if (number_of_particles(det)>0) then + cycle + endif + if (number_of_holes(det)>0) then + cycle + endif + endif + + if (do_ddci) then + logical, external :: is_a_two_holes_two_particles + if (is_a_two_holes_two_particles(det)) then + cycle + endif + endif + + if (do_only_1h1p) then + logical, external :: is_a_1h1p + if (.not.is_a_1h1p(det)) cycle + endif + + + Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) + + sum_e_pert = 0d0 + + do istate=1,N_states + delta_E = E0(istate) - Hii + E_shift + alpha_h_psi = mat(istate, p1, p2) + val = alpha_h_psi + alpha_h_psi + tmp = dsqrt(delta_E * delta_E + val * val) + if (delta_E < 0.d0) then + tmp = -tmp + endif + e_pert = 0.5d0 * (tmp - delta_E) + coef(istate) = e_pert / alpha_h_psi + pt2(istate) = pt2(istate) + e_pert + variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi + norm(istate) = norm(istate) + coef(istate) * coef(istate) + + if (weight_selection /= 5) then + ! Energy selection + sum_e_pert = sum_e_pert + e_pert * selection_weight(istate) + else + ! Variance selection + sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate) + endif + end do + + call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection,n_det_connection,nkeys,keys,values,sze_buff) + + if(sum_e_pert <= buf%mini) then + call add_to_selection_buffer(buf, det, sum_e_pert) + end if + end do + end do +end + + diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 19e53dcc..71442538 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -1,9 +1,5 @@ -use bitmasks -BEGIN_PROVIDER [logical , pert_2rdm ] - implicit none - pert_2rdm = .False. -END_PROVIDER +use bitmasks BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ] implicit none @@ -768,148 +764,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d end do end - -subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf, psi_det_connection, psi_coef_connection, n_det_connection) - use bitmasks - use selection_types - implicit none - - integer, intent(in) :: n_det_connection - double precision, intent(in) :: psi_coef_connection(n_det_connection,N_states) - integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) - integer, intent(in) :: i_generator, sp, h1, h2 - double precision, intent(in) :: mat(N_states, mo_num, mo_num) - logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num) - double precision, intent(in) :: fock_diag_tmp(mo_num) - double precision, intent(in) :: E0(N_states) - double precision, intent(inout) :: pt2(N_states) - double precision, intent(inout) :: variance(N_states) - double precision, intent(inout) :: norm(N_states) - type(selection_buffer), intent(inout) :: buf - logical :: ok - integer :: s1, s2, p1, p2, ib, j, istate - integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states) - double precision, external :: diag_H_mat_elem_fock - double precision :: E_shift - - logical, external :: detEq - double precision, allocatable :: values(:) - integer, allocatable :: keys(:,:) - integer :: nkeys - integer :: sze_buff - sze_buff = 5 * mo_num ** 2 - allocate(keys(4,sze_buff),values(sze_buff)) - nkeys = 0 - if(sp == 3) then - s1 = 1 - s2 = 2 - else - s1 = sp - s2 = sp - end if - call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) - E_shift = 0.d0 - - if (h0_type == 'SOP') then - j = det_to_occ_pattern(i_generator) - E_shift = psi_det_Hii(i_generator) - psi_occ_pattern_Hii(j) - endif - - do p1=1,mo_num - if(bannedOrb(p1, s1)) cycle - ib = 1 - if(sp /= 3) ib = p1+1 - - do p2=ib,mo_num - -! ----- -! /!\ Generating only single excited determinants doesn't work because a -! determinant generated by a single excitation may be doubly excited wrt -! to a determinant of the future. In that case, the determinant will be -! detected as already generated when generating in the future with a -! double excitation. -! -! if (.not.do_singles) then -! if ((h1 == p1) .or. (h2 == p2)) then -! cycle -! endif -! endif -! -! if (.not.do_doubles) then -! if ((h1 /= p1).and.(h2 /= p2)) then -! cycle -! endif -! endif -! ----- - - if(bannedOrb(p2, s2)) cycle - if(banned(p1,p2)) cycle - - - if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle - call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) - - if (do_only_cas) then - integer, external :: number_of_holes, number_of_particles - if (number_of_particles(det)>0) then - cycle - endif - if (number_of_holes(det)>0) then - cycle - endif - endif - - if (do_ddci) then - logical, external :: is_a_two_holes_two_particles - if (is_a_two_holes_two_particles(det)) then - cycle - endif - endif - - if (do_only_1h1p) then - logical, external :: is_a_1h1p - if (.not.is_a_1h1p(det)) cycle - endif - - - Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) - - sum_e_pert = 0d0 - - do istate=1,N_states - delta_E = E0(istate) - Hii + E_shift - alpha_h_psi = mat(istate, p1, p2) - val = alpha_h_psi + alpha_h_psi - tmp = dsqrt(delta_E * delta_E + val * val) - if (delta_E < 0.d0) then - tmp = -tmp - endif - e_pert = 0.5d0 * (tmp - delta_E) - coef(istate) = e_pert / alpha_h_psi - pt2(istate) = pt2(istate) + e_pert - variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi - norm(istate) = norm(istate) + coef(istate) * coef(istate) - - if (weight_selection /= 5) then - ! Energy selection - sum_e_pert = sum_e_pert + e_pert * selection_weight(istate) - else - ! Variance selection - sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate) - endif - end do - - call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection,n_det_connection,nkeys,keys,values,sze_buff) - - if(sum_e_pert <= buf%mini) then - call add_to_selection_buffer(buf, det, sum_e_pert) - end if - end do - end do -end - - subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting) use bitmasks implicit none diff --git a/src/cipsi/update_2rdm.irp.f b/src/cipsi/update_2rdm.irp.f index 07a746d4..3d6409af 100644 --- a/src/cipsi/update_2rdm.irp.f +++ b/src/cipsi/update_2rdm.irp.f @@ -9,5 +9,57 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio double precision, intent(in) :: psi_coef_connection(n_det_connection, N_states) integer, intent(inout) :: keys(4,sze_buff),sze_buff double precision, intent(inout) :: values(sze_buff) + integer :: i + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase, contrib + do i = 1, n_det_connection + call get_excitation(det,psi_det_connection(1,1,i),exc,degree,phase,N_int) + if(degree.gt.2)cycle + contrib = 0.d0 + do j = 1, N_states + contrib += state_average_weight(j) * psi_coef_connection(i,j) * phase * coef(j) + enddo + ! case of single excitations + if(degree == 1)then + if (nkeys+ 2 * elec_alpha_num .ge. sze_buff)then + call update_rdms(nkeys,keys,values,sze_buff) + nkeys = 0 + endif + call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff) + else + ! case of double excitations + if (nkeys+ 4 .ge. sze_buff)then + call update_rdms(nkeys,keys,values,sze_buff) + nkeys = 0 + endif + call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) + endif + enddo + +end + +subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff) + implicit none + integer, intent(in) :: nkeys,sze_buff + integer(bit_kind), intent(in) :: det1(N_int,2) + integer(bit_kind), intent(in) :: det2(N_int,2) + integer,intent(in) :: exc(0:2,2,2) + double precision,intent(in) :: phase, contrib + integer, intent(inout) :: nkeys, keys(4,sze_buff) + double precision, intent(inout):: values(sze_buff) + + + +end + +subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) + implicit none + integer, intent(in) :: nkeys,sze_buff + integer,intent(in) :: exc(0:2,2,2) + double precision,intent(in) :: phase, contrib + integer, intent(inout) :: nkeys, keys(4,sze_buff) + double precision, intent(inout):: values(sze_buff) + end diff --git a/src/two_body_rdm/orb_range_routines_openmp.irp.f b/src/two_body_rdm/orb_range_routines_openmp.irp.f index 97a8ce8a..269f789c 100644 --- a/src/two_body_rdm/orb_range_routines_openmp.irp.f +++ b/src/two_body_rdm/orb_range_routines_openmp.irp.f @@ -166,7 +166,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis !$OMP psi_bilinear_matrix_transp_order, N_st, & !$OMP psi_bilinear_matrix_order_transp_reverse, & !$OMP psi_bilinear_matrix_columns_loc, & - !$OMP psi_bilinear_matrix_transp_rows_loc,norb, & + !$OMP psi_bilinear_matrix_transp_rows_loc, & !$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, state_weights, dim1, & !$OMP ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask) & !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, c_2, & @@ -348,13 +348,13 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis enddo if(alpha_beta.or.spin_trace.or.alpha_alpha)then ! increment the alpha/beta part for single excitations - if (nkeys+ 2 * norb .ge. size(values)) then + if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) nkeys = 0 endif call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) ! increment the alpha/alpha part for single excitations - if (nkeys+4 * norb .ge. size(values)) then + if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) nkeys = 0 endif @@ -381,7 +381,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis c_2(l) = u_t(l,k_a) c_average += c_1(l) * c_2(l) * state_weights(l) enddo - if (nkeys+4 .ge. size(values)) then + if (nkeys+4 .ge. sze_buff) then call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) nkeys = 0 endif @@ -452,13 +452,13 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis enddo if(alpha_beta.or.spin_trace.or.beta_beta)then ! increment the alpha/beta part for single excitations - if (nkeys+2 * norb .ge. size(values)) then + if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) nkeys = 0 endif call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) ! increment the beta /beta part for single excitations - if (nkeys+4 * norb .ge. size(values)) then + if (nkeys+4 * elec_alpha_num .ge. sze_buff) then call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) nkeys = 0 endif @@ -484,7 +484,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis c_2(l) = u_t(l,k_a) c_average += c_1(l) * c_2(l) * state_weights(l) enddo - if (nkeys+4 .ge. size(values)) then + if (nkeys+4 .ge. sze_buff) then call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) nkeys = 0 endif