Fixed MP2 and MP2 wf

plugins/Hartree_Fock/Fock_matrix.irp.f

@@ -73,10 +73,6 @@
      enddo
      
    endif
-   ! Introduce level shift here
-   do i = elec_alpha_num+1, mo_tot_num
-     Fock_matrix_mo(i,i) += level_shift
-   enddo
 
    do i = 1, mo_tot_num
      Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) 

plugins/Hartree_Fock/diagonalize_fock.irp.f

@@ -11,55 +11,55 @@
    double precision, allocatable  :: work(:), F(:,:), S(:,:)
    
 
-   if (mo_tot_num == ao_num) then
-       ! Solve H.C = E.S.C in AO basis set
-
-       allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
-       do j=1,ao_num
-         do i=1,ao_num
-           S(i,j) = ao_overlap(i,j)
-           F(i,j) = Fock_matrix_ao(i,j)
-         enddo
-       enddo
-
-       n = ao_num
-       lwork = 1+6*n + 2*n*n
-       liwork = 3 + 5*n
-       
-       allocate(work(lwork), iwork(liwork) )
-
-       lwork = -1
-       liwork = -1
-
-       call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
-         diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
-
-       if (info /= 0) then
-         print *,  irp_here//' failed : ', info
-         stop 1
-       endif
-       lwork = int(work(1))
-       liwork = iwork(1)
-       deallocate(work,iwork)
-       allocate(work(lwork), iwork(liwork) )
-
-       call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
-         diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
-
-       if (info /= 0) then
-         print *,  irp_here//' failed : ', info
-         stop 1
-       endif
-       do j=1,mo_tot_num
-         do i=1,ao_num
-           eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
-         enddo
-       enddo
-
-       deallocate(work, iwork, F, S)
-
-  else 
-
+!   if (mo_tot_num == ao_num) then
+!       ! Solve H.C = E.S.C in AO basis set
+!
+!       allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
+!       do j=1,ao_num
+!         do i=1,ao_num
+!           S(i,j) = ao_overlap(i,j)
+!           F(i,j) = Fock_matrix_ao(i,j)
+!         enddo
+!       enddo
+!
+!       n = ao_num
+!       lwork = 1+6*n + 2*n*n
+!       liwork = 3 + 5*n
+!       
+!       allocate(work(lwork), iwork(liwork) )
+!
+!       lwork = -1
+!       liwork = -1
+!
+!       call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
+!         diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
+!
+!       if (info /= 0) then
+!         print *,  irp_here//' failed : ', info
+!         stop 1
+!       endif
+!       lwork = int(work(1))
+!       liwork = iwork(1)
+!       deallocate(work,iwork)
+!       allocate(work(lwork), iwork(liwork) )
+!
+!       call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
+!         diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
+!
+!       if (info /= 0) then
+!         print *,  irp_here//' failed : ', info
+!         stop 1
+!       endif
+!       do j=1,mo_tot_num
+!         do i=1,ao_num
+!           eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
+!         enddo
+!       enddo
+!
+!       deallocate(work, iwork, F, S)
+!
+!  else 
+!
        ! Solve H.C = E.C in MO basis set
        
        allocate( F(mo_tot_num_align,mo_tot_num) )
@@ -69,6 +69,12 @@
          enddo
        enddo
        
+
+       ! Insert level shift here
+       do i = elec_alpha_num+1, mo_tot_num
+         Fock_matrix_mo(i,i) += level_shift
+       enddo
+
        n = mo_tot_num
        lwork = 1+6*n + 2*n*n
        liwork = 3 + 5*n
@@ -105,7 +111,12 @@
          0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1))
        deallocate(work, iwork, F)
 
-  endif
+       ! Remove level shift 
+       do i = elec_alpha_num+1, mo_tot_num
+         Fock_matrix_mo(i,i) -= level_shift
+       enddo
+
+!  endif
 
 END_PROVIDER
  

plugins/MP2/mp2_wf.irp.f

@@ -12,8 +12,7 @@ program mp2_wf
   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
+  selection_criterion_factor = 0.d0
   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

plugins/Perturbation/pt2_equations.irp.f

@@ -135,7 +135,6 @@ subroutine pt2_moller_plesset ($arguments)
   endif
 
   if (delta_e /= 0.d0) then
-  !  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
     call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
     h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
   else