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

adding S2 in qsGGW and qsUGW

This commit is contained in:
Pierre-Francois Loos 2023-11-20 12:03:36 +01:00
parent dd164400fb
commit 9388a0d290
6 changed files with 92 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
src/GW/print_qsGGW.f90
View File

@ -1,4 +1,4 @@
subroutine print_qsGGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
subroutine print_qsGGW(nBas,nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,S,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
! Print information for the generalized version of qsGW
@ -8,6 +8,7 @@ subroutine print_qsGGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas2
integer,intent(in) :: nO
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
@ -19,21 +20,34 @@ subroutine print_qsGGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
double precision,intent(in) :: EcRPA
double precision,intent(in) :: Conv
double precision,intent(in) :: thresh
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: c(nBas)
double precision,intent(in) :: SigC(nBas,nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: eHF(nBas2)
double precision,intent(in) :: eGW(nBas2)
double precision,intent(in) :: c(nBas2,nBas2)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: SigC(nBas2,nBas2)
double precision,intent(in) :: Z(nBas2)
double precision,intent(in) :: EqsGW
double precision,intent(in) :: dipole(ncart)
! Local variables
logical :: dump_orb = .false.
integer :: i,j
integer :: p,ixyz,HOMO,LUMO
double precision :: Gap
double precision,external :: trace_matrix
double precision :: Sx,Sy,Sz
double precision :: SmSp,SpSm,Sz2,S2
double precision,allocatable :: Ca(:,:)
double precision,allocatable :: Cb(:,:)
double precision,allocatable :: Paa(:,:)
double precision,allocatable :: Pab(:,:)
double precision,allocatable :: Pba(:,:)
double precision,allocatable :: Pbb(:,:)
! Output variables
! HOMO and LUMO
@ -42,7 +56,59 @@ subroutine print_qsGGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
LUMO = HOMO + 1
Gap = eGW(LUMO)-eGW(HOMO)
! Compute energies
! Density matrices
allocate(Paa(nO,nO),Pab(nO,nO),Pba(nO,nO),Pbb(nO,nO))
allocate(Ca(nBas,nO),Cb(nBas,nO))
Ca(:,:) = C( 1:nBas ,1:nO)
Cb(:,:) = C(nBas+1:nBas2,1:nO)
Paa = matmul(transpose(Ca),matmul(S,Ca))
Pab = matmul(transpose(Ca),matmul(S,Cb))
Pba = matmul(transpose(Cb),matmul(S,Ca))
Pbb = matmul(transpose(Cb),matmul(S,Cb))
! Compute components of S = (Sx,Sy,Sz)
Sx = 0.5d0*(trace_matrix(nO,Pab) + trace_matrix(nO,Pba))
Sy = 0.5d0*(trace_matrix(nO,Pab) - trace_matrix(nO,Pba))
Sz = 0.5d0*(trace_matrix(nO,Paa) - trace_matrix(nO,Pbb))
! Compute <S^2> = <Sx^2> + <Sy^2> + <Sz^2>
SpSm = 0d0
do i=1,nO
do j=1,nO
SpSm = SpSm + Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i)
end do
end do
SpSm = trace_matrix(nO,Paa) + SpSm
SmSp = 0d0
do i=1,nO
do j=1,nO
SmSp = SmSp + Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i)
end do
end do
SmSp = trace_matrix(nO,Pbb) + SmSp
Sz2 = 0d0
do i=1,nO
do j=1,nO
Sz2 = Sz2 + (Paa(i,i) - Pbb(i,i))*(Paa(j,j) - Pbb(j,j)) - (Paa(i,j) - Pbb(i,j))**2
end do
end do
Sz2 = 0.25d0*(dble(nO) + Sz2)
! Compute <S^2> from Sz^2, S^+S^- and S^-S^+
S2 = Sz2 + 0.5d0*(SpSm + SmSp)
call print_GHF_spin(nBas,nBas2,nO,C,S)
! Dump results

27
src/GW/print_qsUGW.f90
View File

@ -1,4 +1,4 @@
subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,SigC,Z,dipole)
subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,Ov,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,SigC,Z,dipole)
! Print one-electron energies and other stuff for qsUGW
@ -22,7 +22,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov,ENuc,ET,EV,EJ,Ex,
double precision,intent(in) :: thresh
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: eGW(nBas,nspin)
double precision,intent(in) :: cGW(nBas,nBas,nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: Ov(nBas,nBas)
double precision,intent(in) :: SigC(nBas,nBas,nspin)
double precision,intent(in) :: Z(nBas,nspin)
@ -36,8 +36,8 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov,ENuc,ET,EV,EJ,Ex,
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
double precision :: S_exact,S2_exact
double precision :: S,S2
double precision :: Sz
double precision :: Sx2,Sy2,Sz2
double precision,external :: trace_matrix
! HOMO and LUMO
@ -54,11 +54,10 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov,ENuc,ET,EV,EJ,Ex,
end if
end do
S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0)
S2 = S2_exact + nO(2) - sum(matmul(transpose(cGW(:,1:nO(1),1)),matmul(Ov,cGW(:,1:nO(2),2)))**2)
S_exact = 0.5d0*dble(nO(1) - nO(2))
S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2)
Sz = 0.5d0*dble(nO(1) - nO(2))
Sx2 = 0.25d0*dble(nO(1) - nO(2)) + 0.5d0*nO(2) - 0.5d0*sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
Sy2 = 0.25d0*dble(nO(1) - nO(2)) + 0.5d0*nO(2) - 0.5d0*sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
Sz2 = 0.25d0*dble(nO(1) - nO(2))**2
! Dump results
@ -147,10 +146,8 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov,ENuc,ET,EV,EJ,Ex,
write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au'
write(*,'(A40,1X,F16.10,A3)') ' qsUGW energy = ',ENuc + EqsGW,' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,F13.6)') ' S (exact) = ',2d0*S_exact + 1d0
write(*,'(A40,F13.6)') ' S = ',2d0*S + 1d0
write(*,'(A40,F13.6)') ' <S**2> (exact) = ',S2_exact
write(*,'(A40,F13.6)') ' <S**2> = ',S2
write(*,'(A40,1X,F10.6)') ' <Sz> = ',Sz
write(*,'(A40,1X,F10.6)') ' <S^2> = ',Sx2+Sy2+Sz2
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A45)') ' Dipole moment (Debye) '
write(*,'(19X,4A10)') 'X','Y','Z','Tot.'
@ -164,12 +161,12 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov,ENuc,ET,EV,EJ,Ex,
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'qsUGW spin-up orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,cGW(:,:,1))
call matout(nBas,nBas,c(:,:,1))
write(*,*)
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'qsUGW spin-down orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,cGW(:,:,2))
call matout(nBas,nBas,c(:,:,2))
write(*,*)
end if

2
src/GW/qsGGW.f90
View File

@ -359,7 +359,7 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Print results
call dipole_moment(nBas2,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsGGW(nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole)
call print_qsGGW(nBas,nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,Ov,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole)
enddo
!------------------------------------------------------------------------

2
src/HF/GHF.f90
View File

@ -288,7 +288,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Compute final GHF energy
call print_GHF(nBas,nBas2,nO,eHF,C,P,Ov,ENuc,ET,EV,EJ,EK,EGHF,dipole)
call print_GHF(nBas,nBas2,nO,eHF,C,Ov,ENuc,ET,EV,EJ,EK,EGHF,dipole)
! Print test values

10
src/HF/print_GHF.f90
View File

@ -1,4 +1,4 @@
subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
subroutine print_GHF(nBas,nBas2,nO,eHF,C,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
! Print one-electron energies and other stuff for GHF
@ -14,7 +14,6 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
double precision,intent(in) :: eHF(nBas2)
double precision,intent(in) :: C(nBas2,nBas2)
double precision,intent(in) :: P(nBas2,nBas2)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET
@ -26,10 +25,11 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
! Local variables
logical :: dump_orb = .false.
integer :: i,j
integer :: ixyz
integer :: mu,nu
integer :: HOMO
integer :: LUMO
double precision :: Gap
@ -45,8 +45,6 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
double precision,external :: trace_matrix
logical :: dump_orb = .false.
! HOMO and LUMO
HOMO = nO
@ -103,7 +101,7 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
S2 = Sz2 + 0.5d0*(SpSm + SmSp)
call print_GHF_spin(nBas, nBas2, nO, C, S)
call print_GHF_spin(nBas,nBas2,nO,C,S)
! Dump results

8
src/HF/print_UHF.f90
View File

@ -27,7 +27,7 @@ subroutine print_UHF(nBas,nO,Ov,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
double precision :: Sx,Sy,Sz
double precision :: Sz
double precision :: Sx2,Sy2,Sz2
logical :: dump_orb = .false.
@ -51,7 +51,6 @@ subroutine print_UHF(nBas,nO,Ov,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
end do
Sz = 0.5d0*dble(nO(1) - nO(2))
! print*,Sz*(Sz+1d0) + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
Sx2 = 0.25d0*dble(nO(1) - nO(2)) + 0.5d0*nO(2) - 0.5d0*sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
Sy2 = 0.25d0*dble(nO(1) - nO(2)) + 0.5d0*nO(2) - 0.5d0*sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
Sz2 = 0.25d0*dble(nO(1) - nO(2))**2
@ -96,12 +95,7 @@ subroutine print_UHF(nBas,nO,Ov,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
write(*,'(A40,1X,F16.6,A3)') ' UHF LUMO b energy = ',LUMO(2)*HatoeV,' eV'
write(*,'(A40,1X,F16.6,A3)') ' UHF HOMOb-LUMOb gap = ',Gap(2)*HatoeV,' eV'
write(*,'(A60)') '---------------------------------------------'
write(*,'(A40,1X,F10.6)') ' <Sx> = ',Sx
write(*,'(A40,1X,F10.6)') ' <Sy> = ',Sy
write(*,'(A40,1X,F10.6)') ' <Sz> = ',Sz
write(*,'(A40,1X,F10.6)') ' <Sx^2> = ',Sx2
write(*,'(A40,1X,F10.6)') ' <Sy^2> = ',Sy2
write(*,'(A40,1X,F10.6)') ' <Sz^2> = ',Sz2
write(*,'(A40,1X,F10.6)') ' <S^2> = ',Sx2+Sy2+Sz2
write(*,'(A60)') '---------------------------------------------'
write(*,'(A45)') ' Dipole moment (Debye) '