From 77a54456e9d7246ad4f27848328fdc24a36b144f Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 9 Nov 2023 17:38:30 +0100 Subject: [PATCH] fix dipole bug in qsGW --- src/GW/qsRGW.f90 | 10 +++------- src/GW/qsUGW.f90 | 23 +++++++++-------------- src/QuAcK/GQuAcK.f90 | 4 +--- src/QuAcK/QuAcK.f90 | 2 +- 4 files changed, 14 insertions(+), 25 deletions(-) diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index f1d1bc7..99bf0c0 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -167,7 +167,7 @@ subroutine qsRGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do ! AO to MO transformation of two-electron integrals do ixyz=1,ncart - call AOtoMO_transform(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + call AOtoMO_transform(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) end do call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO) @@ -202,13 +202,13 @@ subroutine qsRGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + if(nSCF > 1) Conv = maxval(abs(error)) + ! DIIS extrapolation if(max_diis > 1) then - n_diis = min(n_diis+1,max_diis) call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) - end if ! Diagonalize Hamiltonian in AO basis @@ -223,10 +223,6 @@ subroutine qsRGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - ! Save quasiparticles energy for next cycle - - if(nSCF > 1) Conv = maxval(abs(error)) - !------------------------------------------------------------------------ ! Compute total energy !------------------------------------------------------------------------ diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index 2c4b901..4c62060 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -1,7 +1,6 @@ -subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, & - dBSE,dTDA,spin_conserved,spin_flip,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & - nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa, & - dipole_int_bb,PHF,cHF,eHF) +subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & + eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & + dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -181,8 +180,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, !-------------------------------------------------- do ixyz=1,ncart - 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)) + call AOtoMO_transform(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO_transform(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block @@ -247,14 +246,10 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! 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,err_diis(:,1:n_diis,is), & - F_diis(:,1:n_diis,is),err(:,:,is),F(:,:,is)) - end do - else - n_diis = 0 - end if + do is=1,nspin + 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 ! Transform Fock matrix in orthogonal basis diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 6fb893b..773b41b 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -2,7 +2,7 @@ subroutine GQuAcK(doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA, doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & - spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & + TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) @@ -41,8 +41,6 @@ subroutine GQuAcK(doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA, logical,intent(in) :: reg_MP - logical,intent(in) :: spin_conserved - logical,intent(in) :: spin_flip logical,intent(in) :: TDA integer,intent(in) :: maxSCF_GF,max_diis_GF diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 647a99a..81b179d 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -219,7 +219,7 @@ program QuAcK doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & - spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & + TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)