LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-14 01:05:30 +02:00
Code Releases Activity

Merge branch 'abd' into master

Browse Source
This commit is contained in:
AbdAmmar 2023-11-16 15:52:06 +01:00 committed by GitHub
commit 3f7c1f3982
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 32 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
src/HF/print_GHF.f90
View File

@ -41,10 +41,10 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
double precision,allocatable :: Ca(:,:) double precision,allocatable :: Ca(:,:)
double precision,allocatable :: Cb(:,:) double precision,allocatable :: Cb(:,:)
double precision,allocatable :: Paa(:,:), Saa(:,:) double precision,allocatable :: Paa(:,:)
double precision,allocatable :: Pab(:,:), Sab(:,:) double precision,allocatable :: Pab(:,:)
double precision,allocatable :: Pba(:,:), Sba(:,:) double precision,allocatable :: Pba(:,:)
double precision,allocatable :: Pbb(:,:), Sbb(:,:) double precision,allocatable :: Pbb(:,:)
double precision,allocatable :: tmp(:,:) double precision,allocatable :: tmp(:,:)
double precision,allocatable :: Mx(:,:) double precision,allocatable :: Mx(:,:)
@ -80,6 +80,7 @@ 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))
@ -88,6 +89,7 @@ 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
@ -124,6 +126,7 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
- 0.25d0*trace_matrix(nO,matmul(Paa-Pbb,Paa-Pbb)) - 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 = Sz*(Sz+1d0) + trace_matrix(nO,Pbb) + 0.25d0*trace_matrix(nO,Paa+Pbb)
do i=1,nO do i=1,nO
do j=1,nO do j=1,nO
S2 = S2 - 0.25d0*(Paa(i,j) - Pbb(i,j))**2 & S2 = S2 - 0.25d0*(Paa(i,j) - Pbb(i,j))**2 &
@ -139,12 +142,9 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
! Sx2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Pab+Pba)**2 & ! 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)) ! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) + matmul(Pab,Pab))
! Sx2 = trace_matrix(
! Sy2 = 0.25d0*trace_matrix(nO,Paa+Pbb) - 0.25d0*trace_matrix(nO,Pab-Pba)**2 & ! 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)) ! - 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 & ! 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(Paa,Paa) + matmul(Pbb,Pbb)) &
! + 0.25d0*trace_matrix(nO,matmul(Pab,Pba) + matmul(Pba,Pab)) ! + 0.25d0*trace_matrix(nO,matmul(Pab,Pba) + matmul(Pba,Pab))
@ -159,77 +159,33 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
! end do ! end do
! print*,'<S^2> = ',S2 ! print*,'<S^2> = ',S2
! 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
! allocate(Saa(nBas,nBas),Sab(nBas,nBas),Sba(nBas,nBas),Sbb(nBas,nBas))
! allocate(tmp(nBas,nBas))
! Saa = Ca x Sao x Ca.T ! na = 0.d0
! call dgemm("N", "N", nBas, nBas, nBas, 1.d0, Ca, size(Ca, 1), Sao, size(Sao, 1), 0.d0, tmp, size(tmp, 1)) ! nb = 0.d0
! call dgemm("N", "T", nBas, nBas, nBas, 1.d0, tmp, size(tmp, 1), Ca, size(Ca, 1), 0.d0, Saa, size(Saa, 1)) ! do i = 1, nO
! na = na + Paa(i,i)
! Sab = Ca x Sao x Cb.T ! nb = nb + Pbb(i,i)
! call dgemm("N", "N", nBas, nBas, nBas, 1.d0, Ca, size(Ca, 1), Sao, size(Sao, 1), 0.d0, tmp, size(tmp, 1)) ! enddo
! call dgemm("N", "T", nBas, nBas, nBas, 1.d0, tmp, size(tmp, 1), Cb, size(Cb, 1), 0.d0, Sab, size(Sab, 1)) ! Sz = 0.5d0 * (na - nb)
!
! Sba = Cb x Sao x Ca.T ! nonco_z = dble(nO)
! = Sab.T ! do j = 1, nO
! Sba = transpose(Sab) ! do i = 1, nO
! nonco_z = nonco_z - (Paa(i,j) - Pbb(i,j))**2
! Sbb = Cb x Sao x Cb.T ! enddo
! call dgemm("N", "N", nBas, nBas, nBas, 1.d0, Cb, size(Cb, 1), Sao, size(Sao, 1), 0.d0, tmp, size(tmp, 1)) ! enddo
! call dgemm("N", "T", nBas, nBas, nBas, 1.d0, tmp, size(tmp, 1), Cb, size(Cb, 1), 0.d0, Sbb, size(Sbb, 1)) ! nonco_z = 0.25d0 * nonco_z
!
! deallocate(tmp) ! Sz2 = Sz*Sz + nonco_z
!
! TODO ! contam_ghf = 0.d0
! nO = nb of electrons ? ! do j = 1, nO
! na = 0.d0 ! do i = 1, nO
! nb = 0.d0 ! contam_ghf = contam_ghf - (Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i))
! do i = 1, nO ! enddo
! na = na + Saa(i,i) ! enddo
! nb = nb + Sbb(i,i) ! S2 = Sz * (Sz + 1.d0) + nonco_z + contam_ghf
! 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
! If Na > Nb
!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
! contam_ghf = 0.d0
! do j = 1, nO
! do i = 1, nO
! contam_ghf = contam_ghf - (Sab(i,i)*Sba(j,j) - Sab(i,j)*Sba(j,i))
! enddo
! enddo
! S2 = Sz * (Sz + 1.d0) + nonco_z + contam_ghf