From 84e01bfcb675c5f03be8d26f8e91f2e336d49780 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 11 Jun 2020 15:36:16 +0200 Subject: [PATCH] -w --- input/basis | 76 ++----------------- input/molecule | 5 +- input/molecule.xyz | 5 +- input/options | 2 +- src/QuAcK/AOtoMO_integral_transform.f90 | 1 + .../Bethe_Salpeter_AB_matrix_dynamic.f90 | 69 ++++++++++++----- .../Bethe_Salpeter_ZAB_matrix_dynamic.f90 | 64 +++++++++++----- .../Bethe_Salpeter_dynamic_perturbation.f90 | 46 ++++++----- 8 files changed, 131 insertions(+), 137 deletions(-) diff --git a/input/basis b/input/basis index bbe0bfe..0a75001 100644 --- a/input/basis +++ b/input/basis @@ -1,71 +1,7 @@ -1 9 -S 8 - 1 9046.0000000 0.0007000 - 2 1357.0000000 0.0053890 - 3 309.3000000 0.0274060 - 4 87.7300000 0.1032070 - 5 28.5600000 0.2787230 - 6 10.2100000 0.4485400 - 7 3.8380000 0.2782380 - 8 0.7466000 0.0154400 -S 8 - 1 9046.0000000 -0.0001530 - 2 1357.0000000 -0.0012080 - 3 309.3000000 -0.0059920 - 4 87.7300000 -0.0245440 - 5 28.5600000 -0.0674590 - 6 10.2100000 -0.1580780 - 7 3.8380000 -0.1218310 - 8 0.7466000 0.5490030 +1 2 +S 3 + 1 38.4216340 0.0237660 + 2 5.7780300 0.1546790 + 3 1.2417740 0.4696300 S 1 - 1 0.2248000 1.0000000 -S 1 - 1 0.0612400 1.0000000 -P 3 - 1 13.5500000 0.0399190 - 2 2.9170000 0.2171690 - 3 0.7973000 0.5103190 -P 1 - 1 0.2185000 1.0000000 -P 1 - 1 0.0561100 1.0000000 -D 1 - 1 0.8170000 1.0000000 -D 1 - 1 0.2300000 1.0000000 -2 9 -S 8 - 1 9046.0000000 0.0007000 - 2 1357.0000000 0.0053890 - 3 309.3000000 0.0274060 - 4 87.7300000 0.1032070 - 5 28.5600000 0.2787230 - 6 10.2100000 0.4485400 - 7 3.8380000 0.2782380 - 8 0.7466000 0.0154400 -S 8 - 1 9046.0000000 -0.0001530 - 2 1357.0000000 -0.0012080 - 3 309.3000000 -0.0059920 - 4 87.7300000 -0.0245440 - 5 28.5600000 -0.0674590 - 6 10.2100000 -0.1580780 - 7 3.8380000 -0.1218310 - 8 0.7466000 0.5490030 -S 1 - 1 0.2248000 1.0000000 -S 1 - 1 0.0612400 1.0000000 -P 3 - 1 13.5500000 0.0399190 - 2 2.9170000 0.2171690 - 3 0.7973000 0.5103190 -P 1 - 1 0.2185000 1.0000000 -P 1 - 1 0.0561100 1.0000000 -D 1 - 1 0.8170000 1.0000000 -D 1 - 1 0.2300000 1.0000000 - + 1 0.2979640 1.0000000 diff --git a/input/molecule b/input/molecule index 76ebcdf..c78e87e 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 7 7 0 0 + 1 1 1 0 0 # Znuc x y z - N 0. 0. -1.04008632 - N 0. 0. +1.04008632 + He 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index e1773f0..797b5fc 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - N 0.0000000000 0.0000000000 -0.5503900175 - N 0.0000000000 0.0000000000 0.5503900175 + He 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/options b/input/options index dd6fae5..ae799d8 100644 --- a/input/options +++ b/input/options @@ -9,7 +9,7 @@ # GF: maxSCF thresh DIIS n_diis lin renorm BSE TDA eta 256 0.00001 T 5 T 3 T T 0.00367493 # GW/GT: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA_W TDA G0W GW0 lin eta - 256 0.00001 T 5 F F T F F F F T 0.00367493 + 256 0.00001 T 5 F F T T T F F T 0.00367493 # ACFDT: AC Kx XBS F F T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/src/QuAcK/AOtoMO_integral_transform.f90 b/src/QuAcK/AOtoMO_integral_transform.f90 index 6ab35bb..8246fe7 100644 --- a/src/QuAcK/AOtoMO_integral_transform.f90 +++ b/src/QuAcK/AOtoMO_integral_transform.f90 @@ -73,6 +73,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis) do nu=1,nBas ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j)*scr(i,nu,k,l) enddo + print*,i,k,j,l,ERI_MO_basis(i,j,k,l) enddo enddo enddo diff --git a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 index 7e0fd20..25e7dea 100644 --- a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 @@ -1,4 +1,5 @@ -subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,A_dyn,B_dyn) +subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, & + Ap,Am,Bp,Bm) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -18,18 +19,26 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O ! Local variables integer :: maxS - double precision :: chi_A,chi_B,eps,eps_A,eps_B + double precision :: chi_A,chi_B + double precision :: chi_Ap,chi_Bp,eps_Ap,eps_Bp + double precision :: chi_Am,chi_Bm,eps_Am,eps_Bm integer :: i,j,a,b,ia,jb,kc ! Output variables - double precision,intent(out) :: A_dyn(nS,nS) - double precision,intent(out) :: B_dyn(nS,nS) + double precision,intent(out) :: Ap(nS,nS) + double precision,intent(out) :: Am(nS,nS) + + double precision,intent(out) :: Bp(nS,nS) + double precision,intent(out) :: Bm(nS,nS) ! Initialization - A_dyn(:,:) = 0d0 - B_dyn(:,:) = 0d0 + Ap(:,:) = 0d0 + Am(:,:) = 0d0 + + Bp(:,:) = 0d0 + Bm(:,:) = 0d0 ! Number of poles taken into account @@ -56,31 +65,51 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O enddo - A_dyn(ia,jb) = A_dyn(ia,jb) - 4d0*lambda*chi_A - B_dyn(ia,jb) = B_dyn(ia,jb) - 4d0*lambda*chi_B + Ap(ia,jb) = Ap(ia,jb) - 4d0*lambda*chi_A + Am(ia,jb) = Am(ia,jb) - 4d0*lambda*chi_A - chi_A = 0d0 - chi_B = 0d0 + Bp(ia,jb) = Bp(ia,jb) - 4d0*lambda*chi_B + Bm(ia,jb) = Bm(ia,jb) - 4d0*lambda*chi_B + + chi_Ap = 0d0 + chi_Am = 0d0 + + chi_Bp = 0d0 + chi_Bm = 0d0 do kc=1,maxS - eps_A = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*eps_A/(eps_A**2 + eta**2) + eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) + chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) - eps_A = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*eps_A/(eps_A**2 + eta**2) + eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) + chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) - eps_B = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*eps_B/(eps_B**2 + eta**2) + eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) - eps_B = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*eps_B/(eps_B**2 + eta**2) + eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) + + eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) + chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) + + eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) + chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) + + eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) + + eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) enddo - A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi_A + Ap(ia,jb) = Ap(ia,jb) - 2d0*lambda*chi_Ap + Am(ia,jb) = Am(ia,jb) - 2d0*lambda*chi_Am - B_dyn(ia,jb) = B_dyn(ia,jb) - 2d0*lambda*chi_B + Bp(ia,jb) = Bp(ia,jb) - 2d0*lambda*chi_Bp + Bm(ia,jb) = Bm(ia,jb) - 2d0*lambda*chi_Bm enddo enddo diff --git a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 index 26fa9db..9ed4416 100644 --- a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 @@ -1,6 +1,7 @@ -subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,ZA,ZB) +subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, & + ZAp,ZAm,ZBp,ZBm) -! Compute the dynamic part of the Bethe-Salpeter equation matrices +! Compute the dynamic part of the renormalization for the Bethe-Salpeter equation matrices implicit none include 'parameters.h' @@ -18,19 +19,25 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, ! Local variables integer :: maxS - double precision :: chi_A,chi_B - double precision :: eps_A,eps_B + double precision :: chi_Ap,chi_Bp,eps_Ap,eps_Bp + double precision :: chi_Am,chi_Bm,eps_Am,eps_Bm integer :: i,j,a,b,ia,jb,kc ! Output variables - double precision,intent(out) :: ZA(nS,nS) - double precision,intent(out) :: ZB(nS,nS) + double precision,intent(out) :: ZAp(nS,nS) + double precision,intent(out) :: ZAm(nS,nS) + + double precision,intent(out) :: ZBp(nS,nS) + double precision,intent(out) :: ZBm(nS,nS) ! Initialization - ZA(:,:) = 0d0 - ZB(:,:) = 0d0 + ZAp(:,:) = 0d0 + ZAm(:,:) = 0d0 + + ZBp(:,:) = 0d0 + ZBm(:,:) = 0d0 ! Number of poles taken into account @@ -47,28 +54,45 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, do b=nO+1,nBas-nR jb = jb + 1 - chi_A = 0d0 - chi_B = 0d0 + chi_Ap = 0d0 + chi_Am = 0d0 + + chi_Bp = 0d0 + chi_Bm = 0d0 do kc=1,maxS - eps_A = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*(eps_A**2 - eta**2)/(eps_A**2 + eta**2)**2 + eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) + chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 - eps_A = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*(eps_A**2 - eta**2)/(eps_A**2 + eta**2)**2 + eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) + chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 - eps_B = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*(eps_B**2 - eta**2)/(eps_B**2 + eta**2)**2 + eps_Am = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 - eps_B = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*(eps_B**2 - eta**2)/(eps_B**2 + eta**2)**2 + eps_Am = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 + + eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) + chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 + + eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) + chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 + + eps_Bm = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 + + eps_Bm = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 enddo - ZA(ia,jb) = ZA(ia,jb) + 2d0*lambda*chi_A + ZAp(ia,jb) = ZAp(ia,jb) + 2d0*lambda*chi_Ap + ZAm(ia,jb) = ZAm(ia,jb) + 2d0*lambda*chi_Am - ZB(ia,jb) = ZB(ia,jb) + 2d0*lambda*chi_B + ZBp(ia,jb) = ZBp(ia,jb) + 2d0*lambda*chi_Bp + ZBm(ia,jb) = ZBm(ia,jb) + 2d0*lambda*chi_Bm enddo enddo diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index 07d1fc2..9113519 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 @@ -36,17 +36,23 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O 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 :: ZAp_dyn(:,:) - double precision,allocatable :: B_dyn(:,:) - double precision,allocatable :: ZB_dyn(:,:) + double precision,allocatable :: Bp_dyn(:,:) + double precision,allocatable :: ZBp_dyn(:,:) + + double precision,allocatable :: Am_dyn(:,:) + double precision,allocatable :: ZAm_dyn(:,:) + + double precision,allocatable :: Bm_dyn(:,:) + double precision,allocatable :: ZBm_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),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS)) ! Print main components of transition vectors @@ -74,37 +80,37 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O ! Resonant part of the BSE correction for dynamical TDA call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & - A_dyn(:,:)) + Ap_dyn(:,:)) ! Renormalization factor of the resonant parts for dynamical TDA call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & - ZA_dyn(:,:)) + ZAp_dyn(:,:)) - 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 ! Resonant and anti-resonant part of the BSE correction call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & - A_dyn(:,:),B_dyn(:,:)) + Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) ! Renormalization factor of the resonant and anti-resonant parts call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & - ZA_dyn(:,:),ZB_dyn(:,:)) + ZAp_dyn(:,:),ZAm_dyn(:,:),ZBp_dyn(:,:),ZBm_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(:))) + ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) & + - dot_product(Y(:),matmul(ZAm_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(ZBp_dyn(:,:),Y(:))) & + - dot_product(Y(:),matmul(ZBm_dyn(:,:),X(:))) - 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(:))) + OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & + - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & + - dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) end if