LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-29 16:34:50 +02:00
Code Issues Releases Wiki Activity

Cleaned selection.

Browse Source
This commit is contained in:
Anthony Scemama 2014-05-26 13:09:32 +02:00
parent 583e8859d0
commit 37a639c62c
15 changed files with 129 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
scripts/generate_h_apply.py
View File

@ -35,7 +35,7 @@ class H_apply(object):
s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(SHARED) & s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, & !$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, &
!$OMP occ_particle,occ_hole,j_a,k_a,other_spin, & !$OMP occ_particle,occ_hole,j_a,k_a,other_spin, &
!$OMP hole_save,ispin,jj,l_a,hij_elec,hij_tab, & !$OMP hole_save,ispin,jj,l_a, &
!$OMP accu,i_a,hole_tmp,particle_tmp,occ_particle_tmp, & !$OMP accu,i_a,hole_tmp,particle_tmp,occ_particle_tmp, &
!$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,& !$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_1,N_elec_in_key_part_2, &
@ -121,9 +121,9 @@ class H_apply(object):
sum_e_2_pert = sum_e_2_pert_in sum_e_2_pert = sum_e_2_pert_in
sum_norm_pert = sum_norm_pert_in sum_norm_pert = sum_norm_pert_in
sum_H_pert_diag = sum_H_pert_diag_in sum_H_pert_diag = sum_H_pert_diag_in
PROVIDE reference_energy psi_ref_coef psi_ref PROVIDE CI_electronic_energy psi_ref_coef psi_ref
""" """
self.data["keys_work"] += """ self.data["keys_work"] = """
call perturb_buffer_%s(keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & 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) sum_norm_pert,sum_H_pert_diag,N_st,Nint)
"""%(pert,) """%(pert,)
@ -141,7 +141,7 @@ class H_apply(object):
double precision, intent(inout):: pt2(N_st) double precision, intent(inout):: pt2(N_st)
double precision, intent(inout):: norm_pert(N_st) double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag double precision, intent(inout):: H_pert_diag
PROVIDE reference_energy N_det_generators key_pattern_not_in_ref PROVIDE CI_electronic_energy N_det_generators key_pattern_not_in_ref
pt2 = 0.d0 pt2 = 0.d0
norm_pert = 0.d0 norm_pert = 0.d0
H_pert_diag = 0.d0 H_pert_diag = 0.d0

22
src/Dets/H_apply_template.f
View File

@ -11,7 +11,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
$declarations $declarations
integer(bit_kind),intent(in) :: key_in(N_int,2) integer(bit_kind),intent(in) :: key_in(N_int,2)
integer(bit_kind),allocatable :: keys_out(:,:,:) integer(bit_kind),allocatable :: keys_out(:,:,:)
double precision, allocatable :: hij_tab(:)
integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2) integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2)
integer(bit_kind), intent(in) :: hole_2(N_int,2), particl_2(N_int,2) integer(bit_kind), intent(in) :: hole_2(N_int,2), particl_2(N_int,2)
integer(bit_kind) :: hole_save(N_int,2) integer(bit_kind) :: hole_save(N_int,2)
@ -26,7 +25,7 @@ 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_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 :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
double precision :: hij_elec, mo_bielec_integral, thresh double precision :: mo_bielec_integral, thresh
integer, allocatable :: ia_ja_pairs(:,:,:) integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem double precision :: diag_H_mat_elem
integer :: iproc integer :: iproc
@ -38,7 +37,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
iproc = 0 iproc = 0
$omp_parallel $omp_parallel
!$ iproc = omp_get_thread_num() !$ iproc = omp_get_thread_num()
allocate (keys_out(N_int,2,size_max),hij_tab(size_max)) allocate (keys_out(N_int,2,size_max))
!print*,'key_in !!' !print*,'key_in !!'
!call print_key(key_in) !call print_key(key_in)
@ -88,7 +87,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
integer(bit_kind) :: test(N_int,2) integer(bit_kind) :: test(N_int,2)
double precision :: accu double precision :: accu
accu = 0.d0 accu = 0.d0
hij_elec = 0.d0
do ispin=1,2 do ispin=1,2
other_spin = iand(ispin,1)+1 other_spin = iand(ispin,1)+1
$omp_do $omp_do
@ -142,14 +140,11 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
k = ishft(j_b-1,-bit_kind_shift)+1 k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1 l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,other_spin) = ibset(key(k,other_spin),l) key(k,other_spin) = ibset(key(k,other_spin),l)
! call i_H_j(key,key_in,N_int,hij_elec)
! if(dabs(hij_elec)>=thresh) then
key_idx += 1 key_idx += 1
do k=1,N_int do k=1,N_int
keys_out(k,1,key_idx) = key(k,1) keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2) keys_out(k,2,key_idx) = key(k,2)
enddo enddo
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max) ASSERT (key_idx <= size_max)
if (key_idx == size_max) then if (key_idx == size_max) then
$keys_work $keys_work
@ -181,14 +176,11 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
key(k,ispin) = ibset(key(k,ispin),l) key(k,ispin) = ibset(key(k,ispin),l)
!! a^((+)_j_b(ispin) a_i_b(ispin) : mono exc :: orb(i_b,ispin) --> orb(j_b,ispin) !! a^((+)_j_b(ispin) a_i_b(ispin) : mono exc :: orb(i_b,ispin) --> orb(j_b,ispin)
! call i_H_j(key,key_in,N_int,hij_elec)
! if(dabs(hij_elec)>=thresh) then
key_idx += 1 key_idx += 1
do k=1,N_int do k=1,N_int
keys_out(k,1,key_idx) = key(k,1) keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2) keys_out(k,2,key_idx) = key(k,2)
enddo enddo
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max) ASSERT (key_idx <= size_max)
if (key_idx == size_max) then if (key_idx == size_max) then
$keys_work $keys_work
@ -201,7 +193,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
$omp_enddo $omp_enddo
enddo ! ispin enddo ! ispin
$keys_work $keys_work
deallocate (keys_out,hij_tab,ia_ja_pairs) deallocate (keys_out,ia_ja_pairs)
$omp_end_parallel $omp_end_parallel
$finalization $finalization
@ -220,7 +212,6 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
$declarations $declarations
integer(bit_kind),intent(in) :: key_in(N_int,2) integer(bit_kind),intent(in) :: key_in(N_int,2)
integer(bit_kind),allocatable :: keys_out(:,:,:) integer(bit_kind),allocatable :: keys_out(:,:,:)
double precision, allocatable :: hij_tab(:)
integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2) integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2)
integer(bit_kind) :: hole_2(N_int,2), particl_2(N_int,2) integer(bit_kind) :: hole_2(N_int,2), particl_2(N_int,2)
integer(bit_kind) :: hole_save(N_int,2) integer(bit_kind) :: hole_save(N_int,2)
@ -235,7 +226,7 @@ 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_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 :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
double precision :: hij_elec, thresh double precision :: thresh
integer, allocatable :: ia_ja_pairs(:,:,:) integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem double precision :: diag_H_mat_elem
integer :: iproc integer :: iproc
@ -247,7 +238,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
iproc = 0 iproc = 0
$omp_parallel $omp_parallel
!$ iproc = omp_get_thread_num() !$ iproc = omp_get_thread_num()
allocate (keys_out(N_int,2,size_max),hij_tab(size_max)) allocate (keys_out(N_int,2,size_max))
!!!! First couple hole particle !!!! First couple hole particle
do j = 1, N_int do j = 1, N_int
hole(j,1) = iand(hole_1(j,1),key_in(j,1)) hole(j,1) = iand(hole_1(j,1),key_in(j,1))
@ -300,7 +291,6 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
keys_out(k,1,key_idx) = hole(k,1) keys_out(k,1,key_idx) = hole(k,1)
keys_out(k,2,key_idx) = hole(k,2) keys_out(k,2,key_idx) = hole(k,2)
enddo enddo
hij_tab(key_idx) = hij_elec
if (key_idx == size_max) then if (key_idx == size_max) then
$keys_work $keys_work
key_idx = 0 key_idx = 0
@ -309,7 +299,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
$omp_enddo $omp_enddo
enddo ! ispin enddo ! ispin
$keys_work $keys_work
deallocate (keys_out,hij_tab,ia_ja_pairs) deallocate (keys_out,ia_ja_pairs)
$omp_end_parallel $omp_end_parallel
$finalization $finalization

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

@ -200,7 +200,8 @@ logical function det_is_not_or_may_be_in_ref(key,Nint)
! If true, det is not in ref ! If true, det is not in ref
! If false, det may be in ref ! If false, det may be in ref
END_DOC END_DOC
integer(bit_kind), intent(in) :: key(Nint,2), Nint integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key(Nint,2)
integer(bit_kind) :: key_int integer(bit_kind) :: key_int
integer*1 :: key_short(bit_kind) integer*1 :: key_short(bit_kind)
!DIR$ ATTRIBUTES ALIGN : 32 :: key_short !DIR$ ATTRIBUTES ALIGN : 32 :: key_short

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

@ -72,7 +72,6 @@ END_PROVIDER
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ integer, N_det_reference ] BEGIN_PROVIDER [ integer, N_det_reference ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC

6
src/Dets/diagonalize_CI.irp.f
View File

@ -8,6 +8,10 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ]
else else
diag_algorithm = "Lapack" diag_algorithm = "Lapack"
endif endif
if (N_det < N_states) then
diag_algorithm = "Lapack"
endif
END_PROVIDER END_PROVIDER
@ -18,7 +22,7 @@ BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ]
END_DOC END_DOC
integer :: j integer :: j
do j=1,min(N_states,N_det) do j=1,N_states
CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion
enddo enddo

4
src/Makefile.config.example
View File

@ -1,6 +1,6 @@
OPENMP =1 OPENMP =1
PROFILE =0 PROFILE =0
DEBUG = 0 DEBUG = 1
IRPF90_FLAGS+= --align=32 IRPF90_FLAGS+= --align=32
FC = ifort -g FC = ifort -g
@ -25,5 +25,7 @@ endif
ifeq ($(DEBUG),1) ifeq ($(DEBUG),1)
FC += -C -traceback -fpe0 FC += -C -traceback -fpe0
FCFLAGS+= -axSSE2
IRPF90_FLAGS += -a
#FCFLAGS =-O0 #FCFLAGS =-O0
endif endif

23
src/Makefile.config.gfortran Normal file
View File

@ -0,0 +1,23 @@
OPENMP =1
PROFILE =0
DEBUG = 0
IRPF90_FLAGS=
FC = gfortran -ffree-line-length-none
FCFLAGS=
FCFLAGS+= -O2
MKL=-L /opt/intel/composerxe/mkl/lib/intel64/ -lmkl_gf_lp64 -lmkl_core -lmkl_gnu_thread -fopenmp
ifeq ($(PROFILE),1)
FC += -p -g
CXX += -pg
endif
ifeq ($(OPENMP),1)
FC += -fopenmp
endif
ifeq ($(DEBUG),1)
FCFLAGS = -fcheck=all
#FCFLAGS =-O0
endif

29
src/Makefile.config.ifort Normal file
View File

@ -0,0 +1,29 @@
OPENMP =1
PROFILE =0
DEBUG = 0
IRPF90_FLAGS+= --align=32
FC = ifort -g
FCFLAGS=
FCFLAGS+= -axAVX,SSE4.3
FCFLAGS+= -O2
FCFLAGS+= -ip
FCFLAGS+= -opt-prefetch
FCFLAGS+= -ftz
MKL=-mkl=parallel
ifeq ($(PROFILE),1)
FC += -p -g
CXX += -pg
endif
ifeq ($(OPENMP),1)
FC += -openmp
CXX += -fopenmp
endif
ifeq ($(DEBUG),1)
FC += -C -traceback -fpe0
IRPF90_FLAGS += -a
#FCFLAGS =-O0
endif

28
src/Perturbation/README.rst
View File

@ -102,7 +102,7 @@ Documentation
e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> ) e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
.br .br
`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L35>`_ `pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L38>`_
compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
.br .br
for the various n_st states. for the various n_st states.
@ -112,7 +112,7 @@ Documentation
c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert> c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
.br .br
`pt2_epstein_nesbet_2x2_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L119>`_ `pt2_epstein_nesbet_2x2_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L129>`_
compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
.br .br
for the various n_st states. for the various n_st states.
@ -123,7 +123,7 @@ Documentation
.br .br
that can be repeated by repeating all the double excitations that can be repeated by repeating all the double excitations
.br .br
: you repeat all the correlation energy already taken into account in reference_energy(1) : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
.br .br
that could be repeated to this determinant. that could be repeated to this determinant.
.br .br
@ -134,7 +134,7 @@ Documentation
c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert> c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
.br .br
`pt2_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L69>`_ `pt2_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L75>`_
compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
.br .br
for the various n_st states, for the various n_st states,
@ -145,7 +145,7 @@ Documentation
.br .br
that can be repeated by repeating all the double excitations that can be repeated by repeating all the double excitations
.br .br
: you repeat all the correlation energy already taken into account in reference_energy(1) : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
.br .br
that could be repeated to this determinant. that could be repeated to this determinant.
.br .br
@ -156,7 +156,7 @@ Documentation
e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - (<det_pert|H|det_pert> ) ) e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - (<det_pert|H|det_pert> ) )
.br .br
`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L171>`_ `pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L185>`_
compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
.br .br
for the various n_st states, for the various n_st states,
@ -167,7 +167,7 @@ Documentation
.br .br
that can be repeated by repeating all the double excitations that can be repeated by repeating all the double excitations
.br .br
: you repeat all the correlation energy already taken into account in reference_energy(1) : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
.br .br
that could be repeated to this determinant. that could be repeated to this determinant.
.br .br
@ -197,24 +197,18 @@ Documentation
`selection_criterion_min <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L68>`_ `selection_criterion_min <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L68>`_
Threshold to select determinants. Set by selection routines. Threshold to select determinants. Set by selection routines.
`diagonalize <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L18>`_ `n_det_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L11>`_
Undocumented Undocumented
`n_det_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L36>`_ `psi_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L16>`_
Undocumented
`psi_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L41>`_
On what we apply <i|H|psi> for perturbation. If selection, it may be 0.9 of the norm. On what we apply <i|H|psi> for perturbation. If selection, it may be 0.9 of the norm.
`psi_ref_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L42>`_ `psi_ref_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L17>`_
On what we apply <i|H|psi> for perturbation. If selection, it may be 0.9 of the norm. On what we apply <i|H|psi> for perturbation. If selection, it may be 0.9 of the norm.
`psi_ref_size <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L32>`_ `psi_ref_size <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L7>`_
Undocumented Undocumented
`reference_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L6>`_
Reference energy
Needed Modules Needed Modules

7
src/Perturbation/cisd_selection.irp.f
View File

@ -13,13 +13,12 @@ program cisd
pt2 = 1.d0 pt2 = 1.d0
do while (maxval(abs(pt2(1:N_st))) > 1.d-6) do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
E_old = reference_energy(1) E_old = CI_energy(1)
call H_apply_cisd_selection(pt2, norm_pert, H_pert_diag, N_st) call H_apply_cisd_selection(pt2, norm_pert, H_pert_diag, N_st)
print *, 'N_det = ', N_det print *, 'N_det = ', N_det
print *, 'N_states = ', N_states print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2 print *, 'PT2 = ', pt2
print *, 'E = ', E_old+nuclear_repulsion print *, 'E = ', E_old
print *, 'E+PT2 = ', E_old+pt2+nuclear_repulsion print *, 'E+PT2 = ', E_old+pt2
enddo enddo
call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD)
end end

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

@ -19,16 +19,19 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
integer :: i,j integer :: i,j
double precision :: diag_H_mat_elem double precision :: diag_H_mat_elem
double precision, parameter :: eps = tiny(1.d0)
ASSERT (Nint == N_int) ASSERT (Nint == N_int)
ASSERT (Nint > 0) 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) 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) H_pert_diag = diag_H_mat_elem(det_pert,Nint)
print *, H_pert_diag if (dabs(CI_electronic_energy(i) - H_pert_diag) > 1.d-6) then
do i =1,n_st do i =1,n_st
c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag + eps) c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - H_pert_diag)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i) e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
enddo enddo
else
c_pert = 0.d0
e_2_pert = 0.d0
endif
end end
@ -53,15 +56,18 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
integer :: i,j integer :: i,j
double precision :: diag_H_mat_elem,delta_e double precision :: diag_H_mat_elem,delta_e
double precision, parameter :: eps = tiny(1.d0)
ASSERT (Nint == N_int) ASSERT (Nint == N_int)
ASSERT (Nint > 0) ASSERT (Nint > 0)
call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,N_det_ref,psi_ref_size,n_st,i_H_psi_array) 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) H_pert_diag = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st do i =1,n_st
delta_e = H_pert_diag - reference_energy(i) delta_e = H_pert_diag - CI_electronic_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))) 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)+eps) if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
else
c_pert(i) = 0.d0
endif
enddo enddo
end end
@ -86,7 +92,7 @@ subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
! !
! that can be repeated by repeating all the double excitations ! that can be repeated by repeating all the double excitations
! !
! : you repeat all the correlation energy already taken into account in reference_energy(1) ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
! !
! that could be repeated to this determinant. ! that could be repeated to this determinant.
! !
@ -111,8 +117,12 @@ subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
accu_e_corr = accu_e_corr / psi_ref_coef(1,1) accu_e_corr = accu_e_corr / psi_ref_coef(1,1)
H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,n_st do i =1,n_st
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) e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
if (dabs(CI_electronic_energy(i) - H_pert_diag) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - H_pert_diag)
else
c_pert(i) = 0.d0
endif
enddo enddo
end end
@ -136,7 +146,7 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,
! !
! that can be repeated by repeating all the double excitations ! that can be repeated by repeating all the double excitations
! !
! : you repeat all the correlation energy already taken into account in reference_energy(1) ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
! !
! that could be repeated to this determinant. ! that could be repeated to this determinant.
! !
@ -161,9 +171,13 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,
accu_e_corr = accu_e_corr / psi_ref_coef(1,1) accu_e_corr = accu_e_corr / psi_ref_coef(1,1)
H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,n_st do i =1,n_st
delta_e = H_pert_diag - reference_energy(i) delta_e = H_pert_diag - CI_electronic_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))) 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) if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
else
c_pert(i) = 0.d0
endif
enddo enddo
end end
@ -188,7 +202,7 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
! !
! that can be repeated by repeating all the double excitations ! that can be repeated by repeating all the double excitations
! !
! : you repeat all the correlation energy already taken into account in reference_energy(1) ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
! !
! that could be repeated to this determinant. ! that could be repeated to this determinant.
! !
@ -218,10 +232,14 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
accu_e_corr = accu_e_corr / psi_ref_coef(1,1) accu_e_corr = accu_e_corr / psi_ref_coef(1,1)
H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
c_pert(1) = 1.d0/psi_ref_coef(1,1) * i_H_psi_array(1) / (reference_energy(i) - H_pert_diag) c_pert(1) = 1.d0/psi_ref_coef(1,1) * i_H_psi_array(1) / (CI_electronic_energy(i) - H_pert_diag)
e_2_pert(1) = c_pert(i) * h0j e_2_pert(1) = c_pert(i) * h0j
do i =2,n_st do i =2,n_st
c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag) if (dabs(CI_electronic_energy(i) - H_pert_diag) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - H_pert_diag)
else
c_pert(i) = 0.d0
endif
e_2_pert(i) = c_pert(i) * i_H_psi_array(i) e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
enddo enddo

4
src/Perturbation/h_apply.py
View File

@ -3,8 +3,8 @@
from generate_h_apply import * from generate_h_apply import *
s = H_apply("cisd_pt2",openmp=True) s = H_apply("cisd_pt2",openmp=True)
s.set_perturbation("epstein_nesbet_2x2") #s.set_perturbation("epstein_nesbet_2x2")
#s.set_perturbation("epstein_nesbet") s.set_perturbation("epstein_nesbet")
#s["keys_work"] += """ #s["keys_work"] += """
#call fill_H_apply_buffer_cisd(key_idx,keys_out,N_int) #call fill_H_apply_buffer_cisd(key_idx,keys_out,N_int)
#""" #"""

5
src/Perturbation/perturbation_test.irp.f
View File

@ -10,8 +10,7 @@ program cisd
call H_apply_cisd_pt2(pt2, norm_pert, H_pert_diag, N_st) call H_apply_cisd_pt2(pt2, norm_pert, H_pert_diag, N_st)
print *, 'N_det = ', N_det print *, 'N_det = ', N_det
print *, 'N_states = ', N_states print *, 'pt2 = ', pt2(1)
print *, 'pt2 = ', pt2 print *, 'E = ', CI_energy(1)+pt2(1)
print *, 'E = ', reference_energy+pt2+nuclear_repulsion
return return
end end

25
src/Perturbation/temporary_stuff.irp.f
View File

@ -3,31 +3,6 @@ BEGIN_SHELL [ /bin/bash ]
./h_apply.py ./h_apply.py
END_SHELL END_SHELL
BEGIN_PROVIDER [ double precision, reference_energy, (N_states) ]
implicit none
BEGIN_DOC
! Reference energy
END_DOC
integer :: i
call diagonalize(psi_det,psi_coef,reference_energy,size(psi_coef,1),N_det,N_states,N_int,output_CISD)
SOFT_TOUCH psi_det psi_coef
END_PROVIDER
subroutine diagonalize(psi_det_in,psi_coef_in,eigvalues_out,psi_coef_dim,Ndet,N_st,Nint,iunit)
use bitmasks
implicit none
integer,intent(in) :: Nint,Ndet,N_st,psi_coef_dim, iunit
integer(bit_kind), intent(in) :: psi_det_in(Nint,2,Ndet)
double precision, intent(in) :: psi_coef_in(Ndet,N_st)
double precision, intent(out) :: eigvalues_out(N_st)
ASSERT (Nint == N_int)
call davidson_diag(psi_det_in,psi_coef_in,eigvalues_out,psi_coef_dim,Ndet,N_st,Nint,iunit)
end
BEGIN_PROVIDER [ integer, psi_ref_size ] BEGIN_PROVIDER [ integer, psi_ref_size ]
implicit none implicit none

8
src/Utils/map_module.f90
View File

@ -317,7 +317,7 @@ subroutine map_shrink(map,thr)
use map_module use map_module
implicit none implicit none
type (map_type), intent(inout) :: map type (map_type), intent(inout) :: map
type (map_type), intent(in) :: thr real(integral_kind), intent(in) :: thr
integer(map_size_kind) :: i integer(map_size_kind) :: i
integer(map_size_kind) :: icount integer(map_size_kind) :: icount
@ -371,7 +371,7 @@ subroutine map_update(map, key, value, sze, thr)
local_map%n_elements = map%map(idx_cache)%n_elements local_map%n_elements = map%map(idx_cache)%n_elements
do do
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1_8, local_map%n_elements) call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1, local_map%n_elements)
if (idx > 0_8) then if (idx > 0_8) then
local_map%value(idx) = local_map%value(idx) + value(i) local_map%value(idx) = local_map%value(idx) + value(i)
else else
@ -458,7 +458,7 @@ end
subroutine map_get(map, key, value) subroutine map_get(map, key, value)
use map_module use map_module
implicit none implicit none
type (map_type), intent(in) :: map type (map_type), intent(inout) :: map
integer(key_kind), intent(in) :: key integer(key_kind), intent(in) :: key
real(integral_kind), intent(out) :: value real(integral_kind), intent(out) :: value
integer(map_size_kind) :: idx_cache integer(map_size_kind) :: idx_cache
@ -476,7 +476,7 @@ subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx)
integer(key_kind), intent(in) :: key integer(key_kind), intent(in) :: key
integer(cache_map_size_kind), intent(in) :: ibegin, iend integer(cache_map_size_kind), intent(in) :: ibegin, iend
real(integral_kind), intent(out) :: value real(integral_kind), intent(out) :: value
integer(cache_map_size_kind), intent(in) :: idx integer(cache_map_size_kind), intent(inout) :: idx
call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend) call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend)
if (idx > 0) then if (idx > 0) then