From f38c07a4441cb97d42ba5b5ed087bd741916e4d5 Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Thu, 23 Apr 2020 23:13:44 +0200
Subject: [PATCH] dynamic screening in BSE

---
 src/QuAcK/Bethe_Salpeter_Z_matrix_dynamic.f90 | 66 +++++++++++++++++++
 1 file changed, 66 insertions(+)
 create mode 100644 src/QuAcK/Bethe_Salpeter_Z_matrix_dynamic.f90

diff --git a/src/QuAcK/Bethe_Salpeter_Z_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_Z_matrix_dynamic.f90
new file mode 100644
index 0000000..157446f
--- /dev/null
+++ b/src/QuAcK/Bethe_Salpeter_Z_matrix_dynamic.f90
@@ -0,0 +1,66 @@
+subroutine Bethe_Salpeter_Z_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,Z_dyn)
+
+! Compute the dynamic part of the Bethe-Salpeter equation matrices
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: eGW(nBas)
+  double precision,intent(in)   :: OmRPA(nS)
+  double precision,intent(in)   :: OmBSE
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+  
+! Local variables
+
+  integer                       :: maxS
+  double precision              :: chi
+  double precision              :: eps
+  integer                       :: i,j,a,b,ia,jb,kc
+
+! Output variables
+
+  double precision,intent(out)  :: Z_dyn(nS,nS)
+
+! Initialization
+
+  Z_dyn(:,:) = 0d0
+
+! Number of poles taken into account 
+
+  maxS = nS
+
+! Build dynamic A matrix
+
+  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
+ 
+          chi = 0d0
+          do kc=1,maxS
+
+            eps = (OmBSE - OmRPA(kc) - (eGW(a) - eGW(i)))**2 + eta**2
+            chi = chi + rho(i,j,kc)*rho(a,b,kc)*((OmBSE - OmRPA(kc) - (eGW(a) - eGW(i)))/eps)**2
+
+            eps = (OmBSE - OmRPA(kc) - (eGW(b) - eGW(j)))**2 + eta**2
+            chi = chi + rho(i,j,kc)*rho(a,b,kc)*((OmBSE - OmRPA(kc) - (eGW(b) - eGW(j)))/eps)**2
+
+          enddo
+
+          Z_dyn(ia,jb) = Z_dyn(ia,jb) + 2d0*lambda*chi
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+end subroutine Bethe_Salpeter_Z_matrix_dynamic