diff --git a/src/casscf/test_pert_2rdm.irp.f b/src/casscf/test_pert_2rdm.irp.f new file mode 100644 index 00000000..89d7bc8a --- /dev/null +++ b/src/casscf/test_pert_2rdm.irp.f @@ -0,0 +1,29 @@ +program test_pert_2rdm + implicit none + read_wf = .True. + touch read_wf + call get_pert_2rdm + integer :: i,j,k,l,ii,jj,kk,ll + double precision :: accu , get_two_e_integral, integral + accu = 0.d0 + print*,'n_orb_pert_rdm = ',n_orb_pert_rdm + do ii = 1, n_orb_pert_rdm + i = list_orb_pert_rdm(ii) + do jj = 1, n_orb_pert_rdm + j = list_orb_pert_rdm(jj) + do kk = 1, n_orb_pert_rdm + k= list_orb_pert_rdm(kk) + do ll = 1, n_orb_pert_rdm + l = list_orb_pert_rdm(ll) + 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 +! endif + accu += pert_2rdm_provider(ii,jj,kk,ll) * integral + enddo + enddo + enddo + enddo + print*,'accu = ',accu +end diff --git a/src/cipsi/lock_2rdm.irp.f b/src/cipsi/lock_2rdm.irp.f new file mode 100644 index 00000000..e69de29b diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f index 9cf8fba7..85bea747 100644 --- a/src/cipsi/pert_rdm_providers.irp.f +++ b/src/cipsi/pert_rdm_providers.irp.f @@ -1,5 +1,12 @@ use bitmasks +use omp_lib + +BEGIN_PROVIDER [ integer(omp_lock_kind), pert_2rdm_lock] + use f77_zmq + implicit none + call omp_init_lock(pert_2rdm_lock) +END_PROVIDER BEGIN_PROVIDER [logical , pert_2rdm ] implicit none @@ -29,13 +36,13 @@ BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert 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) +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_reverse, 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) + double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection) 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) @@ -136,7 +143,9 @@ 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 + call get_excitation_degree(det,HF_bitmask,degree,N_int) + if(degree == 2)cycle do istate=1,N_states delta_E = E0(istate) - Hii + E_shift alpha_h_psi = mat(istate, p1, p2) @@ -147,6 +156,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 pt2(istate) = pt2(istate) + e_pert variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi norm(istate) = norm(istate) + coef(istate) * coef(istate) @@ -154,19 +164,20 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo 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) + call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,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 + call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) end diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 71442538..248876ef 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -249,7 +249,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d integer,allocatable :: tmp_array(:) integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) logical, allocatable :: banned(:,:,:), bannedOrb(:,:) - double precision, allocatable :: coef_fullminilist(:,:) + double precision, allocatable :: coef_fullminilist_rev(:,:) double precision, allocatable :: mat(:,:,:) @@ -549,9 +549,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d allocate (fullminilist (N_int, 2, fullinteresting(0)), & minilist (N_int, 2, interesting(0)) ) if(pert_2rdm)then - allocate(coef_fullminilist(fullinteresting(0),N_states)) + allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) do i=1,fullinteresting(0) - coef_fullminilist(i,:) = psi_coef_sorted(fullinteresting(i),:) + do j = 1, N_states + coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) + enddo enddo endif do i=1,fullinteresting(0) @@ -608,7 +610,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d if(.not.pert_2rdm)then call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf) else - call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf,fullminilist, coef_fullminilist, fullinteresting(0)) + call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) endif end if enddo @@ -616,7 +618,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d enddo deallocate(fullminilist,minilist) if(pert_2rdm)then - deallocate(coef_fullminilist) + deallocate(coef_fullminilist_rev) endif enddo enddo diff --git a/src/cipsi/update_2rdm.irp.f b/src/cipsi/update_2rdm.irp.f index 7ae42ea8..260c48fd 100644 --- a/src/cipsi/update_2rdm.irp.f +++ b/src/cipsi/update_2rdm.irp.f @@ -1,12 +1,12 @@ use bitmasks -subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection,n_det_connection,nkeys,keys,values,sze_buff) +subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff) implicit none integer, intent(in) :: n_det_connection,sze_buff double precision, intent(in) :: coef(N_states) integer(bit_kind), intent(in) :: det(N_int,2) integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) - double precision, intent(in) :: psi_coef_connection(n_det_connection, N_states) + double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection) integer, intent(inout) :: keys(4,sze_buff),nkeys double precision, intent(inout) :: values(sze_buff) integer :: i,j @@ -18,24 +18,26 @@ 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(i,j) * phase * coef(j) + contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i) * 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) + if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then + call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) 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) + !! case of double excitations + ! if (nkeys + 4 .ge. sze_buff)then + ! call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) + ! nkeys = 0 + ! endif + ! 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 end @@ -49,7 +51,81 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v integer, intent(inout) :: nkeys, keys(4,sze_buff) double precision, intent(inout):: values(sze_buff) + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2),ispin,other_spin + integer :: h1,h2,p1,p2,i + call bitstring_to_list_ab(det1, occ, n_occ_ab, N_int) + + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + p1 = exc(1,2,1) + ispin = 1 + other_spin = 2 + else + ! Mono beta + h1 = exc(1,1,2) + p1 = exc(1,2,2) + ispin = 2 + other_spin = 1 + endif + 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 + p1 = list_orb_reverse_pert_rdm(p1) + !update the alpha/beta part + do i = 1, n_occ_ab(other_spin) + h2 = occ(i,other_spin) + 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) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + 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) + +! 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) = h2 +! keys(2,nkeys) = h1 +! keys(3,nkeys) = p1 +! keys(4,nkeys) = h2 +!enddo end @@ -60,21 +136,88 @@ subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_ double precision,intent(in) :: phase, contrib integer, intent(inout) :: nkeys, keys(4,sze_buff) double precision, intent(inout):: values(sze_buff) + integer :: h1,h2,p1,p2 + if (exc(0,1,1) == 1) then + ! Double alpha/beta + h1 = exc(1,1,1) + h2 = exc(1,1,2) + p1 = exc(1,2,1) + p2 = exc(1,2,2) + ! check if the orbitals involved are within the orbital range + if(list_orb_reverse_pert_rdm(h1).lt.0)return + h1 = list_orb_reverse_pert_rdm(h1) + if(list_orb_reverse_pert_rdm(h2).lt.0)return + h2 = list_orb_reverse_pert_rdm(h2) + if(list_orb_reverse_pert_rdm(p1).lt.0)return + p1 = list_orb_reverse_pert_rdm(p1) + if(list_orb_reverse_pert_rdm(p2).lt.0)return + p2 = list_orb_reverse_pert_rdm(p2) + nkeys += 1 + values(nkeys) = 0.5d0 * contrib * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + nkeys += 1 + values(nkeys) = 0.5d0 * contrib * phase + keys(1,nkeys) = p1 + keys(2,nkeys) = p2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + else + if (exc(0,1,1) == 2) then + ! Double alpha/alpha + h1 = exc(1,1,1) + h2 = exc(2,1,1) + p1 = exc(1,2,1) + p2 = exc(2,2,1) + else if (exc(0,1,2) == 2) then + ! Double beta + h1 = exc(1,1,2) + h2 = exc(2,1,2) + p1 = exc(1,2,2) + p2 = exc(2,2,2) + endif + ! check if the orbitals involved are within the orbital range + if(list_orb_reverse_pert_rdm(h1).lt.0)return + h1 = list_orb_reverse_pert_rdm(h1) + if(list_orb_reverse_pert_rdm(h2).lt.0)return + h2 = list_orb_reverse_pert_rdm(h2) + if(list_orb_reverse_pert_rdm(p1).lt.0)return + p1 = list_orb_reverse_pert_rdm(p1) + if(list_orb_reverse_pert_rdm(p2).lt.0)return + p2 = list_orb_reverse_pert_rdm(p2) + nkeys += 1 + values(nkeys) = 0.5d0 * contrib * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + values(nkeys) = - 0.5d0 * contrib * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * contrib * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + 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) = p2 + endif + end -subroutine update_rdms(nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: nkeys, keys(4,sze_buff),sze_buff - double precision, intent(in) :: values(sze_buff) - integer :: i,h1,h2,p1,p2 - do i = 1, nkeys - h1 = keys(1,i) - h2 = keys(2,i) - p1 = keys(3,i) - p2 = keys(4,i) - pert_2rdm_provider(h1,h2,p1,p2) += values(i) - enddo -end