From 8ba6e368fae93cc87ef7138c81efe10f3ac7bf42 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Mon, 13 Nov 2023 09:13:58 +0100 Subject: [PATCH] added Sz, Sz^2 and S^2 for GHF --- src/HF/print_GHF.f90 | 133 ++++++++++++++++++++++++++++++++----------- 1 file changed, 100 insertions(+), 33 deletions(-) diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index 8fde002..3588fae 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -1,4 +1,4 @@ -subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) +subroutine print_GHF(nBas,nBas2,nO,e,Sao,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) ! Print one-electron energies and other stuff for GHF @@ -11,6 +11,9 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) integer,intent(in) :: nBas2 integer,intent(in) :: nO double precision,intent(in) :: e(nBas2) + ! TODO + ! add AO overlap as input + double precision,intent(in) :: Sao(nBas,nBas) double precision,intent(in) :: C(nBas2,nBas2) double precision,intent(in) :: P(nBas2,nBas2) double precision,intent(in) :: ENuc @@ -23,19 +26,22 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) ! Local variables - integer :: ixyz + integer :: i, j, ixyz integer :: mu,nu integer :: HOMO integer :: LUMO double precision :: Gap - double precision :: Sx2,Sy2,Sz2,S2 + double precision :: Sz,Sx2,Sy2,Sz2,S2 + double precision :: na, nb + double precision :: nonco_z, contam_uhf, xy_perp double precision,allocatable :: Ca(:,:) double precision,allocatable :: Cb(:,:) - double precision,allocatable :: Paa(:,:) - double precision,allocatable :: Pab(:,:) - double precision,allocatable :: Pba(:,:) - double precision,allocatable :: Pbb(:,:) + double precision,allocatable :: Paa(:,:), Saa(:,:) + double precision,allocatable :: Pab(:,:), Sab(:,:) + double precision,allocatable :: Pba(:,:), Sba(:,:) + double precision,allocatable :: Pbb(:,:), Sbb(:,:) + double precision,allocatable :: tmp(:,:) double precision,external :: trace_matrix @@ -56,36 +62,97 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) - allocate(Ca(nBas,nBas2),Cb(nBas,nBas2)) + ! TODO + ! check C size + allocate(Ca(nBas,nBas), Cb(nBas,nBas)) + do i = 1, nBas + do j = 1, nBas + Ca(j,i) = C(j, i) + Cb(j,i) = C(j,nBas+i) + enddo + enddo - Ca(:,:) = C( 1:nBas ,1:nBas2) - Cb(:,:) = C(nBas+1:nBas2,1:nBas2) + allocate(Saa(nBas,nBas),Sab(nBas,nBas),Sba(nBas,nBas),Sbb(nBas,nBas)) + allocate(tmp(nBas,nBas)) + + ! Saa = Ca x Sao x Ca.T + call dgemm("N", "N", nBas, nBas, nBas, 1.d0, Ca, size(Ca, 1), Sao, size(Sao, 1), 0.d0, tmp, size(tmp, 1)) + call dgemm("N", "T", nBas, nBas, nBas, 1.d0, tmp, size(tmp, 1), Ca, size(Ca, 1), 0.d0, Saa, size(Saa, 1)) + + ! Sab = Ca x Sao x Cb.T + call dgemm("N", "N", nBas, nBas, nBas, 1.d0, Ca, size(Ca, 1), Sao, size(Sao, 1), 0.d0, tmp, size(tmp, 1)) + call dgemm("N", "T", nBas, nBas, nBas, 1.d0, tmp, size(tmp, 1), Cb, size(Cb, 1), 0.d0, Sab, size(Sab, 1)) + + ! Sba = Cb x Sao x Ca.T + ! = Sab.T + Sba = transpose(Sab) + + ! Sbb = Cb x Sao x Cb.T + call dgemm("N", "N", nBas, nBas, nBas, 1.d0, Cb, size(Cb, 1), Sao, size(Sao, 1), 0.d0, tmp, size(tmp, 1)) + call dgemm("N", "T", nBas, nBas, nBas, 1.d0, tmp, size(tmp, 1), Cb, size(Cb, 1), 0.d0, Sbb, size(Sbb, 1)) + + deallocate(tmp) + + ! TODO + ! nO = nb of electrons ? + na = 0.d0 + nb = 0.d0 + do i = 1, nO + na = na + Saa(i,i) + nb = nb + Sbb(i,i) + enddo + + nonco_z = dble(nO) + do j = 1, nO + do i = 1, nO + nonco_z = nonco_z - (Saa(i,j) - Sbb(i,j))**2 + enddo + enddo + nonco_z = 0.25d0 * nonco_z + + Sz = 0.5d0 * (na - nb) + Sz2 = Sz*Sz + nonco_z + + contam_uhf = nb + do j = 1, nO + do i = 1, nO + contam_uhf = contam_uhf - (Sab(i,j) - Sba(j,i)) + enddo + enddo + + xy_perp = 0.d0 + do i = 1, nO + xy_perp = xy_perp + (Sba(i,i))**2 + enddo + + S2 = Sz * (Sz + 1.d0) + nonco_z + contam_uhf + xy_perp + ! Compute expectation values of S^2 (WRONG!) - Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2 - do mu=1,nBas - do nu=1,nBas - Sx2 = Sx2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) + Pab(mu,nu)*Pab(nu,mu)) - end do - end do - - Sy2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) - 0.25d0*trace_matrix(nBas,Pab+Pba)**2 - do mu=1,nBas - do nu=1,nBas - Sy2 = Sy2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu)) - end do - end do - - Sz2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab-Pba)**2 - do mu=1,nBas - do nu=1,nBas - Sz2 = Sz2 - 0.25d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu)) - Sz2 = Sz2 + 0.25d0*(Pab(mu,nu)*Pba(nu,mu) - Pba(mu,nu)*Pab(nu,mu)) - end do - end do - - S2 = Sx2 + Sy2 + Sz2 +! Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2 +! do mu=1,nBas +! do nu=1,nBas +! Sx2 = Sx2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) + Pab(mu,nu)*Pab(nu,mu)) +! end do +! end do +! +! Sy2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) - 0.25d0*trace_matrix(nBas,Pab+Pba)**2 +! do mu=1,nBas +! do nu=1,nBas +! Sy2 = Sy2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu)) +! end do +! end do +! +! Sz2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab-Pba)**2 +! do mu=1,nBas +! do nu=1,nBas +! Sz2 = Sz2 - 0.25d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu)) +! Sz2 = Sz2 + 0.25d0*(Pab(mu,nu)*Pba(nu,mu) - Pba(mu,nu)*Pab(nu,mu)) +! end do +! end do +! +! S2 = Sx2 + Sy2 + Sz2 ! Dump results