From 6da5afae6d3b2fda443426d2bf9ee2fad58509d1 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 17 Jun 2020 14:32:37 +0200 Subject: [PATCH] BSE2 iterative --- input/basis | 69 +++++++++++++++++-- input/molecule | 5 +- input/molecule.xyz | 5 +- src/QuAcK/BSE2.f90 | 29 +++++++- src/QuAcK/BSE2_dynamic_perturbation.f90 | 39 ++++++----- src/QuAcK/Bethe_Salpeter.f90 | 12 ++-- .../Bethe_Salpeter_dynamic_perturbation.f90 | 3 +- ...alpeter_dynamic_perturbation_iterative.f90 | 3 +- 8 files changed, 125 insertions(+), 40 deletions(-) diff --git a/input/basis b/input/basis index 0a75001..80356bb 100644 --- a/input/basis +++ b/input/basis @@ -1,7 +1,64 @@ -1 2 -S 3 - 1 38.4216340 0.0237660 - 2 5.7780300 0.1546790 - 3 1.2417740 0.4696300 +1 6 +S 9 +1 9.046000E+03 7.000000E-04 +2 1.357000E+03 5.389000E-03 +3 3.093000E+02 2.740600E-02 +4 8.773000E+01 1.032070E-01 +5 2.856000E+01 2.787230E-01 +6 1.021000E+01 4.485400E-01 +7 3.838000E+00 2.782380E-01 +8 7.466000E-01 1.544000E-02 +9 2.248000E-01 -2.864000E-03 +S 9 +1 9.046000E+03 -1.530000E-04 +2 1.357000E+03 -1.208000E-03 +3 3.093000E+02 -5.992000E-03 +4 8.773000E+01 -2.454400E-02 +5 2.856000E+01 -6.745900E-02 +6 1.021000E+01 -1.580780E-01 +7 3.838000E+00 -1.218310E-01 +8 7.466000E-01 5.490030E-01 +9 2.248000E-01 5.788150E-01 S 1 - 1 0.2979640 1.0000000 +1 2.248000E-01 1.000000E+00 +P 4 +1 1.355000E+01 3.991900E-02 +2 2.917000E+00 2.171690E-01 +3 7.973000E-01 5.103190E-01 +4 2.185000E-01 4.622140E-01 +P 1 +1 2.185000E-01 1.000000E+00 +D 1 +1 8.170000E-01 1.0000000 +2 6 +S 9 +1 9.046000E+03 7.000000E-04 +2 1.357000E+03 5.389000E-03 +3 3.093000E+02 2.740600E-02 +4 8.773000E+01 1.032070E-01 +5 2.856000E+01 2.787230E-01 +6 1.021000E+01 4.485400E-01 +7 3.838000E+00 2.782380E-01 +8 7.466000E-01 1.544000E-02 +9 2.248000E-01 -2.864000E-03 +S 9 +1 9.046000E+03 -1.530000E-04 +2 1.357000E+03 -1.208000E-03 +3 3.093000E+02 -5.992000E-03 +4 8.773000E+01 -2.454400E-02 +5 2.856000E+01 -6.745900E-02 +6 1.021000E+01 -1.580780E-01 +7 3.838000E+00 -1.218310E-01 +8 7.466000E-01 5.490030E-01 +9 2.248000E-01 5.788150E-01 +S 1 +1 2.248000E-01 1.000000E+00 +P 4 +1 1.355000E+01 3.991900E-02 +2 2.917000E+00 2.171690E-01 +3 7.973000E-01 5.103190E-01 +4 2.185000E-01 4.622140E-01 +P 1 +1 2.185000E-01 1.000000E+00 +D 1 +1 8.170000E-01 1.0000000 diff --git a/input/molecule b/input/molecule index c78e87e..76ebcdf 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,5 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 2 7 7 0 0 # Znuc x y z - He 0.0 0.0 0.0 + N 0. 0. -1.04008632 + N 0. 0. +1.04008632 diff --git a/input/molecule.xyz b/input/molecule.xyz index 797b5fc..e1773f0 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,4 @@ - 1 + 2 - He 0.0000000000 0.0000000000 0.0000000000 + N 0.0000000000 0.0000000000 -0.5503900175 + N 0.0000000000 0.0000000000 0.5503900175 diff --git a/src/QuAcK/BSE2.f90 b/src/QuAcK/BSE2.f90 index c698e19..cff9aef 100644 --- a/src/QuAcK/BSE2.f90 +++ b/src/QuAcK/BSE2.f90 @@ -58,10 +58,22 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) ! Compute dynamic correction for BSE via perturbation theory - if(dBSE) & + + if(dBSE) then + + if(evDyn) then + + call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & + ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + else + call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + end if + + end if + end if !------------------- @@ -81,9 +93,20 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & ! Compute dynamic correction for BSE via perturbation theory - if(dBSE) & - call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & + if(dBSE) then + + if(evDyn) then + + call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & + ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + else + + call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + end if + + end if end if diff --git a/src/QuAcK/BSE2_dynamic_perturbation.f90 b/src/QuAcK/BSE2_dynamic_perturbation.f90 index 0efb9bc..703684f 100644 --- a/src/QuAcK/BSE2_dynamic_perturbation.f90 +++ b/src/QuAcK/BSE2_dynamic_perturbation.f90 @@ -35,17 +35,19 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF, double precision,allocatable :: X(:) double precision,allocatable :: Y(:) - double precision,allocatable :: A_dyn(:,:) - double precision,allocatable :: ZA_dyn(:,:) + double precision,allocatable :: Ap_dyn(:,:) + double precision,allocatable :: Am_dyn(:,:) + double precision,allocatable :: ZAp_dyn(:,:) + double precision,allocatable :: ZAm_dyn(:,:) double precision,allocatable :: B_dyn(:,:) double precision,allocatable :: ZB_dyn(:,:) ! Memory allocation - allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),A_dyn(nS,nS),ZA_dyn(nS,nS)) + allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) - if(.not.dTDA) allocate(B_dyn(nS,nS),ZB_dyn(nS,nS)) + if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),B_dyn(nS,nS),ZB_dyn(nS,nS)) ! Print main components of transition vectors @@ -68,28 +70,31 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF, ! Resonant part of the BSE correction for dynamical TDA - call BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI(:,:,:,:),eGF(:),OmBSE(ia),A_dyn(:,:),ZA_dyn(:,:)) + call BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),Ap_dyn,ZAp_dyn) if(dTDA) then - ZDyn(ia) = dot_product(X(:),matmul(ZA_dyn(:,:),X(:))) - OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) + ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) else - ! Anti-resonant part of the BSE correction (frequency independent) + ! Second part of the resonant and anti-resonant part of the BSE correction (frequency independent) - call BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI(:,:,:,:),eGF(:),B_dyn(:,:),ZB_dyn(:,:)) + call BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),Am_dyn,ZAm_dyn) + ZAm_dyn(:,:) = - ZAm_dyn(:,:) - ZDyn(ia) = dot_product(X(:),matmul(ZA_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(ZA_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(ZB_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(ZB_dyn(:,:),X(:))) + call BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn,ZB_dyn) - OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(B_dyn(:,:),X(:))) + ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & + - dot_product(Y,matmul(ZAm_dyn,Y)) & + + dot_product(X,matmul(ZB_dyn,Y)) & + - dot_product(Y,matmul(ZB_dyn,X)) + + OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & + - dot_product(Y,matmul(Am_dyn,Y)) & + + dot_product(X,matmul(B_dyn,Y)) & + - dot_product(Y,matmul(B_dyn,X)) end if diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index 9e4b26a..88ad8fe 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -83,15 +83,15 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) else if(W_BSE) then - call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) end if @@ -137,15 +137,15 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) else if(W_BSE) then - call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) end if diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index 371f3b2..71028b9 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) +subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) ! Compute dynamical effects via perturbation theory for BSE @@ -7,7 +7,6 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS, ! Input variables - logical,intent(in) :: TDA logical,intent(in) :: dTDA double precision,intent(in) :: eta integer,intent(in) :: nBas diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 index 0f058f4..9435c09 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) +subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) ! Compute self-consistently the dynamical effects via perturbation theory for BSE @@ -7,7 +7,6 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO ! Input variables - logical,intent(in) :: TDA logical,intent(in) :: dTDA double precision,intent(in) :: eta integer,intent(in) :: nBas