LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00
Code Releases Activity

clean up S^2 GHF

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-11-16 19:38:32 +01:00
parent efd001b7cf
commit 0aabc1a5f4
1 changed files with 9 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

57
src/HF/print_GHF.f90
View File

@ -80,7 +80,6 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
Pba = matmul(transpose(Cb),matmul(S,Ca)) Pba = matmul(transpose(Cb),matmul(S,Ca))
Pbb = matmul(transpose(Cb),matmul(S,Cb)) Pbb = matmul(transpose(Cb),matmul(S,Cb))
! Compute components of S = (Sx,Sy,Sz) ! Compute components of S = (Sx,Sy,Sz)
Sx = 0.5d0*(trace_matrix(nO,Pab) + trace_matrix(nO,Pba)) Sx = 0.5d0*(trace_matrix(nO,Pab) + trace_matrix(nO,Pba))
@ -89,17 +88,14 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
! Compute <S^2> = <Sx^2> + <Sy^2> + <Sz^2> ! Compute <S^2> = <Sx^2> + <Sy^2> + <Sz^2>
SpSm = 0d0 SpSm = 0d0
do i=1,nO do i=1,nO
do j=1,nO do j=1,nO
SpSm = SpSm + Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i) SpSm = SpSm + Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i)
end do end do
end do end do
SpSm = 0.5d0*(trace_matrix(nO,Paa) + SpSm) SpSm = trace_matrix(nO,Paa) + SpSm
print*,'<S^+S^-> = ',SpSm
! Sx2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Pab+Pba)**2 &
! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) + matmul(Pab,Pab))
SmSp = 0d0 SmSp = 0d0
do i=1,nO do i=1,nO
@ -107,10 +103,8 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
SmSp = SmSp + Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i) SmSp = SmSp + Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i)
end do end do
end do end do
SmSp = 0.5d0*(trace_matrix(nO,Pbb) + SmSp) SmSp = trace_matrix(nO,Pbb) + SmSp
print*,'<S^-S^+> = ',SmSp
! Sy2 = 0.25d0*trace_matrix(nO,Paa+Pbb) - 0.25d0*trace_matrix(nO,Pab-Pba)**2 &
! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) - matmul(Pab,Pab))
Sz2 = 0d0 Sz2 = 0d0
do i=1,nO do i=1,nO
@ -121,45 +115,17 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
Sz2 = 0.25d0*(dble(nO) + Sz2) Sz2 = 0.25d0*(dble(nO) + Sz2)
print*,'<Sz^2> = ',Sz2 print*,'<Sz^2> = ',Sz2
Sz2 = 0.25d0*trace_matrix(nO,Paa+Pbb) & ! Compute <S^2> from Sz^2, S^+S^- and S^-S^+
- 0.25d0*trace_matrix(nO,Paa-Pbb)**2 &
- 0.25d0*trace_matrix(nO,matmul(Paa-Pbb,Paa-Pbb))
S2 = Sz*(Sz+1d0) + trace_matrix(nO,Pbb) + 0.25d0*trace_matrix(nO,Paa+Pbb) S2 = Sz2 + 0.5d0*(SpSm + SmSp)
do i=1,nO
do j=1,nO
S2 = S2 - 0.25d0*(Paa(i,j) - Pbb(i,j))**2 &
+ (Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i))
end do
end do
print*,'<S^2> = ',S2 print*,'<S^2> = ',S2
! Compute <Sx^2> and <Sy^2> from <S^2>, <Sz^2>, <S^+S^-> and <S^-S^+>
Sx2 = 0.5d0*(S2 - Sz2 + 0.5d0*(SmSp + SpSm)) Sx2 = 0.5d0*(S2 - Sz2 + 0.5d0*(SmSp + SpSm))
print*,'<Sx^2> = ',Sx2 print*,'<Sx^2> = ',Sx2
Sy2 = 0.5d0*(S2 - Sz2 - 0.5d0*(SmSp + SpSm)) Sy2 = 0.5d0*(S2 - Sz2 - 0.5d0*(SmSp + SpSm))
print*,'<Sy^2> = ',Sy2 print*,'<Sy^2> = ',Sy2
! Sx2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Pab+Pba)**2 &
! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) + matmul(Pab,Pab))
! Sy2 = 0.25d0*trace_matrix(nO,Paa+Pbb) - 0.25d0*trace_matrix(nO,Pab-Pba)**2 &
! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) - matmul(Pab,Pab))
! Sz2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Paa-Pbb)**2 &
! - 0.25d0*trace_matrix(nO,matmul(Paa,Paa) + matmul(Pbb,Pbb)) &
! + 0.25d0*trace_matrix(nO,matmul(Pab,Pba) + matmul(Pba,Pab))
! S2 = Sz*(Sz+1d0) + trace_matrix(nO,Pbb) + 0.25d0*trace_matrix(nO,Paa+Pbb)
! do i=1,nO
! do j=1,nO
! S2 = S2 - 0.25d0*(Paa(i,j) - Pbb(i,j))**2 &
! + (Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i))
! end do
! end do
! print*,'<S^2> = ',S2
! na = 0.d0 ! na = 0.d0
! nb = 0.d0 ! nb = 0.d0
@ -187,9 +153,6 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
! enddo ! enddo
! S2 = Sz * (Sz + 1.d0) + nonco_z + contam_ghf ! S2 = Sz * (Sz + 1.d0) + nonco_z + contam_ghf
! deallocate(Paa,Pab,Pba,Pbb) ! deallocate(Paa,Pab,Pba,Pbb)
! Check collinearity and coplanarity ! Check collinearity and coplanarity
@ -255,10 +218,8 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
write(*,'(A33,1X,F16.6)') ' <Sy> = ',Sy write(*,'(A33,1X,F16.6)') ' <Sy> = ',Sy
write(*,'(A33,1X,F16.6)') ' <Sz> = ',Sz write(*,'(A33,1X,F16.6)') ' <Sz> = ',Sz
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.6)') ' <Sx**2> = ',Sx2 write(*,'(A33,1X,F16.6)') ' <Sx^2+Sy^2> = ',S2 - Sz2
write(*,'(A33,1X,F16.6)') ' <Sy**2> = ',Sy2
write(*,'(A33,1X,F16.6)') ' <Sz**2> = ',Sz2 write(*,'(A33,1X,F16.6)') ' <Sz**2> = ',Sz2
write(*,'(A33,1X,F16.6)') ' <S**2> = ',Sx2+Sy2+Sz2
write(*,'(A33,1X,F16.6)') ' <S**2> = ',S2 write(*,'(A33,1X,F16.6)') ' <S**2> = ',S2
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A36)') ' Dipole moment (Debye) ' write(*,'(A36)') ' Dipole moment (Debye) '