From 6d4199f989529dd3b95f517fecd5f6f76636361d Mon Sep 17 00:00:00 2001
From: EnzoMonino <enzo.monino@irsamc.ups-tlse.fr>
Date: Mon, 10 Jul 2023 13:53:39 +0200
Subject: [PATCH] add optimization in GF

---
 input/methods                          |  6 ++--
 src/GF/GF2_phBSE2_dynamic_kernel_A.f90 | 37 +++++++++++++------------
 src/GF/GF2_phBSE2_dynamic_kernel_B.f90 | 38 ++++++++++++++------------
 src/GF/GF2_phBSE2_static_kernel_A.f90  | 37 ++++++++++++++-----------
 src/GF/GF2_phBSE2_static_kernel_B.f90  | 38 ++++++++++++++------------
 src/GW/Bethe_Salpeter_A_matrix.f90     | 19 +++++++------
 6 files changed, 97 insertions(+), 78 deletions(-)

diff --git a/input/methods b/input/methods
index d9f3262..237541d 100644
--- a/input/methods
+++ b/input/methods
@@ -11,9 +11,9 @@
 # RPA* RPAx* crRPA ppRPA 
   F    F     F      F
 # G0F2* evGF2* qsGF2* G0F3 evGF3
-  F     F      F      F    F
+  T     F      F      F    F
 # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
-  T     F     F     F        F      F
+  F     F     F     F        F      F
 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
-  T      F      F      F      F      F
+  F      F      F      F      F      F
 # * unrestricted version available
diff --git a/src/GF/GF2_phBSE2_dynamic_kernel_A.f90 b/src/GF/GF2_phBSE2_dynamic_kernel_A.f90
index 31a01b7..63eb315 100644
--- a/src/GF/GF2_phBSE2_dynamic_kernel_A.f90
+++ b/src/GF/GF2_phBSE2_dynamic_kernel_A.f90
@@ -36,15 +36,15 @@ subroutine GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,
 
   if(ispin == 1) then
 
-    ia = 0
-    do i=nC+1,nO
-      do a=nO+1,nBas-nR
-        ia = ia + 1
+    jb = 0
+    do j=nC+1,nO
+      do b=nO+1,nBas-nR
+        jb = (b-nO) + (j-1)*(nBas-nO)
  
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
+        ia = 0
+        do i=nC+1,nO
+          do a=nO+1,nBas-nR
+            ia = (a-nO) + (i-1)*(nBas-nO) 
   
             do k=nC+1,nO
               do c=nO+1,nBas-nR
@@ -97,6 +97,7 @@ subroutine GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,
 
       end do
     end do
+!$omp end parallel do
 
   end if
 
@@ -104,15 +105,16 @@ subroutine GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,
 
   if(ispin == 2) then
 
-    ia = 0
-    do i=nC+1,nO
-      do a=nO+1,nBas-nR
-        ia = ia + 1
- 
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
+    jb = 0
+!$omp parallel do default(private) shared(A_dyn,ZA_dyn,ERI,OmBSE,num,dem,eGF,nO,nBas,eta,nC,nR)
+    do j=nC+1,nO
+      do b=nO+1,nBas-nR
+        jb = (b-nO) + (j-1)*(nBas-nO)
+
+        ia = 0
+        do i=nC+1,nO
+          do a=nO+1,nBas-nR
+            ia = (a-nO) + (i-1)*(nBas-nO)
   
             do k=nC+1,nO
               do c=nO+1,nBas-nR
@@ -161,6 +163,7 @@ subroutine GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,
 
       end do
     end do
+!$omp end parallel do
 
   end if
 
diff --git a/src/GF/GF2_phBSE2_dynamic_kernel_B.f90 b/src/GF/GF2_phBSE2_dynamic_kernel_B.f90
index 2ced77e..f63a640 100644
--- a/src/GF/GF2_phBSE2_dynamic_kernel_B.f90
+++ b/src/GF/GF2_phBSE2_dynamic_kernel_B.f90
@@ -33,15 +33,16 @@ subroutine GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,
 
   if(ispin == 1) then
 
-    ia = 0
-    do i=nC+1,nO
-      do a=nO+1,nBas-nR
-        ia = ia + 1
+    jb = 0
+!$omp parallel do default(private) shared(B_dyn,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
+    do j=nC+1,nO
+      do b=nO+1,nBas-nR
+        jb = (b-nO) + (j-1)*(nBas-nO)
  
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
+        ia = 0
+        do i=nC+1,nO
+          do a=nO+1,nBas-nR
+            ia = (a-nO) + (i-1)*(nBas-nO)
   
             do k=nC+1,nO
               do c=nO+1,nBas-nR
@@ -90,6 +91,7 @@ subroutine GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,
 
       end do
     end do
+!$omp end parallel do
 
   end if
 
@@ -97,16 +99,17 @@ subroutine GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,
 
   if(ispin == 2) then
 
-    ia = 0
-    do i=nC+1,nO
-      do a=nO+1,nBas-nR
-        ia = ia + 1
+    jb = 0
+!$omp parallel do default(private) shared(B_dyn,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
+    do j=nC+1,nO
+      do b=nO+1,nBas-nR
+        jb = (b-nO) + (j-1)*(nBas-nO)
  
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
-  
+        ia = 0
+        do i=nC+1,nO
+          do a=nO+1,nBas-nR
+            ia = (a-nO) + (i-1)*(nBas-nO) 
+
             do k=nC+1,nO
               do c=nO+1,nBas-nR
      
@@ -150,6 +153,7 @@ subroutine GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,
 
       end do
     end do
+!$omp end parallel do
 
   end if
 
diff --git a/src/GF/GF2_phBSE2_static_kernel_A.f90 b/src/GF/GF2_phBSE2_static_kernel_A.f90
index ce8fede..932ceed 100644
--- a/src/GF/GF2_phBSE2_static_kernel_A.f90
+++ b/src/GF/GF2_phBSE2_static_kernel_A.f90
@@ -33,15 +33,17 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e
 
   if(ispin == 1) then
 
-    ia = 0
-    do i=nC+1,nO
-      do a=nO+1,nBas-nR
-        ia = ia + 1
+    jb = 0
+!$omp parallel do default(private) shared(A_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
+    do j=nC+1,nO
+      do b=nO+1,nBas-nR
+        jb = (b-nO) + (j-1)*(nBas-nO)
  
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
+        ia = 0
+        do i=nC+1,nO
+          do a=nO+1,nBas-nR
+            ia = (a-nO) + (i-1)*(nBas-nO)
+            
   
             do k=nC+1,nO
               do c=nO+1,nBas-nR
@@ -90,6 +92,7 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e
 
       end do
     end do
+!$omp end parallel do
 
   end if
 
@@ -97,15 +100,16 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e
 
   if(ispin == 2) then
 
-    ia = 0
-    do i=nC+1,nO
-      do a=nO+1,nBas-nR
-        ia = ia + 1
+    jb = 0
+!$omp parallel do default(private) shared(A_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
+    do j=nC+1,nO
+      do b=nO+1,nBas-nR
+        jb = (b-nO) + (j-1)*(nBas-nO)
  
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
+        ia = 0
+        do i=nC+1,nO
+          do a=nO+1,nBas-nR
+            ia = (a-nO) + (i-1)*(nBas-nO)
   
             do k=nC+1,nO
               do c=nO+1,nBas-nR
@@ -150,6 +154,7 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e
 
       end do
     end do
+!$omp end parallel do
 
   end if
 
diff --git a/src/GF/GF2_phBSE2_static_kernel_B.f90 b/src/GF/GF2_phBSE2_static_kernel_B.f90
index 011949c..fdf8a53 100644
--- a/src/GF/GF2_phBSE2_static_kernel_B.f90
+++ b/src/GF/GF2_phBSE2_static_kernel_B.f90
@@ -33,15 +33,16 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e
 
   if(ispin == 1) then
 
-    ia = 0
-    do i=nC+1,nO
-      do a=nO+1,nBas-nR
-        ia = ia + 1
- 
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
+    jb = 0
+!$omp parallel do default(private) shared(B_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
+    do j=nC+1,nO
+      do b=nO+1,nBas-nR
+        jb = (b-nO) + (j-1)*(nBas-nO) 
+
+        ia = 0
+        do i=nC+1,nO
+          do a=nO+1,nBas-nR
+            ia = (a-nO) + (i-1)*(nBas-nO) 
   
             do k=nC+1,nO
               do c=nO+1,nBas-nR
@@ -90,6 +91,7 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e
 
       end do
     end do
+!$omp end parallel do
 
   end if
 
@@ -97,15 +99,16 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e
 
   if(ispin == 2) then
 
-    ia = 0
-    do i=nC+1,nO
-      do a=nO+1,nBas-nR
-        ia = ia + 1
+    jb = 0
+!$omp parallel do default(private) shared(B_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
+    do j=nC+1,nO
+      do b=nO+1,nBas-nR
+        jb = (b-nO) + (j-1)*(nBas-nO)
  
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
+        ia = 0
+        do i=nC+1,nO
+          do a=nO+1,nBas-nR
+            ia = (a-nO) + (i-1)*(nBas-nO)
   
             do k=nC+1,nO
               do c=nO+1,nBas-nR
@@ -150,6 +153,7 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e
 
       end do
     end do
+!$omp end parallel do
 
   end if
 
diff --git a/src/GW/Bethe_Salpeter_A_matrix.f90 b/src/GW/Bethe_Salpeter_A_matrix.f90
index 28dc945..38de0f2 100644
--- a/src/GW/Bethe_Salpeter_A_matrix.f90
+++ b/src/GW/Bethe_Salpeter_A_matrix.f90
@@ -24,14 +24,16 @@ subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
 
   double precision,intent(out)  :: A_lr(nS,nS)
 
-  ia = 0
-  do i=nC+1,nO
-    do a=nO+1,nBas-nR
-      ia = ia + 1
-      jb = 0
-      do j=nC+1,nO
-        do b=nO+1,nBas-nR
-          jb = jb + 1
+  jb = 0
+!$omp parallel do default(private) shared(A_lr,ERI,Omega,rho,nO,nBas,nS,chi,eps,eta,nC,nR,lambda)
+  do j=nC+1,nO
+    do b=nO+1,nBas-nR
+      jb = (b-nO) + (j-1)*(nBas-nO)
+
+      ia = 0
+      do i=nC+1,nO
+        do a=nO+1,nBas-nR
+          ia = (a-nO) + (i-1)*(nBas-nO) 
  
           chi = 0d0
           do kc=1,nS
@@ -46,5 +48,6 @@ subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
       enddo
     enddo
   enddo
+!$omp end parallel do
 
 end subroutine Bethe_Salpeter_A_matrix