LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-13 17:43:55 +01:00
Code Issues Releases Wiki Activity

Bug in map_integrals.irp.d fixed

Browse Source
This commit is contained in:
Manu 2014-08-21 11:14:30 +02:00
parent 8a966ed732
commit 1229781220
7 changed files with 87 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/BiInts/map_integrals.irp.f
View File

@ -98,11 +98,12 @@ subroutine get_ao_bielec_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_
PROVIDE ao_bielec_integrals_in_map
thresh = ao_integrals_threshold
non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
non_zero_int = 0
do i=1,sze
!DIR$ FORCEINLINE

4
src/CISD_SC2_selected/H_apply.irp.f
View File

@ -6,5 +6,9 @@ from perturbation import perturbations
s = H_apply("PT2",SingleRef=True)
s.set_perturbation("epstein_nesbet_sc2_projected")
print s
s = H_apply("PT2_en_sc2",SingleRef=True)
s.set_perturbation("epstein_nesbet_sc2")
print s
END_SHELL

7
src/CISD_SC2_selected/cisd_sc2_selection.irp.f
View File

@ -11,9 +11,9 @@ program cisd_sc2_selected
allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st),E_old(N_st))
pt2 = 1.d0
perturbation = "epstein_nesbet_sc2"
perturbation = "epstein_nesbet_sc2_projected"
E_old(1) = HF_energy
davidson_threshold = 1.d-6
davidson_threshold = 1.d-8
if (N_det > n_det_max_cisd_sc2) then
call diagonalize_CI_SC2
call save_wavefunction
@ -80,7 +80,8 @@ program cisd_sc2_selected
integer :: imax
print *, 'PT2(SC2) = ', pt2(i)
print *, 'E(SC2) = ', CI_SC2_energy(i)
print *, 'E_before(SC2)+PT2(SC2) = ', (CI_SC2_energy(i)+pt2(i))
print *, 'E_before(SC2)+PT2(SC2) = ', CI_SC2_energy(i)+pt2(i)
print *, 'E_before(SC2)+PT2(SC2)_new = ', CI_SC2_energy(i)+pt2(i)*(1.d0+norm_pert)
print*,'greater coeficient of the state : ',dabs(psi_coef(imax,i))
call get_excitation_degree(ref_bitmask,psi_det(1,1,imax),degree,N_int)

1
src/Dets/diagonalize_CI_SC2.irp.f
View File

@ -33,6 +33,7 @@ END_PROVIDER
do j=1,N_states
do i=1,N_det
! CI_SC2_eigenvectors(i,j) = psi_coef(i,j)
CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j)
enddo
CI_SC2_electronic_energy(j) = CI_electronic_energy(j)

55
src/Dets/slater_rules.irp.f
View File

@ -571,6 +571,61 @@ subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx
end
subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat)
use bitmasks
BEGIN_DOC
! <key|H|psi> for the various Nstate
!
! returns in addition
!
! the array of the index of the non connected determinants to key1
!
! in order to know what double excitation can be repeated on key1
!
! idx_repeat(0) is the number of determinants that can be used
!
! to repeat the excitations
END_DOC
implicit none
integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
double precision, intent(in) :: coef(Ndet_max,Nstate)
double precision, intent(out) :: i_H_psi_array(Nstate)
integer , intent(out) :: idx_repeat(0:Ndet)
integer :: i, ii,j
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet)
ASSERT (Nint > 0)
ASSERT (N_int == Nint)
ASSERT (Nstate > 0)
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
i_H_psi_array = 0.d0
call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat)
print*,'--------'
do ii=1,idx(0)
print*,'--'
i = idx(ii)
!DEC$ FORCEINLINE
call i_H_j(keys(1,1,i),key,Nint,hij)
if (i==1)then
print*,'i==1 !!'
endif
print*,coef(i,1) * hij,coef(i,1),hij
do j = 1, Nstate
i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
enddo
print*,i_H_psi_array(1)
enddo
print*,'------'
end
subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)
use bitmasks

15
src/Perturbation/pert_sc2.irp.f
View File

@ -39,14 +39,14 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
double precision :: tmp
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
!$IVDEP
do i = 1, idx_repeat(0)
accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
enddo
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
@ -54,10 +54,6 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
do i =2,N_st
H_pert_diag(i) = h
! if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then
! c_pert(i) = -1.d0
! e_2_pert(i) = -2.d0
! else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
@ -68,17 +64,16 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
enddo
degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) + &
popcnt(xor( ref_bitmask(1,2), det_pert(1,2)))
popcnt(xor( ref_bitmask(1,2), det_pert(1,2)))
!DEC$ NOUNROLL
do l=2,Nint
degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) + &
popcnt(xor( ref_bitmask(l,2), det_pert(l,2)))
popcnt(xor( ref_bitmask(l,2), det_pert(l,2)))
enddo
if(degree==4)then
! <psi|delta_H|psi>
call i_H_j(ref_bitmask,det_pert,Nint,h0j)
H_pert_diag(1) = e_2_pert(1)
e_2_pert_fonda = H_pert_diag(1)
e_2_pert_fonda = H_pert_diag(1)
do i = 1, N_st
do j = 1, idx_repeat(0)
e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i)

18
src/Selectors_full/e_corr_selectors.irp.f
View File

@ -25,7 +25,13 @@ use bitmasks
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, coef_hf_selector]
&BEGIN_PROVIDER[double precision, inv_selectors_coef_hf]
&BEGIN_PROVIDER[double precision, inv_selectors_coef_hf_squared]
&BEGIN_PROVIDER[double precision, E_corr_per_selectors, (N_det_selectors)]
&BEGIN_PROVIDER[double precision, i_H_HF_per_selectors, (N_det_selectors)]
&BEGIN_PROVIDER[double precision, Delta_E_per_selector, (N_det_selectors)]
&BEGIN_PROVIDER[double precision, E_corr_double_only ]
&BEGIN_PROVIDER[double precision, E_corr_second_order ]
implicit none
BEGIN_DOC
! energy of correlation per determinant respect to the Hartree Fock determinant
@ -39,25 +45,33 @@ END_PROVIDER
! coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants
END_DOC
integer :: i,degree
double precision :: hij,inv_selectors_coef_hf
double precision :: hij,diag_H_mat_elem
E_corr_double_only = 0.d0
E_corr_second_order = 0.d0
do i = 1, N_det_selectors
if(exc_degree_per_selectors(i)==2)then
call i_H_j(ref_bitmask,psi_selectors(1,1,i),N_int,hij)
i_H_HF_per_selectors(i) = hij
E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij
E_corr_double_only += E_corr_per_selectors(i)
E_corr_second_order += hij * hij /(ref_bitmask_energy - diag_H_mat_elem(psi_selectors(1,1,i),N_int))
elseif(exc_degree_per_selectors(i) == 0)then
coef_hf_selector = psi_selectors_coef(i,1)
E_corr_per_selectors(i) = -1000.d0
Delta_E_per_selector(i) = 0.d0
else
E_corr_per_selectors(i) = -1000.d0
endif
enddo
if (dabs(coef_hf_selector) > 1.d-8) then
inv_selectors_coef_hf = 1.d0/coef_hf_selector
inv_selectors_coef_hf_squared = inv_selectors_coef_hf * inv_selectors_coef_hf
else
inv_selectors_coef_hf = 0.d0
inv_selectors_coef_hf_squared = 0.d0
endif
do i = 1,n_double_selectors
E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf
enddo
E_corr_double_only = E_corr_double_only * inv_selectors_coef_hf
END_PROVIDER