From a77093577efe0473b2b184e8b31f855834a1081c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 21 Nov 2018 22:54:05 +0100 Subject: [PATCH] Introduced SOP-PT --- ocaml/qptypes_generator.ml | 3 ++ src/Determinants/determinants.irp.f | 16 ++++++ src/Determinants/occ_pattern.irp.f | 81 ++++++++++++++++++++++++++--- src/FCI/PT2.irp.f | 9 ---- src/FCI/energy.irp.f | 2 + src/FCI/selection.irp.f | 16 ++++-- src/Iterations/print_summary.irp.f | 5 +- src/Perturbation/EZFIO.cfg | 2 +- 8 files changed, 113 insertions(+), 21 deletions(-) diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index 48eb6085..06aa0ab6 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -224,14 +224,17 @@ end = struct | EN | Barycentric | Variance + | SOP [@@deriving sexp] let to_string = function | EN -> \"EN\" | Variance -> \"Variance\" | Barycentric -> \"Barycentric\" + | SOP -> \"SOP\" let of_string s = match (String.lowercase_ascii s) with + | \"sop\" -> SOP | \"en\" -> EN | \"variance\" -> Variance | \"barycentric\" -> Barycentric diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index e08bc793..af2da4b0 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -884,3 +884,19 @@ subroutine apply_hole(det, s1, h1, res, ok, Nint) ok = .true. end subroutine + + + +BEGIN_PROVIDER [ double precision, psi_det_Hii, (N_det) ] + implicit none + BEGIN_DOC + ! for all determinants. + END_DOC + integer :: i,j + double precision, external :: diag_H_mat_elem + do i=1,N_det + psi_det_Hii(i) = diag_H_mat_elem(psi_det(1,1,i),N_int) + enddo +END_PROVIDER + + diff --git a/src/Determinants/occ_pattern.irp.f b/src/Determinants/occ_pattern.irp.f index d6d2668b..c3cf9c47 100644 --- a/src/Determinants/occ_pattern.irp.f +++ b/src/Determinants/occ_pattern.irp.f @@ -1,5 +1,5 @@ use bitmasks -subroutine det_to_occ_pattern(d,o,Nint) +subroutine occ_pattern_of_det(d,o,Nint) use bitmasks implicit none BEGIN_DOC @@ -181,6 +181,7 @@ end ! array of the occ_pattern present in the wf ! psi_occ_pattern(:,1,j) = jth occ_pattern of the wave function : represent all the single occupations ! psi_occ_pattern(:,2,j) = jth occ_pattern of the wave function : represent all the double occupations + ! The occ patterns are sorted by occ_pattern_search_key END_DOC integer :: i,j,k @@ -304,9 +305,18 @@ BEGIN_PROVIDER [ integer, det_to_occ_pattern, (N_det) ] BEGIN_DOC ! Returns the index of the occupation pattern for each determinant END_DOC - integer :: i,j,k + integer :: i,j,k,r,l + integer*8 :: key integer(bit_kind) :: occ(N_int,2) logical :: found + integer*8, allocatable :: bit_tmp(:) + integer*8, external :: occ_pattern_search_key + + allocate(bit_tmp(N_occ_pattern)) + do i=1,N_occ_pattern + bit_tmp(i) = occ_pattern_search_key(psi_occ_pattern(1,1,i),N_int) + enddo + !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE(i,k,j,found,occ) do i=1,N_det @@ -314,24 +324,83 @@ BEGIN_PROVIDER [ integer, det_to_occ_pattern, (N_det) ] occ(k,1) = ieor(psi_det(k,1,i),psi_det(k,2,i)) occ(k,2) = iand(psi_det(k,1,i),psi_det(k,2,i)) enddo - do j=1,N_occ_pattern + + key = occ_pattern_search_key(psi_occ_pattern(1,1,i),N_int) + + ! Binary search + l = 1 + r = N_occ_pattern + do while(r-l > 1) + j = shiftr(r+l,1) + if (bit_tmp(j) < key) then + l = j + else + r = j + endif + enddo + + do while (bit_tmp(l) == key) + if (l == 0) then + print *, '1 bug in ', irp_here + stop -1 + endif found = .True. do k=1,N_int - if ( (occ(k,1) /= psi_occ_pattern(k,1,j)) & - .or. (occ(k,2) /= psi_occ_pattern(k,2,j)) ) then + if ( (occ(k,1) /= psi_occ_pattern(k,1,l)) & + .or. (occ(k,2) /= psi_occ_pattern(k,2,l)) ) then found = .False. exit endif enddo if (found) then - det_to_occ_pattern(i) = j + det_to_occ_pattern(i) = l + r=1 exit endif + l = l-1 + enddo + + do while (bit_tmp(r) == key) + if (r > N_occ_pattern) then + print *, '2 bug in ', irp_here + stop -1 + endif + found = .True. + do k=1,N_int + if ( (occ(k,1) /= psi_occ_pattern(k,1,r)) & + .or. (occ(k,2) /= psi_occ_pattern(k,2,r)) ) then + found = .False. + exit + endif + enddo + if (found) then + det_to_occ_pattern(i) = r + exit + endif + r = r+1 enddo enddo !$OMP END PARALLEL DO END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_occ_pattern_Hii, (N_occ_pattern) ] + implicit none + BEGIN_DOC + ! where |I> is an occupation pattern. This is the minimum Hii of + ! all , where the |i> are the determinants if oI> + END_DOC + integer :: j, i + + psi_occ_pattern_Hii(:) = huge(1.d0) + do i=1,N_det + j = det_to_occ_pattern(i) + psi_occ_pattern_Hii(j) = min(psi_occ_pattern_Hii(j), psi_det_Hii(i)) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, weight_occ_pattern, (N_occ_pattern,N_states) ] implicit none BEGIN_DOC diff --git a/src/FCI/PT2.irp.f b/src/FCI/PT2.irp.f index 4795e0e2..c0c2e01d 100644 --- a/src/FCI/PT2.irp.f +++ b/src/FCI/PT2.irp.f @@ -33,15 +33,6 @@ subroutine run call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm) - print *, 'State ', k - do k=1,N_states - print *, 'N_det = ', N_det - print *, 'PT2 = ', pt2(k) - print *, 'rPT2 = ', rpt2(k) - print *, 'E = ', E_CI_before(k) - print *, 'E+PT2 = ', E_CI_before(k)+pt2(k), ' +/- ', error(k) - print *, '-----' - enddo call ezfio_set_fci_energy(E_CI_before) call ezfio_set_fci_energy_pt2(E_CI_before+pt2) end diff --git a/src/FCI/energy.irp.f b/src/FCI/energy.irp.f index 46117256..b7ba42bb 100644 --- a/src/FCI/energy.irp.f +++ b/src/FCI/energy.irp.f @@ -18,6 +18,8 @@ BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ] pt2_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states) else if (h0_type == "Variance") then pt2_E0_denominator(1:N_states) = psi_energy(1:N_states) !1.d0-nuclear_repulsion + else if (h0_type == "SOP") then + pt2_E0_denominator(1:N_states) = psi_energy(1:N_states) else print *, h0_type, ' not implemented' stop diff --git a/src/FCI/selection.irp.f b/src/FCI/selection.irp.f index e79f2059..4bffdb39 100644 --- a/src/FCI/selection.irp.f +++ b/src/FCI/selection.irp.f @@ -658,6 +658,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d 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 double precision, external :: diag_H_mat_elem_fock + double precision :: E_shift logical, external :: detEq @@ -670,7 +671,13 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d 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_tot_num if(bannedOrb(p1, s1)) cycle ib = 1 @@ -683,10 +690,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) 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 + 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) @@ -699,10 +707,10 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi norm(istate) = norm(istate) + coef * coef - if (h0_type /= "Variance") then - sum_e_pert = sum_e_pert + e_pert * state_average_weight(istate) - else + if (h0_type == "Variance") then sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * state_average_weight(istate) + else + sum_e_pert = sum_e_pert + e_pert * state_average_weight(istate) endif end do diff --git a/src/Iterations/print_summary.irp.f b/src/Iterations/print_summary.irp.f index f785fb2a..42a27786 100644 --- a/src/Iterations/print_summary.irp.f +++ b/src/Iterations/print_summary.irp.f @@ -60,15 +60,18 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_) if (s2_eig) then print *, 'N_sop = ', N_occ_pattern endif + print *, '' do k=1, N_states_p - print*,'State ',k + print*,'* State ',k print *, 'Variance = ', variance_(k) print *, 'PT norm = ', dsqrt(norm_(k)) print *, 'PT2 = ', pt2_(k) + print *, 'rPT2 = ', pt2_(k)*f(k) print *, 'E = ', e_(k) print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_(k), ' +/- ', error_(k) print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_(k)*f(k), ' +/- ', error_(k)*f(k) + print *, '' enddo print *, '-----' diff --git a/src/Perturbation/EZFIO.cfg b/src/Perturbation/EZFIO.cfg index e4c40d4a..f92b1f40 100644 --- a/src/Perturbation/EZFIO.cfg +++ b/src/Perturbation/EZFIO.cfg @@ -26,7 +26,7 @@ default: 1.00 [h0_type] type: Perturbation -doc: Type of zeroth-order Hamiltonian [ EN | Barycentric | Variance ] +doc: Type of zeroth-order Hamiltonian [ EN | Barycentric | Variance | SOP ] interface: ezfio,provider,ocaml default: EN