diff --git a/src/cipsi/EZFIO.cfg b/src/cipsi/EZFIO.cfg index 19b45ac1..7fcf19eb 100644 --- a/src/cipsi/EZFIO.cfg +++ b/src/cipsi/EZFIO.cfg @@ -1,9 +1,3 @@ -[pert_2rdm] -type: logical -doc: If true, computes the one- and two-body rdms with perturbation theory -interface: ezfio,provider,ocaml -default: False - [save_wf_after_selection] type: logical doc: If true, saves the wave function after the selection, before the diagonalization diff --git a/src/cipsi/NEED b/src/cipsi/NEED index bfbc559a..85d01f79 100644 --- a/src/cipsi/NEED +++ b/src/cipsi/NEED @@ -2,5 +2,4 @@ perturbation zmq mpi iterations -two_body_rdm csf diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f deleted file mode 100644 index eca8decc..00000000 --- a/src/cipsi/pert_rdm_providers.irp.f +++ /dev/null @@ -1,183 +0,0 @@ - -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 [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 - -BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)] - implicit none - pert_2rdm_provider = 0.d0 - -END_PROVIDER - -subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, 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_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) - 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) - type(pt2_type), intent(inout) :: pt2_data - type(selection_buffer), intent(inout) :: buf - logical :: ok - integer :: s1, s2, p1, p2, ib, j, istate, jstate - 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 == 'CFG') then - j = det_to_configuration(i_generator) - E_shift = psi_det_Hii(i_generator) - psi_configuration_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 - 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) - 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 - print*,e_pert,coef,alpha_h_psi - pt2_data % pt2(istate) += e_pert - pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi - enddo - - do istate=1,N_states - alpha_h_psi = mat(istate, p1, p2) - e_pert = coef(istate) * alpha_h_psi - do jstate=1,N_states - pt2_data % overlap(jstate,jstate) = coef(istate) * coef(jstate) - enddo - - 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_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/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index b366a268..3594aaf2 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -133,7 +133,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted PROVIDE psi_det_hii selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max if (h0_type == 'CFG') then PROVIDE psi_configuration_hii det_to_configuration diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index eda9642c..f1ec6ff6 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -464,14 +464,14 @@ 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_rev(N_states,fullinteresting(0))) - do i=1,fullinteresting(0) - do j = 1, N_states - coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) - enddo - enddo - endif +! if(pert_2rdm)then +! allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) +! do i=1,fullinteresting(0) +! do j = 1, N_states +! coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) +! enddo +! enddo +! endif do i=1,fullinteresting(0) do k=1,N_int @@ -531,19 +531,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) - if(.not.pert_2rdm)then +! if(.not.pert_2rdm)then call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) - else - call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) - endif +! else +! call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) +! endif end if enddo if(s1 /= s2) monoBdo = .false. enddo deallocate(fullminilist,minilist) - if(pert_2rdm)then - deallocate(coef_fullminilist_rev) - endif +! if(pert_2rdm)then +! deallocate(coef_fullminilist_rev) +! endif enddo enddo deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) diff --git a/src/cipsi/update_2rdm.irp.f b/src/cipsi/update_2rdm.irp.f deleted file mode 100644 index 260c48fd..00000000 --- a/src/cipsi/update_2rdm.irp.f +++ /dev/null @@ -1,223 +0,0 @@ -use bitmasks - -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_reverse(N_states,n_det_connection) - integer, intent(inout) :: keys(4,sze_buff),nkeys - double precision, intent(inout) :: values(sze_buff) - integer :: i,j - 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_reverse(j,i) * phase * coef(j) - enddo - ! case of single excitations - if(degree == 1)then - 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_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 - -subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: 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) - - 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 - -subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: 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) - 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 - - diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 58630709..1bfe87c0 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -22,7 +22,7 @@ subroutine ZMQ_selection(N_in, pt2_data) PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection') diff --git a/src/davidson_dressed/diagonalize_ci.irp.f b/src/davidson_dressed/diagonalize_ci.irp.f index 7619532a..709ee0e6 100644 --- a/src/davidson_dressed/diagonalize_ci.irp.f +++ b/src/davidson_dressed/diagonalize_ci.irp.f @@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] enddo do j=1,min(N_det,N_states) write(st,'(I4)') j - call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st)) + call write_double(6,CI_energy_dressed(j),'Energy dressed of state '//trim(st)) call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st)) enddo diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f index a1d7a5c1..21f6dc21 100644 --- a/src/dft_one_e/mu_erf_dft.irp.f +++ b/src/dft_one_e/mu_erf_dft.irp.f @@ -58,6 +58,7 @@ BEGIN_PROVIDER [double precision, mu_of_r_dft_average] r(3) = final_grid_points(3,i) call dm_dft_alpha_beta_at_r(r,dm_a,dm_b) rho = dm_a + dm_b + if(mu_of_r_dft(i).gt.1.d+3)cycle mu_of_r_dft_average += rho * mu_of_r_dft(i) * final_weight_at_r_vector(i) enddo mu_of_r_dft_average = mu_of_r_dft_average / dble(elec_alpha_num + elec_beta_num) diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 641ee209..8e6285e2 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -98,7 +98,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s enddo endif - call print_energy_components() +! call print_energy_components() end subroutine diff --git a/src/tools/NEED b/src/tools/NEED index c07c9109..0f4e17b0 100644 --- a/src/tools/NEED +++ b/src/tools/NEED @@ -2,3 +2,4 @@ fci mo_two_e_erf_ints aux_quantities hartree_fock +two_body_rdm diff --git a/src/davidson/print_e_components.irp.f b/src/two_body_rdm/print_e_components.irp.f similarity index 100% rename from src/davidson/print_e_components.irp.f rename to src/two_body_rdm/print_e_components.irp.f