Merge pull request #142 from eginer/master

added the All_singles module

plugins/All_singles/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Generators_restart Perturbation Properties Selectors_no_sorted Utils

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.

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

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

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"] = """
     """

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
+

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)