From f89d54c59b6014adc79b91685ead72a6277e9d99 Mon Sep 17 00:00:00 2001 From: pfloos Date: Tue, 28 Nov 2023 10:40:15 +0100 Subject: [PATCH] correct big bug in HF search --- src/GW/RGW.f90 | 18 ++++++------- src/HF/GHF_search.f90 | 44 ++++++++++++------------------- src/HF/RHF_search.f90 | 38 ++++++++++----------------- src/HF/UHF_search.f90 | 55 ++++++++++++++------------------------- src/HF/print_GHF_spin.f90 | 12 ++++----- src/HF/print_RHF.f90 | 2 +- src/LR/phLR.f90 | 8 ++++++ src/QuAcK/GQuAcK.f90 | 5 ++-- src/QuAcK/RQuAcK.f90 | 2 +- src/QuAcK/UQuAcK.f90 | 7 ++--- 10 files changed, 82 insertions(+), 109 deletions(-) diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index f836401..c137e74 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -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 diff --git a/src/HF/GHF_search.f90 b/src/HF/GHF_search.f90 index 796202a..aa7362f 100644 --- a/src/HF/GHF_search.f90 +++ b/src/HF/GHF_search.f90 @@ -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. diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index a60f347..d44e555 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -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. diff --git a/src/HF/UHF_search.f90 b/src/HF/UHF_search.f90 index 708db73..6b16f62 100644 --- a/src/HF/UHF_search.f90 +++ b/src/HF/UHF_search.f90 @@ -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. diff --git a/src/HF/print_GHF_spin.f90 b/src/HF/print_GHF_spin.f90 index ef2a824..ea436c4 100644 --- a/src/HF/print_GHF_spin.f90 +++ b/src/HF/print_GHF_spin.f90 @@ -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(*,*) diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 3a06d31..89b93ec 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -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 diff --git a/src/LR/phLR.f90 b/src/LR/phLR.f90 index 227582e..ca49344 100644 --- a/src/LR/phLR.f90 +++ b/src/LR/phLR.f90 @@ -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 diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index aaadb3a..6060939 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -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 diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index dc44b33..5651d2f 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -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 diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 index 228da29..9bf04f4 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -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)