From 893ee6b65488a28fae8c4111ff871b8d0ccef3c1 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner LCT Date: Sat, 26 Oct 2019 17:05:16 +0200 Subject: [PATCH] pert 2rdm work for the Li atom with HF --- src/casscf/test_pert_2rdm.irp.f | 2 +- src/cipsi/pert_rdm_providers.irp.f | 10 ++-- src/cipsi/update_2rdm.irp.f | 74 ++++++++++++++++-------------- 3 files changed, 48 insertions(+), 38 deletions(-) diff --git a/src/casscf/test_pert_2rdm.irp.f b/src/casscf/test_pert_2rdm.irp.f index 890e8cbe..4e049f7d 100644 --- a/src/casscf/test_pert_2rdm.irp.f +++ b/src/casscf/test_pert_2rdm.irp.f @@ -24,7 +24,7 @@ program test_pert_2rdm integral = get_two_e_integral(i,j,k,l,mo_integrals_map) ! if(dabs(pert_2rdm_provider(ii,jj,kk,ll) * integral).gt.1.d-12)then ! print*,i,j,k,l -! print*,pert_2rdm_provider(ii,jj,kk,ll) * integral,pert_2rdm_provider(ii,jj,kk,ll), pert_2rdm_provider(ii,jj,kk,ll), integral +! print*,pert_2rdm_provider(ii,jj,kk,ll) , integral, pert_2rdm_provider(ii,jj,kk,ll)* integral ! endif accu += pert_2rdm_provider(ii,jj,kk,ll) * integral enddo diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f index 5a18e683..83265fba 100644 --- a/src/cipsi/pert_rdm_providers.irp.f +++ b/src/cipsi/pert_rdm_providers.irp.f @@ -174,9 +174,13 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) sum_e_pert = 0d0 - integer :: degree + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase + call get_excitation(det,HF_bitmask,exc,degree,phase,N_int) call get_excitation_degree(det,HF_bitmask,degree,N_int) - if(degree == 2)cycle + if(degree == 1)cycle +! if (exc(0,1,1) .ne. 1) cycle !only double alpha/beta do istate=1,N_states delta_E = E0(istate) - Hii + E_shift alpha_h_psi = mat(istate, p1, p2) @@ -187,7 +191,7 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo endif e_pert = 0.5d0 * (tmp - delta_E) coef(istate) = e_pert / alpha_h_psi - print*,e_pert,coef,alpha_h_psi +! print*,coef,alpha_h_psi,e_pert pt2(istate) = pt2(istate) + e_pert variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi norm(istate) = norm(istate) + coef(istate) * coef(istate) diff --git a/src/cipsi/update_2rdm.irp.f b/src/cipsi/update_2rdm.irp.f index 900be8e7..5dd832fd 100644 --- a/src/cipsi/update_2rdm.irp.f +++ b/src/cipsi/update_2rdm.irp.f @@ -18,7 +18,7 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio if(degree.gt.2)cycle contrib = 0.d0 do j = 1, N_states - contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i) * phase * coef(j) + contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i) * coef(j) enddo ! case of single excitations if(degree == 1)then @@ -36,8 +36,8 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) endif enddo -!call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) -!nkeys = 0 + call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) + nkeys = 0 end @@ -69,7 +69,7 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v ispin = 2 other_spin = 1 endif - print*,'h1,p1,s1',h1,p1,ispin +!print*,'h1,p1,s1',h1,p1,ispin if(list_orb_reverse_pert_rdm(h1).lt.0)return h1 = list_orb_reverse_pert_rdm(h1) if(list_orb_reverse_pert_rdm(p1).lt.0)return @@ -94,39 +94,41 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v keys(4,nkeys) = p1 enddo !update the same spin part -!do i = 1, n_occ_ab(ispin) -! h2 = occ(i,ispin) -! if(list_orb_reverse_pert_rdm(h2).lt.0)return -! h2 = list_orb_reverse_pert_rdm(h2) + do i = 1, n_occ_ab(ispin) + h2 = occ(i,ispin) + if(list_orb_reverse_pert_rdm(h2).lt.0)return + h2 = list_orb_reverse_pert_rdm(h2) -! nkeys += 1 -! values(nkeys) = 0.5d0 * contrib * phase -! keys(1,nkeys) = h1 -! keys(2,nkeys) = h2 -! keys(3,nkeys) = p1 -! keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = 0.5d0 * contrib * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 -! nkeys += 1 -! values(nkeys) = - 0.5d0 * contrib * phase -! keys(1,nkeys) = h1 -! keys(2,nkeys) = h2 -! keys(3,nkeys) = h2 -! keys(4,nkeys) = p1 -! -! nkeys += 1 -! values(nkeys) = 0.5d0 * contrib * phase -! keys(1,nkeys) = h2 -! keys(2,nkeys) = h1 -! keys(3,nkeys) = h2 -! keys(4,nkeys) = p1 + nkeys += 1 + values(nkeys) = - 0.5d0 * contrib * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * contrib * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 -! nkeys += 1 -! values(nkeys) = - 0.5d0 * contrib * phase -! keys(1,nkeys) = h2 -! keys(2,nkeys) = h1 -! keys(3,nkeys) = p1 -! keys(4,nkeys) = h2 -!enddo + nkeys += 1 + values(nkeys) = - 0.5d0 * contrib * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) + nkeys = 0 end @@ -160,12 +162,14 @@ subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_ keys(2,nkeys) = h2 keys(3,nkeys) = p1 keys(4,nkeys) = p2 +! print*,contrib nkeys += 1 values(nkeys) = 0.5d0 * contrib * phase keys(1,nkeys) = p1 keys(2,nkeys) = p2 keys(3,nkeys) = h1 keys(4,nkeys) = h2 +! print*,contrib else if (exc(0,1,1) == 2) then @@ -218,6 +222,8 @@ subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_ keys(3,nkeys) = p1 keys(4,nkeys) = p2 endif + call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) + nkeys = 0 end