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

added Sz, Sz^2 and S^2 for GHF

Browse Source
This commit is contained in:
Abdallah Ammar 2023-11-13 09:13:58 +01:00
parent ae21a778c1
commit 8ba6e368fa
1 changed files with 100 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

133
src/HF/print_GHF.f90
View File

@ -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 ! 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) :: nBas2
integer,intent(in) :: nO integer,intent(in) :: nO
double precision,intent(in) :: e(nBas2) 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) :: C(nBas2,nBas2)
double precision,intent(in) :: P(nBas2,nBas2) double precision,intent(in) :: P(nBas2,nBas2)
double precision,intent(in) :: ENuc 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 ! Local variables
integer :: ixyz integer :: i, j, ixyz
integer :: mu,nu integer :: mu,nu
integer :: HOMO integer :: HOMO
integer :: LUMO integer :: LUMO
double precision :: Gap 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 :: Ca(:,:)
double precision,allocatable :: Cb(:,:) double precision,allocatable :: Cb(:,:)
double precision,allocatable :: Paa(:,:) double precision,allocatable :: Paa(:,:), Saa(:,:)
double precision,allocatable :: Pab(:,:) double precision,allocatable :: Pab(:,:), Sab(:,:)
double precision,allocatable :: Pba(:,:) double precision,allocatable :: Pba(:,:), Sba(:,:)
double precision,allocatable :: Pbb(:,:) double precision,allocatable :: Pbb(:,:), Sbb(:,:)
double precision,allocatable :: tmp(:,:)
double precision,external :: trace_matrix 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 ) Pba(:,:) = P(nBas+1:nBas2, 1:nBas )
Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) 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) allocate(Saa(nBas,nBas),Sab(nBas,nBas),Sba(nBas,nBas),Sbb(nBas,nBas))
Cb(:,:) = C(nBas+1:nBas2,1:nBas2) 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!) ! Compute expectation values of S^2 (WRONG!)
Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2 ! Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2
do mu=1,nBas ! do mu=1,nBas
do nu=1,nBas ! do nu=1,nBas
Sx2 = Sx2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) + Pab(mu,nu)*Pab(nu,mu)) ! Sx2 = Sx2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) + Pab(mu,nu)*Pab(nu,mu))
end do ! end do
end do ! end do
!
Sy2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) - 0.25d0*trace_matrix(nBas,Pab+Pba)**2 ! Sy2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) - 0.25d0*trace_matrix(nBas,Pab+Pba)**2
do mu=1,nBas ! do mu=1,nBas
do nu=1,nBas ! do nu=1,nBas
Sy2 = Sy2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu)) ! Sy2 = Sy2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu))
end do ! end do
end do ! end do
!
Sz2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab-Pba)**2 ! Sz2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab-Pba)**2
do mu=1,nBas ! do mu=1,nBas
do nu=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*(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)) ! Sz2 = Sz2 + 0.25d0*(Pab(mu,nu)*Pba(nu,mu) - Pba(mu,nu)*Pab(nu,mu))
end do ! end do
end do ! end do
!
S2 = Sx2 + Sy2 + Sz2 ! S2 = Sx2 + Sy2 + Sz2
! Dump results ! Dump results