Added MP2 wave function

configure

@@ -33,8 +33,8 @@ import sys
 from os.path import join
 
 if not any(i in ["--production", "--development"] for i in sys.argv):
-    print __doc__
-    sys.exit()
+    sys.argv += ["--development"]
+
 if len(sys.argv) != 3:
     print __doc__
     sys.exit()

plugins/MP2/H_apply.irp.f

@@ -6,5 +6,9 @@ from perturbation import perturbations
 s = H_apply("mp2")
 s.set_perturbation("Moller_plesset")
 print s
+
+s = H_apply("mp2_selection")
+s.set_selection_pt2("Moller_plesset")
+print s
 END_SHELL
 

plugins/MP2/mp2_wf.irp.f

@@ -0,0 +1,31 @@
+program mp2_wf
+  implicit none
+  BEGIN_DOC
+! Save the MP2 wave function
+  END_DOC
+  integer                        :: i,k
+
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, iter
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st))
+  
+  pt2 = 1.d0
+  selection_criterion = 1.e-12
+  selection_criterion_min = 1.e-12
+  TOUCH selection_criterion_min selection_criterion selection_criterion_factor
+  call H_apply_mp2_selection(pt2, norm_pert, H_pert_diag,  N_st)
+  psi_det = psi_det_sorted
+  psi_coef = psi_coef_sorted
+  print*,'N_det = ',N_det
+  print*,'-----'
+  print *,  'PT2                          = ', pt2(1)
+  print *,  'E                            = ', HF_energy
+  print *,  'E_before     +PT2            = ', HF_energy+pt2(1)
+  N_det = min(N_det,N_det_max)
+  touch N_det psi_det psi_coef
+  call save_wavefunction
+  call ezfio_set_mp2_energy(HF_energy+pt2(1))
+  deallocate(pt2,norm_pert,H_pert_diag)
+end

plugins/Perturbation/Moller_plesset.irp.f

@@ -26,10 +26,18 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
   ASSERT (Nint == N_int)
   ASSERT (Nint > 0)
   call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint)
-  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-  delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
-           (Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
-  delta_e = 1.d0/delta_e
+  if (degree == 2) then
+    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+    delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
+            (Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
+    delta_e = 1.d0/delta_e
+  else if (degree == 1) then
+    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+    delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1) 
+    delta_e = 1.d0/delta_e
+  else
+    delta_e = 0.d0
+  endif
 
   call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
   h = diag_H_mat_elem(det_pert,Nint)

scripts/generate_h_apply.py

@@ -201,7 +201,7 @@ class H_apply(object):
       """
       self.data["size_max"] = "256" 
       self.data["initialization"] = """
-      PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
+      PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
       """
       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, &
@@ -219,7 +219,7 @@ class H_apply(object):
   double precision, intent(inout):: norm_pert(N_st) 
   double precision, intent(inout):: H_pert_diag(N_st)
   double precision               :: delta_pt2(N_st), norm_psi(N_st), pt2_old(N_st)
-  PROVIDE CI_electronic_energy N_det_generators 
+  PROVIDE N_det_generators 
   do k=1,N_st
     pt2(k) = 0.d0
     norm_pert(k) = 0.d0

src/Determinants/H_apply.template.f

@@ -10,6 +10,9 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl
   integer                               :: highest, p1,p2,sp,ni,i,mi,nt,ns
   
   integer(bit_kind), intent(in)          :: key_prev(N_int, 2, *)
+  PROVIDE N_int
+  PROVIDE N_det
+ 
   $declarations