LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

routines for exponential of matrices

This commit is contained in:
Pierre-Francois Loos 2023-11-07 10:25:52 +01:00
parent 8554099673
commit 66f6890082
5 changed files with 85 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF GHF ROHF
F T F F
T F F F
# MP2 MP3
F F
# CCD pCCD DCD CCSD CCSD(T)

6
src/GW/GGW.f90
View File

@ -105,9 +105,9 @@ subroutine GGW(doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT, &
if(doqsGW) then
call wall_time(start_GW)
! call qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
! eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
! dipole_int_AO,dipole_int,PHF,cHF,epsHF)
! call qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
! eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
! dipole_int_AO,dipole_int,PHF,cHF,epsHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW

2
src/GW/qsUGW.f90
View File

@ -57,8 +57,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Local variables
logical :: doGWPT = .false.
integer :: nSCF
integer :: nBasSq
integer :: ispin

20
src/HF/RHF_search.f90
View File

@ -48,6 +48,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: AB(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: R(:,:)
double precision,allocatable :: ExpR(:,:)
integer :: eig
double precision :: kick,step
@ -72,7 +74,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
!-------------------!
nS = (nO - nC)*(nV - nR)
allocate(ERI_MO(nBas,nBas,nBas,nBas),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS))
allocate(ERI_MO(nBas,nBas,nBas,nBas),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS), &
R(nBas,nBas),ExpR(nBas,nBas))
!------------------!
! Search algorithm !
@ -124,7 +127,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
call diagonalize_matrix(nS,AB,Om)
Om(:) = 0.5d0*Om(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real RHF -> Real RHF |'
write(*,*)'-------------------------------------------------------------'
@ -168,6 +171,19 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
end do
end do
! R(:,:) = 0d0
! ia = 0
! do i=nC+1,nO
! do a=nO+1,nBas-nR
! ia = ia + 1
! R(a,i) = +AB(ia,eig)
! R(i,a) = -AB(ia,eig)
! end do
! end do
! call matrix_exponential(nBas,R,ExpR)
! c = matmul(c,ExpR)
else
write(*,'(1X,A40,1X)') 'Well done, RHF solution is stable!'

63
src/utils/utils.f90
View File

@ -44,6 +44,69 @@ function KroneckerDelta(i,j) result(delta)
end function
!------------------------------------------------------------------------
subroutine diagonal_matrix(N,D,A)
! Construct diagonal matrix A from vector D
implicit none
integer,intent(in) :: N
double precision,intent(in) :: D(N)
double precision,intent(out) :: A(N,N)
integer :: i
A(:,:) = 0d0
do i=1,N
A(i,i) = D(i)
end do
end subroutine
!------------------------------------------------------------------------
subroutine matrix_exponential(N,A,ExpA)
! Compute Exp(A)
implicit none
integer,intent(in) :: N
double precision,intent(in) :: A(N,N)
double precision,allocatable :: W(:,:)
double precision,allocatable :: tau(:)
double precision,allocatable :: t(:,:)
double precision,intent(out) :: ExpA(N,N)
! Memory allocation
allocate(W(N,N),tau(N),t(N,N))
! Initialize
ExpA(:,:) = 0d0
! Diagonalize
W(:,:) = - matmul(A,A)
call diagonalize_matrix(N,W,tau)
tau(:) = sqrt(tau(:))
! Construct cos part
call diagonal_matrix(N,cos(tau),t)
t(:,:) = matmul(t,transpose(W))
ExpA(:,:) = ExpA(:,:) + matmul(W,t)
! Construct sin part
call diagonal_matrix(N,sin(tau)/tau,t)
t(:,:) = matmul(t,transpose(W))
t(:,:) = matmul(t,A)
ExpA(:,:) = ExpA(:,:) + matmul(W,t)
end subroutine
!------------------------------------------------------------------------
subroutine matout(m,n,A)