From 94cc75309efc9f0f0a5045f5fbdc36781bdf2704 Mon Sep 17 00:00:00 2001 From: pfloos Date: Wed, 8 Nov 2023 22:58:55 +0100 Subject: [PATCH] AO <-> MO modifs --- src/AOtoMO/AOtoMO_transform.f90 | 7 ++- src/AOtoMO/AOtoMO_transform_GHF.f90 | 4 +- src/AOtoMO/MOtoAO_transform.f90 | 10 +-- src/GF/qsGF2.f90 | 10 ++- src/GF/qsUGF2.f90 | 10 ++- src/GT/qsGTeh.f90 | 10 ++- src/GT/qsGTpp.f90 | 10 ++- src/GT/qsUGTpp.f90 | 10 ++- src/GW/SRG_qsGW.f90 | 36 +++++------ src/GW/print_qsGGW.f90 | 13 ++-- src/GW/print_qsGW.f90 | 13 ++-- src/GW/print_qsUGW.f90 | 25 ++++---- src/GW/qsGGW.f90 | 94 ++++++++++++++--------------- src/GW/qsGW.f90 | 24 +++----- src/GW/qsUGW.f90 | 70 +++++++++------------ src/HF/GHF.f90 | 6 +- src/HF/ROHF_fock_matrix.f90 | 4 +- src/HF/UHF.f90 | 10 +-- src/HF/print_GHF.f90 | 8 +-- src/HF/print_RHF.f90 | 2 +- src/HF/print_UHF.f90 | 2 +- src/QuAcK/RQuAcK.f90 | 3 +- src/QuAcK/UQuAcK.f90 | 7 +-- 23 files changed, 183 insertions(+), 205 deletions(-) diff --git a/src/AOtoMO/AOtoMO_transform.f90 b/src/AOtoMO/AOtoMO_transform.f90 index 410d3ab..5c4b710 100644 --- a/src/AOtoMO/AOtoMO_transform.f90 +++ b/src/AOtoMO/AOtoMO_transform.f90 @@ -1,4 +1,4 @@ -subroutine AOtoMO_transform(nBas,c,A) +subroutine AOtoMO_transform(nBas,c,A,B) ! Perform AO to MO transformation of a matrix A for given coefficients c @@ -8,11 +8,12 @@ subroutine AOtoMO_transform(nBas,c,A) integer,intent(in) :: nBas double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: A(nBas,nBas) ! Output variables - double precision,intent(inout):: A(nBas,nBas) + double precision,intent(out) :: B(nBas,nBas) - A = matmul(transpose(c),matmul(A,c)) + B = matmul(transpose(c),matmul(A,c)) end subroutine diff --git a/src/AOtoMO/AOtoMO_transform_GHF.f90 b/src/AOtoMO/AOtoMO_transform_GHF.f90 index ac022f7..67fc4be 100644 --- a/src/AOtoMO/AOtoMO_transform_GHF.f90 +++ b/src/AOtoMO/AOtoMO_transform_GHF.f90 @@ -10,11 +10,11 @@ subroutine AOtoMO_transform_GHF(nBas,nBas2,Ca,Cb,A,B) integer,intent(in) :: nBas2 double precision,intent(in) :: Ca(nBas,nBas2) double precision,intent(in) :: Cb(nBas,nBas2) - double precision,intent(inout):: A(nBas,nBas) + double precision,intent(in) :: A(nBas,nBas) ! Output variables - double precision,intent(inout):: B(nBas2,nBas2) + double precision,intent(out) :: B(nBas2,nBas2) B = matmul(transpose(Ca),matmul(A,Ca)) & + matmul(transpose(Cb),matmul(A,Cb)) diff --git a/src/AOtoMO/MOtoAO_transform.f90 b/src/AOtoMO/MOtoAO_transform.f90 index 005c84f..923f55e 100644 --- a/src/AOtoMO/MOtoAO_transform.f90 +++ b/src/AOtoMO/MOtoAO_transform.f90 @@ -1,4 +1,4 @@ -subroutine MOtoAO_transform(nBas,S,c,A) +subroutine MOtoAO_transform(nBas,S,c,B,A) ! Perform MO to AO transformation of a matrix A for a given metric S ! and coefficients c @@ -8,7 +8,9 @@ subroutine MOtoAO_transform(nBas,S,c,A) ! Input variables integer,intent(in) :: nBas - double precision,intent(in) :: S(nBas,nBas),c(nBas,nBas) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: B(nBas,nBas) ! Local variables @@ -16,12 +18,12 @@ subroutine MOtoAO_transform(nBas,S,c,A) ! Output variables - double precision,intent(inout):: A(nBas,nBas) + double precision,intent(out) :: A(nBas,nBas) ! Memory allocation allocate(Sc(nBas,nBas)) Sc = matmul(S,c) - A = matmul(Sc,matmul(A,transpose(Sc))) + A = matmul(Sc,matmul(B,transpose(Sc))) end subroutine diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index 793435a..76e0542 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -73,7 +73,6 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr double precision,allocatable :: K(:,:) double precision,allocatable :: SigC(:,:) double precision,allocatable :: SigCp(:,:) - double precision,allocatable :: SigCm(:,:) double precision,allocatable :: Z(:) double precision,allocatable :: error(:,:) @@ -104,7 +103,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr ! Memory allocation allocate(eGF(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & + J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -157,10 +156,9 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr ! Make correlation self-energy Hermitian and transform it back to AO basis - SigCp = 0.5d0*(SigC + transpose(SigC)) - SigCm = 0.5d0*(SigC - transpose(SigC)) + SigC = 0.5d0*(SigC + transpose(SigC)) - call MOtoAO_transform(nBas,S,c,SigCp) + call MOtoAO_transform(nBas,S,c,SigC,SigCp) ! Solve the quasi-particle equation @@ -253,7 +251,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,error,error_diis,F_diis) ! Perform BSE calculation diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 2d24d3b..d53b753 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -84,7 +84,6 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f double precision,allocatable :: K(:,:,:) double precision,allocatable :: SigC(:,:,:) double precision,allocatable :: SigCp(:,:,:) - double precision,allocatable :: SigCm(:,:,:) double precision,allocatable :: Z(:,:) double precision,allocatable :: error(:,:,:) @@ -120,7 +119,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f allocate(eGF2(nBas,nspin),eOld(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), & Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), & - SigCm(nBas,nBas,nspin),Z(nBas,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), & + Z(nBas,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), & F_diis(nBasSq,max_diis,nspin)) ! Initialization @@ -191,12 +190,11 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f ! Make correlation self-energy Hermitian and transform it back to AO basis do is=1,nspin - SigCp(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is))) - SigCm(:,:,is) = 0.5d0*(SigC(:,:,is) - transpose(SigC(:,:,is))) + SigC(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is))) end do do is=1,nspin - call MOtoAO_transform(nBas,S,c(:,:,is),SigCp(:,:,is)) + call MOtoAO_transform(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) end do ! Solve the quasi-particle equation @@ -326,7 +324,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f ! Deallocate memory - deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis) + deallocate(cp,P,F,Fp,J,K,SigC,SigCp,Z,error,error_diis,F_diis) ! Perform BSE calculation diff --git a/src/GT/qsGTeh.f90 b/src/GT/qsGTeh.f90 index 8cd3381..5c7690b 100644 --- a/src/GT/qsGTeh.f90 +++ b/src/GT/qsGTeh.f90 @@ -93,7 +93,6 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d double precision,allocatable :: K(:,:) double precision,allocatable :: Sig(:,:) double precision,allocatable :: Sigp(:,:) - double precision,allocatable :: Sigm(:,:) double precision,allocatable :: Z(:) double precision,allocatable :: error(:,:) @@ -131,7 +130,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! Memory allocation allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Sigm(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & + J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -189,10 +188,9 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! Make correlation self-energy Hermitian and transform it back to AO basis - Sigp = 0.5d0*(Sig + transpose(Sig)) - Sigm = 0.5d0*(Sig - transpose(Sig)) + Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO_transform(nBas,S,c,Sigp) + call MOtoAO_transform(nBas,S,c,Sig,Sigp) ! Solve the quasi-particle equation @@ -278,7 +276,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Sigm,Z,Om,XpY,XmY,rhoL,rhoR,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,error,error_diis,F_diis) ! Perform BSE calculation diff --git a/src/GT/qsGTpp.f90 b/src/GT/qsGTpp.f90 index 67775fa..7c6a481 100644 --- a/src/GT/qsGTpp.f90 +++ b/src/GT/qsGTpp.f90 @@ -92,7 +92,6 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T double precision,allocatable :: K(:,:) double precision,allocatable :: Sig(:,:) double precision,allocatable :: Sigp(:,:) - double precision,allocatable :: Sigm(:,:) double precision,allocatable :: Z(:) double precision,allocatable :: error(:,:) @@ -138,7 +137,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T ! Memory allocation allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Sigm(nBas,nBas),Z(nBas), & + J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & @@ -233,10 +232,9 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T ! Make correlation self-energy Hermitian and transform it back to AO basis - Sigp = 0.5d0*(Sig + transpose(Sig)) - Sigm = 0.5d0*(Sig - transpose(Sig)) + Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO_transform(nBas,S,c,Sigp) + call MOtoAO_transform(nBas,S,c,Sig,Sigp) ! Solve the quasi-particle equation @@ -322,7 +320,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Sigm,Z,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis) ! Perform BSE calculation diff --git a/src/GT/qsUGTpp.f90 b/src/GT/qsUGTpp.f90 index 84f3ad3..4f35488 100644 --- a/src/GT/qsUGTpp.f90 +++ b/src/GT/qsUGTpp.f90 @@ -91,7 +91,6 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & double precision,allocatable :: K(:,:,:) double precision,allocatable :: SigT(:,:,:) double precision,allocatable :: SigTp(:,:,:) - double precision,allocatable :: SigTm(:,:,:) double precision,allocatable :: Z(:,:) double precision,allocatable :: eGT(:,:) double precision,allocatable :: eOld(:,:) @@ -131,7 +130,7 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Memory allocation - allocate(SigT(nBas,nbas,nspin),SigTp(nBas,nbas,nspin),SigTm(nBas,nbas,nspin), & + allocate(SigT(nBas,nbas,nspin),SigTp(nBas,nbas,nspin), & Z(nBas,nspin),eGT(nBas,nspin),eOld(nBas,nspin), & error_diis(nBas,max_diis,nspin),e_diis(nBasSq,max_diis,nspin), & F_diis(nBasSq,max_diis,nspin),error(nBas,nBas,nspin),& @@ -274,12 +273,11 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ! Make correlation self-energy Hermitian and transform it back to AO basis do ispin=1,nspin - SigTp(:,:,ispin) = 0.5d0*(SigT(:,:,ispin) + transpose(SigT(:,:,ispin))) - SigTm(:,:,ispin) = 0.5d0*(SigT(:,:,ispin) - transpose(SigT(:,:,ispin))) + SigT(:,:,ispin) = 0.5d0*(SigT(:,:,ispin) + transpose(SigT(:,:,ispin))) end do do ispin=1,nspin - call MOtoAO_transform(nBas,S,c(:,:,ispin),SigTp(:,:,ispin)) + call MOtoAO_transform(nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin)) end do ! Solve the quasi-particle equation @@ -409,6 +407,6 @@ write(*,*) 'EcGM', EcGM(1) Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb) - deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,SigTm,Z,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,Z,error,error_diis,F_diis) end subroutine diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index 94eb373..5be6bfd 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -77,10 +77,10 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE double precision,allocatable :: F_diis(:,:) double precision,allocatable :: Aph(:,:) double precision,allocatable :: Bph(:,:) - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:) double precision,allocatable :: c(:,:) double precision,allocatable :: cp(:,:) double precision,allocatable :: eGW(:) @@ -91,6 +91,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE double precision,allocatable :: J(:,:) double precision,allocatable :: K(:,:) double precision,allocatable :: SigC(:,:) + double precision,allocatable :: SigCp(:,:) double precision,allocatable :: Z(:) double precision,allocatable :: error(:,:) @@ -127,9 +128,9 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE ! Memory allocation - allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),Z(nBas),Aph(nS,nS),Bph(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & - rho_RPA(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & + J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas),Aph(nS,nS),Bph(nS,nS), & + Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -169,9 +170,8 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE call wall_time(tao1) - dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:) do ixyz=1,ncart - call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) end do call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO) @@ -187,40 +187,39 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) - call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) call wall_time(tlr2) tlr = tlr + tlr2 -tlr1 - if(print_W) call print_excitation_energies('phRPA@SRG-qsGW',ispin,nS,OmRPA) + if(print_W) call print_excitation_energies('phRPA@SRG-qsGW',ispin,nS,Om) ! Compute correlation part of the self-energy call wall_time(tex1) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) + call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho) call wall_time(tex2) tex=tex+tex2-tex1 call wall_time(tsrg1) - call SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z) + call SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) call wall_time(tsrg2) tsrg = tsrg + tsrg2 -tsrg1 - ! Make correlation self-energy Hermitian and transform it back to AO basis call wall_time(tmo1) - call MOtoAO_transform(nBas,S,c,SigC) + call MOtoAO_transform(nBas,S,c,SigC,SigCp) call wall_time(tmo2) tmo = tmo + tmo2 - tmo1 ! Solve the quasi-particle equation - F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigC(:,:) + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:) ! Compute commutator and convergence criteria @@ -241,7 +240,8 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE cp(:,:) = Fp(:,:) call diagonalize_matrix(nBas,cp,eGW) c = matmul(X,cp) - SigC = matmul(transpose(c),matmul(SigC,c)) + + call AOtoMO_transform(nBas,c,SigCp,SigC) ! Compute new density matrix in the AO basis @@ -309,7 +309,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,SigC,Z,Om,XpY,XmY,rho,error,error_diis,F_diis) ! Perform BSE calculation diff --git a/src/GW/print_qsGGW.f90 b/src/GW/print_qsGGW.f90 index 2ef25ad..b244d5b 100644 --- a/src/GW/print_qsGGW.f90 +++ b/src/GW/print_qsGGW.f90 @@ -29,6 +29,7 @@ subroutine print_qsGGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E ! Local variables + logical :: dump_orb = .false. integer :: p,ixyz,HOMO,LUMO double precision :: Gap double precision,external :: trace_matrix @@ -103,11 +104,13 @@ subroutine print_qsGGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E write(*,'(A50)') '-----------------------------------------' write(*,*) - write(*,'(A50)') '---------------------------------------' - write(*,'(A32)') ' qsGGW orbital coefficients' - write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) - write(*,*) + if(dump_orb) then + write(*,'(A50)') '---------------------------------------' + write(*,'(A32)') ' qsGGW orbital coefficients' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,nBas,c) + write(*,*) + end if write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGGW orbital energies' write(*,'(A50)') '---------------------------------------' diff --git a/src/GW/print_qsGW.f90 b/src/GW/print_qsGW.f90 index bf673c6..3d7f7bb 100644 --- a/src/GW/print_qsGW.f90 +++ b/src/GW/print_qsGW.f90 @@ -29,6 +29,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex ! Local variables + logical :: dump_orb = .false. integer :: p,ixyz,HOMO,LUMO double precision :: Gap double precision,external :: trace_matrix @@ -103,11 +104,13 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex write(*,'(A50)') '-----------------------------------------' write(*,*) - write(*,'(A50)') '---------------------------------------' - write(*,'(A32)') ' qsGW MO coefficients' - write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) - write(*,*) + if(dump_orb) then + write(*,'(A50)') '---------------------------------------' + write(*,'(A32)') ' qsGW MO coefficients' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,nBas,c) + write(*,*) + end if write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGW MO energies' write(*,'(A50)') '---------------------------------------' diff --git a/src/GW/print_qsUGW.f90 b/src/GW/print_qsUGW.f90 index 1e65b2b..4523dc2 100644 --- a/src/GW/print_qsUGW.f90 +++ b/src/GW/print_qsUGW.f90 @@ -31,6 +31,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov, & ! Local variables + logical :: dump_orb = .false. integer :: p integer :: ispin,ixyz double precision :: HOMO(nspin) @@ -158,17 +159,19 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov, & ! Print orbitals - write(*,'(A50)') '-----------------------------------------' - write(*,'(A50)') 'qsUGW spin-up orbital coefficients ' - write(*,'(A50)') '-----------------------------------------' - call matout(nBas,nBas,cGW(:,:,1)) - write(*,*) - write(*,'(A50)') '-----------------------------------------' - write(*,'(A50)') 'qsUGW spin-down orbital coefficients ' - write(*,'(A50)') '-----------------------------------------' - call matout(nBas,nBas,cGW(:,:,2)) - write(*,*) + if(dump_orb) then + write(*,'(A50)') '-----------------------------------------' + write(*,'(A50)') 'qsUGW spin-up orbital coefficients ' + write(*,'(A50)') '-----------------------------------------' + call matout(nBas,nBas,cGW(:,:,1)) + write(*,*) + write(*,'(A50)') '-----------------------------------------' + write(*,'(A50)') 'qsUGW spin-down orbital coefficients ' + write(*,'(A50)') '-----------------------------------------' + call matout(nBas,nBas,cGW(:,:,2)) + write(*,*) + end if - endif + end if end subroutine diff --git a/src/GW/qsGGW.f90 b/src/GW/qsGGW.f90 index b00273e..2404dd2 100644 --- a/src/GW/qsGGW.f90 +++ b/src/GW/qsGGW.f90 @@ -65,7 +65,7 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do double precision :: EK,EKaaaa,EKabba,EKbaab,EKbbbb double precision :: EqsGW double precision :: EcRPA - double precision :: EcBSE(nspin) + double precision :: EcBSE double precision :: EcGM double precision :: Conv double precision :: rcond @@ -93,7 +93,6 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do double precision,allocatable :: C(:,:) double precision,allocatable :: Cp(:,:) double precision,allocatable :: eGW(:) - double precision,allocatable :: eOld(:) double precision,allocatable :: P(:,:) double precision,allocatable :: F(:,:) double precision,allocatable :: H(:,:) @@ -101,15 +100,16 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do double precision,allocatable :: X(:,:) double precision,allocatable :: Fp(:,:) double precision,allocatable :: SigC(:,:) + double precision,allocatable :: SigCp(:,:) double precision,allocatable :: Z(:) double precision,allocatable :: err(:,:) ! Hello world write(*,*) - write(*,*)'************************************************' - write(*,*)'| Self-consistent qsGW calculation |' - write(*,*)'************************************************' + write(*,*)'********************************' + write(*,*)'| Generalized qsGW Calculation |' + write(*,*)'********************************' write(*,*) ! Warning @@ -138,17 +138,15 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do ! Memory allocation - - allocate(P(nBas2,nBas2),Jaa(nBas,nBas),Jbb(nBas,nBas), & + allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),P(nBas2,nBas2),Jaa(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),C(nBas2,nBas2),Cp(nBas2,nBas2), & H(nBas2,nBas2),S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2), & - err_diis(nBas2Sq,max_diis),F_diis(nBas2Sq,max_diis)) - - allocate(eGW(nBas2),eOld(nBas2),SigC(nBas2,nBas2),Z(nBas2), & - Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas2,nBas2,nS)) + err_diis(nBas2Sq,max_diis),F_diis(nBas2Sq,max_diis), & + eGW(nBas2),SigC(nBas2,nBas2),SigCp(nBas,nBas),Z(nBas2),Aph(nS,nS),Bph(nS,nS), & + Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas2,nBas2,nS)) ! Initialization @@ -158,7 +156,6 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do Conv = 1d0 P(:,:) = PHF(:,:) eGW(:) = eHF(:) - eOld(:) = eHF(:) c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 err_diis(:,:) = 0d0 @@ -182,6 +179,15 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do H( 1:nBas , 1:nBas ) = Hc(1:nBas,1:nBas) H(nBas+1:nBas2,nBas+1:nBas2) = Hc(1:nBas,1:nBas) +! Construct super density matrix + + P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO))) + + Paa(:,:) = P( 1:nBas , 1:nBas ) + Pab(:,:) = P( 1:nBas ,nBas+1:nBas2) + Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) + Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) + !------------------------------------------------------------------------ ! Main loop !------------------------------------------------------------------------ @@ -204,9 +210,23 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do call exchange_matrix_AO_basis(nBas,Pab,ERI_AO,Kba) call exchange_matrix_AO_basis(nBas,Pbb,ERI_AO,Kbb) +! Build individual Fock matrices + + Faa(:,:) = Hc(:,:) + Jaa(:,:) + Jbb(:,:) + Kaa(:,:) + Fab(:,:) = + Kab(:,:) + Fba(:,:) = + Kba(:,:) + Fbb(:,:) = Hc(:,:) + Jbb(:,:) + Jaa(:,:) + Kbb(:,:) + +! Build super Fock matrix + + F( 1:nBas , 1:nBas ) = Faa(1:nBas,1:nBas) + F( 1:nBas ,nBas+1:nBas2) = Fab(1:nBas,1:nBas) + F(nBas+1:nBas2, 1:nBas ) = Fba(1:nBas,1:nBas) + F(nBas+1:nBas2,nBas+1:nBas2) = Fbb(1:nBas,1:nBas) + ! AO to MO transformation of two-electron integrals - allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),ERI_tmp(nBas2,nBas2,nBas2,nBas2)) + allocate(ERI_tmp(nBas2,nBas2,nBas2,nBas2)) Ca(:,:) = C(1:nBas,1:nBas2) Cb(:,:) = C(nBas+1:nBas2,1:nBas2) @@ -227,7 +247,7 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Cb,Cb,Cb,ERI_AO,ERI_tmp) ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:) - deallocate(Ca,Cb,ERI_tmp) + deallocate(ERI_tmp) ! Compute linear response @@ -249,30 +269,18 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do SigC = 0.5d0*(SigC + transpose(SigC)) -! call MOtoAO_transform(nBas2,S,C,SigC) - -! Build individual Fock matrices - - Faa(:,:) = Hc(:,:) + Jaa(:,:) + Jbb(:,:) + Kaa(:,:) - Fab(:,:) = + Kab(:,:) - Fba(:,:) = + Kba(:,:) - Fbb(:,:) = Hc(:,:) + Jbb(:,:) + Jaa(:,:) + Kbb(:,:) - -! Build super Fock matrix - - F( 1:nBas , 1:nBas ) = Faa(1:nBas,1:nBas) - F( 1:nBas ,nBas+1:nBas2) = Fab(1:nBas,1:nBas) - F(nBas+1:nBas2, 1:nBas ) = Fba(1:nBas,1:nBas) - F(nBas+1:nBas2,nBas+1:nBas2) = Fbb(1:nBas,1:nBas) + call MOtoAO_transform_GHF(nBas2,nBas,S,Ca,Cb,SigC,SigCp) ! ... and add self-energy - F(:,:) = F(:,:) + SigC(:,:) + F(:,:) = F(:,:) + SigCp(:,:) ! Compute commutator and convergence criteria err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + if(nSCF > 1) Conv = maxval(abs(err)) + ! DIIS extrapolation if(max_diis > 1) then @@ -280,6 +288,10 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do call DIIS_extrapolation(rcond,nBas2Sq,nBas2Sq,n_diis,err_diis,F_diis,err,F) end if +! Transform Fock matrix in orthogonal basis + + Fp(:,:) = matmul(transpose(X),matmul(F,X)) + ! Diagonalize Fock matrix to get eigenvectors and eigenvalues Cp(:,:) = Fp(:,:) @@ -289,7 +301,7 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do C(:,:) = matmul(X,Cp) - SigC = matmul(transpose(c),matmul(SigC,c)) + call AOtoMO_transform_GHF(nBas,nBas2,Ca,Cb,SigCp,SigC) ! Form super density matrix @@ -302,11 +314,6 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) - ! Save quasiparticles energy for next cycle - - Conv = maxval(abs(err)) - eOld(:) = eGW(:) - !------------------------------------------------------------------------ ! Compute total energy !------------------------------------------------------------------------ @@ -350,7 +357,7 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do ! Print results call dipole_moment(nBas2,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsGW(nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole) + call print_qsGGW(nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole) enddo !------------------------------------------------------------------------ @@ -377,19 +384,10 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do call GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas2,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(2) - - end if - write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1) - write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW correlation energy =',EcBSE,' au' + write(*,'(2X,A50,F20.10,A4)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index 4f6e2ef..4a44a45 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -85,7 +85,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop double precision,allocatable :: c(:,:) double precision,allocatable :: cp(:,:) double precision,allocatable :: eGW(:) - double precision,allocatable :: eOld(:) double precision,allocatable :: P(:,:) double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) @@ -93,7 +92,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop double precision,allocatable :: K(:,:) double precision,allocatable :: SigC(:,:) double precision,allocatable :: SigCp(:,:) - double precision,allocatable :: SigCm(:,:) double precision,allocatable :: Z(:) double precision,allocatable :: error(:,:) @@ -130,9 +128,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop ! Memory allocation - allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & - Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & + allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & + J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), & + Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -143,7 +141,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop Conv = 1d0 P(:,:) = PHF(:,:) eGW(:) = eHF(:) - eOld(:) = eHF(:) c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 @@ -169,9 +166,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop ! AO to MO transformation of two-electron integrals - dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:) do ixyz=1,ncart - call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) end do call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO) @@ -194,10 +190,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop ! Make correlation self-energy Hermitian and transform it back to AO basis - SigCp = 0.5d0*(SigC + transpose(SigC)) - SigCm = 0.5d0*(SigC - transpose(SigC)) + SigC = 0.5d0*(SigC + transpose(SigC)) - call MOtoAO_transform(nBas,S,c,SigCp) + call MOtoAO_transform(nBas,S,c,SigC,SigCp) ! Solve the quasi-particle equation @@ -222,7 +217,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop cp(:,:) = Fp(:,:) call diagonalize_matrix(nBas,cp,eGW) c = matmul(X,cp) - SigCp = matmul(transpose(c),matmul(SigCp,c)) + call AOtoMO_transform(nBas,c,SigCp,SigC) ! Compute new density matrix in the AO basis @@ -230,8 +225,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop ! Save quasiparticles energy for next cycle - Conv = maxval(abs(error)) - eOld(:) = eGW(:) + if(nSCF > 1) Conv = maxval(abs(error)) !------------------------------------------------------------------------ ! Compute total energy @@ -283,7 +277,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,Om,XpY,XmY,rho,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,error,error_diis,F_diis) ! Perform BSE calculation diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index 340f818..328451e 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -69,7 +69,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, double precision :: ET(nspin) double precision :: EV(nspin) double precision :: EJ(nsp) - double precision :: Ex(nspin) + double precision :: EK(nspin) double precision :: EcRPA double precision :: EcGM(nspin) double precision :: EqsGW @@ -78,7 +78,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, double precision :: Conv double precision :: rcond(nspin) double precision,external :: trace_matrix - double precision,allocatable :: error_diis(:,:,:) + double precision,allocatable :: err_diis(:,:,:) double precision,allocatable :: F_diis(:,:,:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) @@ -87,7 +87,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, double precision,allocatable :: c(:,:,:) double precision,allocatable :: cp(:,:,:) double precision,allocatable :: eGW(:,:) - double precision,allocatable :: eOld(:,:) double precision,allocatable :: P(:,:,:) double precision,allocatable :: F(:,:,:) double precision,allocatable :: Fp(:,:,:) @@ -95,9 +94,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, double precision,allocatable :: K(:,:,:) double precision,allocatable :: SigC(:,:,:) double precision,allocatable :: SigCp(:,:,:) - double precision,allocatable :: SigCm(:,:,:) double precision,allocatable :: Z(:,:) - double precision,allocatable :: error(:,:,:) + double precision,allocatable :: err(:,:,:) ! Hello world @@ -136,25 +134,24 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(eGW(nBas,nspin),eOld(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), & - Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), & - SigCm(nBas,nBas,nspin),Z(nBas,nspin),Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc), & - rho(nBas,nBas,nS_sc,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), & + allocate(eGW(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), & + Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), & + Z(nBas,nspin),Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc), & + rho(nBas,nBas,nS_sc,nspin),err(nBas,nBas,nspin),err_diis(nBasSq,max_diis,nspin), & F_diis(nBasSq,max_diis,nspin)) ! Initialization - nSCF = -1 - n_diis = 0 - ispin = 1 - Conv = 1d0 - P(:,:,:) = PHF(:,:,:) - eGW(:,:) = eHF(:,:) - eOld(:,:) = eHF(:,:) - c(:,:,:) = cHF(:,:,:) - F_diis(:,:,:) = 0d0 - error_diis(:,:,:) = 0d0 - rcond(:) = 0d0 + nSCF = -1 + n_diis = 0 + ispin = 1 + Conv = 1d0 + P(:,:,:) = PHF(:,:,:) + eGW(:,:) = eHF(:,:) + c(:,:,:) = cHF(:,:,:) + F_diis(:,:,:) = 0d0 + err_diis(:,:,:) = 0d0 + rcond(:) = 0d0 !------------------------------------------------------------------------ ! Main loop @@ -181,13 +178,10 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, !-------------------------------------------------- ! AO to MO transformation of two-electron integrals !-------------------------------------------------- - - dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:) - dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:) do ixyz=1,ncart - call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz)) - call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block @@ -228,12 +222,11 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Make correlation self-energy Hermitian and transform it back to AO basis do is=1,nspin - SigCp(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is))) - SigCm(:,:,is) = 0.5d0*(SigC(:,:,is) - transpose(SigC(:,:,is))) + SigC(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is))) end do do is=1,nspin - call MOtoAO_transform(nBas,S,c(:,:,is),SigCp(:,:,is)) + call MOtoAO_transform(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) end do ! Solve the quasi-particle equation @@ -245,18 +238,18 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Check convergence do is=1,nspin - error(:,:,is) = matmul(F(:,:,is),matmul(P(:,:,is),S(:,:))) - matmul(matmul(S(:,:),P(:,:,is)),F(:,:,is)) + err(:,:,is) = matmul(F(:,:,is),matmul(P(:,:,is),S(:,:))) - matmul(matmul(S(:,:),P(:,:,is)),F(:,:,is)) end do - if(nSCF > 1) conv = maxval(abs(error(:,:,:))) + if(nSCF > 1) Conv = maxval(abs(err(:,:,:))) ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) if(minval(rcond(:)) > 1d-7) then do is=1,nspin - if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,error_diis(:,1:n_diis,is), & - F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is)) + if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,is), & + F_diis(:,1:n_diis,is),err(:,:,is),F(:,:,is)) end do else n_diis = 0 @@ -293,11 +286,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is))) end do - ! Save quasiparticles energy for next cycle - - Conv = maxval(abs(eGW(:,:) - eOld(:,:))) - eOld(:,:) = eGW(:,:) - !------------------------------------------------------------------------ ! Compute total energy !------------------------------------------------------------------------ @@ -324,19 +312,19 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Exchange energy do is=1,nspin - Ex(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is))) + EK(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is))) end do ! Total energy - EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(EK(:)) !------------------------------------------------------------------------ ! Print results !------------------------------------------------------------------------ call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,S,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,SigCp,Z,dipole) + call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,S,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,SigCp,Z,dipole) enddo !------------------------------------------------------------------------ @@ -359,7 +347,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Deallocate memory - deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,Om,XpY,XmY,rho,error,error_diis,F_diis) + deallocate(cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis) ! Perform BSE calculation diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90 index bef529e..00967f8 100644 --- a/src/HF/GHF.f90 +++ b/src/HF/GHF.f90 @@ -72,9 +72,9 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Hello world write(*,*) - write(*,*)'***********************************************' - write(*,*)'* Generalized Hartree-Fock calculation *' - write(*,*)'***********************************************' + write(*,*)'****************************************' + write(*,*)'* Generalized Hartree-Fock Calculation *' + write(*,*)'****************************************' write(*,*) ! Useful stuff diff --git a/src/HF/ROHF_fock_matrix.f90 b/src/HF/ROHF_fock_matrix.f90 index 4b35a96..ef2c4c7 100644 --- a/src/HF/ROHF_fock_matrix.f90 +++ b/src/HF/ROHF_fock_matrix.f90 @@ -59,8 +59,8 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,Fa,Fb,F) ! Block-by-block Fock matrix - call AOtoMO_transform(nBas,c,Fa) - call AOtoMO_transform(nBas,c,Fb) + Fa = matmul(transpose(c),matmul(Fa,c)) + Fb = matmul(transpose(c),matmul(Fb,c)) F(1:nC, 1:nC ) = aC*Fa(1:nC, 1:nC ) + bC*Fb(1:nC, 1:nC ) F(1:nC, nC+1:nC+nO ) = Fb(1:nC, nC+1:nC+nO ) diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 9899b9d..058f662 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -35,7 +35,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, integer :: nSCF integer :: nBasSq integer :: n_diis - double precision :: conv + double precision :: Conv double precision :: rcond(nspin) double precision :: ET(nspin) double precision :: EV(nspin) @@ -90,7 +90,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Initialization nSCF = 0 - conv = 1d0 + Conv = 1d0 n_diis = 0 F_diis(:,:,:) = 0d0 @@ -106,7 +106,7 @@ subroutine UHF(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 @@ -136,7 +136,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, err(:,:,ispin) = matmul(F(:,:,ispin),matmul(P(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),P(:,:,ispin)),F(:,:,ispin)) end do - if(nSCF > 1) conv = maxval(abs(err(:,:,:))) + if(nSCF > 1) Conv = maxval(abs(err(:,:,:))) ! DIIS extrapolation @@ -224,7 +224,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! 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,'|' + '|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',Conv,'|' end do write(*,*)'----------------------------------------------------------' diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index b190509..d296f39 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -61,7 +61,7 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) Ca(:,:) = C( 1:nBas ,1:nBas2) Cb(:,:) = C(nBas+1:nBas2,1:nBas2) -! Compute expectation values of S^2 +! Compute expectation values of S^2 (WRONG!) Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2 do mu=1,nBas @@ -109,8 +109,8 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) write(*,'(A32,1X,F16.6,A3)') ' GHF LUMO energy: ',e(LUMO)*HaToeV,' eV' write(*,'(A32,1X,F16.6,A3)') ' GHF HOMO-LUMO gap : ',Gap*HaToeV,' eV' write(*,'(A50)') '-----------------------------------------' - write(*,'(A32,1X,F16.6)') ' :',S2 - write(*,'(A50)') '-----------------------------------------' +! write(*,'(A32,1X,F16.6)') ' :',S2 +! write(*,'(A50)') '-----------------------------------------' write(*,'(A35)') ' Dipole moment (Debye) ' write(*,'(10X,4A10)') 'X','Y','Z','Tot.' write(*,'(10X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD @@ -124,8 +124,8 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) write(*,'(A32)') ' GHF orbital coefficients' write(*,'(A50)') '---------------------------------------' call matout(nBas2,nBas2,c) + write(*,*) end if - write(*,*) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' GHF orbital energies' write(*,'(A50)') '---------------------------------------' diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 7808164..dd7e9a9 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -76,8 +76,8 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) write(*,'(A32)') ' RHF orbital coefficients' write(*,'(A50)') '---------------------------------------' call matout(nBas,nBas,cHF) + write(*,*) end if - write(*,*) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' RHF orbital energies' write(*,'(A50)') '---------------------------------------' diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index d21e1f0..38339e2 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -119,8 +119,8 @@ subroutine print_UHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) write(*,'(A50)') 'UHF spin-down orbital coefficients ' write(*,'(A50)') '-----------------------------------------' call matout(nBas,nBas,c(:,:,2)) + write(*,*) end if - write(*,*) write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' UHF spin-up orbital energies ' write(*,'(A50)') '---------------------------------------' diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index e42afbb..57ce1ef 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -150,9 +150,8 @@ subroutine RQuAcK(doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do ! Read and transform dipole-related integrals - dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:) do ixyz=1,ncart - call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) end do ! 4-index transform diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 index bdc3d7b..4903872 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -156,13 +156,10 @@ subroutine UQuAcK(doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,do write(*,*) ! Read and transform dipole-related integrals - - dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:) - dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:) do ixyz=1,ncart - call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz)) - call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block