diff --git a/src/CI/CIS.f90 b/src/CI/CIS.f90 index b4bf0fa..a4380b3 100644 --- a/src/CI/CIS.f90 +++ b/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 diff --git a/src/GF/G0F3.f90 b/src/GF/G0F3.f90 index ad972cb..90e67a6 100644 --- a/src/GF/G0F3.f90 +++ b/src/GF/G0F3.f90 @@ -430,4 +430,4 @@ call print_G0F3(nBas,nO,e0,Z,eGF3) -end subroutine G0F3 +end subroutine diff --git a/src/GF/SigmaC_GF2.f90 b/src/GF/SigmaC_GF2.f90 index 8f099fb..0fa97aa 100644 --- a/src/GF/SigmaC_GF2.f90 +++ b/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) - 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 i=nC+1,nO + 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 b=nO+1,nBas-nR + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR - 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) + 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) - end do end do end do + end do -end function SigmaC_GF2 +end function diff --git a/src/GF/UG0F2.f90 b/src/GF/UG0F2.f90 index 8c81da8..327d8e4 100644 --- a/src/GF/UG0F2.f90 +++ b/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 diff --git a/src/GF/dSigmaC_GF2.f90 b/src/GF/dSigmaC_GF2.f90 index bf7a74f..e0c3109 100644 --- a/src/GF/dSigmaC_GF2.f90 +++ b/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) - 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 + 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 + 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 diff --git a/src/GF/evGF3.f90 b/src/GF/evGF3.f90 index 8c80c80..8451ee5 100644 --- a/src/GF/evGF3.f90 +++ b/src/GF/evGF3.f90 @@ -491,4 +491,4 @@ endif -end subroutine evGF3 +end subroutine diff --git a/src/GF/evUGF2.f90 b/src/GF/evUGF2.f90 index 7d8c635..6849320 100644 --- a/src/GF/evUGF2.f90 +++ b/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 diff --git a/src/GF/print_G0F2.f90 b/src/GF/print_G0F2.f90 index 87b545b..beb688a 100644 --- a/src/GF/print_G0F2.f90 +++ b/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 diff --git a/src/GF/print_G0F3.f90 b/src/GF/print_G0F3.f90 index 09da736..12cbba8 100644 --- a/src/GF/print_G0F3.f90 +++ b/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 diff --git a/src/GF/print_UG0F2.f90 b/src/GF/print_UG0F2.f90 index 1495401..7a19acd 100644 --- a/src/GF/print_UG0F2.f90 +++ b/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 diff --git a/src/GF/print_evGF2.f90 b/src/GF/print_evGF2.f90 index 07c8b42..cbb0059 100644 --- a/src/GF/print_evGF2.f90 +++ b/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 diff --git a/src/GF/print_evGF3.f90 b/src/GF/print_evGF3.f90 index 6747058..c8cb2ef 100644 --- a/src/GF/print_evGF3.f90 +++ b/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 diff --git a/src/GF/print_evUGF2.f90 b/src/GF/print_evUGF2.f90 index f04f6f3..59a054a 100644 --- a/src/GF/print_evUGF2.f90 +++ b/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 diff --git a/src/GF/print_qsGF2.f90 b/src/GF/print_qsGF2.f90 index 1462e70..bae114b 100644 --- a/src/GF/print_qsGF2.f90 +++ b/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 diff --git a/src/GF/print_qsUGF2.f90 b/src/GF/print_qsUGF2.f90 index f87a9d5..caf3065 100644 --- a/src/GF/print_qsUGF2.f90 +++ b/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 diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 7a53c1d..11e8e8c 100644 --- a/src/GF/qsUGF2.f90 +++ b/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 diff --git a/src/GF/regularized_self_energy_GF2.f90 b/src/GF/regularized_self_energy_GF2.f90 index 53facf8..c8df16d 100644 --- a/src/GF/regularized_self_energy_GF2.f90 +++ b/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 diff --git a/src/GT/SigmaC_GT.f90 b/src/GT/SigmaC_GT.f90 index 770aee0..82a1cda 100644 --- a/src/GT/SigmaC_GT.f90 +++ b/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 - eps = w + e(i) - Omega1(cd) - SigmaC_GT = SigmaC_GT + rho1(p,i,cd)**2*eps/(eps**2 + eta**2) - enddo + do cd=1,nVV + eps = w + e(i) - Omega1(cd) + SigmaC_GT = SigmaC_GT + rho1(p,i,cd)**2*eps/(eps**2 + eta**2) + enddo enddo - write (*,*) SigmaC_GT + !---------------------------------------------- ! Virtual part of the T-matrix self-energy !---------------------------------------------- + do a=nO+1,nBas-nR - do kl=1,nOO - eps = w + e(a) - Omega2(kl) - SigmaC_GT = SigmaC_GT + rho2(p,a,kl)**2*eps/(eps**2 + eta**2) - enddo + do kl=1,nOO + eps = w + e(a) - Omega2(kl) + SigmaC_GT = SigmaC_GT + rho2(p,a,kl)**2*eps/(eps**2 + eta**2) + enddo enddo - write (*,*) SigmaC_GT -end function SigmaC_GT +end function diff --git a/src/GT/dSigmaC_GT.f90 b/src/GT/dSigmaC_GT.f90 index 0e73b11..ee59d99 100644 --- a/src/GT/dSigmaC_GT.f90 +++ b/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 - eps = w + e(i) - Omega1(cd) - dSigmaC_GT = dSigmaC_GT - rho1(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + 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 + + enddo enddo + !---------------------------------------------- ! Virtual part of the T-matrix self-energy +! !---------------------------------------------- + do a=nO+1,nBas-nR - do kl=1,nOO - eps = w + e(a) - Omega2(kl) - dSigmaC_GT = dSigmaC_GT - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + do kl=1,nOO + eps = w + e(a) - Omega2(kl) + dSigmaC_GT = dSigmaC_GT - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo enddo -end function dSigmaC_GT +end function diff --git a/src/GT/soG0T0.f90 b/src/GT/soG0T0.f90 index 540ed47..c749389 100644 --- a/src/GT/soG0T0.f90 +++ b/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, & - Omega1,X1,Y1,Omega2,X2,Y2,EcRPA) + call ppLR(ispin,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,1d0,seHF,sERI,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) diff --git a/src/HF/RHF_stability.f90 b/src/HF/RHF_stability.f90 index bade757..2b2da03 100644 --- a/src/HF/RHF_stability.f90 +++ b/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 linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B) + call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A) + 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 linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B) + call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A) + 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 diff --git a/src/LR/linear_response_BSE.f90 b/src/LR/linear_response_BSE.f90 index 46d2b03..5aa3835 100644 --- a/src/LR/linear_response_BSE.f90 +++ b/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 diff --git a/src/LR/linear_response_pp_BSE.f90 b/src/LR/linear_response_pp_BSE.f90 index 1584ba4..f336a4a 100644 --- a/src/LR/linear_response_pp_BSE.f90 +++ b/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 linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C) + 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 diff --git a/src/LR/phLR.f90 b/src/LR/phLR.f90 index 20d7eca..4b0e48e 100644 --- a/src/LR/phLR.f90 +++ b/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 diff --git a/src/LR/linear_response_A_matrix.f90 b/src/LR/phLR_A.f90 similarity index 59% rename from src/LR/linear_response_A_matrix.f90 rename to src/LR/phLR_A.f90 index efc6e81..8f9e8ee 100644 --- a/src/LR/linear_response_A_matrix.f90 +++ b/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) & - + 2d0*lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a) + 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) 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) & - - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a) + Aph(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + - (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) & - + lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a) + 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) 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 diff --git a/src/LR/linear_response_B_matrix.f90 b/src/LR/phLR_B.f90 similarity index 75% rename from src/LR/linear_response_B_matrix.f90 rename to src/LR/phLR_B.f90 index 5f15159..6e5258d 100644 --- a/src/LR/linear_response_B_matrix.f90 +++ b/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 diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 4eb5be5..d6c189d 100644 --- a/src/LR/ppLR.f90 +++ b/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 linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C) + 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 diff --git a/src/LR/linear_response_B_pp.f90 b/src/LR/ppLR_B.f90 similarity index 80% rename from src/LR/linear_response_B_pp.f90 rename to src/LR/ppLR_B.f90 index 9cc5860..0bff801 100644 --- a/src/LR/linear_response_B_pp.f90 +++ b/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 diff --git a/src/LR/linear_response_C_pp.f90 b/src/LR/ppLR_C.f90 similarity index 70% rename from src/LR/linear_response_C_pp.f90 rename to src/LR/ppLR_C.f90 index 452975e..56f686e 100644 --- a/src/LR/linear_response_C_pp.f90 +++ b/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) & - + 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) = + (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))) 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) & - + lambda*(ERI(a,b,c,d) - ERI(a,b,d,c)) + 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)) 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) & - + lambda*ERI(a,b,c,d) + Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,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 diff --git a/src/LR/linear_response_C_pp_od.f90 b/src/LR/ppLR_C_od.f90 similarity index 78% rename from src/LR/linear_response_C_pp_od.f90 rename to src/LR/ppLR_C_od.f90 index 82b22eb..48606d0 100644 --- a/src/LR/linear_response_C_pp_od.f90 +++ b/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 diff --git a/src/LR/linear_response_D_pp.f90 b/src/LR/ppLR_D.f90 similarity index 69% rename from src/LR/linear_response_D_pp.f90 rename to src/LR/ppLR_D.f90 index d6a799b..4d19ea5 100644 --- a/src/LR/linear_response_D_pp.f90 +++ b/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) & - + 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) = - (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))) 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) & - + lambda*(ERI(i,j,k,l) - ERI(i,j,l,k)) + 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)) 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) & - + lambda*ERI(i,j,k,l) + Dpp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,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 diff --git a/src/LR/linear_response_D_pp_od.f90 b/src/LR/ppLR_D_od.f90 similarity index 77% rename from src/LR/linear_response_D_pp_od.f90 rename to src/LR/ppLR_D_od.f90 index d9c902d..1d1b904 100644 --- a/src/LR/linear_response_D_pp_od.f90 +++ b/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 diff --git a/src/RPA/ACFDT_pp_correlation_energy.f90 b/src/RPA/ACFDT_pp_correlation_energy.f90 index 28991f8..92f8ec1 100644 --- a/src/RPA/ACFDT_pp_correlation_energy.f90 +++ b/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 linear_response_C_pp_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_B (ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,B) + call ppLR_C_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,C) + call ppLR_D_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,D) ! Compute Tr(K x P_lambda)