diff --git a/src/GW/print_qsGGW.f90 b/src/GW/print_qsGGW.f90 index 0eb3405..0215ffe 100644 --- a/src/GW/print_qsGGW.f90 +++ b/src/GW/print_qsGGW.f90 @@ -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 = + + + + 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 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 diff --git a/src/GW/print_qsUGW.f90 b/src/GW/print_qsUGW.f90 index 8975144..dc0c61c 100644 --- a/src/GW/print_qsUGW.f90 +++ b/src/GW/print_qsUGW.f90 @@ -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)') ' (exact) = ',S2_exact - write(*,'(A40,F13.6)') ' = ',S2 + write(*,'(A40,1X,F10.6)') ' = ',Sz + write(*,'(A40,1X,F10.6)') ' = ',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 diff --git a/src/GW/qsGGW.f90 b/src/GW/qsGGW.f90 index 544dffb..cfe3fa7 100644 --- a/src/GW/qsGGW.f90 +++ b/src/GW/qsGGW.f90 @@ -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 !------------------------------------------------------------------------ diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90 index 56b8058..60fd9d2 100644 --- a/src/HF/GHF.f90 +++ b/src/HF/GHF.f90 @@ -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 diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index 77cf31f..4845bdc 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -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 diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index 27fe49a..1cd3ff7 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -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 - write(*,'(A40,1X,F10.6)') ' = ',Sy write(*,'(A40,1X,F10.6)') ' = ',Sz - write(*,'(A40,1X,F10.6)') ' = ',Sx2 - write(*,'(A40,1X,F10.6)') ' = ',Sy2 - write(*,'(A40,1X,F10.6)') ' = ',Sz2 write(*,'(A40,1X,F10.6)') ' = ',Sx2+Sy2+Sz2 write(*,'(A60)') '---------------------------------------------' write(*,'(A45)') ' Dipole moment (Debye) '