rename LR matrices

src/CI/CIS.f90

@@ -45,7 +45,7 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
   if(singlet) then
 
     ispin = 1
-    call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
+    call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
  
     if(dump_matrix) then
       print*,'CIS matrix (singlet state)'
@@ -73,7 +73,7 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
   if(triplet) then
 
     ispin = 2
-    call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
+    call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
  
     if(dump_matrix) then
       print*,'CIS matrix (triplet state)'
@@ -98,4 +98,4 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
 
   endif
 
-end subroutine CIS
+end subroutine 

src/GF/G0F3.f90

@@ -430,4 +430,4 @@
 
   call print_G0F3(nBas,nO,e0,Z,eGF3)
 
-end subroutine G0F3
+end subroutine 

src/GF/SigmaC_GF2.f90

@@ -22,26 +22,29 @@ double precision function SigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
   SigmaC_GF2 = 0d0
 
   ! Occupied part of the correlation self-energy
-  do i=nC+1,nO
-     do j=nC+1,nO
-        do a=nO+1,nBas-nR
 
-           eps = w + eHF(a) - eHF(i) - eHF(j)
+  do i=nC+1,nO
-           SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*eps/(eps**2 + eta**2)
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+
+        eps = w + eHF(a) - eHF(i) - eHF(j)
+        SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*eps/(eps**2 + eta**2)
 
         end do
      end do
   end do
+
   ! Virtual part of the correlation self-energy
+
   do i=nC+1,nO
-     do a=nO+1,nBas-nR
+    do a=nO+1,nBas-nR
-        do b=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
 
-          eps = w + eHF(i) - eHF(a) - eHF(b)
+        eps = w + eHF(i) - eHF(a) - eHF(b)
-          SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*eps/(eps**2 + eta**2)
+        SigmaC_GF2 = SigmaC_GF2 + (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*eps/(eps**2 + eta**2)
 
-        end do
       end do
     end do
+  end do
 
-end function SigmaC_GF2
+end function 

src/GF/UG0F2.f90

@@ -128,4 +128,4 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,
 
   end if
 
-end subroutine UG0F2
+end subroutine 

src/GF/dSigmaC_GF2.f90

@@ -24,26 +24,29 @@ double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
   dSigmaC_GF2 = 0d0
 
   ! Occupied part of the correlation self-energy
-  do i=nC+1,nO
-     do j=nC+1,nO
-        do a=nO+1,nBas-nR
 
-           eps = w + eHF(a) - eHF(i) - eHF(j)
-           dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
-
-        end do
-     end do
-  end do
-  ! Virtual part of the correlation self-energy
   do i=nC+1,nO
+    do j=nC+1,nO
       do a=nO+1,nBas-nR
-        do b=nO+1,nBas-nR
 
-          eps = w + eHF(i) - eHF(a) - eHF(b)
+        eps = w + eHF(a) - eHF(i) - eHF(j)
-          dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+        dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
-        end do
       end do
     end do
+  end do
 
-end function dSigmaC_GF2
+  ! Virtual part of the correlation self-energy
+
+  do i=nC+1,nO
+    do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+
+        eps = w + eHF(i) - eHF(a) - eHF(b)
+        dSigmaC_GF2 = dSigmaC_GF2 - (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+
+      end do
+    end do
+  end do
+
+end function 

src/GF/evGF3.f90

@@ -491,4 +491,4 @@
 
   endif
 
-end subroutine evGF3
+end subroutine 

src/GF/evUGF2.f90

@@ -201,4 +201,4 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
 
   endif
 
-end subroutine evUGF2
+end subroutine 

src/GF/print_G0F2.f90

@@ -50,4 +50,4 @@ subroutine print_G0F2(nBas,nO,eHF,Sig,eGF2,Z,ENuc,ERHF,Ec)
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
-end subroutine print_G0F2
+end subroutine 

src/GF/print_G0F3.f90

@@ -38,4 +38,4 @@ subroutine print_G0F3(nBas,nO,eHF,Z,eGF3)
   write(*,*)'-------------------------------------------------------------'
   write(*,*)
 
-end subroutine print_G0F3
+end subroutine 

src/GF/print_UG0F2.f90

@@ -68,6 +68,4 @@ subroutine print_UG0F2(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec)
               -------------------------------------------------'
   write(*,*)
 
-end subroutine print_UG0F2
+end subroutine 
-
-

src/GF/print_evGF2.f90

@@ -56,4 +56,4 @@ subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2,ENuc,ERHF,Ec)
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
-end subroutine print_evGF2
+end subroutine 

src/GF/print_evGF3.f90

@@ -41,4 +41,4 @@ subroutine print_evGF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
   write(*,*)'-------------------------------------------------------------'
   write(*,*)
 
-end subroutine print_evGF3
+end subroutine 

src/GF/print_evUGF2.f90

@@ -78,4 +78,4 @@ subroutine print_evUGF2(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec)
               -------------------------------------------------'
   write(*,*)
 
-end subroutine print_evUGF2
+end subroutine 

src/GF/print_qsGF2.f90

@@ -116,4 +116,4 @@ subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,P,T,V,J,K,F,SigC,Z, &
   endif
 
 
-end subroutine print_qsGF2
+end subroutine 

src/GF/print_qsUGF2.f90

@@ -175,4 +175,4 @@ subroutine print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,cGF2,PGF2,Ov,T,V,J,K,
 
   endif
 
-end subroutine print_qsUGF2
+end subroutine 

src/GF/qsUGF2.f90

@@ -337,4 +337,4 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
 
   end if
 
-end subroutine qsUGF2
+end subroutine 

src/GF/regularized_self_energy_GF2.f90

@@ -87,4 +87,4 @@ subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC
 
   Z(:) = 1d0/(1d0 - Z(:))
 
-end subroutine regularized_self_energy_GF2
+end subroutine

src/GT/SigmaC_GT.f90

@@ -34,22 +34,23 @@ double precision function SigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rh
 !----------------------------------------------
 ! Occupied part of the T-matrix self-energy 
 !----------------------------------------------
+
   do i=nC+1,nO
-     do cd=1,nVV
+    do cd=1,nVV
-        eps = w + e(i) - Omega1(cd)
+      eps = w + e(i) - Omega1(cd)
-        SigmaC_GT = SigmaC_GT + rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
+      SigmaC_GT = SigmaC_GT + rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
-     enddo
+    enddo
   enddo
-  write (*,*) SigmaC_GT
+
 !----------------------------------------------
 ! Virtual part of the T-matrix self-energy
 !----------------------------------------------
+
   do a=nO+1,nBas-nR
-     do kl=1,nOO
+    do kl=1,nOO
-        eps = w + e(a) - Omega2(kl)
+      eps = w + e(a) - Omega2(kl)
-        SigmaC_GT = SigmaC_GT + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
+      SigmaC_GT = SigmaC_GT + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
-     enddo
+    enddo
   enddo
-  write (*,*) SigmaC_GT
      
-end function SigmaC_GT
+end function 

src/GT/dSigmaC_GT.f90

@@ -35,20 +35,26 @@ double precision function dSigmaC_GT(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,r
 !----------------------------------------------
 ! Occupied part of the T-matrix self-energy 
 !----------------------------------------------
+
   do i=nC+1,nO
-     do cd=1,nVV
+    do cd=1,nVV
-        eps = w + e(i) - Omega1(cd)
+
-        dSigmaC_GT = dSigmaC_GT - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+      eps = w + e(i) - Omega1(cd)
-     enddo
+      dSigmaC_GT = dSigmaC_GT - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+
+    enddo
   enddo
+
 !----------------------------------------------
 ! Virtual part of the T-matrix self-energy
+!
 !----------------------------------------------
+
   do a=nO+1,nBas-nR
-     do kl=1,nOO
+    do kl=1,nOO
-        eps = w + e(a) - Omega2(kl)
+      eps = w + e(a) - Omega2(kl)
-        dSigmaC_GT = dSigmaC_GT - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+      dSigmaC_GT = dSigmaC_GT - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
-     enddo
+    enddo
   enddo
      
-end function dSigmaC_GT
+end function 

src/GT/soG0T0.f90

@@ -83,8 +83,7 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 
 ! Compute linear response
 
-  call linear_response_pp(ispin,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,1d0,seHF,sERI, & 
+  call ppLR(ispin,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,1d0,seHF,sERI,Omega1,X1,Y1,Omega2,X2,Y2,EcRPA)
-                          Omega1,X1,Y1,Omega2,X2,Y2,EcRPA)
 
   call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1)
   call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2)

src/HF/RHF_stability.f90

@@ -38,8 +38,8 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
  
   ispin = 1
 
-  call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
+  call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
-  call linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
+  call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
 
   AB(:,:) = A(:,:) + B(:,:)
 
@@ -111,8 +111,8 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
  
   ispin = 2
 
-  call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
+  call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
-  call linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
+  call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
 
   AB(:,:) = A(:,:) + B(:,:)
 
@@ -144,4 +144,4 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
   write(*,*)'-------------------------------------------------------------'
   write(*,*)
     
-end subroutine RHF_stability
+end subroutine 

src/LR/linear_response_BSE.f90

@@ -48,7 +48,7 @@ subroutine linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda
 
 ! Build A and B matrices 
 
-  call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
+  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
 
   if(BSE) A(:,:) = A(:,:) - A_BSE(:,:)
 
@@ -64,7 +64,7 @@ subroutine linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda
 
   else
 
-    call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
+    call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
 
     if(BSE) B(:,:) = B(:,:) - B_BSE(:,:)
 
@@ -104,4 +104,4 @@ subroutine linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda
 
     Ec = 0.5d0*(sum(Omega) - trace_matrix(nS,A))
 
-end subroutine linear_response_BSE
+end subroutine 

src/LR/linear_response_pp_BSE.f90

@@ -64,8 +64,8 @@ subroutine linear_response_pp_BSE(ispin,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,lambda,
 
 ! Build B, C and D matrices for the pp channel
 
-  call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C)
+  call ppLR_C(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C)
-  call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D)
+  call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D)
 
   if(BSE) then
 
@@ -86,7 +86,7 @@ subroutine linear_response_pp_BSE(ispin,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,lambda,
 
   else
 
-    call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B)
+    call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B)
     if(BSE) B(:,:) = B(:,:) - WB(:,:)
 
   ! Diagonal blocks 
@@ -118,4 +118,4 @@ subroutine linear_response_pp_BSE(ispin,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,lambda,
   if(abs(EcBSE - EcBSE1) > 1d-6 .or. abs(EcBSE - EcBSE2) > 1d-6) & 
     print*,'!!! Issue in pp-BSE linear reponse calculation BSE1 != BSE2 !!!'
 
-end subroutine linear_response_pp_BSE
+end subroutine 

src/LR/phLR.f90

@@ -48,7 +48,7 @@ subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY
 
 ! Build A and B matrices 
 
-  call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
+  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
 
 ! Tamm-Dancoff approximation
 
@@ -62,7 +62,7 @@ subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY
 
   else
 
-    call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
+    call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
 
     ! Build A + B and A - B matrices 
 

src/LR/phLR_A.f90

@@ -1,6 +1,6 @@
-subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A_lr)
+subroutine phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
 
-! Compute linear response
+! Compute resonant block of the ph channel
 
   implicit none
   include 'parameters.h'
@@ -8,7 +8,13 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
 ! Input variables
 
   logical,intent(in)            :: dRPA
-  integer,intent(in)            :: ispin,nBas,nC,nO,nV,nR,nS
+  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
   double precision,intent(in)   :: lambda
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas) 
@@ -22,7 +28,7 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
 
 ! Output variables
 
-  double precision,intent(out)  :: A_lr(nS,nS)
+  double precision,intent(out)  :: Aph(nS,nS)
 
 ! Direct RPA
 
@@ -42,8 +48,8 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
           do b=nO+1,nBas-nR
             jb = jb + 1
  
-            A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+            Aph(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
-                        + 2d0*lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
+                       + 2d0*lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
 
           end  do
         end  do
@@ -65,8 +71,8 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
           do b=nO+1,nBas-nR
             jb = jb + 1
  
-            A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+            Aph(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
-                        - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
+                       - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
 
           end  do
         end  do
@@ -88,8 +94,8 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
           do b=nO+1,nBas-nR
             jb = jb + 1
  
-            A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+            Aph(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
-                        + lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
+                       + lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
 
           end  do
         end  do
@@ -98,4 +104,4 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
 
   end if
 
-end subroutine linear_response_A_matrix
+end subroutine 

src/LR/phLR_B.f90

@@ -1,6 +1,6 @@
-subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_lr)
+subroutine phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
-! Compute linear response
+! Compute the coupling block of the ph channel
 
   implicit none
   include 'parameters.h'
@@ -20,7 +20,7 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
 
 ! Output variables
 
-  double precision,intent(out)  :: B_lr(nS,nS)
+  double precision,intent(out)  :: Bph(nS,nS)
 
 ! Direct RPA
 
@@ -40,7 +40,7 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
           do b=nO+1,nBas-nR
             jb = jb + 1
  
-            B_lr(ia,jb) = 2d0*lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
+            Bph(ia,jb) = 2d0*lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
  
           end do
         end do
@@ -62,7 +62,7 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
           do b=nO+1,nBas-nR
             jb = jb + 1
  
-            B_lr(ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
+            Bph(ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
  
           end do
         end do
@@ -84,7 +84,7 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
           do b=nO+1,nBas-nR
             jb = jb + 1
  
-            B_lr(ia,jb) = lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
+            Bph(ia,jb) = lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
  
           end do
         end do
@@ -93,4 +93,4 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
 
   end if
 
-end subroutine linear_response_B_matrix
+end subroutine 

src/LR/ppLR.f90

@@ -61,8 +61,8 @@ subroutine ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Omega1,X1,Y1,Ome
 
 ! Build B, C and D matrices for the pp channel
 
-  call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C)
+  call ppLR_C(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C)
-  call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D)
+  call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D)
 
   if(TDA) then
 
@@ -76,7 +76,7 @@ subroutine ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Omega1,X1,Y1,Ome
 
   else
 
-    call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B)
+    call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B)
 
   ! Diagonal blocks 
 

src/LR/ppLR_B.f90

@@ -1,4 +1,4 @@
-subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp)
+subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
 
 ! Compute the B matrix of the pp channel
 
@@ -25,7 +25,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp)
 
 ! Output variables
 
-  double precision,intent(out)  :: B_pp(nVV,nOO)
+  double precision,intent(out)  :: Bpp(nVV,nOO)
 
 ! Build B matrix for the singlet manifold
  
@@ -40,7 +40,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp)
           do j=i,nO
             ij = ij + 1
  
-            B_pp(ab,ij) = lambda*(ERI(a,b,i,j) + ERI(a,b,j,i))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
+            Bpp(ab,ij) = lambda*(ERI(a,b,i,j) + ERI(a,b,j,i))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
  
           end do
         end do
@@ -62,7 +62,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp)
          do j=i+1,nO
             ij = ij + 1
  
-            B_pp(ab,ij) = lambda*(ERI(a,b,i,j) - ERI(a,b,j,i))
+            Bpp(ab,ij) = lambda*(ERI(a,b,i,j) - ERI(a,b,j,i))
  
           end do
         end do
@@ -84,7 +84,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp)
          do j=nC+1,nO
             ij = ij + 1
  
-            B_pp(ab,ij) = lambda*ERI(a,b,i,j)
+            Bpp(ab,ij) = lambda*ERI(a,b,i,j)
  
           end do
         end do
@@ -93,4 +93,4 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp)
 
   end if
 
-end subroutine linear_response_B_pp
+end subroutine 

src/LR/ppLR_C.f90

@@ -1,4 +1,4 @@
-subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp)
+subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
 
 ! Compute the C matrix of the pp channel
 
@@ -13,7 +13,6 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
-  integer,intent(in)            :: nOO
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: lambda
   double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
@@ -27,7 +26,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp
 
 ! Output variables
 
-  double precision,intent(out)  :: C_pp(nVV,nVV)
+  double precision,intent(out)  :: Cpp(nVV,nVV)
 
 ! Define the chemical potential
  
@@ -47,8 +46,8 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp
          do d=c,nBas-nR
             cd = cd + 1
  
-            C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+            Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
-                          + lambda*(ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
+                         + lambda*(ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
  
           end do
         end do
@@ -70,8 +69,8 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp
          do d=c+1,nBas-nR
             cd = cd + 1
  
-            C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
+            Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
-                          + lambda*(ERI(a,b,c,d) - ERI(a,b,d,c))
+                         + lambda*(ERI(a,b,c,d) - ERI(a,b,d,c))
  
           end do
         end do
@@ -93,8 +92,8 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp
          do d=nO+1,nBas-nR
             cd = cd + 1
  
-            C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
+            Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
-                          + lambda*ERI(a,b,c,d)
+                         + lambda*ERI(a,b,c,d)
  
           end do
         end do
@@ -103,4 +102,4 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp
 
   end if
 
-end subroutine linear_response_C_pp
+end subroutine 

src/LR/ppLR_C_od.f90

@@ -1,4 +1,4 @@
-subroutine linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,C_pp)
+subroutine ppLR_C_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Cpp)
 
 ! Compute the C matrix of the pp channel (without diagonal term)
 
@@ -20,7 +20,7 @@ subroutine linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,C_p
 
 ! Output variables
 
-  double precision,intent(out)  :: C_pp(nVV,nVV)
+  double precision,intent(out)  :: Cpp(nVV,nVV)
  
 ! Build C matrix for the singlet manifold
 
@@ -35,7 +35,7 @@ subroutine linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,C_p
          do d=c,nBas-nR
             cd = cd + 1
  
-            C_pp(ab,cd) = lambda*(ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
+            Cpp(ab,cd) = lambda*(ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
  
           end do
         end do
@@ -57,7 +57,7 @@ subroutine linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,C_p
          do d=c+1,nBas-nR
             cd = cd + 1
  
-            C_pp(ab,cd) = lambda*(ERI(a,b,c,d) - ERI(a,b,d,c))
+            Cpp(ab,cd) = lambda*(ERI(a,b,c,d) - ERI(a,b,d,c))
  
           end do
         end do
@@ -79,7 +79,7 @@ subroutine linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,C_p
          do d=nO+1,nBas-nR
             cd = cd + 1
  
-            C_pp(ab,cd) = lambda*ERI(a,b,c,d)
+            Cpp(ab,cd) = lambda*ERI(a,b,c,d)
  
           end do
         end do
@@ -88,4 +88,4 @@ subroutine linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,C_p
 
   end if
 
-end subroutine linear_response_C_pp_od
+end subroutine 

src/LR/ppLR_D.f90

@@ -1,4 +1,4 @@
-subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp)
+subroutine ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
 
 ! Compute the D matrix of the pp channel
 
@@ -14,7 +14,6 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
   integer,intent(in)            :: nOO
-  integer,intent(in)            :: nVV
   double precision,intent(in)   :: lambda
   double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
   
@@ -27,7 +26,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp
 
 ! Output variables
 
-  double precision,intent(out)  :: D_pp(nOO,nOO)
+  double precision,intent(out)  :: Dpp(nOO,nOO)
 
 ! Define the chemical potential
 
@@ -47,8 +46,8 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp
          do l=k,nO
             kl = kl + 1
  
-            D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
+            Dpp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
-                          + lambda*(ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
+                         + lambda*(ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
  
           end do
         end do
@@ -70,8 +69,8 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp
          do l=k+1,nO
             kl = kl + 1
  
-            D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
+            Dpp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
-                          + lambda*(ERI(i,j,k,l) - ERI(i,j,l,k))
+                         + lambda*(ERI(i,j,k,l) - ERI(i,j,l,k))
  
           end do
         end do
@@ -93,8 +92,8 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp
          do l=nC+1,nO
             kl = kl + 1
  
-            D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
+            Dpp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
-                          + lambda*ERI(i,j,k,l)
+                         + lambda*ERI(i,j,k,l)
  
           end do
         end do
@@ -103,4 +102,4 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp
 
   end if
 
-end subroutine linear_response_D_pp
+end subroutine 

src/LR/ppLR_D_od.f90

@@ -1,4 +1,4 @@
-subroutine linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,D_pp)
+subroutine ppLR_D_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Dpp)
 
 ! Compute the D matrix of the pp channel (without the diagonal term)
 
@@ -20,7 +20,7 @@ subroutine linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,D_p
 
 ! Output variables
 
-  double precision,intent(out)  :: D_pp(nOO,nOO)
+  double precision,intent(out)  :: Dpp(nOO,nOO)
 
 ! Build the D matrix for the singlet manifold
 
@@ -35,7 +35,7 @@ subroutine linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,D_p
          do l=k,nO
             kl = kl + 1
  
-            D_pp(ij,kl) = lambda*(ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
+            Dpp(ij,kl) = lambda*(ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
  
           end do
         end do
@@ -57,7 +57,7 @@ subroutine linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,D_p
          do l=k+1,nO
             kl = kl + 1
  
-            D_pp(ij,kl) = lambda*(ERI(i,j,k,l) - ERI(i,j,l,k))
+            Dpp(ij,kl) = lambda*(ERI(i,j,k,l) - ERI(i,j,l,k))
  
           end do
         end do
@@ -79,7 +79,7 @@ subroutine linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,D_p
          do l=nC+1,nO
             kl = kl + 1
  
-            D_pp(ij,kl) = lambda*ERI(i,j,k,l)
+            Dpp(ij,kl) = lambda*ERI(i,j,k,l)
  
           end do
         end do
@@ -88,4 +88,4 @@ subroutine linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,D_p
 
   end if
 
-end subroutine linear_response_D_pp_od
+end subroutine 

src/RPA/ACFDT_pp_correlation_energy.f90

@@ -35,9 +35,9 @@ subroutine ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,
 
 ! Build pp matrices
 
-  call linear_response_B_pp   (ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,B)
+  call ppLR_B   (ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,B)
-  call linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,C)
+  call ppLR_C_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,C)
-  call linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,D)
+  call ppLR_D_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,D)
 
 ! Compute Tr(K x P_lambda)