diff --git a/src/GW/GW_ppBSE_dynamic_kernel_D.f90 b/src/GW/GW_ppBSE_dynamic_kernel_D.f90
new file mode 100644
index 0000000..c8433bf
--- /dev/null
+++ b/src/GW/GW_ppBSE_dynamic_kernel_D.f90
@@ -0,0 +1,140 @@
+subroutine GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eGW,Om,rho,OmBSE,KD_dyn,ZD_dyn)
+
+! Compute the dynamic part of the Bethe-Salpeter equation matrices
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: ispin
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  integer,intent(in)            :: nOO
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: eGW(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+  double precision,intent(in)   :: OmBSE
+  
+! Local variables
+
+  double precision              :: dem,num
+  integer                       :: m
+  integer                       :: i,j,k,l
+  integer                       :: ij,kl
+
+! Output variables
+
+  double precision,intent(out)  :: KD_dyn(nOO,nOO)
+  double precision,intent(out)  :: ZD_dyn(nOO,nOO)
+
+! Initialization
+
+  KD_dyn(:,:) = 0d0
+  ZD_dyn(:,:) = 0d0
+
+! Build dynamic A matrix
+
+  if(ispin == 1) then
+
+    ij = 0
+    do i=nC+1,nO
+      do j=i,nO
+        ij = ij + 1
+ 
+        kl = 0
+        do k=nC+1,nO
+          do l=k,nO
+            kl = kl + 1
+  
+            do m=1,nS
+
+              dem = - OmBSE + eGW(k) - Om(m) + eGW(j)
+              num = rho(i,k,m)*rho(j,l,m)
+
+              KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
+              ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+              dem = - OmBSE + eGW(k) - Om(m) + eGW(i)
+              num = rho(j,k,m)*rho(i,l,m)
+
+              KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2)
+              ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+              dem = - OmBSE + eGW(l) - Om(m) + eGW(i)
+              num = rho(i,k,m)*rho(j,l,m)
+
+              KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
+              ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+              dem = - OmBSE + eGW(l) - Om(m) + eGW(j)
+              num = rho(j,k,m)*rho(i,l,m)
+
+              KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2)
+              ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+            end do
+
+          end do
+        end do
+
+      end do
+    end do
+
+  end if
+
+ if(ispin == 2) then
+
+    ij = 0
+    do i=nC+1,nO
+      do j=i+1,nO
+        ij = ij + 1
+ 
+        kl = 0
+        do k=nC+1,nO
+          do l=k+1,nO
+            kl = kl + 1
+  
+            do m=1,nS
+
+              dem = - OmBSE + eGW(k) - Om(m) + eGW(j)
+              num = rho(i,k,m)*rho(j,l,m)
+
+              KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
+              ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+              dem = - OmBSE + eGW(k) - Om(m) + eGW(i)
+              num = rho(j,k,m)*rho(i,l,m)
+
+              KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2)
+              ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+              dem = - OmBSE + eGW(l) - Om(m) + eGW(i)
+              num = rho(i,k,m)*rho(j,l,m)
+
+              KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
+              ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+              dem = - OmBSE + eGW(l) - Om(m) + eGW(j)
+              num = rho(j,k,m)*rho(i,l,m)
+
+              KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2)
+              ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+            end do
+
+          end do
+        end do
+
+      end do
+    end do
+
+  end if
+
+end subroutine