From 90508131e30c0b885962730330ff02c8fb81f8af Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 26 Oct 2023 09:35:48 +0200 Subject: [PATCH] working on GHF --- src/AOtoMO/Coulomb_matrix_AO_basis.f90 | 34 ----- src/AOtoMO/Coulomb_matrix_MO_basis.f90 | 26 ---- src/AOtoMO/Hartree_matrix_AO_basis.f90 | 7 +- src/AOtoMO/Hartree_matrix_MO_basis.f90 | 6 +- src/GF/qsGF2.f90 | 6 +- src/GF/qsUGF2.f90 | 6 +- src/GT/qsGTeh.f90 | 6 +- src/GT/qsGTpp.f90 | 6 +- src/GT/qsUGTpp.f90 | 6 +- src/GW/SRG_qsGW.f90 | 6 +- src/GW/qsGW.f90 | 6 +- src/GW/qsUGW.f90 | 6 +- src/HF/GHF.f90 | 179 ++++++++++++++----------- src/HF/RHF.f90 | 2 +- src/HF/RMOM.f90 | 2 +- src/HF/ROHF.f90 | 6 +- src/HF/UHF.f90 | 6 +- src/HF/mo_fock_exchange_potential.f90 | 41 ------ 18 files changed, 140 insertions(+), 217 deletions(-) delete mode 100644 src/AOtoMO/Coulomb_matrix_AO_basis.f90 delete mode 100644 src/AOtoMO/Coulomb_matrix_MO_basis.f90 delete mode 100644 src/HF/mo_fock_exchange_potential.f90 diff --git a/src/AOtoMO/Coulomb_matrix_AO_basis.f90 b/src/AOtoMO/Coulomb_matrix_AO_basis.f90 deleted file mode 100644 index 7268acb..0000000 --- a/src/AOtoMO/Coulomb_matrix_AO_basis.f90 +++ /dev/null @@ -1,34 +0,0 @@ -subroutine Coulomb_matrix_AO_basis(nBas,P,G,J) - -! Compute Coulomb matrix in the AO basis - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas - double precision,intent(in) :: P(nBas,nBas) - double precision,intent(in) :: G(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: mu,nu,la,si - -! Output variables - - double precision,intent(out) :: J(nBas,nBas) - - J = 0d0 - do si=1,nBas - do nu=1,nBas - do la=1,nBas - do mu=1,nBas - J(mu,nu) = J(mu,nu) + P(la,si)*G(mu,la,nu,si) - enddo - enddo - enddo - enddo - - -end subroutine diff --git a/src/AOtoMO/Coulomb_matrix_MO_basis.f90 b/src/AOtoMO/Coulomb_matrix_MO_basis.f90 deleted file mode 100644 index d64bfe1..0000000 --- a/src/AOtoMO/Coulomb_matrix_MO_basis.f90 +++ /dev/null @@ -1,26 +0,0 @@ -subroutine Coulomb_matrix_MO_basis(nBas,c,P,G,J) - -! Compute Coulomb matrix in the MO basis - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas - double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas) - double precision,intent(in) :: G(nBas,nBas,nBas,nBas) - -! Output variables - - double precision,intent(out) :: J(nBas,nBas) - -! Compute Hartree Hamiltonian in the AO basis - - call Coulomb_matrix_AO_basis(nBas,P,G,J) - -! Transform Coulomb matrix in the MO basis - - J = matmul(transpose(c),matmul(J,c)) - -end subroutine diff --git a/src/AOtoMO/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 index 0d0e5a7..87ef787 100644 --- a/src/AOtoMO/Hartree_matrix_AO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_AO_basis.f90 @@ -1,4 +1,4 @@ -subroutine Hartree_matrix_AO_basis(nBas,P,Hc,G,H) +subroutine Hartree_matrix_AO_basis(nBas,P,G,H) ! Compute Hartree matrix in the AO basis @@ -9,7 +9,7 @@ subroutine Hartree_matrix_AO_basis(nBas,P,Hc,G,H) integer,intent(in) :: nBas double precision,intent(in) :: P(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas),G(nBas,nBas,nBas,nBas) + double precision,intent(in) :: G(nBas,nBas,nBas,nBas) ! Local variables @@ -19,7 +19,8 @@ subroutine Hartree_matrix_AO_basis(nBas,P,Hc,G,H) double precision,intent(out) :: H(nBas,nBas) - H = Hc + H(:,:) = 0d0 + do mu=1,nBas do nu=1,nBas do la=1,nBas diff --git a/src/AOtoMO/Hartree_matrix_MO_basis.f90 b/src/AOtoMO/Hartree_matrix_MO_basis.f90 index 829c0e9..c1496b6 100644 --- a/src/AOtoMO/Hartree_matrix_MO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_MO_basis.f90 @@ -1,4 +1,4 @@ -subroutine Hartree_matrix_MO_basis(nBas,c,P,Hc,G,H) +subroutine Hartree_matrix_MO_basis(nBas,c,P,G,H) ! Compute Hartree matrix in the MO basis @@ -9,7 +9,7 @@ subroutine Hartree_matrix_MO_basis(nBas,c,P,Hc,G,H) integer,intent(in) :: nBas double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas),G(nBas,nBas,nBas,nBas) + double precision,intent(in) :: G(nBas,nBas,nBas,nBas) ! Output variables @@ -17,7 +17,7 @@ subroutine Hartree_matrix_MO_basis(nBas,c,P,Hc,G,H) ! Compute Hartree matrix in the AO basis - call Hartree_matrix_AO_basis(nBas,P,Hc,G,H) + call Hartree_matrix_AO_basis(nBas,P,G,H) ! Transform Hartree matrix in the MO basis diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index a0bd1c9..793435a 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -131,9 +131,9 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr nSCF = nSCF + 1 - ! Buid Coulomb matrix + ! Buid Hartree matrix - call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) ! Compute exchange part of the self-energy @@ -208,7 +208,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr EV = trace_matrix(nBas,matmul(P,V)) - ! Coulomb energy + ! Hartree energy EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 6bbfe7c..2d24d3b 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -146,10 +146,10 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f nSCF = nSCF + 1 - ! Buid Coulomb matrix + ! Buid Hartree matrix do is=1,nspin - call Coulomb_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),J(:,:,is)) + call Hartree_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),J(:,:,is)) end do ! Compute exchange part of the self-energy @@ -277,7 +277,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f EV(is) = trace_matrix(nBas,matmul(P(:,:,is),V(:,:))) end do - ! Coulomb energy + ! Hartree energy EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) EJ(2) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) & diff --git a/src/GT/qsGTeh.f90 b/src/GT/qsGTeh.f90 index 416a18e..8cd3381 100644 --- a/src/GT/qsGTeh.f90 +++ b/src/GT/qsGTeh.f90 @@ -158,9 +158,9 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d nSCF = nSCF + 1 - ! Buid Coulomb matrix + ! Buid Hartree matrix - call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) ! Compute exchange part of the self-energy @@ -240,7 +240,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d EV = trace_matrix(nBas,matmul(P,V)) - ! Coulomb energy + ! Hartree energy EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) diff --git a/src/GT/qsGTpp.f90 b/src/GT/qsGTpp.f90 index 8a86dea..67775fa 100644 --- a/src/GT/qsGTpp.f90 +++ b/src/GT/qsGTpp.f90 @@ -172,9 +172,9 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T nSCF = nSCF + 1 - ! Buid Coulomb matrix + ! Buid Hartree matrix - call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) ! Compute exchange part of the self-energy @@ -284,7 +284,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T EV = trace_matrix(nBas,matmul(P,V)) - ! Coulomb energy + ! Hartree energy EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) diff --git a/src/GT/qsUGTpp.f90 b/src/GT/qsUGTpp.f90 index 4595bf5..84f3ad3 100644 --- a/src/GT/qsUGTpp.f90 +++ b/src/GT/qsUGTpp.f90 @@ -172,9 +172,9 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & nSCF = nSCF + 1 -! Buid Coulomb matrix +! Buid Hartree matrix do ispin=1,nspin - call Coulomb_matrix_AO_basis(nBas,P(:,:,ispin),ERI_AO(:,:,:,:), & + call Hartree_matrix_AO_basis(nBas,P(:,:,ispin),ERI_AO(:,:,:,:), & J(:,:,ispin)) end do @@ -362,7 +362,7 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) end do -! Coulomb energy +! Hartree energy EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) EJ(2) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) & diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index c3bd4ac..3d4a7b7 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -154,9 +154,9 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE nSCF = nSCF + 1 - ! Buid Coulomb matrix + ! Buid Hartree matrix call wall_time(t1) - call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) ! Compute exchange part of the self-energy @@ -259,7 +259,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE EV = trace_matrix(nBas,matmul(P,V)) - ! Coulomb energy + ! Hartree energy EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index 37c60d8..652d972 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -158,9 +158,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop nSCF = nSCF + 1 - ! Buid Coulomb matrix + ! Buid Hartree matrix - call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) ! Compute exchange part of the self-energy @@ -239,7 +239,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop EV = trace_matrix(nBas,matmul(P,V)) - ! Coulomb energy + ! Hartree energy EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index e63e561..9c43f4c 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -167,10 +167,10 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, nSCF = nSCF + 1 - ! Buid Coulomb matrix + ! Buid Hartree matrix do is=1,nspin - call Coulomb_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),J(:,:,is)) + call Hartree_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),J(:,:,is)) end do ! Compute exchange part of the self-energy @@ -307,7 +307,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, EV(is) = trace_matrix(nBas,matmul(P(:,:,is),V(:,:))) end do - ! Coulomb energy + ! Hartree energy EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) EJ(2) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) & diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90 index e065dd9..abe46fa 100644 --- a/src/HF/GHF.f90 +++ b/src/HF/GHF.f90 @@ -1,5 +1,5 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P) + nBas,nBas2,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P) ! Perform unrestricted Hartree-Fock calculation @@ -11,10 +11,11 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, integer,intent(in) :: maxSCF integer,intent(in) :: max_diis integer,intent(in) :: guess_type - logical,intent(in) :: mix + double precision,intent(in) :: mix double precision,intent(in) :: level_shift double precision,intent(in) :: thresh integer,intent(in) :: nBas + integer,intent(in) :: nBas2 integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -34,9 +35,12 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, integer :: nSCF integer :: nBasSq + integer :: nBas2Sq + integer :: nOa + integer :: nOb integer :: n_diis - double precision :: conv - double precision :: rcond(nspin) + double precision :: Conv + double precision :: rcond double precision :: ET(nspin) double precision :: EV(nspin) double precision :: EJ(nsp) @@ -47,9 +51,10 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, double precision,allocatable :: Jaa(:,:),Jab(:,:),Jba(:,:),Jbb(:,:) double precision,allocatable :: Kaa(:,:),Kab(:,:),Kba(:,:),Kbb(:,:) double precision,allocatable :: Faa(:,:),Fab(:,:),Fba(:,:),Fbb(:,:) - double precision,allocatable :: Cp(:,:,:,:) + double precision,allocatable :: Paa(:,:),Pab(:,:),Pba(:,:),Pbb(:,:) double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) + double precision,allocatable :: Cp(:,:) double precision,allocatable :: err(:,:) double precision,allocatable :: err_diis(:,:) double precision,allocatable :: F_diis(:,:) @@ -60,9 +65,9 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Output variables double precision,intent(out) :: EHF - double precision,intent(out) :: e(nBas) - double precision,intent(out) :: c(nBas,nBas) - double precision,intent(out) :: P(nBas,nBas) + double precision,intent(out) :: e(nBas2) + double precision,intent(out) :: C(nBas2,nBas2) + double precision,intent(out) :: P(nBas2,nBas2) ! Hello world @@ -75,28 +80,44 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Useful stuff nBasSq = nBas*nBas + nBas2Sq = nBas2*nBas2 + + nOa = nO(1) + nOb = nO(2) ! Memory allocation - allocate(Caa(nBas,nBas),Jab(F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), & - K(nBas,nBas,nspin),err(nBas,nBas,nspin),cp(nBas,nBas,nspin), & - err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) + allocate(Caa(nBas,nBas),Cab(nBas,nBas),Cba(nBas,nBas),Cbb(nBas,nBas), & + Jaa(nBas,nBas),Jab(nBas,nBas),Jba(nBas,nBas),Jbb(nBas,nBas), & + Kaa(nBas,nBas),Kab(nBas,nBas),Kba(nBas,nBas),Kbb(nBas,nBas), & + Faa(nBas,nBas),Fab(nBas,nBas),Fba(nBas,nBas),Fbb(nBas,nBas), & + Paa(nBas,nBas),Pab(nBas,nBas),Pba(nBas,nBas),Pbb(nBas,nBas), & + F(nBas2,nBas2),Fp(nBas2,nBas2),Cp(nBas2,nBas2),err(nBas2,nBas2), & + err_diis(nBas2Sq,max_diis),F_diis(nBas2Sq,max_diis)) ! Guess coefficients and demsity matrices - do ispin=1,nspin - call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin)) - P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) - end do +! do ispin=1,nspin +! call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin)) +! P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) +! end do ! Initialization nSCF = 0 - conv = 1d0 + Conv = 1d0 - n_diis = 0 - F_diis(:,:,:) = 0d0 - err_diis(:,:,:) = 0d0 + n_diis = 0 + F_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 + +! Construct super overlap matrix + +! TO DO + +! Construct super orthogonalization matrix + +! TO DO !------------------------------------------------------------------------ ! Main SCF loop @@ -108,88 +129,90 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, '|','#','|','E(UHF)','|','Ex(UHF)','|','Conv','|' write(*,*)'----------------------------------------------------------' - do while(conv > thresh .and. nSCF < maxSCF) + do while(Conv > thresh .and. nSCF < maxSCF) ! Increment nSCF = nSCF + 1 -! Build Coulomb repulsion +! Build individual Hartree matrices - do ispin=1,nspin - call Coulomb_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) - end do + call Hartree_matrix_AO_basis(nBas,Paa,ERI,Jaa) + call Hartree_matrix_AO_basis(nBas,Pab,ERI,Jab) + call Hartree_matrix_AO_basis(nBas,Pba,ERI,Jba) + call Hartree_matrix_AO_basis(nBas,Pbb,ERI,Jbb) -! Compute exchange potential +! Compute individual exchange matrices - do ispin=1,nspin - call exchange_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),K(:,:,ispin)) - end do + call exchange_matrix_AO_basis(nBas,Paa,ERI,Kaa) + call exchange_matrix_AO_basis(nBas,Pab,ERI,Kab) + call exchange_matrix_AO_basis(nBas,Pba,ERI,Kba) + call exchange_matrix_AO_basis(nBas,Pbb,ERI,Kbb) -! Build Fock operator +! Build individual Fock matrices - do ispin=1,nspin - F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin) - end do + Faa(:,:) = Hc(:,:) + Jaa(:,:) + Jab(:,:) + Kaa(:,:) + Fab(:,:) = Hc(:,:) + Jab(:,:) + Jba(:,:) + Kab(:,:) + Fba(:,:) = Hc(:,:) + Jba(:,:) + Jab(:,:) + Kba(:,:) + Fbb(:,:) = Hc(:,:) + Jbb(:,:) + Jba(:,:) + Kbb(:,:) + +! Build super Fock matrix + + F( 1:nBas , 1:nBas ) = Faa(1:nBas,1:nBas) + F(nBas+1:nBas2, 1:nBas ) = Fab(1:nBas,1:nBas) + F( 1:nBas ,nBas+1:nBas2) = Fba(1:nBas,1:nBas) + F(nBas+1:nBas2,nBas+1:nBas2) = Fbb(1:nBas,1:nBas) ! Check convergence - do ispin=1,nspin - err(:,:,ispin) = matmul(F(:,:,ispin),matmul(P(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),P(:,:,ispin)),F(:,:,ispin)) - end do +! err(:,:) = matmul(F(:,:),matmul(P(:,:),S(:,:))) - matmul(matmul(S(:,:),P(:,:)),F(:,:)) - if(nSCF > 1) conv = maxval(abs(err(:,:,:))) +! if(nSCF > 1) conv = maxval(abs(err(:,:))) ! DIIS extrapolation - if(max_diis > 1) then +! if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - do ispin=1,nspin - if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,ispin), & - F_diis(:,1:n_diis,ispin),err(:,:,ispin),F(:,:,ispin)) - end do + call DIIS_extrapolation(rcond,nBas2Sq,nBas2Sq,n_diis,err_diis(:,1:n_diis),F_diis(:,1:n_diis),err,F) - end if +! end if ! Level-shifting if(level_shift > 0d0 .and. Conv > thresh) then - do ispin=1,nspin - call level_shifting(level_shift,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin)) - end do + call level_shifting(level_shift,nBas,nOa+nOb,S,C,F) end if ! Transform Fock matrix in orthogonal basis - do ispin=1,nspin - Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:))) - end do + Fp(:,:) = matmul(transpose(X),matmul(F,X)) ! Diagonalize Fock matrix to get eigenvectors and eigenvalues - cp(:,:,:) = Fp(:,:,:) - do ispin=1,nspin - call diagonalize_matrix(nBas,cp(:,:,ispin),e(:,ispin)) - end do + Cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas2,Cp,e) ! Back-transform eigenvectors in non-orthogonal basis - do ispin=1,nspin - c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin)) - end do + C(:,:) = matmul(X,Cp) + +! Form individual coefficient matrices + +! TO DO ! Mix guess for UHF solution in singlet states - if(mix .and. nSCF == 1) call mix_guess(nBas,nO,c) +! if(nSCF == 1) call mix_guess(nBas,nO,mix,c) -! Compute density matrix +! Compute individual density matrices - do ispin=1,nspin - P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) - end do + Paa(:,:) = matmul(Caa(:,1:nOa),transpose(Caa(:,1:nOa))) + Pab(:,:) = matmul(Cab(:,1:nOb),transpose(Cab(:,1:nOb))) + Pba(:,:) = matmul(Cba(:,1:nOa),transpose(Cba(:,1:nOa))) + Pbb(:,:) = matmul(Cbb(:,1:nOb),transpose(Cbb(:,1:nOb))) !------------------------------------------------------------------------ ! Compute UHF energy @@ -197,36 +220,36 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Kinetic energy - do ispin=1,nspin - ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) - end do +! do ispin=1,nspin +! ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) +! end do ! Potential energy - do ispin=1,nspin - EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) - end do +! do ispin=1,nspin +! EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) +! end do -! Coulomb energy +! Hartree energy - EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) - EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) - EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) +! EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) +! EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) +! EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) ! Exchange energy - do ispin=1,nspin - Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin))) - end do +! do ispin=1,nspin +! Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin))) +! end do ! Total energy - EHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) +! EHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) ! Dump results - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',conv,'|' +! write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & +! '|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',conv,'|' end do write(*,*)'----------------------------------------------------------' @@ -250,7 +273,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Compute final UHF energy - call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) +! call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) +! call print_GHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) end subroutine diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 5f61db6..4fae804 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -109,7 +109,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc ! Build Fock matrix - call Coulomb_matrix_AO_basis(nBas,P,ERI,J) + call Hartree_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) diff --git a/src/HF/RMOM.f90 b/src/HF/RMOM.f90 index 0106bf9..9d48dd6 100644 --- a/src/HF/RMOM.f90 +++ b/src/HF/RMOM.f90 @@ -99,7 +99,7 @@ subroutine RMOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P) ! Build Fock matrix - call Coulomb_matrix_AO_basis(nBas,P,ERI,J) + call Hartree_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) F(:,:) = Hc(:,:) + J(:,:) + 0.5*K(:,:) diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index ee24b32..03d6c1a 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -115,10 +115,10 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc nSCF = nSCF + 1 -! Build Coulomb repulsion +! Build Hartree repulsion do ispin=1,nspin - call Coulomb_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) + call Hartree_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) end do ! Compute exchange potential @@ -195,7 +195,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) end do -! Coulomb energy +! Hartree energy EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 87f37e2..8435170 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -112,10 +112,10 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, nSCF = nSCF + 1 -! Build Coulomb repulsion +! Build Hartree repulsion do ispin=1,nspin - call Coulomb_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) + call Hartree_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) end do ! Compute exchange potential @@ -205,7 +205,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) end do -! Coulomb energy +! Hartree energy EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) diff --git a/src/HF/mo_fock_exchange_potential.f90 b/src/HF/mo_fock_exchange_potential.f90 deleted file mode 100644 index c30357e..0000000 --- a/src/HF/mo_fock_exchange_potential.f90 +++ /dev/null @@ -1,41 +0,0 @@ -subroutine mo_fock_exchange_potential(nBas,c,P,ERI,Vx) - -! Compute the exchange potential in the MO basis - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas - double precision,intent(in) :: c(nBas,nBas) - double precision,intent(in) :: P(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: mu,nu - integer :: q - double precision,allocatable :: Fx(:,:) - -! Output variables - - double precision,intent(out) :: Vx(nBas) - -! Compute Vx - - allocate(Fx(nBas,nBas)) - call exchange_matrix_AO_basis(nBas,P,ERI,Fx) - - Vx(:) = 0d0 - do q=1,nBas - do mu=1,nBas - do nu=1,nBas - Vx(q) = Vx(q) + c(mu,q)*Fx(mu,nu)*c(nu,q) - end do - end do - end do - - deallocate(Fx) - -end subroutine