diff --git a/src/LR/oscillator_strength_ph.f90 b/src/LR/oscillator_strength_ph.f90 index 25aa2d1..5591cc4 100644 --- a/src/LR/oscillator_strength_ph.f90 +++ b/src/LR/oscillator_strength_ph.f90 @@ -1,4 +1,4 @@ -subroutine oscillator_strength_ph(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY,XmY,os) +subroutine oscillator_strength_ph(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,XmY,os) ! Compute linear response @@ -15,7 +15,7 @@ subroutine oscillator_strength_ph(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY, integer,intent(in) :: nS integer,intent(in) :: maxS double precision :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) @@ -54,7 +54,7 @@ subroutine oscillator_strength_ph(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY, f(:,:) = sqrt(2d0)*f(:,:) do ia=1,maxS - os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) + os(ia) = 2d0/3d0*Om(ia)*sum(f(ia,:)**2) end do write(*,*) '---------------------------------------------------------------' diff --git a/src/LR/oscillator_strength_pp.f90 b/src/LR/oscillator_strength_pp.f90 index e29ead6..4906bb7 100644 --- a/src/LR/oscillator_strength_pp.f90 +++ b/src/LR/oscillator_strength_pp.f90 @@ -1,4 +1,4 @@ -subroutine oscillator_strength_pp(nBas,nC,nO,nV,nR,nOO,nVV,maxOO,maxVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2,os1,os2) +subroutine oscillator_strength_pp(nBas,nC,nO,nV,nR,nOO,nVV,maxOO,maxVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2,os1,os2) ! Compute linear response @@ -17,10 +17,10 @@ subroutine oscillator_strength_pp(nBas,nC,nO,nV,nR,nOO,nVV,maxOO,maxVV,dipole_in integer,intent(in) :: maxOO integer,intent(in) :: maxVV double precision,intent(in) :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: X1(nVV,nVV) double precision,intent(in) :: Y1(nOO,nVV) - double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: X2(nVV,nOO) double precision,intent(in) :: Y2(nOO,nOO) @@ -73,7 +73,7 @@ subroutine oscillator_strength_pp(nBas,nC,nO,nV,nR,nOO,nVV,maxOO,maxVV,dipole_in f1(:,:) = sqrt(2d0)*f1(:,:) do ab=1,maxVV - os1(ab) = +2d0/3d0*abs(Omega1(ab))*sum(f1(ab,:)**2) + os1(ab) = +2d0/3d0*abs(Om1(ab))*sum(f1(ab,:)**2) end do @@ -101,7 +101,7 @@ subroutine oscillator_strength_pp(nBas,nC,nO,nV,nR,nOO,nVV,maxOO,maxVV,dipole_in f2(:,:) = sqrt(2d0)*f2(:,:) do ij=1,maxOO - os2(ij) = 2d0/3d0*abs(Omega2(ij))*sum(f2(ij,:)**2) + os2(ij) = 2d0/3d0*abs(Om2(ij))*sum(f2(ij,:)**2) end do write(*,*) '---------------------------------------------------------------' diff --git a/src/LR/phLR.f90 b/src/LR/phLR.f90 index 9e24cf9..511d008 100644 --- a/src/LR/phLR.f90 +++ b/src/LR/phLR.f90 @@ -1,4 +1,4 @@ -subroutine phLR(TDA,nS,A,B,EcRPA,Om,XpY,XmY) +subroutine phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) ! Compute linear response @@ -9,13 +9,11 @@ subroutine phLR(TDA,nS,A,B,EcRPA,Om,XpY,XmY) logical,intent(in) :: TDA integer,intent(in) :: nS - double precision,intent(in) :: A(nS,nS) - double precision,intent(in) :: B(nS,nS) + double precision,intent(in) :: Aph(nS,nS) + double precision,intent(in) :: Bph(nS,nS) ! Local variables - integer :: i,j,k - double precision :: trace_matrix double precision,allocatable :: ApB(:,:) double precision,allocatable :: AmB(:,:) @@ -39,7 +37,7 @@ subroutine phLR(TDA,nS,A,B,EcRPA,Om,XpY,XmY) if(TDA) then - XpY(:,:) = A(:,:) + XpY(:,:) = Aph(:,:) call diagonalize_matrix(nS,XpY,Om) XpY(:,:) = transpose(XpY(:,:)) XmY(:,:) = XpY(:,:) @@ -48,8 +46,8 @@ subroutine phLR(TDA,nS,A,B,EcRPA,Om,XpY,XmY) ! Build A + B and A - B matrices - ApB = A + B - AmB = A - B + ApB(:,:) = Aph(:,:) + Bph(:,:) + AmB(:,:) = Aph(:,:) - Bph(:,:) ! Diagonalize linear response matrix @@ -81,6 +79,6 @@ subroutine phLR(TDA,nS,A,B,EcRPA,Om,XpY,XmY) ! Compute the RPA correlation energy - EcRPA = 0.5d0*(sum(Om) - trace_matrix(nS,A)) + EcRPA = 0.5d0*(sum(Om) - trace_matrix(nS,Aph)) end subroutine diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index b921c74..9ce816a 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -1,4 +1,4 @@ -subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA) +subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! Compute the pp channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013) @@ -10,9 +10,9 @@ subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA) logical,intent(in) :: TDA integer,intent(in) :: nOO integer,intent(in) :: nVV - double precision,intent(in) :: B(nVV,nOO) - double precision,intent(in) :: C(nVV,nVV) - double precision,intent(in) :: D(nOO,nOO) + double precision,intent(in) :: Bpp(nVV,nOO) + double precision,intent(in) :: Cpp(nVV,nVV) + double precision,intent(in) :: Dpp(nOO,nOO) ! Local variables @@ -49,25 +49,25 @@ subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA) if(TDA) then - X1(:,:) = +C(:,:) + X1(:,:) = +Cpp(:,:) Y1(:,:) = 0d0 if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1) X2(:,:) = 0d0 - Y2(:,:) = -D(:,:) + Y2(:,:) = -Dpp(:,:) if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2) else ! Diagonal blocks - M( 1:nVV , 1:nVV) = + C(1:nVV,1:nVV) - M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - D(1:nOO,1:nOO) + M( 1:nVV , 1:nVV) = + Cpp(1:nVV,1:nVV) + M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO) ! Off-diagonal blocks - M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO) - M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO)) + M( 1:nVV ,nVV+1:nOO+nVV) = - Bpp(1:nVV,1:nOO) + M(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp(1:nVV,1:nOO)) ! Diagonalize the p-p matrix @@ -81,9 +81,9 @@ subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! Compute the RPA correlation energy - EcRPA = 0.5d0*( sum(Om1) - sum(Om2) - trace_matrix(nVV,C) - trace_matrix(nOO,D) ) - EcRPA1 = +sum(Om1) - trace_matrix(nVV,C) - EcRPA2 = -sum(Om2) - trace_matrix(nOO,D) + EcRPA = 0.5d0*( sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp) ) + EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp) + EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp) if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' diff --git a/src/LR/print_transition_vectors_ph.f90 b/src/LR/print_transition_vectors_ph.f90 index 3723bc7..efb6a1e 100644 --- a/src/LR/print_transition_vectors_ph.f90 +++ b/src/LR/print_transition_vectors_ph.f90 @@ -1,4 +1,4 @@ -subroutine print_transition_vectors_ph(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY) +subroutine print_transition_vectors_ph(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) ! Print transition vectors for linear response calculation @@ -15,7 +15,7 @@ subroutine print_transition_vectors_ph(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_i integer,intent(in) :: nR integer,intent(in) :: nS double precision :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) @@ -37,7 +37,7 @@ subroutine print_transition_vectors_ph(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_i ! Compute oscillator strengths os(:) = 0d0 - if(spin_allowed) call oscillator_strength_ph(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY,XmY,os) + if(spin_allowed) call oscillator_strength_ph(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,XmY,os) ! Print details about excitations @@ -56,7 +56,7 @@ subroutine print_transition_vectors_ph(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_i print*,'-------------------------------------------------------------' write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & - ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2 + ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2 print*,'-------------------------------------------------------------' jb = 0 diff --git a/src/LR/print_transition_vectors_pp.f90 b/src/LR/print_transition_vectors_pp.f90 index a339dcd..dae8d0c 100644 --- a/src/LR/print_transition_vectors_pp.f90 +++ b/src/LR/print_transition_vectors_pp.f90 @@ -1,4 +1,4 @@ -subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) +subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) ! Print transition vectors for p-p linear response calculation @@ -16,10 +16,10 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip integer,intent(in) :: nOO integer,intent(in) :: nVV double precision,intent(in) :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: Om1(nVV) double precision,intent(in) :: X1(nVV,nVV) double precision,intent(in) :: Y1(nOO,nVV) - double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: Om2(nOO) double precision,intent(in) :: X2(nVV,nOO) double precision,intent(in) :: Y2(nOO,nOO) @@ -47,7 +47,7 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip os2(:) = 0d0 if(spin_allowed) & - call oscillator_strength_pp(nBas,nC,nO,nV,nR,nOO,nVV,maxOO,maxVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2,os1,os2) + call oscillator_strength_pp(nBas,nC,nO,nV,nR,nOO,nVV,maxOO,maxVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2,os1,os2) !-----------------------------------------------! ! Print details about excitations for pp sector ! @@ -70,7 +70,7 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip print*,'-------------------------------------------------------------' write(*,'(A20,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & - ' p-p excitation n. ',ab,': ',Omega1(ab)*HaToeV,' eV',' f = ',os1(ab),' = ',S2 + ' p-p excitation n. ',ab,': ',Om1(ab)*HaToeV,' eV',' f = ',os1(ab),' = ',S2 print*,'-------------------------------------------------------------' if(spin_allowed) then @@ -141,7 +141,7 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip print*,'-------------------------------------------------------------' write(*,'(A20,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & - ' h-h excitation n. ',ij,': ',Omega2(ij)*HaToeV,' eV',' f = ',os2(ij),' = ',S2 + ' h-h excitation n. ',ij,': ',Om2(ij)*HaToeV,' eV',' f = ',os2(ij),' = ',S2 print*,'-------------------------------------------------------------' if(spin_allowed) then