PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-04-30 04:04:50 +02:00
Code Releases Activity

rCCD for Petros

Browse Source
This commit is contained in:
Pierre-Francois Loos 2025-01-21 15:01:59 +01:00
parent ad1888b021
commit 7d63abee67
4 changed files with 94 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/CC/EE_EOM_CCD_1h1p.f90
View File

@ -28,7 +28,8 @@ 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(:,:)
integer,allocatable :: order(:)
@ -46,7 +47,7 @@ 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),VL(nS,nS),VR(nS,nS))
allocate(order(nS))
@ -134,17 +135,22 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
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
end subroutine

3
src/CC/rCCD.f90
View File

@ -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

32
src/utils/sort.f90
View File

@ -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

52
src/utils/wrap_lapack.f90
View File

@ -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)!!'