mirror of
https://github.com/pfloos/quack
synced 2025-03-09 18:22:25 +01:00
Merge branch 'master' of https://github.com/pfloos/QuAcK
This commit is contained in:
commit
7488f0c845
@ -28,10 +28,22 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
double precision,allocatable :: Wovvo(:,:,:,:)
|
||||
double precision,allocatable :: H(:,:)
|
||||
double precision,allocatable :: Om(:)
|
||||
double precision,allocatable :: Z(:,:)
|
||||
double precision,allocatable :: VL(:,:)
|
||||
double precision,allocatable :: VR(:,:)
|
||||
double precision,allocatable :: Leom(:,:,:)
|
||||
double precision,allocatable :: Reom(:,:,:)
|
||||
|
||||
integer :: nstate,m
|
||||
double precision :: Ex,tmp
|
||||
|
||||
integer,allocatable :: order(:)
|
||||
|
||||
double precision,allocatable :: rdm1_oo(:,:)
|
||||
double precision,allocatable :: rdm1_vv(:,:)
|
||||
|
||||
double precision,allocatable :: rdm2_oovv(:,:,:,:)
|
||||
double precision,allocatable :: rdm2_ovvo(:,:,:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
@ -46,10 +58,9 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS),Z(nS,nS))
|
||||
allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS))
|
||||
allocate(order(nS))
|
||||
|
||||
|
||||
! Form one-body terms
|
||||
|
||||
do a=1,nV-nR
|
||||
@ -132,19 +143,142 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
|
||||
! Diagonalize EOM Hamiltonian
|
||||
|
||||
allocate(VL(nS,nS),VR(nS,nS))
|
||||
|
||||
if(nS > 0) then
|
||||
|
||||
call diagonalize_general_matrix(nS,H,Om,Z)
|
||||
call diagonalize_general_matrix_LR(nS,H,Om,VL,VR)
|
||||
|
||||
do ia=1,nS
|
||||
order(ia) = ia
|
||||
end do
|
||||
|
||||
call quick_sort(Om,order,nS)
|
||||
call set_order(Z,order,nS,nS)
|
||||
call set_order_LR(VL,VR,order,nS,nS)
|
||||
|
||||
call print_excitation_energies('EE-EOM-CCD','spinorbital',nS,Om)
|
||||
|
||||
! write(*,*) 'Right Eigenvectors'
|
||||
! call matout(nS,nS,VR)
|
||||
|
||||
! call matout(nS,3,VR(:,1:3))
|
||||
|
||||
end if
|
||||
|
||||
allocate(Leom(nO,nV,nS),Reom(nO,nV,nS))
|
||||
|
||||
do m=1,nS
|
||||
ia = 0
|
||||
do i=1,nO
|
||||
do a=1,nV
|
||||
ia = ia + 1
|
||||
Leom(i,a,m) = VL(ia,m)
|
||||
Reom(i,a,m) = VR(ia,m)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
deallocate(VL,VR)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! EOM section
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
allocate(rdm1_oo(nO,nO),rdm1_vv(nV,nV))
|
||||
allocate(rdm2_oovv(nO,nO,nV,nV),rdm2_ovvo(nO,nV,nV,nO))
|
||||
|
||||
nstate = 1
|
||||
|
||||
tmp = 0d0
|
||||
do i=1,nO
|
||||
do a=1,nV
|
||||
tmp = tmp + Leom(i,a,nstate)*Reom(i,a,nstate)
|
||||
end do
|
||||
end do
|
||||
print*,tmp
|
||||
|
||||
rdm1_oo(:,:) = 0d0
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do c=1,nV
|
||||
|
||||
rdm1_oo(i,j) = rdm1_oo(i,j) - Reom(i,c,nstate)*Leom(j,c,nstate)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
rdm1_vv(:,:) = 0d0
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
do k=1,nO
|
||||
|
||||
rdm1_vv(a,b) = rdm1_vv(a,b) + Reom(k,b,nstate)*Leom(k,a,nstate)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
rdm2_ovvo(:,:,:,:) = 0d0
|
||||
do i=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
do j=1,nO
|
||||
|
||||
rdm2_ovvo(i,a,b,j) = Reom(i,b,nstate)*Leom(j,a,nstate)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
rdm2_oovv(:,:,:,:) = 0d0
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
|
||||
do k=1,nO
|
||||
do c=1,nV
|
||||
|
||||
rdm2_oovv(i,j,a,b) = rdm2_oovv(i,j,a,b) &
|
||||
+ Reom(j,b,nstate)*t(k,i,c,a)*Leom(k,c,nstate) &
|
||||
- Reom(i,b,nstate)*t(k,j,c,a)*Leom(k,c,nstate) &
|
||||
- Reom(j,a,nstate)*t(k,i,c,b)*Leom(k,c,nstate) &
|
||||
+ Reom(i,a,nstate)*t(k,j,c,b)*Leom(k,c,nstate)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
Ex = 0d0
|
||||
|
||||
do i=1,nO
|
||||
Ex = Ex + rdm1_oo(i,i)*eO(i)
|
||||
end do
|
||||
|
||||
do a=1,nV
|
||||
Ex = Ex + rdm1_vv(a,a)*eV(a)
|
||||
end do
|
||||
|
||||
do i=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
do j=1,nO
|
||||
|
||||
Ex = Ex + rdm2_ovvo(i,a,b,j)*OVVO(i,a,b,j) + 0.25d0*rdm2_oovv(i,j,a,b)*OOVV(i,j,a,b)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
print*,'Ex = ',Ex
|
||||
print*,'Om = ',Om(nstate)
|
||||
|
||||
|
||||
end subroutine
|
||||
|
||||
@ -220,8 +220,7 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu
|
||||
|
||||
! EE-EOM-CCD (1h1p)
|
||||
|
||||
! if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
|
||||
if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
|
||||
if(dotest) then
|
||||
|
||||
|
||||
@ -60,6 +60,38 @@ recursive subroutine rec_quicksort(x,iorder,isize,first,last,level)
|
||||
end if
|
||||
end subroutine
|
||||
|
||||
subroutine set_order_LR(x,y,iorder,isize,jsize)
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: isize,jsize
|
||||
double precision :: x(isize,jsize)
|
||||
double precision :: y(isize,jsize)
|
||||
double precision,allocatable :: xtmp(:,:)
|
||||
double precision,allocatable :: ytmp(:,:)
|
||||
integer :: iorder(*)
|
||||
integer :: i,j
|
||||
|
||||
allocate(xtmp(isize,jsize),ytmp(isize,jsize))
|
||||
|
||||
do i=1,isize
|
||||
do j=1,jsize
|
||||
xtmp(i,j) = x(i,iorder(j))
|
||||
ytmp(i,j) = y(i,iorder(j))
|
||||
end do
|
||||
end do
|
||||
|
||||
do i=1,isize
|
||||
do j=1,jsize
|
||||
x(i,j) = xtmp(i,j)
|
||||
y(i,j) = ytmp(i,j)
|
||||
end do
|
||||
end do
|
||||
|
||||
deallocate(xtmp,ytmp)
|
||||
|
||||
end subroutine
|
||||
|
||||
subroutine set_order(x,iorder,isize,jsize)
|
||||
|
||||
implicit none
|
||||
|
||||
@ -1,3 +1,53 @@
|
||||
subroutine diagonalize_general_matrix_LR(N,A,WR,VL,VR)
|
||||
|
||||
! Diagonalize a non-symmetric square matrix
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: N
|
||||
double precision,intent(inout):: A(N,N)
|
||||
double precision,intent(out) :: VL(N,N)
|
||||
double precision,intent(out) :: VR(N,N)
|
||||
double precision,intent(out) :: WR(N)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i
|
||||
double precision :: tmp
|
||||
integer :: lwork,info
|
||||
double precision,allocatable :: work(:),WI(:)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(work(1),WI(N))
|
||||
|
||||
lwork = -1
|
||||
call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
|
||||
lwork = int(work(1))
|
||||
|
||||
deallocate(work)
|
||||
|
||||
allocate(work(lwork))
|
||||
|
||||
call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
|
||||
|
||||
do i=1,N
|
||||
tmp = dot_product(vl(:,i),vr(:,i))
|
||||
vl(:,i) = vl(:,i)/tmp
|
||||
end do
|
||||
|
||||
call matout(N,N,matmul(transpose(VL),VR))
|
||||
|
||||
deallocate(work,WI)
|
||||
|
||||
if(info /= 0) then
|
||||
print*,'Problem in diagonalize_general_matrix (dgeev)!!'
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
|
||||
subroutine diagonalize_general_matrix(N,A,WR,VR)
|
||||
|
||||
! Diagonalize a non-symmetric square matrix
|
||||
@ -31,7 +81,7 @@ subroutine diagonalize_general_matrix(N,A,WR,VR)
|
||||
|
||||
call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
|
||||
|
||||
deallocate(work, WI, VL)
|
||||
deallocate(work,WI,VL)
|
||||
|
||||
if(info /= 0) then
|
||||
print*,'Problem in diagonalize_general_matrix (dgeev)!!'
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user