From 95807231fe68fe549ee6527f0c17fe747e6cf60c Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Mon, 17 Jul 2023 14:21:39 +0200
Subject: [PATCH] add MP module

---
 src/MP/MP.f90       | 69 +++++++++++++++++++++++++++++++++++++++++++++
 src/QuAcK/QuAcK.f90 | 49 +++++---------------------------
 2 files changed, 76 insertions(+), 42 deletions(-)
 create mode 100644 src/MP/MP.f90

diff --git a/src/MP/MP.f90 b/src/MP/MP.f90
new file mode 100644
index 0000000..fabe963
--- /dev/null
+++ b/src/MP/MP.f90
@@ -0,0 +1,69 @@
+subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
+
+! Moller-Plesset module
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: doMP2
+  logical,intent(in)            :: doMP3
+  logical,intent(in)            :: unrestricted
+
+  logical,intent(in)            :: regularize
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: epsHF(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+! Output variables
+
+!------------------------------------------------------------------------
+! Compute MP3 energy
+!------------------------------------------------------------------------                               
+
+  if(doMP2) then    
+                    
+    if(unrestricted) then
+                    
+      call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
+                    
+    else
+
+      call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+  
+    end if          
+                    
+  end if
+
+!------------------------------------------------------------------------
+! Compute MP3 energy
+!------------------------------------------------------------------------                               
+                    
+  if(doMP3) then    
+
+    if(unrestricted) then
+
+      write(*,*) 'MP3 NYI for UHF reference'
+      stop
+
+    else
+
+      call MP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+
+    end if
+
+  end if
+
+end subroutine 
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 3c0fdf5..bc7f46e 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -7,7 +7,7 @@ program QuAcK
   logical                       :: doHF,doRHF,doUHF,doRMOM,doUMOM
   logical                       :: dostab
   logical                       :: doKS
-  logical                       :: doMP2,doMP3
+  logical                       :: doMP,doMP2,doMP3
   logical                       :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
   logical                       :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
   logical                       :: doCIS,doCIS_D,doCID,doCISD,doFCI
@@ -229,7 +229,7 @@ program QuAcK
   end if
 
 !------------------------------------------------------------------------
-! KS module
+! Kohn-Sham module
 !------------------------------------------------------------------------
 
   if(doKS) then
@@ -360,23 +360,15 @@ program QuAcK
   end if
 
 !------------------------------------------------------------------------
-! Compute MP2 energy
+! Moller-Plesset module
 !------------------------------------------------------------------------
 
-  if(doMP2) then
+  doMP = doMP2 .or. doMP3
+
+  if(doMP) then
 
     call cpu_time(start_MP)
-
-    if(unrestricted) then
-
-      call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ENuc,EHF,epsHF)
-
-    else
-
-      call MP2(regMP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
-
-    end if
-
+    call MP(doMP2,doMP3,unrestricted,regMP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
     call cpu_time(end_MP)
 
     t_MP = end_MP - start_MP
@@ -385,33 +377,6 @@ program QuAcK
 
   end if
 
-!------------------------------------------------------------------------
-! Compute MP3 energy
-!------------------------------------------------------------------------
-
-  if(doMP3) then
-
-    call cpu_time(start_MP)
-    
-    if(unrestricted) then
-
-      write(*,*) 'MP3 NYI for UHF reference'
-      stop
-
-    else
-
-      call MP3(nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
-
-    end if
-
-    call cpu_time(end_MP)
-
-    t_MP = end_MP - start_MP
-    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP3 = ',t_MP,' seconds'
-    write(*,*)
-
-  end if
-
 !------------------------------------------------------------------------
 ! Perform CCD calculation
 !------------------------------------------------------------------------