mirror of
https://github.com/pfloos/quack
synced 2024-12-23 04:43:42 +01:00
add S2
This commit is contained in:
parent
d738e335ba
commit
866ed94153
180
src/LR/S2_expval.f90
Normal file
180
src/LR/S2_expval.f90
Normal file
@ -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 <S**2> 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(:,:))
|
||||
|
||||
!-------------------------!
|
||||
! <S**2> 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)
|
||||
|
||||
!------------------------------------------!
|
||||
! <S**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
|
||||
|
||||
!------------------------------------------!
|
||||
! <S**2> 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
|
||||
Loading…
Reference in New Issue
Block a user