From 92dd252be00138d1811fbc22493c4bf3c89d3d46 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 16 Feb 2016 11:14:19 +0100 Subject: [PATCH] added the All_singles module --- plugins/All_singles/NEEDED_CHILDREN_MODULES | 1 + plugins/All_singles/README.rst | 12 ++++ plugins/All_singles/all_1h_1p.irp.f | 75 ++++++++++++++++++++ plugins/All_singles/all_singles.irp.f | 76 +++++++++++++++++++++ scripts/generate_h_apply.py | 14 ++++ src/Bitmask/bitmask_cas_routines.irp.f | 44 ++++++++++++ src/Determinants/H_apply.template.f | 5 ++ 7 files changed, 227 insertions(+) create mode 100644 plugins/All_singles/NEEDED_CHILDREN_MODULES create mode 100644 plugins/All_singles/README.rst create mode 100644 plugins/All_singles/all_1h_1p.irp.f create mode 100644 plugins/All_singles/all_singles.irp.f diff --git a/plugins/All_singles/NEEDED_CHILDREN_MODULES b/plugins/All_singles/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..bb97ddb9 --- /dev/null +++ b/plugins/All_singles/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Generators_restart Perturbation Properties Selectors_no_sorted Utils diff --git a/plugins/All_singles/README.rst b/plugins/All_singles/README.rst new file mode 100644 index 00000000..b4b3f517 --- /dev/null +++ b/plugins/All_singles/README.rst @@ -0,0 +1,12 @@ +=========== +All_singles +=========== + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/All_singles/all_1h_1p.irp.f b/plugins/All_singles/all_1h_1p.irp.f new file mode 100644 index 00000000..e6d48749 --- /dev/null +++ b/plugins/All_singles/all_1h_1p.irp.f @@ -0,0 +1,75 @@ +program restart_more_singles + BEGIN_DOC + ! Generates and select single and double excitations of type 1h-1p + ! on the top of a given restart wave function of type CAS + END_DOC + read_wf = .true. + touch read_wf + print*,'ref_bitmask_energy = ',ref_bitmask_energy + call routine + +end +subroutine routine + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:) + integer :: N_st, degree + integer :: n_det_before + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + i = 0 + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + pt2=-1.d0 + E_before = ref_bitmask_energy + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + n_det_before = N_det + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_just_1h_1p(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + E_before = CI_energy + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + call save_wavefunction + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif + + enddo + + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1)) + enddo + endif + + call save_wavefunction + deallocate(pt2,norm_pert) +end diff --git a/plugins/All_singles/all_singles.irp.f b/plugins/All_singles/all_singles.irp.f new file mode 100644 index 00000000..3b5c5cce --- /dev/null +++ b/plugins/All_singles/all_singles.irp.f @@ -0,0 +1,76 @@ +program restart_more_singles + BEGIN_DOC + ! Generates and select single excitations + ! on the top of a given restart wave function + END_DOC + read_wf = .true. + touch read_wf + call routine + +end +subroutine routine + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision,allocatable :: E_before(:) + integer :: n_det_before + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + i = 0 + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + pt2=-1.d0 + E_before = ref_bitmask_energy + pt2_max = 1.d-10 + n_det_max = 200000 + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + n_det_before = N_det + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_just_mono(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy + print*,'pt2 = ',pt2 + print*,'E+PT2 = ',E_before + pt2 + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy + call save_wavefunction + + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif + enddo + + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1)) + enddo + endif + call save_wavefunction + deallocate(pt2,norm_pert,E_before) +end diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 02524c3d..b52e4445 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -24,6 +24,8 @@ skip init_main filter_integrals filter2h2p +filter_only_1h1p_single +filter_only_1h1p_double filterhole filterparticle do_double_excitations @@ -150,6 +152,18 @@ class H_apply(object): self["filterparticle"] = """ if(iand(ibset(0_bit_kind,j_a),hole(k_a,other_spin)).eq.0_bit_kind )cycle """ + + def filter_only_1h1p(self): + self["filter_only_1h1p_single"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1h1p(hole).eq..False.) cycle + """ + self["filter_only_1h1p_double"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1h1p(key).eq..False.) cycle + """ + + def unset_skip(self): self["skip"] = """ """ diff --git a/src/Bitmask/bitmask_cas_routines.irp.f b/src/Bitmask/bitmask_cas_routines.irp.f index 776d4546..4441fb22 100644 --- a/src/Bitmask/bitmask_cas_routines.irp.f +++ b/src/Bitmask/bitmask_cas_routines.irp.f @@ -445,3 +445,47 @@ integer function number_of_particles_verbose(key_in) + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) end +logical function is_a_1h1p(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_1h1p = .False. + if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.1)then + is_a_1h1p = .True. + endif + +end + +logical function is_a_1h(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_1h = .False. + if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.0)then + is_a_1h = .True. + endif + +end + +logical function is_a_1p(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_1p = .False. + if(number_of_holes(key_in).eq.0 .and. number_of_particles(key_in).eq.1)then + is_a_1p = .True. + endif + +end + +logical function is_a_2p(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_2p = .False. + if(number_of_holes(key_in).eq.0 .and. number_of_particles(key_in).eq.2)then + is_a_2p = .True. + endif + +end + diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 48e5d335..3f99f705 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -165,6 +165,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl endif logical :: check_double_excitation + logical :: is_a_1h1p logical :: b_cycle check_double_excitation = .True. iproc = iproc_in @@ -298,6 +299,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,other_spin) = ibset(key(k,other_spin),l) $filter2h2p + $filter_only_1h1p_double key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -348,6 +350,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,ispin) = ibset(key(k,ispin),l) $filter2h2p + $filter_only_1h1p_double key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -418,6 +421,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato integer(bit_kind) :: key_mask(N_int, 2) logical :: check_double_excitation + logical :: is_a_1h1p key_mask(:,:) = 0_bit_kind @@ -490,6 +494,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato $filterparticle hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) $filter2h2p + $filter_only_1h1p_single key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = hole(k,1)