LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-25 22:52:15 +02:00
Code Issues Releases Wiki Activity

Dangerous commit

Browse Source
This commit is contained in:
Anthony Scemama 2014-05-21 16:37:54 +02:00
parent 0ddfe658cc
commit a345c0ce31
12 changed files with 362 additions and 320 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
scripts/check_dependencies.sh
View File

@ -51,7 +51,7 @@ DEPS=$(unique_list $DEPS_LONG)
if [[ ! "$COMMAND_LINE" == "$DEPS" ]]
then
DEPS=$(check_dependencies.sh $DEPS)
DEPS=$(${QPACKAGE_ROOT}/scripts/check_dependencies.sh ${DEPS})
fi
echo $DEPS

38
scripts/generate_h_apply.py
View File

@ -10,7 +10,8 @@ subroutine
parameters
initialization
declarations
keys_work
keys_work_locked
keys_work_unlocked
finalization
""".split()
@ -24,7 +25,10 @@ class H_apply(object):
self.openmp = openmp
if openmp:
s["subroutine"] += "_OpenMP"
self.selection_pt2 = None
self.perturbation = None
#s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) &
s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, &
@ -34,7 +38,7 @@ class H_apply(object):
!$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,&
!$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, &
!$OMP N_elec_in_key_hole_2,ia_ja_pairs) &
!$OMP SHARED(key_in,N_int,elec_num_tab, &
!$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, &
!$OMP hole_1, particl_1, hole_2, particl_2, &
!$OMP lck,thresh,elec_alpha_num)"""
s["omp_init_lock"] = "call omp_init_lock(lck)"
@ -72,6 +76,8 @@ class H_apply(object):
return buffer
def set_perturbation(self,pert):
if self.perturbation is not None:
raise
self.perturbation = pert
if pert is not None:
self.data["parameters"] = ",sum_e_2_pert_in,sum_norm_pert_in,sum_H_pert_diag_in,N_st,Nint"
@ -83,16 +89,17 @@ class H_apply(object):
double precision :: sum_e_2_pert(N_st)
double precision :: sum_norm_pert(N_st)
double precision :: sum_H_pert_diag
double precision :: e_2_pert_buffer(N_st,size_max)
double precision :: coef_pert_buffer(N_st,size_max)
"""
self.data["size_max"] = "256"
self.data["initialization"] = """
E_ref = diag_H_mat_elem(key_in,N_int)
sum_e_2_pert = sum_e_2_pert_in
sum_norm_pert = sum_norm_pert_in
sum_H_pert_diag = sum_H_pert_diag_in
"""
self.data["keys_work"] += """
call perturb_buffer_%s(keys_out,key_idx,sum_e_2_pert, &
self.data["keys_work_unlocked"] += """
call perturb_buffer_%s(keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,Nint)
"""%(pert,)
self.data["finalization"] = """
@ -101,10 +108,29 @@ class H_apply(object):
sum_H_pert_diag_in = sum_H_pert_diag
"""
if self.openmp:
self.data["omp_set_lock"] = ""
self.data["omp_unset_lock"] = ""
self.data["omp_test_lock"] = ".False."
self.data["omp_parallel"] += """&
!$OMP SHARED(N_st) &
!$OMP SHARED(N_st,Nint) PRIVATE(e_2_pert_buffer,coef_pert_buffer) &
!$OMP REDUCTION(+:sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)"""
def set_selection_pt2(self,pert):
if self.selection_pt2 is not None:
raise
self.set_perturbation(pert)
self.selection_pt2 = pert
if pert is not None:
self.data["size_max"] = str(1024*128)
self.data["keys_work_unlocked"] = """
e_2_pert_buffer = 0.d0
coef_pert_buffer = 0.d0
""" + self.data["keys_work_unlocked"]
self.data["keys_work_locked"] = """
call fill_H_apply_buffer_selection(key_idx,keys_out,e_2_pert_buffer,coef_pert_buffer,N_st,N_int)
"""
self.data["omp_test_lock"] = "omp_test_lock(lck)"
self.data["omp_set_lock"] = "call omp_set_lock(lck)"
self.data["omp_unset_lock"] = "call omp_unset_lock(lck)"

22
src/Dets/H_apply.irp.f
View File

@ -33,11 +33,12 @@ subroutine resize_H_apply_buffer_det(new_size)
integer, intent(in) :: new_size
integer(bit_kind), allocatable :: buffer_det(:,:,:)
double precision, allocatable :: buffer_coef(:,:)
double precision, allocatable :: buffer_e2(:,:)
integer :: i,j,k
integer :: Ndet
ASSERT (new_size > 0)
allocate ( buffer_det(N_int,2,new_size), buffer_coef(new_size,N_states) )
allocate ( buffer_det(N_int,2,new_size), buffer_coef(new_size,N_states), buffer_e2(new_size,N_states) )
do i=1,min(new_size,H_apply_buffer_N_det)
do k=1,N_int
@ -48,9 +49,10 @@ subroutine resize_H_apply_buffer_det(new_size)
ASSERT (sum(popcnt(H_apply_buffer_det(:,2,i))) == elec_beta_num )
enddo
do k=1,N_states
do i=1,min(new_size,H_apply_buffer_N_det)
buffer_coef(i,k) = H_apply_buffer_coef(i,k)
enddo
do i=1,min(new_size,H_apply_buffer_N_det)
buffer_coef(i,k) = H_apply_buffer_coef(i,k)
buffer_e2(i,k) = H_apply_buffer_e2(i,k)
enddo
enddo
H_apply_buffer_size = new_size
@ -70,20 +72,23 @@ subroutine resize_H_apply_buffer_det(new_size)
do k=1,N_states
do i=1,H_apply_buffer_N_det
H_apply_buffer_coef(i,k) = buffer_coef(i,k)
H_apply_buffer_e2(i,k) = buffer_e2(i,k)
enddo
enddo
deallocate (buffer_det, buffer_coef)
SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det
deallocate (buffer_det, buffer_coef, buffer_e2)
SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det H_apply_buffer_e2
end
BEGIN_PROVIDER [ integer(bit_kind), H_apply_buffer_det,(N_int,2,H_apply_buffer_size) ]
&BEGIN_PROVIDER [ double precision, H_apply_buffer_coef,(H_apply_buffer_size,N_states) ]
&BEGIN_PROVIDER [ double precision, H_apply_buffer_e2,(H_apply_buffer_size,N_states) ]
&BEGIN_PROVIDER [ integer, H_apply_buffer_N_det ]
implicit none
BEGIN_DOC
! Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines.
! Buffer of determinants/coefficients/perturbative energy for H_apply.
! Uninitialized. Filled by H_apply subroutines.
END_DOC
H_apply_buffer_N_det = 0
@ -148,8 +153,9 @@ subroutine copy_H_apply_buffer_to_wf
psi_coef(i+N_det_old,k) = H_apply_buffer_coef(i,k)
enddo
enddo
H_apply_buffer_N_det = 0
SOFT_TOUCH psi_det psi_coef
SOFT_TOUCH psi_det psi_coef H_apply_buffer_N_det H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_e2
end

66
src/Dets/H_apply_template.f
View File

@ -7,6 +7,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
! particles.
! Assume N_int is already provided.
END_DOC
integer,parameter :: size_max = $size_max
$declarations
integer(omp_lock_kind) :: lck
integer(bit_kind),intent(in) :: key_in(N_int,2)
@ -26,13 +27,12 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
integer,parameter :: size_max = $size_max
double precision :: hij_elec, mo_bielec_integral, thresh
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem
PROVIDE mo_integrals_map ref_bitmask_energy
PROVIDE mo_bielec_integrals_in_map
PROVIDE mo_integrals_map ref_bitmask_energy key_pattern_not_in_ref
PROVIDE mo_bielec_integrals_in_map reference_energy psi_ref_coef psi_ref
$set_i_H_j_threshold
@ -156,8 +156,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$keys_work_unlocked
$omp_set_lock
$keys_work
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
@ -165,7 +166,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
enddo
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work
$keys_work_unlocked
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
@ -204,8 +206,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$keys_work_unlocked
$omp_set_lock
$keys_work
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
@ -213,7 +216,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
enddo
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work
$keys_work_locked
$keys_work_unlocked
$omp_unset_lock
key_idx = 0
endif
@ -222,8 +226,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
enddo ! ii
$omp_enddo
enddo ! ispin
$keys_work_unlocked
$omp_set_lock
$keys_work
$keys_work_locked
$omp_unset_lock
deallocate (keys_out,hij_tab,ia_ja_pairs)
$omp_end_parallel
@ -241,6 +246,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
! particles.
! Assume N_int is already provided.
END_DOC
integer,parameter :: size_max = $size_max
$declarations
integer(omp_lock_kind) :: lck
integer(bit_kind),intent(in) :: key_in(N_int,2)
@ -260,13 +266,12 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
integer,parameter :: size_max = $size_max
double precision :: hij_elec, thresh
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem
PROVIDE mo_integrals_map ref_bitmask_energy
PROVIDE mo_bielec_integrals_in_map
PROVIDE mo_integrals_map ref_bitmask_energy key_pattern_not_in_ref
PROVIDE mo_bielec_integrals_in_map reference_energy psi_ref_coef psi_ref
$set_i_H_j_threshold
@ -311,7 +316,6 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
integer(bit_kind) :: test(N_int,2)
double precision :: accu
accu = 0.d0
hij_elec = 0.d0
do ispin=1,2
other_spin = iand(ispin,1)+1
$omp_do
@ -325,33 +329,33 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
k_a = ishft(j_a-1,-bit_kind_shift)+1
l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
call i_H_j(hole,key_in,N_int,hij_elec)
if(dabs(hij_elec) .ge. thresh)then
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = hole(k,1)
keys_out(k,2,key_idx) = hole(k,2)
enddo
hij_tab(key_idx) = hij_elec
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work
$omp_unset_lock
key_idx = 0
endif
endif
if (key_idx == size_max) then
$omp_set_lock
$keys_work
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = hole(k,1)
keys_out(k,2,key_idx) = hole(k,2)
enddo
hij_tab(key_idx) = hij_elec
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work_unlocked
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
endif
if (key_idx == size_max) then
$keys_work_unlocked
$omp_set_lock
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
enddo ! ii
$omp_enddo
enddo ! ispin
$keys_work_unlocked
$omp_set_lock
$keys_work
$keys_work_locked
$omp_unset_lock
deallocate (keys_out,hij_tab,ia_ja_pairs)
$omp_end_parallel

470
src/Dets/connected_to_ref.irp.f
View File

@ -1,254 +1,256 @@
integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
use bitmasks
implicit none
integer, intent(in) :: Nint, N_past_in, Ndet
integer(bit_kind), intent(in) :: keys(ishft(Nint,-1),2,Ndet)
integer(bit_kind), intent(in) :: key(ishft(Nint,-1),2)
double precision, intent(in) :: thresh
integer :: N_past
integer :: i, l
integer :: degree_x2
logical :: det_is_not_or_may_be_in_ref, t
double precision :: hij_elec
! output : 0 : not connected
! i : connected to determinant i of the past
! -i : is the ith determinant of the refernce wf keys
ASSERT (Nint == N_int)
connected_to_ref = 0
N_past = max(1,N_past_in)
if (Nint == 1) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
if(degree_x2 == 0)then
use bitmasks
implicit none
integer, intent(in) :: Nint, N_past_in, Ndet
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
double precision, intent(in) :: thresh
integer :: N_past
integer :: i, l
integer :: degree_x2
logical :: det_is_not_or_may_be_in_ref, t
double precision :: hij_elec
! output : 0 : not connected
! i : connected to determinant i of the past
! -i : is the ith determinant of the refernce wf keys
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
connected_to_ref = 0
N_past = max(1,N_past_in)
if (Nint == 1) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
if(degree_x2 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
do i=N_past,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)) ) then
cycle
endif
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
do i=N_past,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)) ) then
cycle
endif
connected_to_ref = -i
return
enddo
return
else if (Nint==2) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i))) + &
popcnt(xor( key(2,1), keys(2,1,i))) + &
popcnt(xor( key(2,2), keys(2,2,i)))
if(degree_x2 == 0)then
enddo
return
else if (Nint==2) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i))) + &
popcnt(xor( key(2,1), keys(2,1,i))) + &
popcnt(xor( key(2,2), keys(2,2,i)))
if(degree_x2 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
!DIR$ LOOP COUNT (1000)
do i=N_past+1,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)).or. &
(key(2,1) /= keys(2,1,i)).or. &
(key(2,2) /= keys(2,2,i)) ) then
cycle
endif
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
!DIR$ LOOP COUNT (1000)
do i=N_past+1,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)).or. &
(key(2,1) /= keys(2,1,i)).or. &
(key(2,2) /= keys(2,2,i)) ) then
cycle
endif
connected_to_ref = -i
return
enddo
return
else if (Nint==3) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i))) + &
popcnt(xor( key(2,1), keys(2,1,i))) + &
popcnt(xor( key(2,2), keys(2,2,i))) + &
popcnt(xor( key(3,1), keys(3,1,i))) + &
popcnt(xor( key(3,2), keys(3,2,i)))
if(degree_x2 == 0)then
enddo
return
else if (Nint==3) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i))) + &
popcnt(xor( key(2,1), keys(2,1,i))) + &
popcnt(xor( key(2,2), keys(2,2,i))) + &
popcnt(xor( key(3,1), keys(3,1,i))) + &
popcnt(xor( key(3,2), keys(3,2,i)))
if(degree_x2 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
do i=N_past+1,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)).or. &
(key(2,1) /= keys(2,1,i)).or. &
(key(2,2) /= keys(2,2,i)).or. &
(key(3,1) /= keys(3,1,i)).or. &
(key(3,2) /= keys(3,2,i)) ) then
cycle
endif
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
do i=N_past+1,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)).or. &
(key(2,1) /= keys(2,1,i)).or. &
(key(2,2) /= keys(2,2,i)).or. &
(key(3,1) /= keys(3,1,i)).or. &
(key(3,2) /= keys(3,2,i)) ) then
cycle
endif
connected_to_ref = -i
return
enddo
return
else
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
!DEC$ LOOP COUNT MIN(3)
do l=2,Nint
degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) + &
popcnt(xor( key(l,2), keys(l,2,i)))
enddo
if(degree_x2 == 0)then
connected_to_ref = -i
enddo
return
else
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
!DEC$ LOOP COUNT MIN(3)
do l=2,Nint
degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +&
popcnt(xor( key(l,2), keys(l,2,i)))
enddo
if(degree_x2 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
do i=N_past+1,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)) ) then
cycle
else
connected_to_ref = -1
!DEC$ LOOP COUNT MIN(3)
do l=2,Nint
if ( (key(l,1) /= keys(l,1,i)).or. &
(key(l,2) /= keys(l,2,i)) ) then
connected_to_ref = 0
exit
endif
enddo
if (connected_to_ref /= 0) then
return
endif
endif
enddo
endif
endif
do i=N_past+1,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)) ) then
cycle
else
connected_to_ref = -1
!DEC$ LOOP COUNT MIN(3)
do l=2,Nint
if ( (key(l,1) /= keys(l,1,i)).or. &
(key(l,2) /= keys(l,2,i)) ) then
connected_to_ref = 0
exit
endif
enddo
if (connected_to_ref /= 0) then
return
endif
endif
enddo
endif
end
logical function det_is_not_or_may_be_in_ref(key,Nint)
implicit none
! If true, det is not in ref
! If false, det may be in ref
integer, intent(in) :: key(Nint,2), Nint
integer :: key_int
integer*1 :: key_short(4)
!DIR$ ATTRIBUTES ALIGN : 32 :: key_short
equivalence (key_int,key_short)
integer :: i, ispin
det_is_not_or_may_be_in_ref = .False.
do ispin=1,2
do i=1,Nint
key_int = key(i,ispin)
if ( &
key_pattern_not_in_ref(key_short(1), i, ispin) &
.or.key_pattern_not_in_ref(key_short(2), i, ispin) &
.or.key_pattern_not_in_ref(key_short(3), i, ispin) &
.or.key_pattern_not_in_ref(key_short(4), i, ispin) &
) then
det_is_not_or_may_be_in_ref = .True.
use bitmasks
implicit none
BEGIN_DOC
! If true, det is not in ref
! If false, det may be in ref
END_DOC
integer(bit_kind), intent(in) :: key(Nint,2), Nint
integer(bit_kind) :: key_int
integer*1 :: key_short(bit_kind)
!DIR$ ATTRIBUTES ALIGN : 32 :: key_short
equivalence (key_int,key_short)
integer :: i, ispin, k
det_is_not_or_may_be_in_ref = .False.
do ispin=1,2
do i=1,Nint
key_int = key(i,ispin)
do k=1,bit_kind
det_is_not_or_may_be_in_ref = &
det_is_not_or_may_be_in_ref .or. &
key_pattern_not_in_ref(key_short(k), i, ispin)
enddo
if(det_is_not_or_may_be_in_ref) then
return
endif
endif
enddo
enddo
enddo
end
BEGIN_PROVIDER [ logical, key_pattern_not_in_ref, (-128:127,N_int,2) ]
implicit none
BEGIN_DOC
! Min and max values of the integers of the keys of the reference
END_DOC
integer :: i, j, ispin
integer :: key
integer*1 :: key_short(4)
equivalence (key,key_short)
integer :: idx
key_pattern_not_in_ref = .True.
do j=1,N_det
do ispin=1,2
do i=1,N_int
key = psi_det(i,ispin,j)
key_pattern_not_in_ref( key_short(1), i, ispin ) = .False.
key_pattern_not_in_ref( key_short(2), i, ispin ) = .False.
key_pattern_not_in_ref( key_short(3), i, ispin ) = .False.
key_pattern_not_in_ref( key_short(4), i, ispin ) = .False.
enddo
use bitmasks
implicit none
BEGIN_DOC
! Min and max values of the integers of the keys of the reference
END_DOC
integer :: i, j, ispin
integer(bit_kind) :: key
integer*1 :: key_short(bit_kind)
equivalence (key,key_short)
integer :: idx, k
key_pattern_not_in_ref = .True.
do j=1,N_det
do ispin=1,2
do i=1,N_int
key = psi_det(i,ispin,j)
do k=1,bit_kind
key_pattern_not_in_ref( key_short(k), i, ispin ) = .False.
enddo
enddo
enddo
enddo
enddo
END_PROVIDER

3
src/Dets/filter_connected.irp.f
View File

@ -95,6 +95,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (sze > 0)
l=1
@ -118,8 +119,8 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then

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

@ -492,9 +492,9 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
use bitmasks
implicit none
integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate
integer, intent(in) :: keys(Nint,2,Ndet_max)
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)
integer, intent(in) :: key(Nint,2)
double precision, intent(out) :: i_H_psi_array(Nstate)
integer :: i, ii,j
@ -503,6 +503,11 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
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(keys,key,Nint,Ndet,idx)
do ii=1,idx(0)
@ -512,6 +517,7 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
do j = 1, Nstate
i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
enddo
print *, 'x', coef(i,1), hij, i_H_psi_array(1)
enddo
end

30
src/Hartree_Fock_AOs/README.rst
View File

@ -1,30 +0,0 @@
===================
Hartree-Fock Module
===================
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.

2
src/NEEDED_MODULES
View File

@ -1 +1 @@
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation Selection

20
src/Perturbation/README.rst
View File

@ -82,6 +82,26 @@ Documentation
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L1>`_
compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
.br
for the various n_st states.
.br
c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
.br
e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
.br
`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L33>`_
compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
.br
for the various n_st states.
.br
e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> - E(i) ) - sqrt( ( <det_pert|H|det_pert> - E(i)) ^2 + 4 <psi(i)|H|det_pert>^2 )
.br
c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
.br
Needed Modules

12
src/Perturbation/epstein_nesbet.irp.f
View File

@ -24,7 +24,7 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
c_pert(i) = i_H_psi_array(i) / (E_ref(i) - H_pert_diag)
c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
enddo
@ -53,12 +53,16 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
double precision :: diag_H_mat_elem,delta_e
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array)
print *, 'coefs'
print *, psi_ref_coef(1:N_det_ref,1)
print *, '-----'
call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,N_det_ref,psi_ref_size,n_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
delta_e = H_pert_diag - E_ref(i)
delta_e = H_pert_diag - reference_energy(i)
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
enddo
print *, e_2_pert, delta_e , i_H_psi_array
end

7
src/Perturbation/perturbation_template.f
View File

@ -2,7 +2,7 @@ BEGIN_SHELL [ /usr/bin/env python ]
import perturbation
END_SHELL
subroutine perturb_buffer_$PERT(buffer,buffer_size,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint)
subroutine perturb_buffer_$PERT(buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint)
implicit none
BEGIN_DOC
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
@ -11,7 +11,8 @@ subroutine perturb_buffer_$PERT(buffer,buffer_size,sum_e_2_pert,sum_norm_pert,su
integer, intent(in) :: Nint, N_st, buffer_size
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st),sum_H_pert_diag(N_st)
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag
integer :: i,k, c_ref
integer :: connected_to_ref
@ -33,6 +34,8 @@ subroutine perturb_buffer_$PERT(buffer,buffer_size,sum_e_2_pert,sum_norm_pert,su
c_pert,e_2_pert,H_pert_diag,Nint,n_det_ref,n_st)
do k = 1,N_st
e_2_pert_buffer(k,i) = e_2_pert(k)
coef_pert_buffer(k,i) = c_pert(k)
sum_norm_pert(k) += c_pert(k) * c_pert(k)
sum_e_2_pert(k) += e_2_pert(k)
sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag