From ff8fcfa42bcd9924870c5bf9077b244cd500f1ae Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 28 Nov 2017 17:11:44 +0100
Subject: [PATCH] qdpt test

---
 plugins/Full_CI/H_apply.irp.f            |   2 +
 plugins/Perturbation/pt2_equations.irp.f | 106 +++++++++++++++++++++++
 2 files changed, 108 insertions(+)

diff --git a/plugins/Full_CI/H_apply.irp.f b/plugins/Full_CI/H_apply.irp.f
index 014d0dbd..4f63e29e 100644
--- a/plugins/Full_CI/H_apply.irp.f
+++ b/plugins/Full_CI/H_apply.irp.f
@@ -4,11 +4,13 @@ from generate_h_apply import *
 
 s = H_apply("FCI")
 s.set_selection_pt2("epstein_nesbet_2x2")
+#s.set_selection_pt2("qdpt")
 s.unset_skip()
 print s
 
 s = H_apply("FCI_PT2")
 s.set_perturbation("epstein_nesbet_2x2")
+#s.set_perturbation("qdpt")
 s.unset_skip()
 s.unset_openmp()
 print s
diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f
index 57b3a814..f7ca1006 100644
--- a/plugins/Perturbation/pt2_equations.irp.f
+++ b/plugins/Perturbation/pt2_equations.irp.f
@@ -45,6 +45,58 @@ subroutine pt2_epstein_nesbet ($arguments)
   
 end
 
+subroutine pt2_qdpt ($arguments)
+  use bitmasks
+  implicit none
+  $declarations
+  
+  BEGIN_DOC
+  ! compute the QDPT first order coefficient and second order energetic contribution
+  !
+  ! for the various N_st states.
+  !
+  ! c_pert(i) = <psi(i)|H|det_pert>/( <psi(i)|H|psi(i)> - <det_pert|H|det_pert> )
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem_fock, h, E, diag_H_mat_elem, hij
+  double precision               :: i_H_psi_array(N_st)
+  integer :: degree
+  double precision :: delta_E
+  PROVIDE  selection_criterion
+
+  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)
+  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)
+  c_pert = 0.d0
+  do j=1,N_det_selectors
+    call get_excitation_degree(det_ref, psi_selectors(1,1,j), degree, Nint)
+    if (degree > 2) then
+      E = diag_H_mat_elem(psi_selectors(1,1,j),Nint)
+    else
+      E = diag_H_mat_elem_fock(det_ref,det_ref,fock_diag_tmp,Nint)
+    endif
+    delta_E = E-h
+!    delta_E = electronic_energy(1) - h
+    call i_H_j(psi_selectors(1,1,j),det_pert,Nint,hij)
+    if (dabs(delta_e) > 1.d-3) then
+        do i =1,N_st
+          c_pert(i) += psi_selectors_coef(j,i) * hij / delta_e
+        enddo
+    endif
+  enddo
+  do i =1,N_st
+    e_2_pert(i) = c_pert(i)*i_H_psi_array(i)
+    H_pert_diag(i) = h*c_pert(i)*c_pert(i)
+  enddo
+  
+end
+
 
 subroutine pt2_decontracted ($arguments)
   use bitmasks
@@ -249,6 +301,60 @@ subroutine pt2_moller_plesset ($arguments)
   
 end
 
+subroutine pt2_moller_plesset_general ($arguments)
+  use bitmasks
+  implicit none
+  $declarations
+  
+  BEGIN_DOC
+  ! compute the general Moller-Plesset perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various n_st states.
+  !
+  ! c_pert(i) = <psi(i)|H|det_pert>/(difference of orbital energies) 
+  !
+  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies) 
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem_fock
+  integer                        :: exc(0:2,2,2)
+  integer                        :: degree
+  double precision               :: phase,delta_e,h
+  double precision               :: i_H_psi_array(N_st)
+  integer                        :: h1,h2,p1,p2,s1,s2
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  call get_excitation(det_ref,det_pert,exc,degree,phase,Nint)
+  if (degree == 2) then
+    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+    delta_e = (mo_energy_expval(1,h1,s1,1) - mo_energy_expval(1,p1,s1,2)) + &
+              (mo_energy_expval(1,h2,s2,1) - mo_energy_expval(1,p2,s2,2))
+  else if (degree == 1) then
+    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+    delta_e = mo_energy_expval(1,h1,s1,1) - mo_energy_expval(1,p1,s1,2) 
+  else
+    delta_e = 0.d0
+  endif
+
+  if (dabs(delta_e) > 1.d-10) then
+    delta_e = 1.d0/delta_e
+    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
+    i_H_psi_array(:) = 0.d0
+    h = 0.d0
+  endif
+  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
+  
+end
+
+
 
 subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
   use bitmasks