From 1ba8f49949998e9ffbad1d0a82799451a5f8f711 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 28 May 2014 23:12:00 +0200 Subject: [PATCH] Corrected bugs --- scripts/generate_h_apply.py | 58 ++++++++----- src/CISD/SC2.irp.f | 3 +- src/CISD_selected/H_apply.irp.f | 3 +- src/CISD_selected/README.rst | 6 ++ src/CISD_selected/cisd_selection.irp.f | 12 ++- src/Dets/H_apply.irp.f | 2 +- src/Dets/H_apply_template.f | 72 +++++++++------- src/Dets/connected_to_ref.irp.f | 12 +++ src/Dets/determinants.irp.f | 20 ++--- src/Full_CI/ASSUMPTIONS.rst | 0 src/Full_CI/H_apply.irp.f | 10 +++ src/Full_CI/Makefile | 8 ++ src/Full_CI/NEEDED_MODULES | 1 + src/Full_CI/README.rst | 36 ++++++++ src/Full_CI/full_ci.irp.f | 30 +++++++ src/NEEDED_MODULES | 2 +- src/Perturbation/Moller_plesset.irp.f | 7 +- src/Perturbation/epstein_nesbet.irp.f | 101 ++++++++++++----------- src/Perturbation/perturbation_template.f | 4 +- 19 files changed, 260 insertions(+), 127 deletions(-) create mode 100644 src/Full_CI/ASSUMPTIONS.rst create mode 100644 src/Full_CI/H_apply.irp.f create mode 100644 src/Full_CI/Makefile create mode 100644 src/Full_CI/NEEDED_MODULES create mode 100644 src/Full_CI/README.rst create mode 100644 src/Full_CI/full_ci.irp.f diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 4a5ad695..73848488 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -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: diff --git a/src/CISD/SC2.irp.f b/src/CISD/SC2.irp.f index efd85202..35f19e2b 100644 --- a/src/CISD/SC2.irp.f +++ b/src/CISD/SC2.irp.f @@ -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 diff --git a/src/CISD_selected/H_apply.irp.f b/src/CISD_selected/H_apply.irp.f index 9b761963..1a01b53e 100644 --- a/src/CISD_selected/H_apply.irp.f +++ b/src/CISD_selected/H_apply.irp.f @@ -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 diff --git a/src/CISD_selected/README.rst b/src/CISD_selected/README.rst index 95499a8f..4fed26e5 100644 --- a/src/CISD_selected/README.rst +++ b/src/CISD_selected/README.rst @@ -8,6 +8,12 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +`h_apply_cisd_selection `_ + Undocumented + +`cisd `_ + Undocumented + Needed Modules diff --git a/src/CISD_selected/cisd_selection.irp.f b/src/CISD_selected/cisd_selection.irp.f index 40691065..6d88f61e 100644 --- a/src/CISD_selected/cisd_selection.irp.f +++ b/src/CISD_selected/cisd_selection.irp.f @@ -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 diff --git a/src/Dets/H_apply.irp.f b/src/Dets/H_apply.irp.f index 3a4922fb..be94a7a0 100644 --- a/src/Dets/H_apply.irp.f +++ b/src/Dets/H_apply.irp.f @@ -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 diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index adeff13d..5e310788 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -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)) 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)) 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 diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index 266898da..d16f8ce6 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -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 diff --git a/src/Full_CI/ASSUMPTIONS.rst b/src/Full_CI/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/Full_CI/H_apply.irp.f b/src/Full_CI/H_apply.irp.f new file mode 100644 index 00000000..c9b44cbd --- /dev/null +++ b/src/Full_CI/H_apply.irp.f @@ -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 + diff --git a/src/Full_CI/Makefile b/src/Full_CI/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/Full_CI/Makefile @@ -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 diff --git a/src/Full_CI/NEEDED_MODULES b/src/Full_CI/NEEDED_MODULES new file mode 100644 index 00000000..8a732901 --- /dev/null +++ b/src/Full_CI/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_full Hartree_Fock MonoInts MOs Nuclei Output Selectors_full Utils Perturbation diff --git a/src/Full_CI/README.rst b/src/Full_CI/README.rst new file mode 100644 index 00000000..f63f5a06 --- /dev/null +++ b/src/Full_CI/README.rst @@ -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 `_ +* `BiInts `_ +* `Bitmask `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Generators_full `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Selectors_full `_ +* `Utils `_ +* `Perturbation `_ + diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f new file mode 100644 index 00000000..796dbe6a --- /dev/null +++ b/src/Full_CI/full_ci.irp.f @@ -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 diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 725707f3..6151539c 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -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 diff --git a/src/Perturbation/Moller_plesset.irp.f b/src/Perturbation/Moller_plesset.irp.f index bca48358..7435f70c 100644 --- a/src/Perturbation/Moller_plesset.irp.f +++ b/src/Perturbation/Moller_plesset.irp.f @@ -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 diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f index b609cd02..d6a903c7 100644 --- a/src/Perturbation/epstein_nesbet.irp.f +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -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) = /( E(i) - ) ! @@ -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 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^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 diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f index 55aa7f56..e3033820 100644 --- a/src/Perturbation/perturbation_template.f +++ b/src/Perturbation/perturbation_template.f @@ -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