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

correct big bug in HF search

This commit is contained in:
Pierre-Francois Loos 2023-11-28 10:40:15 +01:00
parent 337a950172
commit f89d54c59b
10 changed files with 82 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/GW/RGW.f90
View File

@ -1,9 +1,9 @@
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF)
ERI_AO,ERI_MO,dipole_int_AO,dipole_int,PHF,cHF,eHF)
! GW module
! Restricted GW module
implicit none
include 'parameters.h'
@ -60,7 +60,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -76,7 +76,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -93,7 +93,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -110,7 +110,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
dipole_int_AO,dipole_int,PHF,cHF,eHF)
call wall_time(end_GW)
@ -128,7 +128,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, &
singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI, &
singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
dipole_int_AO,dipole_int,PHF,cHF,eHF)
call wall_time(end_GW)
@ -145,7 +145,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doufG0W0) then
call wall_time(start_GW)
call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -161,7 +161,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doufGW) then
call wall_time(start_GW)
call ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
call ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW

44
src/HF/GHF_search.f90
View File

@ -1,5 +1,6 @@
subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P)
subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO, &
X,EGHF,e,c,P)
! Search for GHF solutions
@ -29,7 +30,9 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas2,nBas2,nBas2,nBas2)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nBas2,nBas2,ncart)
! Local variables
@ -39,7 +42,6 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
logical :: unstab
integer :: guess
double precision,allocatable :: ERI_MO(:,:,:,:)
double precision,allocatable :: ERI_tmp(:,:,:,:)
double precision,allocatable :: Ca(:,:),Cb(:,:)
integer :: nS
@ -55,13 +57,12 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
double precision,allocatable :: R(:,:)
double precision,allocatable :: ExpR(:,:)
integer :: eig
double precision :: kick,step
integer :: ixyz
! Output variables
double precision,intent(out) :: EHF
double precision,intent(out) :: EGHF
double precision,intent(out) :: e(nBas2)
double precision,intent(inout):: c(nBas2,nBas2)
double precision,intent(out) :: P(nBas2,nBas2)
@ -80,8 +81,7 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
nS = (nO - nC)*(nV - nR)
allocate(ERI_MO(nBas2,nBas2,nBas2,nBas2),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS), &
R(nBas2,nBas2),ExpR(nBas2,nBas2))
allocate(Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS),R(nBas2,nBas2),ExpR(nBas2,nBas2))
!------------------!
! Search algorithm !
@ -99,7 +99,7 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
call wall_time(start_HF)
call GHF(.false.,maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P)
nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EGHF,e,c,P)
call wall_time(end_HF)
t_HF = end_HF - start_HF
@ -120,6 +120,12 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
Ca(:,:) = c(1:nBas,1:nBas2)
Cb(:,:) = c(nBas+1:nBas2,1:nBas2)
! Transform dipole-related integrals
do ixyz=1,ncart
call AOtoMO_GHF(nBas,nBas2,Ca,Cb,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
! 4-index transform
call AOtoMO_ERI_GHF(nBas,nBas2,Ca,Ca,ERI_AO,ERI_tmp)
@ -171,7 +177,7 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
write(*,'(1X,A40,1X)') 'Too bad, GHF solution is unstable!'
write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue:',Om(1),' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(GHF) = ',ENuc + EHF,' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(GHF) = ',ENuc + EGHF,' au'
write(*,*)
write(*,'(1X,A40,1X,A10)') 'Which one would you like to follow?','[Exit:0]'
read(*,*) eig
@ -184,20 +190,6 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
if(eig == 0) return
step = 1d0
! do mu=1,nBas2
! ia = 0
! do i=nC+1,nO
! kick = 0d0
! do a=nO+1,nBas2-nR
! ia = ia + 1
! kick = kick + AB(ia,eig)*c(mu,a)
! end do
! c(mu,i) = c(mu,i) + step*kick
! end do
! end do
R(:,:) = 0d0
ia = 0
do i=nC+1,nO
@ -211,13 +203,11 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
call matrix_exponential(nBas2,R,ExpR)
c = matmul(c,ExpR)
! call matout(nBas2,nBas2,matmul(transpose(ExpR),ExpR))
else
write(*,'(1X,A40,1X)') 'Well done, GHF solution is stable!'
write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(GHF) = ',ENuc + EHF,' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(GHF) = ',ENuc + EGHF,' au'
unstab = .false.

38
src/HF/RHF_search.f90
View File

@ -1,5 +1,6 @@
subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P)
subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO, &
X,ERHF,e,c,P)
! Search for RHF solutions
@ -27,7 +28,9 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart)
! Local variables
@ -37,7 +40,6 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
logical :: unstab
integer :: guess
double precision,allocatable :: ERI_MO(:,:,:,:)
integer :: nS
integer,parameter :: maxS = 20
@ -51,12 +53,12 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
double precision,allocatable :: R(:,:)
double precision,allocatable :: ExpR(:,:)
integer :: ixyz
integer :: eig
double precision :: kick,step
! Output variables
double precision,intent(out) :: EHF
double precision,intent(out) :: ERHF
double precision,intent(out) :: e(nBas)
double precision,intent(inout):: c(nBas,nBas)
double precision,intent(out) :: P(nBas,nBas)
@ -74,8 +76,7 @@ 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), &
R(nBas,nBas),ExpR(nBas,nBas))
allocate(Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS),R(nBas,nBas),ExpR(nBas,nBas))
!------------------!
! Search algorithm !
@ -92,7 +93,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
call wall_time(start_HF)
call RHF(.false.,maxSCF,thresh,max_diis,guess,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P)
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,e,c,P)
call wall_time(end_HF)
t_HF = end_HF - start_HF
@ -107,6 +108,9 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
write(*,*)
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
do ixyz=1,ncart
call AOtoMO(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
call wall_time(end_AOtoMO)
@ -144,7 +148,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
write(*,'(1X,A40,1X)') 'Too bad, RHF solution is unstable!'
write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue:',Om(1),' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(RHF) = ',ENuc + EHF,' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(RHF) = ',ENuc + ERHF,' au'
write(*,*)
write(*,'(1X,A40,1X,A10)') 'Which one would you like to follow?','[Exit:0]'
read(*,*) eig
@ -157,20 +161,6 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
if(eig == 0) return
step = 1d0
! do mu=1,nBas
! ia = 0
! do i=nC+1,nO
! kick = 0d0
! do a=nO+1,nBas-nR
! ia = ia + 1
! kick = kick + AB(ia,eig)*c(mu,a)
! end do
! c(mu,i) = c(mu,i) + step*kick
! end do
! end do
R(:,:) = 0d0
ia = 0
do i=nC+1,nO
@ -188,7 +178,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
write(*,'(1X,A40,1X)') 'Well done, RHF solution is stable!'
write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(RHF) = ',ENuc + EHF,' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(RHF) = ',ENuc + ERHF,' au'
unstab = .false.

55
src/HF/UHF_search.f90
View File

@ -1,5 +1,6 @@
subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P)
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_AO,dipole_int_aa,dipole_int_bb,X,EUHF,e,c,P)
! Search for UHF solutions
@ -28,7 +29,12 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_bb(nBas,nBas,ncart)
! Local variables
@ -38,9 +44,6 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
logical :: unstab
integer :: guess
double precision,allocatable :: ERI_aaaa(:,:,:,:)
double precision,allocatable :: ERI_aabb(:,:,:,:)
double precision,allocatable :: ERI_bbbb(:,:,:,:)
integer :: nS(nspin)
@ -56,12 +59,12 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
double precision,allocatable :: R(:,:)
double precision,allocatable :: ExpR(:,:)
integer :: ixyz
integer :: eig
double precision :: kick,step
! Output variables
double precision,intent(out) :: EHF
double precision,intent(out) :: EUHF
double precision,intent(out) :: e(nBas,nspin)
double precision,intent(inout):: c(nBas,nBas,nspin)
double precision,intent(out) :: P(nBas,nBas,nspin)
@ -84,7 +87,6 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(ERI_aaaa(nBas,nBas,nBas,nBas),ERI_aabb(nBas,nBas,nBas,nBas),ERI_bbbb(nBas,nBas,nBas,nBas))
allocate(Om_sc(nS_sc),A_sc(nS_sc,nS_sc),B_sc(nS_sc,nS_sc),AB_sc(nS_sc,nS_sc),R(nBas,nBas),ExpR(nBas,nBas))
!------------------!
@ -103,7 +105,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
call wall_time(start_HF)
call UHF(.false.,maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P)
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,e,c,P)
call wall_time(end_HF)
t_HF = end_HF - start_HF
@ -119,6 +121,13 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
! Transform dipole-related integrals
do ixyz=1,ncart
call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block
call AOtoMO_ERI_UHF(1,1,nBas,c,ERI_AO,ERI_aaaa)
@ -164,7 +173,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
write(*,'(1X,A40,1X)') 'Too bad, UHF solution is unstable!'
write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue:',Om_sc(1),' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(UHF) = ',ENuc + EHF,' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(UHF) = ',ENuc + EUHF,' au'
write(*,*)
write(*,'(1X,A40,1X,A10)') 'Which one would you like to follow?','[Exit:0]'
read(*,*) eig
@ -177,22 +186,8 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
if(eig == 0) return
step = 1d0
! Spin-up kick
! do mu=1,nBas
! ia = 0
! do i=nC(1)+1,nO(1)
! kick = 0d0
! do a=nO(1)+1,nBas-nR(1)
! ia = ia + 1
! kick = kick + AB_sc(ia,eig)*c(mu,a,1)
! end do
! c(mu,i,1) = c(mu,i,1) + step*kick
! end do
! end do
R(:,:) = 0d0
ia = 0
do i=nC(1)+1,nO(1)
@ -208,18 +203,6 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Spin-down kick
! do mu=1,nBas
! ia = nS_aa
! do i=nC(2)+1,nO(2)
! kick = 0d0
! do a=nO(2)+1,nBas-nR(2)
! ia = ia + 1
! kick = kick + AB_sc(ia,eig)*c(mu,a,2)
! end do
! c(mu,i,2) = c(mu,i,2) + step*kick
! end do
! end do
R(:,:) = 0d0
ia = nS_aa
do i=nC(2)+1,nO(2)
@ -237,7 +220,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
write(*,'(1X,A40,1X)') 'Well done, UHF solution is stable!'
write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om_sc(1),' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(UHF) = ',ENuc + EHF,' au'
write(*,'(1X,A40,1X,F15.10,A3)') 'E(UHF) = ',ENuc + EUHF,' au'
unstab = .false.

12
src/HF/print_GHF_spin.f90
View File

@ -116,8 +116,8 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S)
Sc_yx = Sc_yx - (+0.5d0,0.d0) * (Pab(j,i) + Pba(j,i)) * (0.d0,-0.5d0) * (Pab(i,j) - Pba(i,j))
enddo
enddo
write(*,'(A15,2F10.6)') ' < Sx Sy > = ',Sc_xy
write(*,'(A15,2F10.6)') ' < Sy Sx > = ',Sc_yx
write(*,'(A15,2F10.6)') ' < Sx.Sy > = ',Sc_xy
write(*,'(A15,2F10.6)') ' < Sy.Sx > = ',Sc_yx
Sc_xz = Sc_x * Sc_z
Sc_zx = Sc_x * Sc_z
@ -129,8 +129,8 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S)
Sc_zx = Sc_zx - (+0.5d0,0.d0) * (Pab(j,i) + Pba(j,i)) * (+0.5d0,0.d0) * (Paa(i,j) - Pbb(i,j))
enddo
enddo
write(*,'(A15,2F10.6)') ' < Sx Sz > = ',Sc_xz
write(*,'(A15,2F10.6)') ' < Sz Sx > = ',Sc_zx
write(*,'(A15,2F10.6)') ' < Sx.Sz > = ',Sc_xz
write(*,'(A15,2F10.6)') ' < Sz.Sx > = ',Sc_zx
Sc_yz = Sc_y * Sc_z
Sc_zy = Sc_y * Sc_z
@ -142,8 +142,8 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S)
Sc_zy = Sc_zy - (0.d0,-0.5d0) * (Pab(j,i) - Pba(j,i)) * (+0.5d0,0.d0) * (Paa(i,j) - Pbb(i,j))
enddo
enddo
write(*,'(A15,2F10.6)') ' < Sy Sz > = ',Sc_yz
write(*,'(A15,2F10.6)') ' < Sz Sy > = ', Sc_zy
write(*,'(A15,2F10.6)') ' < Sy.Sz > = ',Sc_yz
write(*,'(A15,2F10.6)') ' < Sz.Sy > = ', Sc_zy
write(*,*)

2
src/HF/print_RHF.f90
View File

@ -27,7 +27,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
double precision :: Gap
double precision :: S,S2
logical :: dump_orb = .false.
logical :: dump_orb = .true.
! HOMO and LUMO

8
src/LR/phLR.f90
View File

@ -60,6 +60,8 @@ subroutine phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call dgemm('N','N',nS,nS,nS,1d0,ApB,size(ApB,1),AmBSq,size(AmBSq,1),0d0,tmp,size(tmp,1))
call dgemm('N','N',nS,nS,nS,1d0,AmBSq,size(AmBSq,1),tmp,size(tmp,1),0d0,Z,size(Z,1))
! Z = matmul(AmBSq,matmul(ApB,AmBSq))
call diagonalize_matrix(nS,Z,Om)
if(minval(Om) < 0d0) &
@ -72,6 +74,12 @@ subroutine phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBIv,size(AmBIv,1),0d0,XmY,size(XmY,1))
call DA(nS,1d0*dsqrt(Om),XmY)
! XpY = matmul(transpose(Z),AmBSq)
! call DA(nS,1d0/sqrt(Om),XpY)
! XmY = matmul(transpose(Z),AmBIv)
! call DA(nS,1d0*sqrt(Om),XmY)
end if

5
src/QuAcK/GQuAcK.f90
View File

@ -139,7 +139,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
! Transform dipole-related integrals
do ixyz=1,ncart
call AOtoMO_GHF(nBas,nBas2,Ca,Cb,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
call AOtoMO_GHF(nBas,nBas2,Ca,Cb,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
! 4-index transform
@ -186,7 +186,8 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
call wall_time(start_stab)
call GHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EGHF,eHF,cHF,PHF)
nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO, &
X,EGHF,eHF,cHF,PHF)
call wall_time(end_stab)
t_stab = end_stab - start_stab

2
src/QuAcK/RQuAcK.f90
View File

@ -188,7 +188,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
call wall_time(start_stab)
call RHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF)
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,X,ERHF,eHF,cHF,PHF)
call wall_time(end_stab)
t_stab = end_stab - start_stab

7
src/QuAcK/UQuAcK.f90
View File

@ -161,8 +161,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
! Read and transform dipole-related integrals
do ixyz=1,ncart
call AOtoMO(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
call AOtoMO(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block
@ -205,7 +205,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
call wall_time(start_stab)
call UHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF)
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO, &
dipole_int_aa,dipole_int_bb,X,EUHF,eHF,cHF,PHF)
call wall_time(end_stab)