LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

rename ULR routines

This commit is contained in:
Pierre-Francois Loos 2023-07-21 13:04:29 +02:00
parent 5660146de5
commit c476b258ff
26 changed files with 322 additions and 344 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/CI/UCIS.f90
View File

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

98
src/GT/UG0T0.f90
View File

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

84
src/GT/evUGT.f90
View File

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

60
src/GT/qsUGT.f90
View File

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

22
src/GW/GW_ppBSE_static_kernel_B.f90
View File

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

22
src/GW/GW_ppBSE_static_kernel_C.f90
View File

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

22
src/GW/GW_ppBSE_static_kernel_D.f90
View File

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

14
src/GW/GW_self_energy.f90
View File

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

10
src/GW/GW_self_energy_diag.f90
View File

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

11
src/GW/SRG_qsGW.f90
View File

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

8
src/GW/UG0W0.f90
View File

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

4
src/GW/evUGW.f90
View File

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

4
src/GW/qsUGW.f90
View File

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

17
src/GW/unrestricted_Bethe_Salpeter.f90
View File

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

12
src/HF/UHF_stability.f90
View File

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

58
src/LR/unrestricted_linear_response.f90 → src/LR/phULR.f90
View File

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

21
src/LR/unrestricted_linear_response_A_matrix.f90 → src/LR/phULR_A.f90
View File

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

21
src/LR/unrestricted_linear_response_B_matrix.f90 → src/LR/phULR_B.f90
View File

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

39
src/LR/unrestricted_linear_response_pp.f90 → src/LR/ppULR.f90
View File

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

13
src/LR/unrestricted_linear_response_B_pp.f90 → src/LR/ppULR_B.f90
View File

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

13
src/LR/unrestricted_linear_response_C_pp.f90 → src/LR/ppULR_C.f90
View File

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

13
src/LR/unrestricted_linear_response_D_pp.f90 → src/LR/ppULR_D.f90
View File

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

8
src/RPA/URPA.f90
View File

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

8
src/RPA/URPAx.f90
View File

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

56
src/RPA/ppURPA.f90
View File

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

22
src/RPA/unrestricted_ACFDT.f90
View File

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