10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-22 20:35:19 +01:00

Corrected bugs

This commit is contained in:
Anthony Scemama 2014-05-28 23:12:00 +02:00
parent 9897ba6c19
commit 1ba8f49949
19 changed files with 260 additions and 127 deletions

View File

@ -16,6 +16,8 @@ keys_work
copy_buffer
finalization
generate_psi_guess
init_thread
deinit_thread
""".split()
class H_apply(object):
@ -42,7 +44,7 @@ class H_apply(object):
!$OMP N_elec_in_key_hole_2,ia_ja_pairs,iproc) &
!$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, &
!$OMP hole_1, particl_1, hole_2, particl_2, &
!$OMP thresh,elec_alpha_num,i_generator)"""
!$OMP elec_alpha_num,i_generator)"""
s["omp_end_parallel"] = "!$OMP END PARALLEL"
s["omp_master"] = "!$OMP MASTER"
s["omp_end_master"] = "!$OMP END MASTER"
@ -54,7 +56,7 @@ class H_apply(object):
s["generate_psi_guess"] = """
! Sort H_jj to find the N_states lowest states
integer :: i,k
integer :: i
integer, allocatable :: iorder(:)
double precision, allocatable :: H_jj(:)
double precision, external :: diag_h_mat_elem
@ -81,9 +83,6 @@ class H_apply(object):
for k in s:
s[k] = ""
s["size_max"] = str(1024*128)
s["set_i_H_j_threshold"] = """
thresh = H_apply_threshold
"""
s["copy_buffer"] = "call copy_h_apply_buffer_to_wf"
self.data = s
@ -109,28 +108,41 @@ class H_apply(object):
integer, intent(in) :: N_st,Nint
double precision, intent(inout) :: sum_e_2_pert_in(N_st)
double precision, intent(inout) :: sum_norm_pert_in(N_st)
double precision, intent(inout) :: sum_H_pert_diag_in
double precision, intent(inout) :: sum_H_pert_diag_in(N_st)
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)
double precision :: sum_H_pert_diag(N_st)
double precision, allocatable :: e_2_pert_buffer(:,:)
double precision, allocatable :: coef_pert_buffer(:,:)
ASSERT (Nint == N_int)
"""
self.data["init_thread"] = """
allocate (e_2_pert_buffer(N_st,size_max), coef_pert_buffer(N_st,size_max))
do k=1,N_st
sum_e_2_pert(k) = 0.d0
sum_norm_pert(k) = 0.d0
sum_H_pert_diag(k) = 0.d0
enddo
"""
self.data["deinit_thread"] = """
!$OMP CRITICAL
do k=1,N_st
sum_e_2_pert_in(k) = sum_e_2_pert_in(k) + sum_e_2_pert(k)
sum_norm_pert_in(k) = sum_norm_pert_in(k) + sum_norm_pert(k)
sum_H_pert_diag_in(k) = sum_H_pert_diag_in(k) + sum_H_pert_diag(k)
enddo
!$OMP END CRITICAL
deallocate (e_2_pert_buffer, coef_pert_buffer)
"""
self.data["size_max"] = "256"
self.data["initialization"] = """
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
PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors
"""
self.data["keys_work"] = """
call perturb_buffer_%s(i_generator,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,N_int)
"""%(pert,)
self.data["finalization"] = """
sum_e_2_pert_in = sum_e_2_pert
sum_norm_pert_in = sum_norm_pert
sum_H_pert_diag_in = sum_H_pert_diag
"""
self.data["copy_buffer"] = ""
self.data["generate_psi_guess"] = ""
@ -140,17 +152,19 @@ class H_apply(object):
self.data["decls_main"] = """ integer, intent(in) :: N_st
double precision, intent(inout):: pt2(N_st)
double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag
double precision, intent(inout):: H_pert_diag(N_st)
PROVIDE CI_electronic_energy N_det_generators key_pattern_not_in_ref
pt2 = 0.d0
norm_pert = 0.d0
H_pert_diag = 0.d0
do k=1,N_st
pt2(k) = 0.d0
norm_pert(k) = 0.d0
H_pert_diag(k) = 0.d0
enddo
"""
if self.openmp:
self.data["omp_parallel"] += """&
!$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)"""
!$OMP SHARED(N_st) PRIVATE(e_2_pert_buffer,coef_pert_buffer) &
!$OMP PRIVATE(sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)"""
def set_selection_pt2(self,pert):
if self.selection_pt2 is not None:

View File

@ -143,7 +143,8 @@ end
subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: key_in(Nint,2),key_1(Nint,2),key_2(Nint,2),Nint
integer(bit_kind), intent(in) :: key_in(Nint,2),key_1(Nint,2),key_2(Nint,2)
integer :: Nint
integer,intent(out) :: i_ok
integer :: ispin,i_hole,k_hole,j_hole,i_particl,k_particl,j_particl,i_trou,degree,exc(0:2,2,2)
double precision :: phase

View File

@ -14,8 +14,7 @@ subroutine H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st
implicit none
character*(64), intent(in) :: perturbation
integer, intent(in) :: N_st
double precision, intent(inout):: pt2(N_st), norm_pert(N_st)
double precision,intent(inout) :: H_pert_diag
double precision, intent(inout):: pt2(N_st), norm_pert(N_st), H_pert_diag(N_st)
BEGIN_SHELL [ /usr/bin/env python ]
from perturbation import perturbations

View File

@ -8,6 +8,12 @@ Documentation
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`h_apply_cisd_selection <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/H_apply.irp.f#L13>`_
Undocumented
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/cisd_selection.irp.f#L1>`_
Undocumented
Needed Modules

View File

@ -3,27 +3,25 @@ program cisd
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:)
double precision :: H_pert_diag, E_old
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, iter
character*(64) :: perturbation
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st))
allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st))
pt2 = 1.d0
perturbation = "epstein_nesbet"
do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
E_old = CI_energy(1)
call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', E_old
print *, 'E+PT2 = ', E_old+pt2
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
if (abort_all) then
exit
endif
enddo
deallocate(pt2,norm_pert)
deallocate(pt2,norm_pert,H_pert_diag)
end

View File

@ -19,7 +19,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ]
! Uninitialized. Filled by H_apply subroutines.
END_DOC
integer :: iproc, sze
sze = 10
sze = 100
if (.not.associated(H_apply_buffer)) then
allocate(H_apply_buffer(0:nproc-1))
iproc = 0

View File

@ -14,31 +14,33 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
integer(bit_kind),allocatable :: keys_out(:,:,:)
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) :: hole_save(N_int,2)
integer(bit_kind) :: key(N_int,2),hole(N_int,2), particle(N_int,2)
integer(bit_kind) :: hole_tmp(N_int,2), particle_tmp(N_int,2)
integer(bit_kind), allocatable :: hole_save(:,:)
integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:)
integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:)
integer :: ii,i,jj,j,k,ispin,l
integer :: occ_particle(N_int*bit_kind_size,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_particle_tmp(N_int*bit_kind_size,2)
integer :: occ_hole_tmp(N_int*bit_kind_size,2)
integer, allocatable :: occ_particle(:,:), occ_hole(:,:)
integer, allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:)
integer :: kk,pp,other_spin,key_idx
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)
double precision :: mo_bielec_integral, thresh
double precision :: mo_bielec_integral
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem
integer :: iproc
$set_i_H_j_threshold
PROVIDE H_apply_threshold
$initialization
iproc = 0
$omp_parallel
!$ iproc = omp_get_thread_num()
allocate (keys_out(N_int,2,size_max))
allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), &
key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),&
particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), &
occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),&
occ_hole_tmp(N_int*bit_kind_size,2))
$init_thread
!print*,'key_in !!'
!call print_key(key_in)
@ -142,7 +144,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
j_b = occ_particle_tmp(jjj,other_spin)
ASSERT (j_b > 0)
ASSERT (j_b <= mo_tot_num)
if(dabs( mo_bielec_integral(j_a,j_b,i_a,i_b))<thresh)cycle
key = hole
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
@ -179,7 +180,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
ASSERT (j_b > 0)
ASSERT (j_b <= mo_tot_num)
if (j_b <= j_a) cycle
if(dabs( mo_bielec_integral(j_a,j_b,i_b,i_a))<thresh)cycle
key = hole
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
@ -205,7 +205,13 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
$omp_enddo
enddo ! ispin
$keys_work
deallocate (keys_out,ia_ja_pairs)
$deinit_thread
deallocate (ia_ja_pairs, &
keys_out, hole_save, &
key,hole, particle, hole_tmp,&
particle_tmp, occ_particle, &
occ_hole, occ_particle_tmp,&
occ_hole_tmp)
$omp_end_parallel
$finalization
abort_here = abort_all
@ -224,34 +230,35 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator $parameters
$declarations
integer ,intent(in) :: i_generator
integer(bit_kind),intent(in) :: key_in(N_int,2)
integer(bit_kind),intent(in) :: hole_1(N_int,2), particl_1(N_int,2)
integer(bit_kind),allocatable :: keys_out(:,:,:)
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_save(N_int,2)
integer(bit_kind) :: key(N_int,2),hole(N_int,2), particle(N_int,2)
integer(bit_kind) :: hole_tmp(N_int,2), particle_tmp(N_int,2)
integer(bit_kind),allocatable :: hole_save(:,:)
integer(bit_kind),allocatable :: key(:,:),hole(:,:), particle(:,:)
integer(bit_kind),allocatable :: hole_tmp(:,:), particle_tmp(:,:)
integer(bit_kind),allocatable :: hole_2(:,:), particl_2(:,:)
integer :: ii,i,jj,j,k,ispin,l
integer :: occ_particle(N_int*bit_kind_size,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_particle_tmp(N_int*bit_kind_size,2)
integer :: occ_hole_tmp(N_int*bit_kind_size,2)
integer,allocatable :: occ_particle(:,:), occ_hole(:,:)
integer,allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:)
integer :: kk,pp,other_spin,key_idx
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)
double precision :: thresh
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem
integer :: iproc
$set_i_H_j_threshold
PROVIDE H_apply_threshold
$initialization
iproc = 0
$omp_parallel
!$ iproc = omp_get_thread_num()
allocate (keys_out(N_int,2,size_max))
allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), &
key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),&
particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), &
occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),&
occ_hole_tmp(N_int*bit_kind_size,2))
$init_thread
!!!! First couple hole particle
do j = 1, N_int
hole(j,1) = iand(hole_1(j,1),key_in(j,1))
@ -312,7 +319,13 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator $parameters
$omp_enddo
enddo ! ispin
$keys_work
deallocate (keys_out,ia_ja_pairs)
$deinit_thread
deallocate (ia_ja_pairs, &
keys_out, hole_save, &
key,hole, particle, hole_tmp,&
particle_tmp, occ_particle, &
occ_hole, occ_particle_tmp,&
occ_hole_tmp)
$omp_end_parallel
$finalization
@ -330,9 +343,8 @@ subroutine $subroutine($params_main)
$decls_main
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_reference psi_generators
integer :: i_generator
integer :: i_generator, k
print *, irp_here
do i_generator=1,N_det_generators
call $subroutine_monoexc(psi_generators(1,1,i_generator), &
generators_bitmask(1,1,s_hole ,i_bitmask_gen), &

View File

@ -32,6 +32,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
endif
if (degree_x2 > 5) then
cycle
else
connected_to_ref = i
return
endif
enddo
@ -64,6 +67,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
endif
if (degree_x2 > 5) then
cycle
else
connected_to_ref = i
return
endif
enddo
@ -101,6 +107,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
endif
if (degree_x2 > 5) then
cycle
else
connected_to_ref = i
return
endif
enddo
@ -140,6 +149,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
endif
if (degree_x2 > 5) then
cycle
else
connected_to_ref = i
return
endif
enddo

View File

@ -57,18 +57,18 @@ END_PROVIDER
ifirst = 1
psi_det = 0_bit_kind
psi_coef = 0.d0
integer :: i
do i=1,N_int
psi_det(i,1,1) = HF_bitmask(i,1)
psi_det(i,2,1) = HF_bitmask(i,2)
enddo
do i=1,N_states
psi_coef(i,i) = 1.d0
enddo
endif
integer :: i
do i=1,N_int
psi_det(i,1,1) = HF_bitmask(i,1)
psi_det(i,2,1) = HF_bitmask(i,2)
enddo
do i=1,N_states
psi_coef(i,i) = 1.d0
enddo
END_PROVIDER

View File

10
src/Full_CI/H_apply.irp.f Normal file
View File

@ -0,0 +1,10 @@
use bitmasks
BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import *
s = H_apply("FCI",openmp=True)
s.set_selection_pt2("epstein_nesbet_2x2")
print s
END_SHELL

8
src/Full_CI/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

View File

@ -0,0 +1 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_full Hartree_Fock MonoInts MOs Nuclei Output Selectors_full Utils Perturbation

36
src/Full_CI/README.rst Normal file
View File

@ -0,0 +1,36 @@
==============
Full_CI Module
==============
Performs a perturbatively selected Full-CI.
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
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>`_
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
* `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>`_
* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_

30
src/Full_CI/full_ci.irp.f Normal file
View File

@ -0,0 +1,30 @@
program cisd
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:)
double precision :: H_pert_diag, E_old
integer :: N_st, degree
character*(64) :: perturbation
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st))
pt2 = 1.d0
diag_algorithm = "Lapack"
do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
print *, '-----'
E_old = CI_energy(1)
call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
if (abort_all) then
exit
endif
enddo
deallocate(pt2,norm_pert)
end

View File

@ -1 +1 @@
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation SingleRefMethod CISD_selected Selectors_full MP2
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation SingleRefMethod CISD_selected Selectors_full MP2 Generators_full Full_CI

View File

@ -3,7 +3,7 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
@ -21,7 +21,7 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
double precision :: diag_H_mat_elem
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase,delta_e
double precision :: phase,delta_e,h
integer :: h1,h2,p1,p2,s1,s2
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
@ -32,8 +32,9 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
delta_e = 1.d0/delta_e
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
H_pert_diag(i) = h
c_pert(i) = i_H_psi_array(i) *delta_e
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
enddo

View File

@ -1,15 +1,15 @@
subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various n_st states.
! for the various N_st states.
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
@ -18,14 +18,15 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem
double precision :: diag_H_mat_elem, h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
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)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,N_st
H_pert_diag(i) = h
if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
c_pert(i) = 0.d0
@ -35,18 +36,18 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
end
subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
!
! for the various n_st states.
! for the various N_st states.
!
! 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 )
!
@ -55,13 +56,14 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,delta_e
double precision :: diag_H_mat_elem,delta_e, h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_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 - CI_electronic_energy(i)
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,N_st
H_pert_diag(i) = h
delta_e = h - 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)))
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
@ -72,19 +74,19 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
end
subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various n_st states,
! for the various N_st states,
!
! but with the correction in the denominator
!
@ -105,40 +107,41 @@ subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,accu_e_corr,hij
double precision :: diag_H_mat_elem,accu_e_corr,hij,h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
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)
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
do i = 1, idx_repeat(0)
call i_H_j(psi_selectors(1,1,idx_repeat(i)),det_pert,Nint,hij)
accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1)
enddo
accu_e_corr = accu_e_corr / psi_selectors_coef(1,1)
H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,n_st
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,N_st
H_pert_diag(i) = h
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)
if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
else
c_pert(i) = 0.d0
endif
enddo
end
subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
!
! for the various n_st states.
! for the various N_st states.
!
! but with the correction in the denominator
!
@ -159,19 +162,20 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,accu_e_corr,hij,delta_e
double precision :: diag_H_mat_elem,accu_e_corr,hij,delta_e,h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
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)
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
do i = 1, idx_repeat(0)
call i_H_j(psi_selectors(1,1,idx_repeat(i)),det_pert,Nint,hij)
accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1)
enddo
accu_e_corr = accu_e_corr / psi_selectors_coef(1,1)
H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,n_st
delta_e = H_pert_diag - CI_electronic_energy(i)
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,N_st
H_pert_diag(i) = h
delta_e = h - 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)))
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
@ -182,19 +186,19 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,
end
subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various n_st states,
! for the various N_st states,
!
! but with the correction in the denominator
!
@ -219,10 +223,10 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j
double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
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)
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
call i_H_j(ref_bitmask,det_pert,Nint,h0j)
do i = 1, idx_repeat(0)
@ -230,13 +234,14 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1)
enddo
accu_e_corr = accu_e_corr / psi_selectors_coef(1,1)
H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
c_pert(1) = 1.d0/psi_selectors_coef(1,1) * i_H_psi_array(1) / (CI_electronic_energy(i) - H_pert_diag)
c_pert(1) = 1.d0/psi_selectors_coef(1,1) * i_H_psi_array(1) / (CI_electronic_energy(i) - h)
e_2_pert(1) = c_pert(i) * h0j
do i =2,n_st
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)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
else
c_pert(i) = 0.d0
endif

View File

@ -13,7 +13,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
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)
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
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
integer :: i,k, c_ref
integer :: connected_to_ref
@ -38,7 +38,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
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
sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag(k)
enddo
enddo