diff --git a/src/LR/S2_expval.f90 b/src/LR/S2_expval.f90 new file mode 100644 index 0000000..e8aed8f --- /dev/null +++ b/src/LR/S2_expval.f90 @@ -0,0 +1,180 @@ +subroutine S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2) + +! Compute for linear response excited states + + implicit none + include 'parameters.h' + +! Input variables + + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + integer,intent(in) :: nSa + integer,intent(in) :: nSb + integer,intent(in) :: nSt + integer,intent(in) :: maxS + double precision,intent(in) :: c(nBas,nBas,nspin) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: XpY(nSt,nSt) + double precision,intent(in) :: XmY(nSt,nSt) + +! Local variables + + integer :: m + integer :: ia,i,a + double precision :: S2_exact + double precision :: S2_gs + double precision,allocatable :: Xa(:,:), Xb(:,:), Ya(:,:), Yb(:,:) + double precision,allocatable :: Xat(:,:),Xbt(:,:),Yat(:,:),Ybt(:,:) + double precision,allocatable :: OO(:,:), OV(:,:), VO(:,:), VV(:,:) + double precision,allocatable :: OOt(:,:),OVt(:,:),VOt(:,:),VVt(:,:) + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: S2(maxS) + +! Memory allocation + + allocate(OO(nO(1)-nC(1),nO(2)-nC(2)), OV(nO(1)-nC(1),nV(2)-nR(2)), VO(nV(1)-nR(1),nO(2)-nC(2)), VV(nV(1)-nR(1),nV(2)-nR(2)), & + OOt(nO(2)-nC(2),nO(1)-nC(1)),OVt(nV(2)-nR(2),nO(1)-nC(1)),VOt(nO(2)-nC(2),nV(1)-nR(1)),VVt(nV(2)-nR(2),nV(1)-nR(1))) + +! Overlap matrix between spin-up and spin-down orbitals + + OO(:,:) = matmul(transpose(c(:,nC(1)+1:nO(1) ,1)),matmul(S,c(:,nC(2)+1:nO(2) ,2))) + OV(:,:) = matmul(transpose(c(:,nC(1)+1:nO(1) ,1)),matmul(S,c(:,nO(2)+1:nBas-nR(2),2))) + VO(:,:) = matmul(transpose(c(:,nO(1)+1:nBas-nR(1),1)),matmul(S,c(:,nC(2)+1:nO(2) ,2))) + VV(:,:) = matmul(transpose(c(:,nO(1)+1:nBas-nR(1),1)),matmul(S,c(:,nO(2)+1:nBas-nR(2),2))) + + OOt(:,:) = transpose(OO(:,:)) + OVt(:,:) = transpose(OV(:,:)) + VOt(:,:) = transpose(VO(:,:)) + VVt(:,:) = transpose(VV(:,:)) + +!-------------------------! +! for ground state ! +!-------------------------! + + S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0) + S2_gs = S2_exact + dble(nO(2)) - sum(OO(:,:)**2) + +!------------------------------------------! +! for spin-conserved-excited states ! +!------------------------------------------! + + if(ispin == 1) then + + allocate(Xa(nO(1)-nC(1),nV(1)-nR(1)), Ya(nO(1)-nC(1),nV(1)-nR(1)), Xb(nO(2)-nC(2),nV(2)-nR(2)), Yb(nO(2)-nC(2),nV(2)-nR(2)), & + Xat(nV(1)-nR(1),nO(1)-nC(1)),Yat(nV(1)-nR(1),nO(1)-nC(1)),Xbt(nV(2)-nR(2),nO(2)-nC(2)),Ybt(nV(2)-nR(2),nO(2)-nC(2))) + + do m=1,maxS + + ia = 0 + do i=nC(1)+1,nO(1) + do a=1,nV(1)-nR(1) + ia = ia + 1 + Xa(i,a) = 0.5d0*(XpY(m,ia) + XmY(m,ia)) + Ya(i,a) = 0.5d0*(XpY(m,ia) - XmY(m,ia)) + end do + end do + + ia = 0 + do i=nC(2)+1,nO(2) + do a=1,nV(2)-nR(2) + ia = ia + 1 + Xb(i,a) = 0.5d0*(XpY(m,nSa+ia) + XmY(m,nSa+ia)) + Yb(i,a) = 0.5d0*(XpY(m,nSa+ia) - XmY(m,nSa+ia)) + end do + end do + + Xat(:,:) = transpose(Xa(:,:)) + Xbt(:,:) = transpose(Xb(:,:)) + Yat(:,:) = transpose(Ya(:,:)) + Ybt(:,:) = transpose(Yb(:,:)) + + S2(m) = S2_gs & + + trace_matrix(nV(1),matmul(Xat,matmul(OO,matmul(OOt,Xa)))) & + + trace_matrix(nV(2),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) & + - trace_matrix(nO(1),matmul(Xa,matmul(VO,matmul(VOt,Xat)))) & + - trace_matrix(nO(2),matmul(Xb,matmul(OVt,matmul(OV,Xbt)))) & + - 2d0*trace_matrix(nO(1),matmul(OO,matmul(Xb,matmul(VVt,Xat)))) & + + - 2d0*trace_matrix(nV(2),matmul(OVt,matmul(Xa,matmul(VO,Yb)))) & + - 2d0*trace_matrix(nV(1),matmul(VO,matmul(Xb,matmul(OVt,Ya)))) & + + - trace_matrix(nV(1),matmul(Yat,matmul(OO,matmul(OOt,Ya)))) & + - trace_matrix(nV(2),matmul(Ybt,matmul(OOt,matmul(OO,Yb)))) & + + trace_matrix(nO(1),matmul(Ya,matmul(VO,matmul(VOt,Yat)))) & + + trace_matrix(nO(2),matmul(Yb,matmul(OVt,matmul(OV,Ybt)))) & + + 2d0*trace_matrix(nO(1),matmul(Ya,matmul(VV,matmul(Ybt,OOt)))) + + end do + + end if + +!------------------------------------------! +! for spin-conserved-excited states ! +!------------------------------------------! + + if(ispin == 2) then + + allocate(Xa(nO(1)-nC(1),nV(2)-nR(2)), Ya(nO(1)-nC(1),nV(2)-nR(2)), Xb(nO(2)-nC(2),nV(1)-nR(1)), Yb(nO(2)-nC(2),nV(1)-nR(1)), & + Xat(nV(2)-nR(2),nO(1)-nC(1)),Yat(nV(2)-nR(2),nO(1)-nC(1)),Xbt(nV(1)-nR(1),nO(2)-nC(2)),Ybt(nV(1)-nR(1),nO(2)-nC(2))) + + do m=1,maxS + + ia = 0 + do i=nC(1)+1,nO(1) + do a=1,nV(2)-nR(2) + ia = ia + 1 + Xa(i,a) = 0.5d0*(XpY(m,ia) + XmY(m,ia)) + Ya(i,a) = 0.5d0*(XpY(m,ia) - XmY(m,ia)) + end do + end do + + ia = 0 + do i=nC(2)+1,nO(2) + do a=1,nV(1)-nR(1) + ia = ia + 1 + Xb(i,a) = 0.5d0*(XpY(m,nSa+ia) + XmY(m,nSa+ia)) + Yb(i,a) = 0.5d0*(XpY(m,nSa+ia) - XmY(m,nSa+ia)) + end do + end do + + Xat(:,:) = transpose(Xa(:,:)) + Xbt(:,:) = transpose(Xb(:,:)) + Yat(:,:) = transpose(Ya(:,:)) + Ybt(:,:) = transpose(Yb(:,:)) + + S2(m) = S2_gs + dble(nO(2) - nO(1)) + 1d0 + + S2(m) = S2(m) & + + + trace_matrix(nV(1),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) & + - trace_matrix(nO(2),matmul(Xb,matmul(VO,matmul(VOt,Xbt)))) & + + trace_matrix(nO(2),matmul(Xb,VO))**2 & + + trace_matrix(nV(2),matmul(Yat,matmul(OO,matmul(OOt,Ya)))) & + + trace_matrix(nO(1),matmul(Ya,matmul(OVt,matmul(OV,Yat)))) & + + trace_matrix(nO(1),matmul(Ya,OVt))**2 & + - 2d0*trace_matrix(nO(2),matmul(Xb,VO))*trace_matrix(nO(1),matmul(Ya,OVt)) & + + + trace_matrix(nV(2),matmul(Xat,matmul(OO,matmul(OOt,Xa)))) & + - trace_matrix(nO(1),matmul(Xa,matmul(OVt,matmul(OV,Xat)))) & + + trace_matrix(nO(1),matmul(Xa,OVt))**2 & + + trace_matrix(nV(1),matmul(Ybt,matmul(OOt,matmul(OO,Yb)))) & + - trace_matrix(nO(2),matmul(Yb,matmul(VO,matmul(VOt,Ybt)))) & + + trace_matrix(nV(1),matmul(Ybt,VOt))**2 & + - 2d0*trace_matrix(nO(1),matmul(Xa,OVt))*trace_matrix(nO(2),matmul(Yb,VO)) + + end do + + end if + +end subroutine