From c476b258ffa29f8dcc31f53737e40e0e3cc2eb55 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 21 Jul 2023 13:04:29 +0200 Subject: [PATCH] rename ULR routines --- src/CI/UCIS.f90 | 6 +- src/GT/UG0T0.f90 | 98 +++++++++---------- src/GT/evUGT.f90 | 84 ++++++++-------- src/GT/qsUGT.f90 | 60 ++++++------ src/GW/GW_ppBSE_static_kernel_B.f90 | 22 ++--- src/GW/GW_ppBSE_static_kernel_C.f90 | 22 ++--- src/GW/GW_ppBSE_static_kernel_D.f90 | 22 ++--- src/GW/GW_self_energy.f90 | 14 +-- src/GW/GW_self_energy_diag.f90 | 10 +- src/GW/SRG_qsGW.f90 | 11 +-- src/GW/UG0W0.f90 | 8 +- src/GW/evUGW.f90 | 4 +- src/GW/qsUGW.f90 | 4 +- src/GW/unrestricted_Bethe_Salpeter.f90 | 17 ++-- src/HF/UHF_stability.f90 | 12 +-- ...stricted_linear_response.f90 => phULR.f90} | 58 ++++++----- ...near_response_A_matrix.f90 => phULR_A.f90} | 21 ++-- ...near_response_B_matrix.f90 => phULR_B.f90} | 21 ++-- ...icted_linear_response_pp.f90 => ppULR.f90} | 39 +++----- ...d_linear_response_B_pp.f90 => ppULR_B.f90} | 13 ++- ...d_linear_response_C_pp.f90 => ppULR_C.f90} | 13 ++- ...d_linear_response_D_pp.f90 => ppULR_D.f90} | 13 ++- src/RPA/URPA.f90 | 8 +- src/RPA/URPAx.f90 | 8 +- src/RPA/ppURPA.f90 | 56 +++++------ src/RPA/unrestricted_ACFDT.f90 | 22 ++--- 26 files changed, 322 insertions(+), 344 deletions(-) rename src/LR/{unrestricted_linear_response.f90 => phULR.f90} (63%) rename src/LR/{unrestricted_linear_response_A_matrix.f90 => phULR_A.f90} (81%) rename src/LR/{unrestricted_linear_response_B_matrix.f90 => phULR_B.f90} (79%) rename src/LR/{unrestricted_linear_response_pp.f90 => ppULR.f90} (64%) rename src/LR/{unrestricted_linear_response_B_pp.f90 => ppULR_B.f90} (81%) rename src/LR/{unrestricted_linear_response_C_pp.f90 => ppULR_C.f90} (82%) rename src/LR/{unrestricted_linear_response_D_pp.f90 => ppULR_D.f90} (81%) diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90 index 6c51ea2..386147a 100644 --- a/src/CI/UCIS.f90 +++ b/src/CI/UCIS.f90 @@ -68,8 +68,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E allocate(A_sc(nS_sc,nS_sc),Omega_sc(nS_sc)) - call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eHF, & - ERI_aaaa,ERI_aabb,ERI_bbbb,A_sc) + call phULR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,A_sc) if(dump_matrix) then print*,'CIS matrix (spin-conserved transitions)' @@ -109,8 +108,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf)) - call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,lambda,eHF, & - ERI_aaaa,ERI_aabb,ERI_bbbb,A_sf) + call phULR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,lambda,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,A_sf) if(dump_matrix) then print*,'CIS matrix (spin-conserved transitions)' diff --git a/src/GT/UG0T0.f90 b/src/GT/UG0T0.f90 index a79f411..5b1fee5 100644 --- a/src/GT/UG0T0.f90 +++ b/src/GT/UG0T0.f90 @@ -53,11 +53,11 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM - double precision,allocatable :: Omega1ab(:),Omega1aa(:),Omega1bb(:) + double precision,allocatable :: Om1ab(:),Om1aa(:),Om1bb(:) double precision,allocatable :: X1ab(:,:),X1aa(:,:),X1bb(:,:) double precision,allocatable :: Y1ab(:,:),Y1aa(:,:),Y1bb(:,:) double precision,allocatable :: rho1ab(:,:,:),rho1aa(:,:,:),rho1bb(:,:,:) - double precision,allocatable :: Omega2ab(:),Omega2aa(:),Omega2bb(:) + double precision,allocatable :: Om2ab(:),Om2aa(:),Om2bb(:) double precision,allocatable :: X2ab(:,:),X2aa(:,:),X2bb(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:) @@ -97,14 +97,14 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & ! Memory allocation - allocate(Omega1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & - Omega2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & + allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & + Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), & - Omega1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & - Omega2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & + Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & + Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & - Omega1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & - Omega2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & + Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & + Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & SigX(nBas,nspin),SigT(nBas,nspin),Z(nBas,nspin)) @@ -118,15 +118,15 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & ! Compute linear response - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPab,nHaa,nHab,nHbb,nHab,1d0,eHF,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & - Omega2ab,X2ab,Y2ab,EcRPA(ispin)) -call matout(nHab,nPab,Y1ab) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPab,nHaa,nHab,nHbb,nHab,1d0,eHF,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, & + Om2ab,X2ab,Y2ab,EcRPA(ispin)) + ! EcRPA(ispin) = 1d0*EcRPA(ispin) - call print_excitation('pp-RPA (N+2)',iblock,nPab,Omega1ab(:)) - call print_excitation('pp-RPA (N-2)',iblock,nHab,Omega2ab(:)) + call print_excitation('pp-RPA (N+2)',iblock,nPab,Om1ab(:)) + call print_excitation('pp-RPA (N-2)',iblock,nHab,Om2ab(:)) !---------------------------------------------- ! alpha-alpha block @@ -137,16 +137,16 @@ call matout(nHab,nPab,Y1ab) ! Compute linear response - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPaa,nHaa,nHab,nHbb,nHaa,1d0,eHF,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & - Omega2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPaa,nHaa,nHab,nHbb,nHaa,1d0,eHF,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, & + Om2aa,X2aa,Y2aa,EcRPA(ispin)) ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) - call print_excitation('pp-RPA (N+2)',iblock,nPaa,Omega1aa(:)) - call print_excitation('pp-RPA (N-2)',iblock,nHaa,Omega2aa(:)) + call print_excitation('pp-RPA (N+2)',iblock,nPaa,Om1aa(:)) + call print_excitation('pp-RPA (N-2)',iblock,nHaa,Om2aa(:)) !---------------------------------------------- ! beta-beta block @@ -157,16 +157,16 @@ call matout(nHab,nPab,Y1ab) ! Compute linear response - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPbb,nHaa,nHab,nHbb,nHbb,1d0,eHF,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & - Omega2bb,X2bb,Y2bb,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPbb,nHaa,nHab,nHbb,nHbb,1d0,eHF,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, & + Om2bb,X2bb,Y2bb,EcRPA(ispin)) ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) - call print_excitation('pp-RPA (N+2)',iblock,nPbb,Omega1bb(:)) - call print_excitation('pp-RPA (N-2)',iblock,nHbb,Omega2bb(:)) + call print_excitation('pp-RPA (N+2)',iblock,nPbb,Om1bb(:)) + call print_excitation('pp-RPA (N-2)',iblock,nHbb,Om2bb(:)) !---------------------------------------------- ! Compute T-matrix version of the self-energy @@ -200,14 +200,14 @@ call matout(nHab,nPab,Y1ab) rho1bb,X2bb,Y2bb,rho2bb) call unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& - nPab,nPbb,eHF,Omega1aa,Omega1ab,Omega1bb,& - rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,& - Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) + nPab,nPbb,eHF,Om1aa,Om1ab,Om1bb,& + rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,& + Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,& - nPaa,nPab,nPbb,eHF,Omega1aa,Omega1ab,& - Omega1bb,rho1aa,rho1ab,rho1bb, & - Omega2aa,Omega2ab,Omega2bb,rho2aa, & + nPaa,nPab,nPbb,eHF,Om1aa,Om1ab,& + Om1bb,rho1aa,rho1ab,rho1bb, & + Om2aa,Om2ab,Om2bb,rho2aa, & rho2ab,rho2bb,Z) @@ -251,20 +251,20 @@ call matout(nHab,nPab,Y1ab) ispin = 1 iblock = 3 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & - Omega2ab,X2ab,Y2ab,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, & + Om2ab,X2ab,Y2ab,EcRPA(ispin)) !alpha-alpha block ispin = 2 iblock = 4 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & - Omega2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, & + Om2aa,X2aa,Y2aa,EcRPA(ispin)) Ecaa = EcRPA(2) @@ -272,10 +272,10 @@ call matout(nHab,nPab,Y1ab) iblock = 7 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & - Omega2bb,X2bb,Y2bb,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, & + Om2bb,X2bb,Y2bb,EcRPA(ispin)) Ecbb = EcRPA(2) EcRPA(2) = Ecaa + Ecbb @@ -286,8 +286,8 @@ call matout(nHab,nPab,Y1ab) ! Free memory - deallocate(Omega1ab,X1ab,Y1ab,Omega2ab,X2ab,Y2ab,rho1ab,rho2ab, & - Omega1aa,X1aa,Y1aa,Omega2aa,X2aa,Y2aa,rho1aa,rho2aa, & - Omega1bb,X1bb,Y1bb,Omega2bb,X2bb,Y2bb,rho1bb,rho2bb) + deallocate(Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab, & + Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & + Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb) -end subroutine UG0T0 +end subroutine diff --git a/src/GT/evUGT.f90 b/src/GT/evUGT.f90 index 78efda5..bea6bcf 100644 --- a/src/GT/evUGT.f90 +++ b/src/GT/evUGT.f90 @@ -59,11 +59,11 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM - double precision,allocatable :: Omega1ab(:),Omega1aa(:),Omega1bb(:) + double precision,allocatable :: Om1ab(:),Om1aa(:),Om1bb(:) double precision,allocatable :: X1ab(:,:),X1aa(:,:),X1bb(:,:) double precision,allocatable :: Y1ab(:,:),Y1aa(:,:),Y1bb(:,:) double precision,allocatable :: rho1ab(:,:,:),rho1aa(:,:,:),rho1bb(:,:,:) - double precision,allocatable :: Omega2ab(:),Omega2aa(:),Omega2bb(:) + double precision,allocatable :: Om2ab(:),Om2aa(:),Om2bb(:) double precision,allocatable :: X2ab(:,:),X2aa(:,:),X2bb(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:) @@ -104,14 +104,14 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Memory allocation - allocate(Omega1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & - Omega2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & + allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & + Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), & - Omega1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & - Omega2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & + Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & + Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & - Omega1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & - Omega2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & + Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & + Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & SigX(nBas,nspin),SigT(nBas,nspin),Z(nBas,nspin), & eGT(nBas,nspin),eOld(nBas,nspin),error_diis(nBas,max_diis,nspin), & @@ -154,10 +154,10 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & - Omega2ab,X2ab,Y2ab,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, & + Om2ab,X2ab,Y2ab,EcRPA(ispin)) ! EcRPA(ispin) = 1d0*EcRPA(ispin) @@ -170,10 +170,10 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & - Omega2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, & + Om2aa,X2aa,Y2aa,EcRPA(ispin)) ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) @@ -187,10 +187,10 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & - Omega2bb,X2bb,Y2bb,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, & + Om2bb,X2bb,Y2bb,EcRPA(ispin)) ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) @@ -227,14 +227,14 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & rho1bb,X2bb,Y2bb,rho2bb) call unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& - nPab,nPbb,eGT,Omega1aa,Omega1ab,Omega1bb,& - rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,& - Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) + nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,& + rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,& + Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,& - nPaa,nPab,nPbb,eGT,Omega1aa,Omega1ab,& - Omega1bb,rho1aa,rho1ab,rho1bb, & - Omega2aa,Omega2ab,Omega2bb,rho2aa, & + nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,& + Om1bb,rho1aa,rho1ab,rho1bb, & + Om2aa,Om2ab,Om2bb,rho2aa, & rho2ab,rho2bb,Z) @@ -261,20 +261,20 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ispin = 1 iblock = 3 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & - Omega2ab,X2ab,Y2ab,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, & + Om2ab,X2ab,Y2ab,EcRPA(ispin)) !alpha-alpha block ispin = 2 iblock = 4 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & - Omega2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, & + Om2aa,X2aa,Y2aa,EcRPA(ispin)) Ecaa = EcRPA(2) @@ -282,10 +282,10 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & iblock = 7 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & - Omega2bb,X2bb,Y2bb,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, & + Om2bb,X2bb,Y2bb,EcRPA(ispin)) Ecbb = EcRPA(2) EcRPA(2) = Ecaa + Ecbb @@ -324,8 +324,8 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Free memory - deallocate(Omega1ab,X1ab,Y1ab,Omega2ab,X2ab,Y2ab,rho1ab,rho2ab, & - Omega1aa,X1aa,Y1aa,Omega2aa,X2aa,Y2aa,rho1aa,rho2aa, & - Omega1bb,X1bb,Y1bb,Omega2bb,X2bb,Y2bb,rho1bb,rho2bb) + deallocate(Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab, & + Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & + Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb) -end subroutine evUGT +end subroutine diff --git a/src/GT/qsUGT.f90 b/src/GT/qsUGT.f90 index 137325e..9508374 100644 --- a/src/GT/qsUGT.f90 +++ b/src/GT/qsUGT.f90 @@ -74,11 +74,11 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM(nspin) - double precision,allocatable :: Omega1ab(:),Omega1aa(:),Omega1bb(:) + double precision,allocatable :: Om1ab(:),Om1aa(:),Om1bb(:) double precision,allocatable :: X1ab(:,:),X1aa(:,:),X1bb(:,:) double precision,allocatable :: Y1ab(:,:),Y1aa(:,:),Y1bb(:,:) double precision,allocatable :: rho1ab(:,:,:),rho1aa(:,:,:),rho1bb(:,:,:) - double precision,allocatable :: Omega2ab(:),Omega2aa(:),Omega2bb(:) + double precision,allocatable :: Om2ab(:),Om2aa(:),Om2bb(:) double precision,allocatable :: X2ab(:,:),X2aa(:,:),X2bb(:,:) double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:) @@ -138,14 +138,14 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), & Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin)) - allocate(Omega1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & - Omega2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & + allocate(Om1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & + Om2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), & - Omega1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & - Omega2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & + Om1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & + Om2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & - Omega1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & - Omega2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & + Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & + Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb)) !Initialization @@ -208,10 +208,10 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & - Omega2ab,X2ab,Y2ab,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab, & + Om2ab,X2ab,Y2ab,EcRPA(ispin)) ! EcRPA(ispin) = 1d0*EcRPA(ispin) @@ -224,10 +224,10 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & - Omega2aa,X2aa,Y2aa,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa, & + Om2aa,X2aa,Y2aa,EcRPA(ispin)) ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) @@ -242,10 +242,10 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Compute linear response - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & - Omega2bb,X2bb,Y2bb,EcRPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb, & + Om2bb,X2bb,Y2bb,EcRPA(ispin)) ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) @@ -286,14 +286,14 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & rho1bb,X2bb,Y2bb,rho2bb) call unrestricted_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& - nPab,nPbb,eGT,Omega1aa,Omega1ab,Omega1bb,& - rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,& - Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) + nPab,nPbb,eGT,Om1aa,Om1ab,Om1bb,& + rho1aa,rho1ab,rho1bb,Om2aa,Om2ab,& + Om2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,& - nPaa,nPab,nPbb,eGT,Omega1aa,Omega1ab,& - Omega1bb,rho1aa,rho1ab,rho1bb, & - Omega2aa,Omega2ab,Omega2bb,rho2aa, & + nPaa,nPab,nPbb,eGT,Om1aa,Om1ab,& + Om1bb,rho1aa,rho1ab,rho1bb, & + Om2aa,Om2ab,Om2bb,rho2aa, & rho2ab,rho2bb,Z) @@ -433,10 +433,10 @@ write(*,*) 'EcGM', EcGM(1) ! Free memory - deallocate(Omega1ab,X1ab,Y1ab,Omega2ab,X2ab,Y2ab,rho1ab,rho2ab, & - Omega1aa,X1aa,Y1aa,Omega2aa,X2aa,Y2aa,rho1aa,rho2aa, & - Omega1bb,X1bb,Y1bb,Omega2bb,X2bb,Y2bb,rho1bb,rho2bb) + deallocate(Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab, & + Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & + Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb) deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,SigTm,Z,error,error_diis,F_diis) -end subroutine qsUGT +end subroutine diff --git a/src/GW/GW_ppBSE_static_kernel_B.f90 b/src/GW/GW_ppBSE_static_kernel_B.f90 index 654cf04..a334e58 100644 --- a/src/GW/GW_ppBSE_static_kernel_B.f90 +++ b/src/GW/GW_ppBSE_static_kernel_B.f90 @@ -1,4 +1,4 @@ -subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega,rho,WB) +subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,WB) ! Compute the VVOO block of the static screening W for the pp-BSE @@ -19,7 +19,7 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -54,9 +54,9 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda chi = 0d0 do m=1,nS - eps = Omega(m)**2 + eta**2 - chi = chi + rho(a,i,m)*rho(b,j,m)*Omega(m)/eps & - + rho(a,j,m)*rho(b,i,m)*Omega(m)/eps + eps = Om(m)**2 + eta**2 + chi = chi + rho(a,i,m)*rho(b,j,m)*Om(m)/eps & + + rho(a,j,m)*rho(b,i,m)*Om(m)/eps enddo WB(ab,ij) = + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) @@ -85,9 +85,9 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda chi = 0d0 do m=1,nS - eps = Omega(m)**2 + eta**2 - chi = chi + rho(a,i,m)*rho(b,j,m)*Omega(m)/eps & - - rho(a,j,m)*rho(b,i,m)*Omega(m)/eps + eps = Om(m)**2 + eta**2 + chi = chi + rho(a,i,m)*rho(b,j,m)*Om(m)/eps & + - rho(a,j,m)*rho(b,i,m)*Om(m)/eps enddo WB(ab,ij) = + 4d0*lambda*chi @@ -116,9 +116,9 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda chi = 0d0 do m=1,nS - eps = Omega(m)**2 + eta**2 - chi = chi + rho(a,i,m)*rho(b,j,m)*Omega(m)/eps & - - rho(a,j,m)*rho(b,i,m)*Omega(m)/eps + eps = Om(m)**2 + eta**2 + chi = chi + rho(a,i,m)*rho(b,j,m)*Om(m)/eps & + - rho(a,j,m)*rho(b,i,m)*Om(m)/eps enddo WB(ab,ij) = + 2d0*lambda*chi diff --git a/src/GW/GW_ppBSE_static_kernel_C.f90 b/src/GW/GW_ppBSE_static_kernel_C.f90 index 42cf8cc..b537cb2 100644 --- a/src/GW/GW_ppBSE_static_kernel_C.f90 +++ b/src/GW/GW_ppBSE_static_kernel_C.f90 @@ -1,4 +1,4 @@ -subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega,rho,WC) +subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,WC) ! Compute the VVVV block of the static screening W for the pp-BSE @@ -19,7 +19,7 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -54,9 +54,9 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda chi = 0d0 do m=1,nS - eps = Omega(m)**2 + eta**2 - chi = chi + rho(a,c,m)*rho(b,d,m)*Omega(m)/eps & - + rho(a,d,m)*rho(b,c,m)*Omega(m)/eps + eps = Om(m)**2 + eta**2 + chi = chi + rho(a,c,m)*rho(b,d,m)*Om(m)/eps & + + rho(a,d,m)*rho(b,c,m)*Om(m)/eps enddo WC(ab,cd) = + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) @@ -85,9 +85,9 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda chi = 0d0 do m=1,nS - eps = Omega(m)**2 + eta**2 - chi = chi + rho(a,c,m)*rho(b,d,m)*Omega(m)/eps & - - rho(a,d,m)*rho(b,c,m)*Omega(m)/eps + eps = Om(m)**2 + eta**2 + chi = chi + rho(a,c,m)*rho(b,d,m)*Om(m)/eps & + - rho(a,d,m)*rho(b,c,m)*Om(m)/eps enddo WC(ab,cd) = + 4d0*lambda*chi @@ -116,9 +116,9 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda chi = 0d0 do m=1,nS - eps = Omega(m)**2 + eta**2 - chi = chi + rho(a,c,m)*rho(b,d,m)*Omega(m)/eps & - - rho(a,d,m)*rho(b,c,m)*Omega(m)/eps + eps = Om(m)**2 + eta**2 + chi = chi + rho(a,c,m)*rho(b,d,m)*Om(m)/eps & + - rho(a,d,m)*rho(b,c,m)*Om(m)/eps enddo WC(ab,cd) = + 2d0*lambda*chi diff --git a/src/GW/GW_ppBSE_static_kernel_D.f90 b/src/GW/GW_ppBSE_static_kernel_D.f90 index e588881..5d3b9c2 100644 --- a/src/GW/GW_ppBSE_static_kernel_D.f90 +++ b/src/GW/GW_ppBSE_static_kernel_D.f90 @@ -1,4 +1,4 @@ -subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega,rho,WD) +subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,WD) ! Compute the OOOO block of the static screening W for the pp-BSE @@ -19,7 +19,7 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -54,9 +54,9 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda chi = 0d0 do m=1,nS - eps = Omega(m)**2 + eta**2 - chi = chi + rho(i,k,m)*rho(j,l,m)*Omega(m)/eps & - + rho(i,l,m)*rho(j,k,m)*Omega(m)/eps + eps = Om(m)**2 + eta**2 + chi = chi + rho(i,k,m)*rho(j,l,m)*Om(m)/eps & + + rho(i,l,m)*rho(j,k,m)*Om(m)/eps enddo WD(ij,kl) = + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) @@ -85,9 +85,9 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda chi = 0d0 do m=1,nS - eps = Omega(m)**2 + eta**2 - chi = chi + rho(i,k,m)*rho(j,l,m)*Omega(m)/eps & - - rho(i,l,m)*rho(j,k,m)*Omega(m)/eps + eps = Om(m)**2 + eta**2 + chi = chi + rho(i,k,m)*rho(j,l,m)*Om(m)/eps & + - rho(i,l,m)*rho(j,k,m)*Om(m)/eps enddo WD(ij,kl) = + 4d0*lambda*chi @@ -116,9 +116,9 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda chi = 0d0 do m=1,nS - eps = Omega(m)**2 + eta**2 - chi = chi + rho(i,k,m)*rho(j,l,m)*Omega(m)/eps & - - rho(i,l,m)*rho(j,k,m)*Omega(m)/eps + eps = Om(m)**2 + eta**2 + chi = chi + rho(i,k,m)*rho(j,l,m)*Om(m)/eps & + - rho(i,l,m)*rho(j,k,m)*Om(m)/eps enddo WD(ij,kl) = + 2d0*lambda*chi diff --git a/src/GW/GW_self_energy.f90 b/src/GW/GW_self_energy.f90 index 9e7f852..9838168 100644 --- a/src/GW/GW_self_energy.f90 +++ b/src/GW/GW_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) +subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) ! Compute correlation part of the self-energy and the renormalization factor @@ -15,7 +15,7 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -43,7 +43,7 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) ! Occupied part of the correlation self-energy !$OMP PARALLEL & -!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & +!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Om) & !$OMP PRIVATE(jb,i,q,p,eps,num) & !$OMP DEFAULT(NONE) !$OMP DO @@ -52,7 +52,7 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) do jb=1,nS do i=nC+1,nO - eps = e(p) - e(i) + Omega(jb) + eps = e(p) - e(i) + Om(jb) num = 2d0*rho(p,i,jb)*rho(q,i,jb) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -67,7 +67,7 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) ! Virtual part of the correlation self-energy !$OMP PARALLEL & -!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & +!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Om) & !$OMP PRIVATE(jb,a,q,p,eps,num) & !$OMP DEFAULT(NONE) !$OMP DO @@ -76,7 +76,7 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) do jb=1,nS do a=nO+1,nBas-nR - eps = e(p) - e(a) - Omega(jb) + eps = e(p) - e(a) - Om(jb) num = 2d0*rho(p,a,jb)*rho(q,a,jb) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -95,7 +95,7 @@ subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) do a=nO+1,nBas-nR do i=nC+1,nO - eps = e(a) - e(i) + Omega(jb) + eps = e(a) - e(i) + Om(jb) num = 4d0*rho(a,i,jb)*rho(a,i,jb) EcGM = EcGM - num*eps/(eps**2 + eta**2) diff --git a/src/GW/GW_self_energy_diag.f90 b/src/GW/GW_self_energy_diag.f90 index c8c0272..dbc02ab 100644 --- a/src/GW/GW_self_energy_diag.f90 +++ b/src/GW/GW_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) +subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) ! Compute diagonal of the correlation part of the self-energy and the renormalization factor @@ -15,7 +15,7 @@ subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -44,7 +44,7 @@ subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) do i=nC+1,nO do jb=1,nS - eps = e(p) - e(i) + Omega(jb) + eps = e(p) - e(i) + Om(jb) num = 2d0*rho(p,i,jb)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -59,7 +59,7 @@ subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) do a=nO+1,nBas-nR do jb=1,nS - eps = e(p) - e(a) - Omega(jb) + eps = e(p) - e(a) - Om(jb) num = 2d0*rho(p,a,jb)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 @@ -75,7 +75,7 @@ subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) do a=nO+1,nBas-nR do jb=1,nS - eps = e(a) - e(i) + Omega(jb) + eps = e(a) - e(i) + Om(jb) num = 4d0*rho(a,i,jb)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index df88aef..338b0c9 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -63,7 +63,6 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE double precision :: EqsGW double precision :: EcRPA double precision :: EcBSE(nspin) - double precision :: EcAC(nspin) double precision :: EcGM double precision :: Conv double precision :: rcond @@ -341,14 +340,14 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE end if - call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) + call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1) - write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 3862679..599a5e9 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -120,8 +120,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,sp ispin = 1 - call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) if(print_W) call print_excitation('RPA@UHF ',5,nS_sc,OmRPA) @@ -177,8 +177,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,sp ! Compute RPA correlation energy - call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) ! Dump results diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index 43fe783..199c2a1 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -146,8 +146,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE ! Compute screening - call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) !----------------------! ! Excitation densities ! diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index 1dfce97..9774122 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -205,8 +205,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE ! Compute linear response - call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) !----------------------! ! Excitation densities ! diff --git a/src/GW/unrestricted_Bethe_Salpeter.f90 b/src/GW/unrestricted_Bethe_Salpeter.f90 index cd50b1e..50b9801 100644 --- a/src/GW/unrestricted_Bethe_Salpeter.f90 +++ b/src/GW/unrestricted_Bethe_Salpeter.f90 @@ -79,8 +79,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_f ! Compute spin-conserved RPA screening - call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) @@ -98,9 +98,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_f ! Compute spin-conserved BSE excitation energies - call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin), & - OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) + call phULR(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin),OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) call print_excitation('BSE@UGW ',5,nS_sc,OmBSE_sc) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & cW,S,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) @@ -133,9 +132,9 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_f ! Compute spin-flip BSE excitation energies - call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin), & - OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) + call phULR(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin), & + OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) call print_excitation('BSE@UGW ',6,nS_sf,OmBSE_sf) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & @@ -154,4 +153,4 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_f end if -end subroutine unrestricted_Bethe_Salpeter +end subroutine diff --git a/src/HF/UHF_stability.f90 b/src/HF/UHF_stability.f90 index 8550855..d1d0ac4 100644 --- a/src/HF/UHF_stability.f90 +++ b/src/HF/UHF_stability.f90 @@ -51,10 +51,8 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) ispin = 1 - call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eHF, & - ERI_aaaa,ERI_aabb,ERI_bbbb,A_sc) - call unrestricted_linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & - ERI_aaaa,ERI_aabb,ERI_bbbb,B_sc) + call phULR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,A_sc) + call phULR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,B_sc) AB_sc(:,:) = A_sc(:,:) + B_sc(:,:) @@ -136,10 +134,8 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) ispin = 2 - call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,1d0,eHF, & - ERI_aaaa,ERI_aabb,ERI_bbbb,A_sf) - call unrestricted_linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,1d0, & - ERI_aaaa,ERI_aabb,ERI_bbbb,B_sf) + call phULR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,A_sf) + call phULR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,B_sf) AB_sf(:,:) = A_sf(:,:) + B_sf(:,:) diff --git a/src/LR/unrestricted_linear_response.f90 b/src/LR/phULR.f90 similarity index 63% rename from src/LR/unrestricted_linear_response.f90 rename to src/LR/phULR.f90 index b940f45..d3bcd82 100644 --- a/src/LR/unrestricted_linear_response.f90 +++ b/src/LR/phULR.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & - e,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,Omega,XpY,XmY) +subroutine phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,Om,XpY,XmY) ! Compute linear response for unrestricted formalism @@ -34,8 +34,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ! Local variables double precision,external :: trace_matrix - double precision,allocatable :: A(:,:) - double precision,allocatable :: B(:,:) + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: ApB(:,:) double precision,allocatable :: AmB(:,:) double precision,allocatable :: AmBSq(:,:) @@ -45,84 +45,82 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ! Output variables double precision,intent(out) :: EcRPA - double precision,intent(out) :: Omega(nSt) + double precision,intent(out) :: Om(nSt) double precision,intent(out) :: XpY(nSt,nSt) double precision,intent(out) :: XmY(nSt,nSt) ! Memory allocation - allocate(A(nSt,nSt),B(nSt,nSt),ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt)) + allocate(Aph(nSt,nSt),Bph(nSt,nSt),ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt)) ! Build A and B matrices - call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,A) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) if(BSE) & call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,A) + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Aph) ! Tamm-Dancoff approximation if(TDA) then - B(:,:) = 0d0 - XpY(:,:) = A(:,:) - call diagonalize_matrix(nSt,XpY,Omega) + Bph(:,:) = 0d0 + XpY(:,:) = Aph(:,:) + call diagonalize_matrix(nSt,XpY,Om) XpY(:,:) = transpose(XpY(:,:)) XmY(:,:) = XpY(:,:) else - call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,B) + call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) if(BSE) & call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,B) + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Bph) ! Build A + B and A - B matrices - ApB = A + B - AmB = A - B + ApB(:,:) = Aph(:,:) + Bph(:,:) + AmB(:,:) = Aph(:,:) - Bph(:,:) ! Diagonalize linear response matrix - call diagonalize_matrix(nSt,AmB,Omega) + call diagonalize_matrix(nSt,AmB,Om) - if(minval(Omega) < 0d0) & + if(minval(Om) < 0d0) & call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') ! do ia=1,nSt - ! if(Omega(ia) < 0d0) Omega(ia) = 0d0 + ! if(Om(ia) < 0d0) Om(ia) = 0d0 ! end do - call ADAt(nSt,AmB,1d0*sqrt(Omega),AmBSq) - call ADAt(nSt,AmB,1d0/sqrt(Omega),AmBIv) + call ADAt(nSt,AmB,1d0*sqrt(Om),AmBSq) + call ADAt(nSt,AmB,1d0/sqrt(Om),AmBIv) Z = matmul(AmBSq,matmul(ApB,AmBSq)) - call diagonalize_matrix(nSt,Z,Omega) + call diagonalize_matrix(nSt,Z,Om) - if(minval(Omega) < 0d0) & + if(minval(Om) < 0d0) & call print_warning('You may have instabilities in linear response: negative excitations!!') ! do ia=1,nSt - ! if(Omega(ia) < 0d0) Omega(ia) = 0d0 + ! if(Om(ia) < 0d0) Om(ia) = 0d0 ! end do - Omega = sqrt(Omega) + Om = sqrt(Om) XpY = matmul(transpose(Z),AmBSq) - call DA(nSt,1d0/sqrt(Omega),XpY) + call DA(nSt,1d0/sqrt(Om),XpY) XmY = matmul(transpose(Z),AmBIv) - call DA(nSt,1d0*sqrt(Omega),XmY) + call DA(nSt,1d0*sqrt(Om),XmY) end if ! Compute the RPA correlation energy - EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nSt,A)) + EcRPA = 0.5d0*(sum(Om) - trace_matrix(nSt,Aph)) -end subroutine unrestricted_linear_response +end subroutine diff --git a/src/LR/unrestricted_linear_response_A_matrix.f90 b/src/LR/phULR_A.f90 similarity index 81% rename from src/LR/unrestricted_linear_response_A_matrix.f90 rename to src/LR/phULR_A.f90 index 83cbf33..197f9f4 100644 --- a/src/LR/unrestricted_linear_response_A_matrix.f90 +++ b/src/LR/phULR_A.f90 @@ -1,5 +1,4 @@ -subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - e,ERI_aaaa,ERI_aabb,ERI_bbbb,A_lr) +subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) ! Compute linear response @@ -33,7 +32,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa ! Output variables - double precision,intent(out) :: A_lr(nSt,nSt) + double precision,intent(out) :: Aph(nSt,nSt) ! Direct RPA @@ -57,7 +56,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - A_lr(ia,jb) = (e(a,1) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + Aph(ia,jb) = (e(a,1) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + lambda*ERI_aaaa(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,b,j,a) end do @@ -76,7 +75,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - A_lr(ia,nSa+jb) = lambda*ERI_aabb(i,b,a,j) + Aph(ia,nSa+jb) = lambda*ERI_aabb(i,b,a,j) end do end do @@ -94,7 +93,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - A_lr(nSa+ia,jb) = lambda*ERI_aabb(b,i,j,a) + Aph(nSa+ia,jb) = lambda*ERI_aabb(b,i,j,a) end do end do @@ -112,7 +111,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - A_lr(nSa+ia,nSa+jb) = (e(a,2) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + Aph(nSa+ia,nSa+jb) = (e(a,2) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + lambda*ERI_bbbb(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,b,j,a) end do @@ -128,7 +127,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa if(ispin == 2) then - A_lr(:,:) = 0d0 + Aph(:,:) = 0d0 ! abab block @@ -140,7 +139,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do j=nC(1)+1,nO(1) do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + Aph(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - (1d0 - delta_dRPA)*lambda*ERI_aabb(i,b,j,a) end do end do @@ -158,7 +157,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + Aph(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - (1d0 - delta_dRPA)*lambda*ERI_aabb(b,i,a,j) end do @@ -169,4 +168,4 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end if -end subroutine unrestricted_linear_response_A_matrix +end subroutine diff --git a/src/LR/unrestricted_linear_response_B_matrix.f90 b/src/LR/phULR_B.f90 similarity index 79% rename from src/LR/unrestricted_linear_response_B_matrix.f90 rename to src/LR/phULR_B.f90 index fef8498..a7bc6a9 100644 --- a/src/LR/unrestricted_linear_response_B_matrix.f90 +++ b/src/LR/phULR_B.f90 @@ -1,5 +1,4 @@ -subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,B_lr) +subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) ! Compute linear response @@ -32,7 +31,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa ! Output variables - double precision,intent(out) :: B_lr(nSt,nSt) + double precision,intent(out) :: Bph(nSt,nSt) ! Direct RPA @@ -56,7 +55,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - B_lr(ia,jb) = lambda*ERI_aaaa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,j,b,a) + Bph(ia,jb) = lambda*ERI_aaaa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,j,b,a) end do end do @@ -74,7 +73,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - B_lr(ia,nSa+jb) = lambda*ERI_aabb(i,j,a,b) + Bph(ia,nSa+jb) = lambda*ERI_aabb(i,j,a,b) end do end do @@ -92,7 +91,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - B_lr(nSa+ia,jb) = lambda*ERI_aabb(j,i,b,a) + Bph(nSa+ia,jb) = lambda*ERI_aabb(j,i,b,a) end do end do @@ -110,7 +109,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - B_lr(nSa+ia,nSa+jb) = lambda*ERI_bbbb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,j,b,a) + Bph(nSa+ia,nSa+jb) = lambda*ERI_bbbb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,j,b,a) end do end do @@ -125,7 +124,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa if(ispin == 2) then - B_lr(:,:) = 0d0 + Bph(:,:) = 0d0 ! abba block @@ -138,7 +137,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - B_lr(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(i,j,b,a) + Bph(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(i,j,b,a) end do end do @@ -156,7 +155,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - B_lr(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(j,i,a,b) + Bph(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(j,i,a,b) end do end do @@ -166,4 +165,4 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end if -end subroutine unrestricted_linear_response_B_matrix +end subroutine diff --git a/src/LR/unrestricted_linear_response_pp.f90 b/src/LR/ppULR.f90 similarity index 64% rename from src/LR/unrestricted_linear_response_pp.f90 rename to src/LR/ppULR.f90 index 72aa75b..49afb79 100644 --- a/src/LR/unrestricted_linear_response_pp.f90 +++ b/src/LR/ppULR.f90 @@ -1,7 +1,5 @@ -subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,& - nPt,nHaa,nHab,nHbb,nHt,lambda,e,ERI_aaaa,& - ERI_aabb,ERI_bbbb,Omega1,X1,Y1,Omega2,X2,Y2,& - EcRPA) +subroutine ppULR(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,nHab,nHbb,nHt,lambda,e,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! Compute linear response for unrestricted formalism @@ -41,14 +39,14 @@ subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab, double precision,allocatable :: D(:,:) double precision,allocatable :: M(:,:) double precision,allocatable :: Z(:,:) - double precision,allocatable :: Omega(:) + double precision,allocatable :: Om(:) ! Output variables - double precision,intent(out) :: Omega1(nPt) + double precision,intent(out) :: Om1(nPt) double precision,intent(out) :: X1(nPt,nPt) double precision,intent(out) :: Y1(nHt,nPt) - double precision,intent(out) :: Omega2(nHt) + double precision,intent(out) :: Om2(nHt) double precision,intent(out) :: X2(nPt,nHt) double precision,intent(out) :: Y2(nHt,nHt) double precision,intent(out) :: EcRPA @@ -58,19 +56,14 @@ subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab, allocate(C(nPt,nPt),B(nPt,nHt),D(nHt,nHt),M(nPt+nHt,nPt+nHt),Z(nPt+nHt,nPt+nHt),& - Omega(nPt+nHt)) + Om(nPt+nHt)) ! Build C, B and D matrices for the pp channel - call unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,& - lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,C) - - call unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,& - nHaa,nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,& - ERI_bbbb,B) - - call unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,& - lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,D) + call ppULR_C(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,C) + call ppULR_B(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb, & + ERI_bbbb,B) + call ppULR_D(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,D) ! Diagonal blocks @@ -84,20 +77,20 @@ subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab, ! Diagonalize the p-h matrix - if(nHt+nPt > 0) call diagonalize_general_matrix(nHt+nPt,M,Omega,Z) + if(nHt+nPt > 0) call diagonalize_general_matrix(nHt+nPt,M,Om,Z) ! Split the various quantities in p-p and h-h parts - call sort_ppRPA(nHt,nPt,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),& + call sort_ppRPA(nHt,nPt,Om(:),Z(:,:),Om1(:),X1(:,:),Y1(:,:),Om2(:),X2(:,:),& Y2(:,:)) ! Compute the RPA correlation energy - EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nPt,C(:,:)) & + EcRPA = 0.5d0*( sum(Om1(:)) - sum(Om2(:)) - trace_matrix(nPt,C(:,:)) & - trace_matrix(nHt,D(:,:)) ) - EcRPA1 = +sum(Omega1(:)) - trace_matrix(nPt,C(:,:)) - EcRPA2 = -sum(Omega2(:)) - trace_matrix(nHt,D(:,:)) + EcRPA1 = +sum(Om1(:)) - trace_matrix(nPt,C(:,:)) + EcRPA2 = -sum(Om2(:)) - trace_matrix(nHt,D(:,:)) if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' -end subroutine unrestricted_linear_response_pp +end subroutine diff --git a/src/LR/unrestricted_linear_response_B_pp.f90 b/src/LR/ppULR_B.f90 similarity index 81% rename from src/LR/unrestricted_linear_response_B_pp.f90 rename to src/LR/ppULR_B.f90 index e1c57bc..6f489c4 100644 --- a/src/LR/unrestricted_linear_response_B_pp.f90 +++ b/src/LR/ppULR_B.f90 @@ -1,5 +1,4 @@ -subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa, & - nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,B_pp) +subroutine ppULR_B(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bpp) ! Compute linear response @@ -36,7 +35,7 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP ! Output variables - double precision,intent(out) :: B_pp(nPt,nHt) + double precision,intent(out) :: Bpp(nPt,nHt) eF = 0d0 @@ -54,7 +53,7 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP do j=nC(2)+1,nO(2) ij = ij + 1 - B_pp(ab,ij) = lambda*ERI_aabb(a,b,i,j) + Bpp(ab,ij) = lambda*ERI_aabb(a,b,i,j) end do end do @@ -77,7 +76,7 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP do j=i+1,nO(1) ij = ij + 1 - B_pp(ab,ij) = lambda*(ERI_aaaa(a,b,i,j) - ERI_aaaa(a,b,j,i)) + Bpp(ab,ij) = lambda*(ERI_aaaa(a,b,i,j) - ERI_aaaa(a,b,j,i)) end do end do @@ -98,7 +97,7 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP do j=i+1,nO(2) ij = ij + 1 - B_pp(ab,ij) = lambda*(ERI_bbbb(a,b,i,j) - ERI_bbbb(a,b,j,i)) + Bpp(ab,ij) = lambda*(ERI_bbbb(a,b,i,j) - ERI_bbbb(a,b,j,i)) end do end do @@ -107,4 +106,4 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP end if -end subroutine unrestricted_linear_response_B_pp +end subroutine diff --git a/src/LR/unrestricted_linear_response_C_pp.f90 b/src/LR/ppULR_C.f90 similarity index 82% rename from src/LR/unrestricted_linear_response_C_pp.f90 rename to src/LR/ppULR_C.f90 index dd14c72..01519bb 100644 --- a/src/LR/unrestricted_linear_response_C_pp.f90 +++ b/src/LR/ppULR_C.f90 @@ -1,5 +1,4 @@ -subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,& - lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,C_pp) +subroutine ppULR_C(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Cpp) ! Compute linear response @@ -33,7 +32,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP ! Output variables - double precision,intent(out) :: C_pp(nPt,nPt) + double precision,intent(out) :: Cpp(nPt,nPt) eF = 0d0 @@ -50,7 +49,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP do c=nO(1)+1,nBas-nR(1) do d=nO(2)+1,nBas-nR(2) cd = cd + 1 - C_pp(ab,cd) = (e(a,1) + e(b,2))*Kronecker_delta(a,c) & + Cpp(ab,cd) = (e(a,1) + e(b,2))*Kronecker_delta(a,c) & *Kronecker_delta(b,d) + lambda*ERI_aabb(a,b,c,d) end do end do @@ -72,7 +71,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP do d=c+1,nBas-nR(1) cd = cd + 1 - C_pp(ab,cd) = (e(a,1) + e(b,1) - eF)*Kronecker_delta(a,c)& + Cpp(ab,cd) = (e(a,1) + e(b,1) - eF)*Kronecker_delta(a,c)& *Kronecker_delta(b,d) + lambda*(ERI_aaaa(a,b,c,d) & - ERI_aaaa(a,b,d,c)) @@ -96,7 +95,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP do d=c+1,nBas-nR(2) cd = cd + 1 - C_pp(ab,cd) = (e(a,2) + e(b,2) - eF)*Kronecker_delta(a,c) & + Cpp(ab,cd) = (e(a,2) + e(b,2) - eF)*Kronecker_delta(a,c) & *Kronecker_delta(b,d) + lambda*(ERI_bbbb(a,b,c,d) & - ERI_bbbb(a,b,d,c)) @@ -107,4 +106,4 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP end if -end subroutine unrestricted_linear_response_C_pp +end subroutine diff --git a/src/LR/unrestricted_linear_response_D_pp.f90 b/src/LR/ppULR_D.f90 similarity index 81% rename from src/LR/unrestricted_linear_response_D_pp.f90 rename to src/LR/ppULR_D.f90 index 28c4c78..4ea0a3d 100644 --- a/src/LR/unrestricted_linear_response_D_pp.f90 +++ b/src/LR/ppULR_D.f90 @@ -1,5 +1,4 @@ -subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,& - lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,D_pp) +subroutine ppULR_D(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Dpp) ! Compute linear response @@ -33,7 +32,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH ! Output variables - double precision,intent(out) :: D_pp(nHt,nHt) + double precision,intent(out) :: Dpp(nHt,nHt) eF = 0d0 @@ -50,7 +49,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH do k=nC(1)+1,nO(1) do l=nC(2)+1,nO(2) kl = kl + 1 - D_pp(ij,kl) = -(e(i,1) + e(j,2))*Kronecker_delta(i,k)& + Dpp(ij,kl) = -(e(i,1) + e(j,2))*Kronecker_delta(i,k)& *Kronecker_delta(j,l) +lambda*ERI_aabb(i,j,k,l) end do end do @@ -72,7 +71,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH do l=k+1,nO(1) kl = kl + 1 - D_pp(ij,kl) = -(e(i,1) + e(j,1) - eF)*Kronecker_delta(i,k)& + Dpp(ij,kl) = -(e(i,1) + e(j,1) - eF)*Kronecker_delta(i,k)& *Kronecker_delta(j,l) + lambda*(ERI_aaaa(i,j,k,l) & - ERI_aaaa(i,j,l,k)) @@ -95,7 +94,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH do l=k+1,nO(2) kl = kl + 1 - D_pp(ij,kl) = -(e(i,2) + e(j,2) - eF)*Kronecker_delta(i,k) & + Dpp(ij,kl) = -(e(i,2) + e(j,2) - eF)*Kronecker_delta(i,k) & *Kronecker_delta(j,l) + lambda*(ERI_bbbb(i,j,k,l) & - ERI_bbbb(i,j,l,k)) @@ -106,4 +105,4 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH end if -end subroutine unrestricted_linear_response_D_pp +end subroutine diff --git a/src/RPA/URPA.f90 b/src/RPA/URPA.f90 index e08fb11..64874dd 100644 --- a/src/RPA/URPA.f90 +++ b/src/RPA/URPA.f90 @@ -84,8 +84,8 @@ subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) - call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) + call phULR(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPA ',5,nS_sc,Omega_sc) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & c,S,Omega_sc,XpY_sc,XmY_sc) @@ -108,8 +108,8 @@ subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) + call phULR(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPA ',6,nS_sf,Omega_sf) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & c,S,Omega_sf,XpY_sf,XmY_sf) diff --git a/src/RPA/URPAx.f90 b/src/RPA/URPAx.f90 index 450e913..5ef5b16 100644 --- a/src/RPA/URPAx.f90 +++ b/src/RPA/URPAx.f90 @@ -85,8 +85,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) - call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) + call phULR(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPAx ',5,nS_sc,Omega_sc) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & c,S,Omega_sc,XpY_sc,XmY_sc) @@ -109,8 +109,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf) + call phULR(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPAx ',6,nS_sf,Omega_sf) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & c,S,Omega_sf,XpY_sf,XmY_sf) diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90 index 752e98d..a8db313 100644 --- a/src/RPA/ppURPA.f90 +++ b/src/RPA/ppURPA.f90 @@ -28,10 +28,10 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH integer :: ispin,iblock integer :: nPaa,nPbb,nPab,nP_sc,nP_sf integer :: nHaa,nHbb,nHab,nH_sc,nH_sf - double precision,allocatable :: Omega1sc(:),Omega1sf(:) + double precision,allocatable :: Om1sc(:),Om1sf(:) double precision,allocatable :: X1sc(:,:),X1sf(:,:) double precision,allocatable :: Y1sc(:,:),Y1sf(:,:) - double precision,allocatable :: Omega2sc(:),Omega2sf(:) + double precision,allocatable :: Om2sc(:),Om2sf(:) double precision,allocatable :: X2sc(:,:),X2sf(:,:) double precision,allocatable :: Y2sc(:,:),Y2sf(:,:) @@ -64,16 +64,16 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ! Memory allocation - allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), & - Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc)) + allocate(Om1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), & + Om2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc)) - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc, & - Omega2sc,X2sc,Y2sc,Ec_ppURPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1sc,X1sc,Y1sc, & + Om2sc,X2sc,Y2sc,Ec_ppURPA(ispin)) - call print_excitation('pp-RPA (N+2)',iblock,nP_sc,Omega1sc) - call print_excitation('pp-RPA (N-2)',iblock,nH_sc,Omega2sc) + call print_excitation('pp-RPA (N+2)',iblock,nP_sc,Om1sc) + call print_excitation('pp-RPA (N-2)',iblock,nH_sc,Om2sc) !alpha-alpha block @@ -90,18 +90,18 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH nH_sf = nHaa - allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & - Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) + allocate(Om1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & + Om2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf, & - Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Om1sf,X1sf,Y1sf, & + Om2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) - call print_excitation('pp-RPA (N+2)',iblock,nP_sf,Omega1sf) - call print_excitation('pp-RPA (N-2)',iblock,nH_sf,Omega2sf) + call print_excitation('pp-RPA (N+2)',iblock,nP_sf,Om1sf) + call print_excitation('pp-RPA (N-2)',iblock,nH_sf,Om2sf) - deallocate(Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf) + deallocate(Om1sf,X1sf,Y1sf,Om2sf,X2sf,Y2sf) !beta-beta block @@ -110,16 +110,16 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH nP_sf = nPbb nH_sf = nHbb - allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & - Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) + allocate(Om1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & + Om2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,& - nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,& - ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,& - Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) + call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,& + nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,& + ERI_aabb,ERI_bbbb,Om1sf,X1sf,Y1sf,& + Om2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) - call print_excitation('pp-RPA (N+2)',iblock,nP_sf,Omega1sf) - call print_excitation('pp-RPA (N-2)',iblock,nH_sf,Omega2sf) + call print_excitation('pp-RPA (N+2)',iblock,nP_sf,Om1sf) + call print_excitation('pp-RPA (N-2)',iblock,nH_sf,Om2sf) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -152,4 +152,4 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ! end if -end subroutine ppURPA +end subroutine diff --git a/src/RPA/unrestricted_ACFDT.f90 b/src/RPA/unrestricted_ACFDT.f90 index bdd9435..80093d2 100644 --- a/src/RPA/unrestricted_ACFDT.f90 +++ b/src/RPA/unrestricted_ACFDT.f90 @@ -93,8 +93,8 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) - call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) ! Spin-conserved manifold @@ -120,14 +120,14 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons if(doXBS) then - call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) end if - call unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Omega_sc,XpY_sc,XmY_sc) + call phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Omega_sc,XpY_sc,XmY_sc) call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, & ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,Ec(iAC,ispin)) @@ -174,14 +174,14 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons if(doXBS) then - call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) end if - call unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Omega_sf,XpY_sf,XmY_sf) + call phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Omega_sf,XpY_sf,XmY_sf) call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, & ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,Ec(iAC,ispin)) @@ -203,4 +203,4 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons end if -end subroutine unrestricted_ACFDT +end subroutine