diff --git a/src/GW/excitation_density.f90 b/src/GW/excitation_density.f90
index cba698e..0baa5fb 100644
--- a/src/GW/excitation_density.f90
+++ b/src/GW/excitation_density.f90
@@ -19,19 +19,25 @@ subroutine excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
   double precision,intent(out)  :: rho(nBas,nBas,nS)
 
   rho(:,:,:) = 0d0   
-
-  jb = 0
-  do j=nC+1,nO
-    do b=nO+1,nBas-nR
-      jb = jb + 1
-      do ia=1,nS
-        do q=nC+1,nBas-nR
-          do p=nC+1,nBas-nR
-            rho(p,q,ia) = rho(p,q,ia) + ERI(p,j,q,b)*XpY(ia,jb)
-          enddo
+  !$OMP PARALLEL &
+  !$OMP SHARED(nC,nBas,nR,nO,nS,rho,ERI,XpY) &
+  !$OMP PRIVATE(q,p,jb,ia) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do q=nC+1,nBas-nR
+     do p=nC+1,nBas-nR
+        jb = 0
+        do j=nC+1,nO
+           do b=nO+1,nBas-nR
+              jb = jb + 1
+              do ia=1,nS
+                 rho(p,q,ia) = rho(p,q,ia) + ERI(p,j,q,b)*XpY(ia,jb)
+              enddo
+           enddo
         enddo
-      enddo
-    enddo
+     enddo
   enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
 
 end subroutine excitation_density