From 881220a97f15e1d267dc7c7a1e42177c09e1af08 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 1 Apr 2019 15:53:39 +0200 Subject: [PATCH] fix DIIS in GW --- input/basis | 100 ++++++++++++++------ input/methods | 2 +- input/molecule | 4 +- input/weight | 60 +++++++----- src/QuAcK/G0W0.f90 | 2 +- src/QuAcK/QuAcK.f90 | 7 +- src/QuAcK/evGW.f90 | 13 ++- src/QuAcK/excitation_density_RI.f90 | 65 +++++++++++++ src/QuAcK/qsGW.f90 | 7 +- src/QuAcK/self_energy_correlation.f90 | 2 +- src/QuAcK/self_energy_correlation_diag.f90 | 4 +- src/{QuAcK => utils}/DIIS_extrapolation.f90 | 0 12 files changed, 197 insertions(+), 69 deletions(-) create mode 100644 src/QuAcK/excitation_density_RI.f90 rename src/{QuAcK => utils}/DIIS_extrapolation.f90 (100%) diff --git a/input/basis b/input/basis index e220e53..1eb213c 100644 --- a/input/basis +++ b/input/basis @@ -1,29 +1,73 @@ -1 6 -S 8 1.00 - 2940.0000000 0.0006800 - 441.2000000 0.0052360 - 100.5000000 0.0266060 - 28.4300000 0.0999930 - 9.1690000 0.2697020 - 3.1960000 0.4514690 - 1.1590000 0.2950740 - 0.1811000 0.0125870 -S 8 1.00 - 2940.0000000 -0.0001230 - 441.2000000 -0.0009660 - 100.5000000 -0.0048310 - 28.4300000 -0.0193140 - 9.1690000 -0.0532800 - 3.1960000 -0.1207230 - 1.1590000 -0.1334350 - 0.1811000 0.5307670 +1 36 S 1 1.00 - 0.0589000 1.0000000 -P 3 1.00 - 3.6190000 0.0291110 - 0.7110000 0.1693650 - 0.1951000 0.5134580 -P 1 1.00 - 0.0601800 1.0000000 -D 1 1.00 - 0.2380000 1.0000000 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 diff --git a/input/methods b/input/methods index 40fd2a7..aaae8e8 100644 --- a/input/methods +++ b/input/methods @@ -9,6 +9,6 @@ # GF2 GF3 F F # G0W0 evGW qsGW - T T F + T T T # MCMP2 F diff --git a/input/molecule b/input/molecule index 6a6f6d1..9bc4ca6 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 2 2 0 0 + 1 1 1 0 0 # Znuc x y z - Be 0.0 0.0 0.0 + X 0.0 0.0 0.0 diff --git a/input/weight b/input/weight index e220e53..35be48a 100644 --- a/input/weight +++ b/input/weight @@ -1,29 +1,39 @@ -1 6 -S 8 1.00 - 2940.0000000 0.0006800 - 441.2000000 0.0052360 - 100.5000000 0.0266060 - 28.4300000 0.0999930 - 9.1690000 0.2697020 - 3.1960000 0.4514690 - 1.1590000 0.2950740 - 0.1811000 0.0125870 -S 8 1.00 - 2940.0000000 -0.0001230 - 441.2000000 -0.0009660 - 100.5000000 -0.0048310 - 28.4300000 -0.0193140 - 9.1690000 -0.0532800 - 3.1960000 -0.1207230 - 1.1590000 -0.1334350 - 0.1811000 0.5307670 +1 10 +S 9 1.00 + 6863.0000000 0.0002360 + 1030.0000000 0.0018260 + 234.7000000 0.0094520 + 66.5600000 0.0379570 + 21.6900000 0.1199650 + 7.7340000 0.2821620 + 2.9160000 0.4274040 + 1.1300000 0.2662780 + 0.1101000 -0.0072750 +S 9 1.00 + 6863.0000000 -0.0000430 + 1030.0000000 -0.0003330 + 234.7000000 -0.0017360 + 66.5600000 -0.0070120 + 21.6900000 -0.0231260 + 7.7340000 -0.0581380 + 2.9160000 -0.1145560 + 1.1300000 -0.1359080 + 0.1101000 0.5774410 S 1 1.00 - 0.0589000 1.0000000 + 0.2577000 1.0000000 +S 1 1.00 + 0.0440900 1.0000000 P 3 1.00 - 3.6190000 0.0291110 - 0.7110000 0.1693650 - 0.1951000 0.5134580 + 7.4360000 0.0107360 + 1.5770000 0.0628540 + 0.4352000 0.2481800 P 1 1.00 - 0.0601800 1.0000000 + 0.1438000 1.0000000 +P 1 1.00 + 0.0499400 1.0000000 D 1 1.00 - 0.2380000 1.0000000 + 0.3480000 1.0000000 +D 1 1.00 + 0.1803000 1.0000000 +F 1 1.00 + 0.3250000 1.0000000 diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index aa122b9..101432f 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -1,5 +1,5 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,P,ERI_AO_basis,ERI_MO_basis,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,P,cHF,eHF,eG0W0) ! Perform G0W0 calculation diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 52c2503..c6b179a 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -405,7 +405,7 @@ program QuAcK call cpu_time(start_G0W0) call G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,PHF,ERI_AO_basis,ERI_MO_basis,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) t_G0W0 = end_G0W0 - start_G0W0 @@ -421,9 +421,8 @@ program QuAcK if(doevGW) then call cpu_time(start_evGW) - call evGW(maxSCF_GW,thresh_GW,n_diis_GW, & - COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eHF) + call evGW(maxSCF_GW,thresh_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize, & + nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0) call cpu_time(end_evGW) t_evGW = end_evGW - start_evGW diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 12bfb51..c7e8ac5 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -20,7 +20,8 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani logical :: dRPA,linear_mixing integer :: ispin,nSCF,n_diis - double precision :: Conv,EcRPA,lambda + double precision :: rcond + double precision :: Conv,EcRPA,EcGM,lambda double precision,allocatable :: error_diis(:,:),e_diis(:,:) double precision,allocatable :: eGW(:),eOld(:),Z(:) double precision,allocatable :: H(:,:),SigmaC(:) @@ -98,13 +99,13 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(G0W) then call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),SigmaC) + Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigmaC) call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) else call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eGW, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),SigmaC) + Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigmaC) call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) endif @@ -139,7 +140,11 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani else n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW) + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW) + +! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 endif diff --git a/src/QuAcK/excitation_density_RI.f90 b/src/QuAcK/excitation_density_RI.f90 new file mode 100644 index 0000000..967f0bb --- /dev/null +++ b/src/QuAcK/excitation_density_RI.f90 @@ -0,0 +1,65 @@ +subroutine excitation_density_RI(nBas,nC,nO,nR,nS,c,G,XpY,rho) + +! Compute excitation densities in the RI approximation + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nR,nS + double precision,intent(in) :: c(nBas,nBas),G(nBas,nBas,nBas,nBas),XpY(nS,nS) + +! Local variables + + double precision,allocatable :: scr(:,:,:) + integer :: mu,nu,la,si,ia,jb,x,y,j,b + +! Output variables + + double precision,intent(out) :: rho(nBas,nBas,nS) + +! Memory allocation + allocate(scr(nBas,nBas,nS)) + + rho(:,:,:) = 0d0 + do nu=1,nBas + do si=1,nBas + do ia=1,nS + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + rho(nu,si,ia) = rho(nu,si,ia) + c(nu,j)*XpY(ia,jb)*c(si,b) + enddo + enddo + enddo + enddo + enddo + + scr(:,:,:) = 0d0 + do mu=1,nBas + do la=1,nBas + do ia=1,nS + do nu=1,nBas + do si=1,nBas + scr(mu,la,ia) = scr(mu,la,ia) + G(mu,nu,la,si)*rho(nu,si,ia) + enddo + enddo + enddo + enddo + enddo + + rho(:,:,:) = 0d0 + do ia=1,nS + do x=nC+1,nBas-nR + do y=nC+1,nBas-nR + do mu=1,nBas + do la=1,nBas + rho(x,y,ia) = rho(x,y,ia) + c(mu,x)*scr(mu,la,ia)*c(la,y) + enddo + enddo + enddo + enddo + enddo + +end subroutine excitation_density_RI diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index acb1544..87da073 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -22,6 +22,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani logical :: dRPA integer :: nSCF,nBasSq,ispin,n_diis double precision :: EcRPA,EcGM,Conv + double precision :: rcond double precision,external :: trace_matrix double precision,allocatable :: error_diis(:,:),F_diis(:,:) double precision,allocatable :: Omega(:,:),XpY(:,:,:),rho(:,:,:,:),rhox(:,:,:,:) @@ -136,7 +137,11 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + +! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 ! Diagonalize Hamiltonian in AO basis diff --git a/src/QuAcK/self_energy_correlation.f90 b/src/QuAcK/self_energy_correlation.f90 index d1d1ad2..ccc3524 100644 --- a/src/QuAcK/self_energy_correlation.f90 +++ b/src/QuAcK/self_energy_correlation.f90 @@ -66,7 +66,7 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, enddo enddo - EcGM=0d0 + EcGM = 0d0 do i=nC+1,nO EcGM = EcGM + SigC(i,i) enddo diff --git a/src/QuAcK/self_energy_correlation_diag.f90 b/src/QuAcK/self_energy_correlation_diag.f90 index 8427acc..8b4679b 100644 --- a/src/QuAcK/self_energy_correlation_diag.f90 +++ b/src/QuAcK/self_energy_correlation_diag.f90 @@ -66,7 +66,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega ! GM correlation energy - EcGM=0d0 + EcGM = 0d0 do i=nC+1,nO EcGM = EcGM + SigC(i) enddo @@ -109,7 +109,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega ! GM correlation energy - EcGM=0d0 + EcGM = 0d0 do i=nC+1,nO do a=nO+1,nBas-nR jb = 0 diff --git a/src/QuAcK/DIIS_extrapolation.f90 b/src/utils/DIIS_extrapolation.f90 similarity index 100% rename from src/QuAcK/DIIS_extrapolation.f90 rename to src/utils/DIIS_extrapolation.f90