mirror of
https://github.com/pfloos/quack
synced 2025-01-03 01:56:09 +01:00
Merge pull request #8 from pfloos/master
impose biorthogonality in GT to be able to generate reference values for benchmarking
This commit is contained in:
AbdAmmar
GitHub
commit
2ca8557cb3
@ -1,5 +1,5 @@
|
||||
subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
|
||||
maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,eHF,cHF)
|
||||
maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc,ERHF,eHF,cHF)
|
||||
|
||||
! Coupled-cluster module
|
||||
|
||||
@ -34,7 +34,8 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: cHF(nBas,nBas)
|
||||
double precision,intent(in) :: Hc(nBas,nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
@ -47,7 +48,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
if(doCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -64,7 +65,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
|
||||
call wall_time(start_CC)
|
||||
call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, &
|
||||
ERI,ENuc,ERHF,eHF)
|
||||
ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -82,7 +83,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
if(doCCSD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -98,7 +99,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
if(dodrCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -114,7 +115,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
if(dorCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -130,7 +131,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
if(docrCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -146,7 +147,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
if(dolCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -162,7 +163,8 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
if(dopCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,eHF,cHF)
|
||||
call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERHF,eHF,cHF)
|
||||
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
|
||||
344
src/CC/pCCD.f90
344
src/CC/pCCD.f90
@ -1,4 +1,4 @@
|
||||
subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,eHF,cHF)
|
||||
subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERHF,eHF,cHF)
|
||||
|
||||
! pair CCD module
|
||||
|
||||
@ -8,25 +8,32 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: maxIt
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
integer,intent(in) :: nBas,nC,nO,nV,nR
|
||||
double precision,intent(in) :: ENuc,ERHF
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: cHF(nBas,nBas)
|
||||
double precision,intent(in) :: Hc(nBas,nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: p,q,r,s,t,u
|
||||
integer :: p,q,r,s,t,u,w
|
||||
integer :: pq,rs
|
||||
integer :: i,j,a,b
|
||||
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
integer :: nItAmp
|
||||
integer :: nItOrb
|
||||
double precision :: CvgAmp
|
||||
double precision :: CvgOrb
|
||||
double precision :: ECC
|
||||
double precision :: EcCC
|
||||
|
||||
@ -56,11 +63,15 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
double precision :: tr_2rdm
|
||||
|
||||
double precision :: E1,E2
|
||||
double precision,allocatable :: c(:,:)
|
||||
double precision,allocatable :: h(:,:)
|
||||
double precision,allocatable :: ERI_MO(:,:,:,:)
|
||||
double precision,allocatable :: grad(:)
|
||||
double precision,allocatable :: tmp(:,:,:,:)
|
||||
double precision,allocatable :: hess(:,:)
|
||||
double precision,allocatable :: eig(:)
|
||||
double precision,allocatable :: hessInv(:,:)
|
||||
double precision,allocatable :: Kap(:,:)
|
||||
double precision,allocatable :: ExpKap(:,:)
|
||||
|
||||
integer :: O,V,N
|
||||
integer :: n_diis
|
||||
@ -74,9 +85,9 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'**************************************'
|
||||
write(*,*)'| pair CCD calculation |'
|
||||
write(*,*)'**************************************'
|
||||
write(*,*)'*******************************'
|
||||
write(*,*)'* Restricted pCCD Calculation *'
|
||||
write(*,*)'*******************************'
|
||||
write(*,*)
|
||||
|
||||
! Useful quantities
|
||||
@ -85,53 +96,76 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
V = nV - nR
|
||||
N = O + V
|
||||
|
||||
! Form energy denominator
|
||||
!------------------------------------!
|
||||
! Star Loop for orbital optimization !
|
||||
!------------------------------------!
|
||||
|
||||
allocate(ERI_MO(N,N,N,N))
|
||||
allocate(c(N,N),h(N,N))
|
||||
allocate(eO(O),eV(V),delta_OV(O,V))
|
||||
allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V))
|
||||
|
||||
c(:,:) = cHF(nC+1:nBas-nR,nC+1:nBas-nR)
|
||||
|
||||
CvgOrb = 1d0
|
||||
nItOrb = 0
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| Orbital Optimization for pCCD |'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(CvgOrb > thresh .and. nItOrb < 1)
|
||||
|
||||
nItOrb = nItOrb + 1
|
||||
|
||||
! Transform integrals
|
||||
|
||||
h = matmul(transpose(c),matmul(Hc(nC+1:nBas-nR,nC+1:nBas-nR),c))
|
||||
call AOtoMO_ERI_RHF(N,c,ERI_AO(nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR),ERI_MO)
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
eO(:) = eHF(nC+1:nO)
|
||||
eV(:) = eHF(nO+1:nBas-nR)
|
||||
|
||||
call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV)
|
||||
do i=1,O
|
||||
do a=1,V
|
||||
delta_OV(i,a) = eV(a) - eO(i)
|
||||
end do
|
||||
end do
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V))
|
||||
! Create integral batches
|
||||
|
||||
do i=1,O
|
||||
do j=1,O
|
||||
OOOO(i,j) = ERI(nC+i,nC+i,nC+j,nC+j)
|
||||
OOOO(i,j) = ERI_MO(i,i,j,j)
|
||||
end do
|
||||
end do
|
||||
|
||||
do i=1,O
|
||||
do a=1,V
|
||||
OOVV(i,a) = ERI(nC+i,nC+i,nO+a,nO+a)
|
||||
OVOV(i,a) = ERI(nC+i,nO+a,nC+i,nO+a)
|
||||
OVVO(i,a) = ERI(nC+i,nO+a,nO+a,nC+i)
|
||||
OOVV(i,a) = ERI_MO(i,i,O+a,O+a)
|
||||
OVOV(i,a) = ERI_MO(i,O+a,i,O+a)
|
||||
OVVO(i,a) = ERI_MO(i,O+a,O+a,i)
|
||||
end do
|
||||
end do
|
||||
|
||||
do a=1,V
|
||||
do b=1,V
|
||||
VVVV(a,b) = ERI(nO+a,nO+a,nO+b,nO+b)
|
||||
VVVV(a,b) = ERI_MO(O+a,O+a,O+b,O+b)
|
||||
end do
|
||||
end do
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(t2(O,V),r2(O,V),yO(O,O),yV(V,V))
|
||||
|
||||
! Memory allocation for DIIS
|
||||
!----------------------------!
|
||||
! Star Loop for t amplitudes !
|
||||
!----------------------------!
|
||||
|
||||
allocate(t2(O,V),r2(O,V),yO(O,O))
|
||||
allocate(err_diis(O*V,max_diis),t2_diis(O*V,max_diis))
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute t ampltiudes
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
CvgAmp = 1d0
|
||||
nItAmp = 0
|
||||
ECC = ERHF
|
||||
EcCC = 0d0
|
||||
|
||||
@ -140,9 +174,6 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
t2_diis(:,:) = 0d0
|
||||
err_diis(:,:) = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| pCCD calculation: t amplitudes |'
|
||||
@ -151,11 +182,11 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
'|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
do while(CvgAmp > thresh .and. nItAmp < maxIt)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
nItAmp = nItAmp + 1
|
||||
|
||||
! Form intermediate array
|
||||
|
||||
@ -163,7 +194,6 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
|
||||
! Compute residual
|
||||
|
||||
|
||||
r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t2(:,:) &
|
||||
- 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t2(:,:))*t2(:,:)
|
||||
|
||||
@ -183,7 +213,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(:,:)))
|
||||
CvgAmp = maxval(abs(r2(:,:)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
@ -191,7 +221,12 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
EcCC = trace_matrix(V,matmul(transpose(OOVV),t2))
|
||||
EcCC = 0d0
|
||||
do i=1,O
|
||||
do a=1,V
|
||||
EcCC = EcCC + OOVV(i,a)*t2(i,a)
|
||||
end do
|
||||
end do
|
||||
|
||||
! Dump results
|
||||
|
||||
@ -207,53 +242,45 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
end if
|
||||
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
|
||||
'|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
|
||||
'|',nItAmp,'|',ECC+ENuc,'|',EcCC,'|',CvgAmp,'|'
|
||||
|
||||
end do
|
||||
write(*,*)'----------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
!---------------------------!
|
||||
! End Loop for t amplitudes !
|
||||
!---------------------------!
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
deallocate(r2,yO)
|
||||
deallocate(err_diis,t2_diis)
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nItAmp == maxIt) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed for t ampitudes '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)'! Convergence failed for t ampitudes !'
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
|
||||
stop
|
||||
|
||||
end if
|
||||
|
||||
! Deallocate memory
|
||||
|
||||
deallocate(err_diis,t2_diis)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(z2(O,V))
|
||||
|
||||
! Memory allocation for DIIS
|
||||
!-----------------------------!
|
||||
! Start Loop for z amplitudes !
|
||||
!-----------------------------!
|
||||
|
||||
allocate(z2(O,V),r2(O,V),yO(O,O),yV(V,V))
|
||||
allocate(err_diis(O*V,max_diis),z2_diis(O*V,max_diis))
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute z ampltiudes
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
CvgAmp = 1d0
|
||||
nItAmp = 0
|
||||
|
||||
n_diis = 0
|
||||
z2_diis(:,:) = 0d0
|
||||
err_diis(:,:) = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| pCCD calculation: z amplitudes |'
|
||||
@ -262,11 +289,11 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
'|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
do while(CvgAmp > thresh .and. nItAmp < maxIt)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
nItAmp = nItAmp + 1
|
||||
|
||||
! Form intermediate array
|
||||
|
||||
@ -296,7 +323,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(:,:)))
|
||||
CvgAmp = maxval(abs(r2(:,:)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
@ -312,45 +339,44 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
end if
|
||||
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
|
||||
'|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
|
||||
'|',nItAmp,'|',ECC+ENuc,'|',EcCC,'|',CvgAmp,'|'
|
||||
|
||||
end do
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
!---------------------------!
|
||||
! End Loop for z ampltiudes !
|
||||
!---------------------------!
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
deallocate(r2,yO,yV)
|
||||
deallocate(err_diis,z2_diis)
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nItAmp == maxIt) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)'! Convergence failed for z ampltiudes !'
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
|
||||
stop
|
||||
|
||||
end if
|
||||
|
||||
! Deallocate memory
|
||||
|
||||
deallocate(err_diis,z2_diis,r2)
|
||||
|
||||
!--------------------------!
|
||||
! Compute density matrices !
|
||||
!--------------------------!
|
||||
!--------------------------!
|
||||
! Compute density matrices !
|
||||
!--------------------------!
|
||||
|
||||
allocate(rdm1(N,N),rdm2(N,N,N,N))
|
||||
allocate(xOO(O,O),xVV(V,V),xOV(O,V))
|
||||
|
||||
xOO(:,:) = matmul(t2,transpose(z2))
|
||||
xVV(:,:) = matmul(transpose(z2),t2)
|
||||
xOV(:,:) = matmul(t2,matmul(transpose(z2),t2))
|
||||
|
||||
! Form 1RDM
|
||||
|
||||
allocate(rdm1(N,N))
|
||||
! Form 1RDM
|
||||
|
||||
rdm1(:,:) = 0d0
|
||||
|
||||
@ -362,7 +388,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
rdm1(O+a,O+a) = 2d0*xVV(a,a)
|
||||
end do
|
||||
|
||||
! Check 1RDM
|
||||
! Check 1RDM
|
||||
|
||||
tr_1rdm = trace_matrix(N,rdm1)
|
||||
write(*,'(A25,F16.10)') ' --> Trace of the 1RDM = ',tr_1rdm
|
||||
@ -371,12 +397,10 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
write(*,*) ' !!! Your 1RDM seems broken !!! '
|
||||
write(*,*)
|
||||
|
||||
! write(*,*) '1RDM is diagonal at the pCCD level:'
|
||||
! call matout(N,N,rdm1)
|
||||
! write(*,*) '1RDM is diagonal at the pCCD level:'
|
||||
! call matout(N,N,rdm1)
|
||||
|
||||
! Form 2RM
|
||||
|
||||
allocate(rdm2(N,N,N,N))
|
||||
! Form 2RM
|
||||
|
||||
rdm2(:,:,:,:) = 0d0
|
||||
|
||||
@ -472,7 +496,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
rdm2(O+a,O+a,O+a,O+a) = 2d0*xVV(a,a)
|
||||
end do
|
||||
|
||||
! Check 2RDM
|
||||
! Check 2RDM
|
||||
|
||||
tr_2rdm = trace_matrix(N**2,rdm2)
|
||||
write(*,'(A25,F16.10)') ' --> Trace of the 2RDM = ',tr_2rdm
|
||||
@ -481,13 +505,13 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
write(*,*) ' !!! Your 2RDM seems broken !!! '
|
||||
write(*,*)
|
||||
|
||||
! write(*,*) '2RDM is not diagonal at the pCCD level:'
|
||||
! call matout(N**2,N**2,rdm2)
|
||||
! write(*,*) '2RDM is not diagonal at the pCCD level:'
|
||||
! call matout(N**2,N**2,rdm2)
|
||||
|
||||
! Compute electronic energy
|
||||
deallocate(xOO,xVV,xOV)
|
||||
deallocate(t2,z2)
|
||||
|
||||
allocate(h(N,N))
|
||||
h = matmul(transpose(cHF),matmul(Hc,cHF))
|
||||
! Compute electronic energy
|
||||
|
||||
E1 = 0d0
|
||||
E2 = 0d0
|
||||
@ -497,7 +521,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
E1 = E1 + rdm1(p,q)*h(p,q)
|
||||
do r=1,N
|
||||
do s=1,N
|
||||
E2 = E2 + rdm2(p,q,r,s)*ERI(p,q,r,s)
|
||||
E2 = E2 + rdm2(p,q,r,s)*ERI_MO(p,q,r,s)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
@ -511,7 +535,9 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
write(*,'(A25,F16.10)') ' Total energy = ',E1 + E2 + ENuc
|
||||
write(*,*)
|
||||
|
||||
! Compute gradient
|
||||
!--------------------------!
|
||||
! Compute orbital gradient !
|
||||
!--------------------------!
|
||||
|
||||
allocate(grad(N**2))
|
||||
|
||||
@ -530,7 +556,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
do r=1,N
|
||||
do s=1,N
|
||||
do t=1,N
|
||||
grad(pq) = grad(pq) + (ERI(r,s,p,t)*rdm2(r,s,q,t) - ERI(q,t,r,s)*rdm2(p,t,r,s))
|
||||
grad(pq) = grad(pq) + (ERI_MO(r,s,p,t)*rdm2(r,s,q,t) - ERI_MO(q,t,r,s)*rdm2(p,t,r,s))
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
@ -540,8 +566,18 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
|
||||
write(*,*) 'Orbital gradient at the pCCD level:'
|
||||
call matout(N,N,grad)
|
||||
write(*,*)
|
||||
|
||||
! Compute Hessian
|
||||
! Check convergence of orbital optimization
|
||||
|
||||
CvgOrb = maxval(abs(grad))
|
||||
write(*,*) ' Iteration',nItOrb,'for pCCD orbital optimization'
|
||||
write(*,*) ' Convergence of orbital gradient = ',CvgOrb
|
||||
write(*,*)
|
||||
|
||||
!-------------------------!
|
||||
! Compute orbital Hessian !
|
||||
!-------------------------!
|
||||
|
||||
allocate(hess(N**2,N**2),tmp(N,N,N,N))
|
||||
|
||||
@ -550,7 +586,6 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
do p=1,N
|
||||
do q=1,N
|
||||
|
||||
rs = 0
|
||||
do r=1,N
|
||||
do s=1,N
|
||||
|
||||
@ -565,9 +600,9 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
end do
|
||||
|
||||
do u=1,N
|
||||
do v=1,N
|
||||
do w=1,N
|
||||
|
||||
tmp(p,q,r,s) = tmp(p,q,r,s) + ERI(u,v,p,r)*rdm2(u,v,q,s) + ERI(q,s,u,v)*rdm2(p,r,u,v)
|
||||
tmp(p,q,r,s) = tmp(p,q,r,s) + ERI_MO(u,w,p,r)*rdm2(u,w,q,s) + ERI_MO(q,s,u,w)*rdm2(p,r,u,w)
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -576,19 +611,19 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
do u=1,N
|
||||
|
||||
tmp(p,q,r,s) = tmp(p,q,r,s) - ( &
|
||||
ERI(s,t,p,u)*rdm2(r,t,q,u) + ERI(t,s,p,u)*rdm2(t,r,q,u) &
|
||||
+ ERI(q,u,r,t)*rdm2(p,u,s,t) + ERI(q,u,t,r)*rdm2(p,u,t,s) )
|
||||
ERI_MO(s,t,p,u)*rdm2(r,t,q,u) + ERI_MO(t,s,p,u)*rdm2(t,r,q,u) &
|
||||
+ ERI_MO(q,u,r,t)*rdm2(p,u,s,t) + ERI_MO(q,u,t,r)*rdm2(p,u,t,s) )
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
do t=1,N
|
||||
do u=1,N
|
||||
do v=1,N
|
||||
do w=1,N
|
||||
|
||||
tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( &
|
||||
Kronecker_delta(q,r)*(ERI(u,v,p,t)*rdm2(u,v,s,t) + ERI(s,t,u,v)*rdm2(p,t,u,v)) &
|
||||
+ Kronecker_delta(p,s)*(ERI(q,t,u,v)*rdm2(r,t,u,v) + ERI(u,v,r,t)*rdm2(u,v,q,t)) )
|
||||
Kronecker_delta(q,r)*(ERI_MO(u,w,p,t)*rdm2(u,w,s,t) + ERI_MO(s,t,u,w)*rdm2(p,t,u,w)) &
|
||||
+ Kronecker_delta(p,s)*(ERI_MO(q,t,u,w)*rdm2(r,t,u,w) + ERI_MO(u,w,r,t)*rdm2(u,w,q,t)) )
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -600,7 +635,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
end do
|
||||
end do
|
||||
|
||||
! Flatten Hessian matrix and add permutations
|
||||
! Flatten Hessian matrix and add permutations
|
||||
|
||||
pq = 0
|
||||
do p=1,N
|
||||
@ -614,7 +649,8 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
|
||||
rs = rs + 1
|
||||
|
||||
hess(pq,rs) = tmp(p,q,r,s) - tmp(q,p,r,s) - tmp(p,q,s,r) + tmp(q,p,s,r)
|
||||
hess(pq,rs) = tmp(p,r,q,s) - tmp(r,p,q,s) - tmp(p,r,s,q) + tmp(r,p,s,q)
|
||||
!! hess(pq,rs) = tmp(p,q,r,s) - tmp(q,p,r,s) - tmp(p,q,s,r) + tmp(q,p,s,r)
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -622,15 +658,81 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
|
||||
end do
|
||||
end do
|
||||
|
||||
call matout(N**2,N**2,hess)
|
||||
deallocate(rdm1,rdm2,tmp)
|
||||
|
||||
deallocate(tmp)
|
||||
allocate(hessInv(N**2,N**2))
|
||||
|
||||
allocate(eig(N**2))
|
||||
call inverse_matrix(N**2,hess,hessInv)
|
||||
|
||||
call diagonalize_matrix(N**2,hess,eig)
|
||||
deallocate(hess)
|
||||
|
||||
call vecout(N**2,eig)
|
||||
allocate(Kap(N,N))
|
||||
|
||||
Kap(:,:) = 0d0
|
||||
|
||||
pq = 0
|
||||
do p=1,N
|
||||
do q=1,N
|
||||
|
||||
pq = pq + 1
|
||||
|
||||
rs = 0
|
||||
do r=1,N
|
||||
do s=1,N
|
||||
|
||||
rs = rs + 1
|
||||
|
||||
Kap(p,q) = Kap(p,q) - hessInv(pq,rs)*grad(rs)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
deallocate(hessInv,grad)
|
||||
|
||||
write(*,*) 'kappa'
|
||||
call matout(N,N,Kap)
|
||||
write(*,*)
|
||||
|
||||
allocate(ExpKap(N,N))
|
||||
call matrix_exponential(N,Kap,ExpKap)
|
||||
deallocate(Kap)
|
||||
|
||||
write(*,*) 'e^kappa'
|
||||
call matout(N,N,ExpKap)
|
||||
write(*,*)
|
||||
|
||||
write(*,*) 'Old orbitals'
|
||||
call matout(N,N,c)
|
||||
write(*,*)
|
||||
|
||||
c = matmul(c,ExpKap)
|
||||
deallocate(ExpKap)
|
||||
|
||||
write(*,*) 'Rotated orbitals'
|
||||
call matout(N,N,c)
|
||||
write(*,*)
|
||||
|
||||
end do
|
||||
|
||||
!-----------------------------------!
|
||||
! End Loop for orbital optimization !
|
||||
!-----------------------------------!
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nItOrb == maxIt) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)'! Convergence failed for orbital optimization !'
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
|
||||
stop
|
||||
|
||||
end if
|
||||
|
||||
! Testing zone
|
||||
|
||||
|
||||
134
src/LR/ppLR.f90
134
src/LR/ppLR.f90
@ -1,93 +1,123 @@
|
||||
subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
|
||||
|
||||
! Solve the pp-RPA linear eigenvalue problem
|
||||
! ---
|
||||
|
||||
subroutine ppLR(TDA, nOO, nVV, Bpp, Cpp, Dpp, Om1, X1, Y1, Om2, X2, Y2, EcRPA)
|
||||
|
||||
!
|
||||
! Solve the pp-RPA linear eigenvalue problem
|
||||
!
|
||||
! right eigen-problem: H R = R w
|
||||
! left eigen-problem: H.T L = L w
|
||||
!
|
||||
! where L.T R = 1
|
||||
!
|
||||
!
|
||||
! (+C +B)
|
||||
! H = ( ) where C = C.T and D = D.T
|
||||
! (-B.T -D)
|
||||
!
|
||||
! (w1 0) (X1 X2) (+X1 +X2)
|
||||
! w = ( ), R = ( ) and L = ( )
|
||||
! (0 w2) (Y1 Y2) (-Y1 -Y2)
|
||||
!
|
||||
!
|
||||
! the normalisation condition reduces to
|
||||
!
|
||||
! X1.T X2 - Y1.T Y2 = 0
|
||||
! X1.T X1 - Y1.T Y1 = 1
|
||||
! X2.T X2 - Y2.T Y2 = 1
|
||||
!
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
logical, intent(in) :: TDA
|
||||
integer, intent(in) :: nOO, nVV
|
||||
double precision, intent(in) :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO)
|
||||
double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV)
|
||||
double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO)
|
||||
double precision, intent(out) :: EcRPA
|
||||
|
||||
logical,intent(in) :: TDA
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
double precision,intent(in) :: Bpp(nVV,nOO)
|
||||
double precision,intent(in) :: Cpp(nVV,nVV)
|
||||
double precision,intent(in) :: Dpp(nOO,nOO)
|
||||
logical :: imp_bio, verbose
|
||||
integer :: i, j, N
|
||||
double precision :: EcRPA1, EcRPA2
|
||||
double precision :: thr_d, thr_nd, thr_deg
|
||||
double precision,allocatable :: M(:,:), Z(:,:), Om(:)
|
||||
|
||||
! Local variables
|
||||
double precision, external :: trace_matrix
|
||||
|
||||
double precision :: trace_matrix
|
||||
double precision :: EcRPA1
|
||||
double precision :: EcRPA2
|
||||
double precision,allocatable :: M(:,:)
|
||||
double precision,allocatable :: Z(:,:)
|
||||
double precision,allocatable :: Om(:)
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Om1(nVV)
|
||||
double precision,intent(out) :: X1(nVV,nVV)
|
||||
double precision,intent(out) :: Y1(nOO,nVV)
|
||||
double precision,intent(out) :: Om2(nOO)
|
||||
double precision,intent(out) :: X2(nVV,nOO)
|
||||
double precision,intent(out) :: Y2(nOO,nOO)
|
||||
double precision,intent(out) :: EcRPA
|
||||
N = nOO + nVV
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Om(nOO+nVV))
|
||||
|
||||
!-------------------------------------------------!
|
||||
! Solve the p-p eigenproblem !
|
||||
!-------------------------------------------------!
|
||||
! !
|
||||
! | C B | | X1 X2 | | w1 0 | | X1 X2 | !
|
||||
! | | | | = | | | | !
|
||||
! | -Bt -D | | Y1 Y2 | | 0 w2 | | Y1 Y2 | !
|
||||
! !
|
||||
!-------------------------------------------------!
|
||||
allocate(M(N,N), Z(N,N), Om(N))
|
||||
|
||||
if(TDA) then
|
||||
|
||||
X1(:,:) = +Cpp(:,:)
|
||||
Y1(:,:) = 0d0
|
||||
if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1)
|
||||
if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1)
|
||||
|
||||
X2(:,:) = 0d0
|
||||
Y2(:,:) = -Dpp(:,:)
|
||||
if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2)
|
||||
if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2)
|
||||
|
||||
else
|
||||
|
||||
! Diagonal blocks
|
||||
|
||||
M( 1:nVV , 1:nVV) = + Cpp(1:nVV,1:nVV)
|
||||
M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO)
|
||||
|
||||
! Off-diagonal blocks
|
||||
|
||||
M( 1:nVV ,nVV+1:nOO+nVV) = - Bpp(1:nVV,1:nOO)
|
||||
M(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp(1:nVV,1:nOO))
|
||||
|
||||
! call matout(nOO,nOO,Dpp)
|
||||
!! Diagonalize the p-p matrix
|
||||
!if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z)
|
||||
!! Split the various quantities in p-p and h-h parts
|
||||
!call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2)
|
||||
|
||||
! Diagonalize the p-p matrix
|
||||
|
||||
if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z)
|
||||
thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1
|
||||
thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1
|
||||
thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not
|
||||
imp_bio = .True. ! impose bi-orthogonality
|
||||
verbose = .False.
|
||||
call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
|
||||
|
||||
! Split the various quantities in p-p and h-h parts
|
||||
do i = 1, nOO
|
||||
Om2(i) = Om(i)
|
||||
do j = 1, nVV
|
||||
X2(j,i) = Z(j,i)
|
||||
enddo
|
||||
do j = 1, nOO
|
||||
Y2(j,i) = Z(nVV+j,i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
|
||||
do i = 1, nVV
|
||||
Om1(i) = Om(nOO+i)
|
||||
do j = 1, nVV
|
||||
X1(j,i) = M(j,nOO+i)
|
||||
enddo
|
||||
do j = 1, nOO
|
||||
Y1(j,i) = M(nVV+j,nOO+i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end if
|
||||
|
||||
! Compute the RPA correlation energy
|
||||
! Compute the RPA correlation energy
|
||||
EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp))
|
||||
EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp)
|
||||
EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp)
|
||||
|
||||
EcRPA = 0.5d0*( sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp) )
|
||||
EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp)
|
||||
EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp)
|
||||
|
||||
if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) &
|
||||
if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then
|
||||
print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'
|
||||
endif
|
||||
|
||||
deallocate(M, Z, Om)
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
|
||||
@ -228,7 +228,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
|
||||
call wall_time(start_CC)
|
||||
call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
|
||||
maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,Hc,ERI_MO,ENuc,ERHF,eHF,cHF)
|
||||
maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc,ERHF,eHF,cHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
|
||||
@ -89,8 +89,8 @@ FIX_ORDER_OF_LIBS=-Wl,--start-group
|
||||
|
||||
if sys.platform in ["linux", "linux2"]:
|
||||
# compiler = compile_gfortran_linux
|
||||
# compiler = compile_ifort_linux
|
||||
compiler = compile_olympe
|
||||
compiler = compile_ifort_linux
|
||||
# compiler = compile_olympe
|
||||
elif sys.platform == "darwin":
|
||||
compiler = compile_gfortran_mac
|
||||
else:
|
||||
|
||||
571
src/utils/non_sym_diag.f90
Normal file
571
src/utils/non_sym_diag.f90
Normal file
@ -0,0 +1,571 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine diagonalize_nonsym_matrix(N, A, L, e_re, thr_d, thr_nd, thr_deg, imp_bio, verbose)
|
||||
|
||||
! Diagonalize a non-symmetric matrix
|
||||
!
|
||||
! Output
|
||||
! right-eigenvectors are saved in A
|
||||
! left-eigenvectors are saved in L
|
||||
! eigenvalues are saved in e = e_re + i e_im
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: N
|
||||
logical, intent(in) :: imp_bio, verbose
|
||||
double precision, intent(in) :: thr_d, thr_nd, thr_deg
|
||||
double precision, intent(inout) :: A(N,N)
|
||||
double precision, intent(out) :: e_re(N), L(N,N)
|
||||
|
||||
integer :: i, j, ii
|
||||
integer :: lwork, info
|
||||
double precision :: accu_d, accu_nd
|
||||
integer, allocatable :: iorder(:), deg_num(:)
|
||||
double precision, allocatable :: Atmp(:,:), Ltmp(:,:), work(:), e_im(:)
|
||||
double precision, allocatable :: S(:,:)
|
||||
|
||||
if(verbose) then
|
||||
print*, ' Starting a non-Hermitian diagonalization ...'
|
||||
print*, ' Good Luck ;)'
|
||||
print*, ' imp_bio = ', imp_bio
|
||||
endif
|
||||
|
||||
! ---
|
||||
! diagonalize
|
||||
|
||||
allocate(Atmp(N,N), e_im(N))
|
||||
Atmp(1:N,1:N) = A(1:N,1:N)
|
||||
|
||||
allocate(work(1))
|
||||
lwork = -1
|
||||
call dgeev('V', 'V', N, Atmp, N, e_re, e_im, L, N, A, N, work, lwork, info)
|
||||
if(info .gt. 0) then
|
||||
print*,'dgeev failed !!', info
|
||||
stop
|
||||
endif
|
||||
|
||||
lwork = max(int(work(1)), 1)
|
||||
deallocate(work)
|
||||
allocate(work(lwork))
|
||||
|
||||
call dgeev('V', 'V', N, Atmp, N, e_re, e_im, L, N, A, N, work, lwork, info)
|
||||
if(info .ne. 0) then
|
||||
print*,'dgeev failed !!', info
|
||||
stop
|
||||
endif
|
||||
|
||||
deallocate(Atmp, WORK)
|
||||
|
||||
|
||||
! ---
|
||||
! check if eigenvalues are real
|
||||
|
||||
i = 1
|
||||
ii = 0
|
||||
do while(i .le. N)
|
||||
if(dabs(e_im(i)) .gt. 1.d-12) then
|
||||
ii = ii + 1
|
||||
if(verbose) then
|
||||
print*, ' Warning: complex eigenvalue !'
|
||||
print*, i, e_re(i), e_im(i)
|
||||
if(dabs(e_im(i)/e_re(i)) .lt. 1.d-6) then
|
||||
print*, ' small enouph to be igored'
|
||||
else
|
||||
print*, ' IMAGINARY PART IS SIGNIFANT !!!'
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
i = i + 1
|
||||
enddo
|
||||
|
||||
if(verbose) then
|
||||
if(ii .eq. 0) print*, ' congratulations :) eigenvalues are real-valued !!'
|
||||
endif
|
||||
|
||||
|
||||
! ---
|
||||
! track & sort the real eigenvalues
|
||||
|
||||
allocate(Atmp(N,N), Ltmp(N,N), iorder(N))
|
||||
|
||||
do i = 1, N
|
||||
iorder(i) = i
|
||||
enddo
|
||||
call quick_sort(e_re, iorder, N)
|
||||
|
||||
Atmp(:,:) = A(:,:)
|
||||
Ltmp(:,:) = L(:,:)
|
||||
do i = 1, N
|
||||
do j = 1, N
|
||||
A(j,i) = Atmp(j,iorder(i))
|
||||
L(j,i) = Ltmp(j,iorder(i))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
deallocate(Atmp, Ltmp, iorder)
|
||||
|
||||
|
||||
|
||||
|
||||
! ---
|
||||
! check bi-orthog
|
||||
|
||||
allocate(S(N,N))
|
||||
call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .false., verbose)
|
||||
|
||||
if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(N))/dble(N) .lt. thr_d)) then
|
||||
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are normalized and bi-orthogonalized'
|
||||
endif
|
||||
|
||||
elseif((accu_nd .lt. thr_nd) .and. (dabs(accu_d - dble(N)) .gt. thr_d)) then
|
||||
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are not normalized but bi-orthogonalized'
|
||||
endif
|
||||
|
||||
call check_biorthog_binormalize(N, N, L, A, thr_d, thr_nd, .true.)
|
||||
call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .true., verbose)
|
||||
|
||||
else
|
||||
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are not normalized neither bi-orthogonalized'
|
||||
endif
|
||||
|
||||
allocate(deg_num(N))
|
||||
call reorder_degen_eigvec(N, thr_deg, deg_num, e_re, L, A)
|
||||
call impose_biorthog_degen_eigvec(N, deg_num, e_re, L, A)
|
||||
deallocate(deg_num)
|
||||
|
||||
call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .false., verbose)
|
||||
if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(N))/dble(N) .lt. thr_d)) then
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are now normalized and bi-orthogonalized'
|
||||
endif
|
||||
elseif((accu_nd .lt. thr_nd) .and. (dabs(accu_d - dble(N)) .gt. thr_d)) then
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are now not normalized but bi-orthogonalized'
|
||||
endif
|
||||
call check_biorthog_binormalize(N, N, L, A, thr_d, thr_nd, .true.)
|
||||
call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .true., verbose)
|
||||
else
|
||||
if(verbose) then
|
||||
print*, ' bi-orthogonalization failed !'
|
||||
endif
|
||||
if(imp_bio) then
|
||||
print*, ' bi-orthogonalization failed !'
|
||||
deallocate(S)
|
||||
stop
|
||||
endif
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
deallocate(S)
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ifnot, verbose)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n, m
|
||||
logical, intent(in) :: stop_ifnot, verbose
|
||||
double precision, intent(in) :: Vl(n,m), Vr(n,m)
|
||||
double precision, intent(in) :: thr_d, thr_nd
|
||||
double precision, intent(out) :: accu_d, accu_nd, S(m,m)
|
||||
|
||||
integer :: i, j
|
||||
double precision, allocatable :: SS(:,:)
|
||||
|
||||
if(verbose) then
|
||||
print *, ' check bi-orthogonality'
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, Vl, size(Vl, 1), Vr, size(Vr, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
accu_d = 0.d0
|
||||
accu_nd = 0.d0
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
if(i==j) then
|
||||
accu_d = accu_d + dabs(S(i,i))
|
||||
else
|
||||
accu_nd = accu_nd + S(j,i) * S(j,i)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
accu_nd = dsqrt(accu_nd) / dble(m)
|
||||
|
||||
if(verbose) then
|
||||
if((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d) then
|
||||
print *, ' non bi-orthogonal vectors !'
|
||||
print *, ' accu_nd = ', accu_nd
|
||||
print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m)
|
||||
else
|
||||
print *, ' vectors are bi-orthogonals'
|
||||
endif
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
if(stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) then
|
||||
print *, ' non bi-orthogonal vectors !'
|
||||
print *, ' accu_nd = ', accu_nd
|
||||
print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m)
|
||||
stop
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n, m
|
||||
logical, intent(in) :: stop_ifnot
|
||||
double precision, intent(in) :: thr_d, thr_nd
|
||||
double precision, intent(inout) :: Vl(n,m), Vr(n,m)
|
||||
|
||||
integer :: i, j
|
||||
double precision :: accu_d, accu_nd, s_tmp
|
||||
double precision, allocatable :: S(:,:)
|
||||
|
||||
! ---
|
||||
|
||||
allocate(S(m,m))
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, Vl, size(Vl, 1), Vr, size(Vr, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
do i = 1, m
|
||||
if(S(i,i) .lt. 0.d0) then
|
||||
do j = 1, n
|
||||
Vl(j,i) = -1.d0 * Vl(j,i)
|
||||
enddo
|
||||
S(i,i) = -S(i,i)
|
||||
endif
|
||||
enddo
|
||||
|
||||
accu_d = 0.d0
|
||||
accu_nd = 0.d0
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
if(i==j) then
|
||||
accu_d = accu_d + S(i,i)
|
||||
else
|
||||
accu_nd = accu_nd + S(j,i) * S(j,i)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
accu_nd = dsqrt(accu_nd) / dble(m)
|
||||
|
||||
! ---
|
||||
|
||||
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then
|
||||
|
||||
do i = 1, m
|
||||
if(S(i,i) <= 0.d0) then
|
||||
print *, ' overap negative'
|
||||
print *, i, S(i,i)
|
||||
exit
|
||||
endif
|
||||
if(dabs(S(i,i) - 1.d0) .gt. thr_d) then
|
||||
s_tmp = 1.d0 / dsqrt(S(i,i))
|
||||
do j = 1, n
|
||||
Vl(j,i) = Vl(j,i) * s_tmp
|
||||
Vr(j,i) = Vr(j,i) * s_tmp
|
||||
enddo
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, Vl, size(Vl, 1), Vr, size(Vr, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
accu_d = 0.d0
|
||||
accu_nd = 0.d0
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
if(i==j) then
|
||||
accu_d = accu_d + S(i,i)
|
||||
else
|
||||
accu_nd = accu_nd + S(j,i) * S(j,i)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
accu_nd = dsqrt(accu_nd) / dble(m)
|
||||
|
||||
deallocate(S)
|
||||
|
||||
! ---
|
||||
|
||||
if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) ) then
|
||||
print *, accu_nd, thr_nd
|
||||
print *, dabs(accu_d-dble(m))/dble(m), thr_d
|
||||
print *, ' biorthog_binormalize failed !'
|
||||
stop
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine reorder_degen_eigvec(n, thr_deg, deg_num, e0, L0, R0)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n
|
||||
double precision, intent(in) :: thr_deg
|
||||
double precision, intent(inout) :: e0(n), L0(n,n), R0(n,n)
|
||||
integer, intent(out) :: deg_num(n)
|
||||
|
||||
logical :: complex_root
|
||||
integer :: i, j, k, m, ii, j_tmp
|
||||
double precision :: ei, ej, de
|
||||
double precision :: accu_d, accu_nd
|
||||
double precision :: e0_tmp, L0_tmp(n), R0_tmp(n)
|
||||
double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:)
|
||||
|
||||
do i = 1, n
|
||||
deg_num(i) = 1
|
||||
enddo
|
||||
|
||||
do i = 1, n-1
|
||||
ei = e0(i)
|
||||
|
||||
! already considered in degen vectors
|
||||
if(deg_num(i) .eq. 0) cycle
|
||||
|
||||
ii = 0
|
||||
do j = i+1, n
|
||||
ej = e0(j)
|
||||
de = dabs(ei - ej)
|
||||
|
||||
if(de .lt. thr_deg) then
|
||||
ii = ii + 1
|
||||
|
||||
j_tmp = i + ii
|
||||
|
||||
deg_num(j_tmp) = 0
|
||||
|
||||
e0_tmp = e0(j_tmp)
|
||||
e0(j_tmp) = e0(j)
|
||||
e0(j) = e0_tmp
|
||||
|
||||
L0_tmp(1:n) = L0(1:n,j_tmp)
|
||||
L0(1:n,j_tmp) = L0(1:n,j)
|
||||
L0(1:n,j) = L0_tmp(1:n)
|
||||
|
||||
R0_tmp(1:n) = R0(1:n,j_tmp)
|
||||
R0(1:n,j_tmp) = R0(1:n,j)
|
||||
R0(1:n,j) = R0_tmp(1:n)
|
||||
endif
|
||||
enddo
|
||||
|
||||
deg_num(i) = ii + 1
|
||||
enddo
|
||||
|
||||
ii = 0
|
||||
do i = 1, n
|
||||
if(deg_num(i) .gt. 1) then
|
||||
ii = ii + 1
|
||||
endif
|
||||
enddo
|
||||
|
||||
if(ii .eq. 0) then
|
||||
print*, ' WARNING: bi-orthogonality is lost but there is no degeneracies'
|
||||
print*, ' rotations may change energy'
|
||||
stop
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n, deg_num(n)
|
||||
double precision, intent(in) :: e0(n)
|
||||
double precision, intent(inout) :: L0(n,n), R0(n,n)
|
||||
|
||||
logical :: complex_root
|
||||
integer :: i, j, k, m
|
||||
double precision :: ei, ej, de
|
||||
double precision :: accu_d, accu_nd
|
||||
double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:)
|
||||
|
||||
!do i = 1, n
|
||||
! if(deg_num(i) .gt. 1) then
|
||||
! print *, ' degen on', i, deg_num(i), e0(i)
|
||||
! endif
|
||||
!enddo
|
||||
|
||||
! ---
|
||||
|
||||
do i = 1, n
|
||||
m = deg_num(i)
|
||||
|
||||
if(m .gt. 1) then
|
||||
|
||||
allocate(L(n,m), R(n,m), S(m,m))
|
||||
|
||||
do j = 1, m
|
||||
L(1:n,j) = L0(1:n,i+j-1)
|
||||
R(1:n,j) = R0(1:n,i+j-1)
|
||||
enddo
|
||||
|
||||
! ---
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, L, size(L, 1), R, size(R, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
accu_nd = 0.d0
|
||||
do j = 1, m
|
||||
do k = 1, m
|
||||
if(j==k) cycle
|
||||
accu_nd = accu_nd + dabs(S(j,k))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if(accu_nd .lt. 1d-12) then
|
||||
deallocate(S, L, R)
|
||||
cycle
|
||||
endif
|
||||
|
||||
call impose_biorthog_svd(n, m, L, R)
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, L, size(L, 1), R, size(R, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
accu_nd = 0.d0
|
||||
do j = 1, m
|
||||
do k = 1, m
|
||||
if(j==k) cycle
|
||||
accu_nd = accu_nd + dabs(S(j,k))
|
||||
enddo
|
||||
enddo
|
||||
if(accu_nd .gt. 1d-12) then
|
||||
print*, ' accu_nd =', accu_nd
|
||||
print*, ' your strategy for degenerates orbitals failed !'
|
||||
print*, m, 'deg on', i
|
||||
stop
|
||||
endif
|
||||
|
||||
deallocate(S)
|
||||
|
||||
! ---
|
||||
|
||||
do j = 1, m
|
||||
L0(1:n,i+j-1) = L(1:n,j)
|
||||
R0(1:n,i+j-1) = R(1:n,j)
|
||||
enddo
|
||||
|
||||
deallocate(L, R)
|
||||
|
||||
endif
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine impose_biorthog_svd(n, m, L, R)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n, m
|
||||
double precision, intent(inout) :: L(n,m), R(n,m)
|
||||
|
||||
integer :: i, j, num_linear_dependencies
|
||||
double precision :: threshold
|
||||
double precision, allocatable :: S(:,:), tmp(:,:)
|
||||
double precision, allocatable :: U(:,:), V(:,:), Vt(:,:), D(:)
|
||||
|
||||
allocate(S(m,m))
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, L, size(L, 1), R, size(R, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
! ---
|
||||
|
||||
allocate(U(m,m), Vt(m,m), D(m))
|
||||
|
||||
call svd(S, m, U, m, D, Vt, m, m, m)
|
||||
|
||||
deallocate(S)
|
||||
|
||||
threshold = 1.d-6
|
||||
num_linear_dependencies = 0
|
||||
do i = 1, m
|
||||
if(abs(D(i)) <= threshold) then
|
||||
D(i) = 0.d0
|
||||
num_linear_dependencies = num_linear_dependencies + 1
|
||||
else
|
||||
D(i) = 1.d0 / dsqrt(D(i))
|
||||
endif
|
||||
enddo
|
||||
if(num_linear_dependencies > 0) then
|
||||
write(*,*) ' linear dependencies = ', num_linear_dependencies
|
||||
write(*,*) ' m = ', m
|
||||
stop
|
||||
endif
|
||||
|
||||
allocate(V(m,m))
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
V(j,i) = Vt(i,j)
|
||||
enddo
|
||||
enddo
|
||||
deallocate(Vt)
|
||||
|
||||
! ---
|
||||
|
||||
! R <-- R x V x D^{-0.5}
|
||||
! L <-- L x U x D^{-0.5}
|
||||
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
V(j,i) = V(j,i) * D(i)
|
||||
U(j,i) = U(j,i) * D(i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
allocate(tmp(n,m))
|
||||
tmp(:,:) = R(:,:)
|
||||
call dgemm( 'N', 'N', n, m, m, 1.d0 &
|
||||
, tmp, size(tmp, 1), V, size(V, 1) &
|
||||
, 0.d0, R, size(R, 1))
|
||||
|
||||
tmp(:,:) = L(:,:)
|
||||
call dgemm( 'N', 'N', n, m, m, 1.d0 &
|
||||
, tmp, size(tmp, 1), U, size(U, 1) &
|
||||
, 0.d0, L, size(L, 1))
|
||||
|
||||
deallocate(tmp, U, V, D)
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
Loading…
Reference in New Issue
Block a user