LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-09-27 03:51:01 +02:00
Code Issues Releases Wiki Activity

Introduced SOP-PT

Browse Source
This commit is contained in:
Anthony Scemama 2018-11-21 22:54:05 +01:00
parent 63aa16cc53
commit a77093577e
8 changed files with 113 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
ocaml/qptypes_generator.ml
View File

@ -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

16
src/Determinants/determinants.irp.f
View File

@ -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
! <i|h|i> 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

81
src/Determinants/occ_pattern.irp.f
View File

@ -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
! <I|h|I> where |I> is an occupation pattern. This is the minimum Hii of
! all <i|h|i>, 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

9
src/FCI/PT2.irp.f
View File

@ -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

2
src/FCI/energy.irp.f
View File

@ -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

16
src/FCI/selection.irp.f
View File

@ -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

5
src/Iterations/print_summary.irp.f
View File

@ -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 *, '-----'

2
src/Perturbation/EZFIO.cfg
View File

@ -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