diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index dfed554..b7c9432 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -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 :: Cb(:,:) - double precision,allocatable :: Paa(:,:), Saa(:,:) - double precision,allocatable :: Pab(:,:), Sab(:,:) - double precision,allocatable :: Pba(:,:), Sba(:,:) - double precision,allocatable :: Pbb(:,:), Sbb(:,:) + double precision,allocatable :: Paa(:,:) + double precision,allocatable :: Pab(:,:) + double precision,allocatable :: Pba(:,:) + double precision,allocatable :: Pbb(:,:) double precision,allocatable :: tmp(:,:) 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)) 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)) @@ -88,100 +89,51 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) ! Compute = + + -! 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)) + 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)) -! Sx2 = trace_matrix( + 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)) -! 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) -! 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)) + 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*,' = ',S2 -! 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*,' = ',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 - 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 - - ! 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 +! na = 0.d0 +! nb = 0.d0 +! do i = 1, nO +! na = na + Paa(i,i) +! nb = nb + Pbb(i,i) +! enddo +! Sz = 0.5d0 * (na - nb) +! +! nonco_z = dble(nO) +! do j = 1, nO +! do i = 1, nO +! nonco_z = nonco_z - (Paa(i,j) - Pbb(i,j))**2 +! enddo +! enddo +! nonco_z = 0.25d0 * nonco_z +! +! Sz2 = Sz*Sz + nonco_z +! +! contam_ghf = 0.d0 +! do j = 1, nO +! do i = 1, nO +! contam_ghf = contam_ghf - (Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i)) +! enddo +! enddo +! S2 = Sz * (Sz + 1.d0) + nonco_z + contam_ghf