PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:42 +01:00
Code Releases Activity

cleanup LR

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-20 22:01:11 +02:00
parent 6df119aaba
commit 752edf29bf
6 changed files with 38 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/LR/oscillator_strength_ph.f90
View File

@ -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(*,*) '---------------------------------------------------------------'

10
src/LR/oscillator_strength_pp.f90
View File

@ -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(*,*) '---------------------------------------------------------------'

16
src/LR/phLR.f90
View File

@ -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

26
src/LR/ppLR.f90
View File

@ -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 !!!'

8
src/LR/print_transition_vectors_ph.f90
View File

@ -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),' <S**2> = ',S2
' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2
print*,'-------------------------------------------------------------'
jb = 0

12
src/LR/print_transition_vectors_pp.f90
View File

@ -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),' <S**2> = ',S2
' p-p excitation n. ',ab,': ',Om1(ab)*HaToeV,' eV',' f = ',os1(ab),' <S**2> = ',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),' <S**2> = ',S2
' h-h excitation n. ',ij,': ',Om2(ij)*HaToeV,' eV',' f = ',os2(ij),' <S**2> = ',S2
print*,'-------------------------------------------------------------'
if(spin_allowed) then