diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90 index cf6f5f6..cace39e 100644 --- a/src/GT/RG0T0pp.f90 +++ b/src/GT/RG0T0pp.f90 @@ -64,6 +64,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d double precision,allocatable :: Z(:) double precision,allocatable :: eGT(:) double precision,allocatable :: eGTlin(:) + integer, allocatable :: supp_data_int(:) + double precision, allocatable :: supp_data_dbl(:) double precision, allocatable :: Om(:), R(:,:) @@ -94,11 +96,11 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Dimensions of the pp-RPA linear reponse matrices - nOOs = nO*(nO + 1)/2 - nVVs = nV*(nV + 1)/2 + !nOOs = nO*(nO + 1)/2 + !nVVs = nV*(nV + 1)/2 - !nOOs = nO*nO - !nVVs = nV*nV + nOOs = nO*nO + nVVs = nV*nV nOOt = nO*(nO - 1)/2 nVVt = nV*(nV - 1)/2 @@ -118,8 +120,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d !---------------------------------------------- ispin = 1 - iblock = 1 - !iblock = 3 + !iblock = 1 + iblock = 3 ! Compute linear response @@ -138,7 +140,19 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d !n_states = nOOs + 5 !n_states_diag = n_states + 4 !allocate(Om(nOOs+nVVs), R(nOOs+nVVs,n_states_diag)) - !call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOs, nVVs, 1.d0, eHF, 0.d0, ERI, Om, R, n_states, n_states_diag, "RPA") + !allocate(supp_data_dbl(1), supp_data_int(1)) + !supp_data_int(1) = 0 + !supp_data_dbl(1) = 0.d0 + !call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOs, nVVs, & + ! 1.d0, & ! lambda + ! eHF(1), & + ! 0.d0, & ! eF + ! ERI(1,1,1,1), & + ! supp_data_int(1), 1, & + ! supp_data_dbl(1), 1, & + ! Om(1), R(1,1), n_states, n_states_diag, "RPA") + !deallocate(supp_data_dbl, supp_data_int) + !deallocate(Om, R) !stop if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s) @@ -169,7 +183,19 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d !n_states = nOOt + 5 !n_states_diag = n_states + 4 !allocate(Om(nOOt+nVVt), R(nOOt+nVVt,n_states_diag)) - !call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOt, nVVt, 1.d0, eHF, 0.d0, ERI, Om, R, n_states, n_states_diag, "RPA") + !allocate(supp_data_dbl(1), supp_data_int(1)) + !supp_data_int(1) = 0 + !supp_data_dbl(1) = 0.d0 + !call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOt, nVVt, & + ! 1.d0, & ! lambda + ! eHF(1), & + ! 0.d0, & ! eF + ! ERI(1,1,1,1), & + ! supp_data_int(1), 1, & + ! supp_data_dbl(1), 1, & + ! Om(1), R(1,1), n_states, n_states_diag, "RPA") + !deallocate(Om, R) + !deallocate(supp_data_dbl) !stop if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t) diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90 index cfe44ed..d652005 100644 --- a/src/GW/RGW_ppBSE.f90 +++ b/src/GW/RGW_ppBSE.f90 @@ -1,4 +1,4 @@ -subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) +subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies at the pp level @@ -15,16 +15,16 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS logical,intent(in) :: triplet double precision,intent(in) :: eta - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS - double precision,intent(in) :: eW(nBas) - double precision,intent(in) :: eGW(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eW(nOrb) + double precision,intent(in) :: eGW(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables @@ -37,6 +37,10 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS integer :: nOO integer :: nVV + integer :: p, q, m + integer :: i_data, supp_data_dbl_size, supp_data_int_size + integer :: n_states, n_states_diag + double precision,allocatable :: Aph(:,:) double precision,allocatable :: Bph(:,:) @@ -62,6 +66,10 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS double precision,allocatable :: KC_sta(:,:) double precision,allocatable :: KD_sta(:,:) + integer, allocatable :: supp_data_int(:) + double precision,allocatable :: supp_data_dbl(:) + double precision,allocatable :: Om(:), R(:,:) + ! Output variables double precision,intent(out) :: EcBSE(nspin) @@ -74,16 +82,16 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS isp_W = 1 EcRPA = 0d0 - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), & Aph(nS,nS),Bph(nS,nS)) - call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph) + call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph) - if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) deallocate(XpY_RPA,XmY_RPA,Aph,Bph) @@ -104,42 +112,100 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS nOO = nO*(nO+1)/2 nVV = nV*(nV+1)/2 - allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & - Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & - Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), & - KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV)) + allocate(Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) ! Compute BSE excitation energies - call RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) - call RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) - if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) + ! --- + ! LAPACK + ! --- + + allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) + allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) + + call RGW_ppBSE_static_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) + call RGW_ppBSE_static_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) + if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) - if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) + call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) + if(.not.TDA) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + deallocate(Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) - ! TODO - !call ppLR_RGW_ppBSE_davidson(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + print*, 'LAPACK:' + print*, Om2 + print*, Om1 - call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + ! --- + + + + ! --- + ! Davidson + ! --- + + n_states = nOO + 5 + n_states_diag = n_states + 4 + allocate(Om(nOO+nVV), R(nOO+nVV,n_states_diag)) + + supp_data_int = 1 + allocate(supp_data_int(supp_data_int_size)) + supp_data_int(1) = nS + + supp_data_dbl_size = nS + nOrb*nOrb*nS + 1 + allocate(supp_data_dbl(supp_data_dbl_size)) + ! scalars + supp_data_dbl(1) = eta + i_data = 1 + ! rho_RPA + do m = 1, nS + do q = 1, nOrb + do p = 1, nOrb + i_data = i_data + 1 + supp_data_dbl(i_data) = rho_RPA(p,q,m) + enddo + enddo + enddo + ! OmRPA + do m = 1, nS + i_data = i_data + 1 + supp_data_dbl(i_data) = OmRPA(m) + enddo + + call ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, & + 1.d0, & ! lambda + eGW(1), & + 0.d0, & ! eF + ERI(1,1,1,1), & + supp_data_int(1), supp_data_int_size, & + supp_data_dbl(1), supp_data_dbl_size, & + Om(1), R(1,1), n_states, n_states_diag, "GW") + + deallocate(Om, R) + deallocate(supp_data_dbl, supp_data_int) + stop + + ! --- + + call ppLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) !----------------------------------------------------! ! Compute the dynamical screening at the ppBSE level ! !----------------------------------------------------! if(dBSE) & - call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & + call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) - deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) + deallocate(Om1,X1,Y1,Om2,X2,Y2) end if !------------------- @@ -159,40 +225,97 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS nOO = nO*(nO-1)/2 nVV = nV*(nV-1)/2 - allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & - Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & - Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), & - KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV)) + allocate(Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) ! Compute BSE excitation energies - call RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) - call RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) - if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) + ! --- + ! LAPACK + ! --- - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) - if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) + allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) + + call RGW_ppBSE_static_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) + call RGW_ppBSE_static_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) + if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) + + call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) + call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) + if(.not.TDA) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + deallocate(Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) + + !print*, 'LAPACK:' + !print*, Om2 + !print*, Om1 + + ! --- + ! Davidson + ! --- + + !n_states = nOO + 5 + !n_states_diag = n_states + 4 + !allocate(Om(nOO+nVV), R(nOO+nVV,n_states_diag)) + + !supp_data_int = 1 + !allocate(supp_data_int(supp_data_int_size)) + !supp_data_int(1) = nS + + !supp_data_dbl_size = nS + nOrb*nOrb*nS + 1 + !allocate(supp_data_dbl(supp_data_dbl_size)) + !! scalars + !supp_data_dbl(1) = eta + !i_data = 1 + !! rho_RPA + !do m = 1, nS + ! do q = 1, nOrb + ! do p = 1, nOrb + ! i_data = i_data + 1 + ! supp_data_dbl(i_data) = rho_RPA(p,q,m) + ! enddo + ! enddo + !enddo + !! OmRPA + !do m = 1, nS + ! i_data = i_data + 1 + ! supp_data_dbl(i_data) = OmRPA(m) + !enddo + + !call ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, & + ! 1.d0, & ! lambda + ! eGW(1), & + ! 0.d0, & ! eF + ! ERI(1,1,1,1), & + ! supp_data_int(1), supp_data_int_size, & + ! supp_data_dbl(1), supp_data_dbl_size, & + ! Om(1), R(1,1), n_states, n_states_diag, "GW") + + !deallocate(Om, R) + !deallocate(supp_data_dbl, supp_data_int) + !stop + + ! --- EcBSE(ispin) = 3d0*EcBSE(ispin) - call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + call ppLR_transition_vectors(.false.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) !----------------------------------------------------! ! Compute the dynamical screening at the ppBSE level ! !----------------------------------------------------! if(dBSE) & - call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & + call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) - deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) + deallocate(Om1,X1,Y1,Om2,X2,Y2) end if diff --git a/src/GW/RGW_ppBSE_static_kernel_B.f90 b/src/GW/RGW_ppBSE_static_kernel_B.f90 index cb6cd6c..c28e479 100644 --- a/src/GW/RGW_ppBSE_static_kernel_B.f90 +++ b/src/GW/RGW_ppBSE_static_kernel_B.f90 @@ -1,4 +1,4 @@ -subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,KB) +subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,KB) ! Compute the VVOO block of the static screening W for the pp-BSE @@ -8,7 +8,7 @@ subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -18,9 +18,9 @@ subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd integer,intent(in) :: nVV double precision,intent(in) :: eta double precision,intent(in) :: lambda - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: rho(nOrb,nOrb,nS) ! Local variables @@ -44,8 +44,8 @@ subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd if(ispin == 1) then ab = 0 - do a=nO+1,nBas-nR - do b=a,nBas-nR + do a=nO+1,nOrb-nR + do b=a,nOrb-nR ab = ab + 1 ij = 0 do i=nC+1,nO @@ -75,8 +75,8 @@ subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd if(ispin == 2) then ab = 0 - do a=nO+1,nBas-nR - do b=a+1,nBas-nR + do a=nO+1,nOrb-nR + do b=a+1,nOrb-nR ab = ab + 1 ij = 0 do i=nC+1,nO diff --git a/src/GW/RGW_ppBSE_static_kernel_C.f90 b/src/GW/RGW_ppBSE_static_kernel_C.f90 index b7636b4..d799fa5 100644 --- a/src/GW/RGW_ppBSE_static_kernel_C.f90 +++ b/src/GW/RGW_ppBSE_static_kernel_C.f90 @@ -1,4 +1,4 @@ -subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI,Om,rho,KC) +subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,lambda,ERI,Om,rho,KC) ! Compute the VVVV block of the static screening W for the pp-BSE @@ -8,7 +8,7 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -17,9 +17,9 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER integer,intent(in) :: nVV double precision,intent(in) :: eta double precision,intent(in) :: lambda - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: rho(nOrb,nOrb,nS) ! Local variables @@ -44,21 +44,21 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER if(ispin == 1) then - a0 = nBas - nR - nO + a0 = nOrb - nR - nO lambda4 = 4.d0 * lambda eta2 = eta * eta - allocate(tmp_m(nBas,nBas,nS)) - allocate(tmp(nBas,nBas,nBas,nBas)) + allocate(tmp_m(nOrb,nOrb,nS)) + allocate(tmp(nOrb,nOrb,nOrb,nOrb)) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(m, c, a, eps) & - !$OMP SHARED(nS, nBas, eta2, Om, rho, tmp_m) + !$OMP SHARED(nS, nOrb, eta2, Om, rho, tmp_m) !$OMP DO do m = 1, nS eps = Om(m) / (Om(m)*Om(m) + eta2) - do c = 1, nBas - do a = 1, nBas + do c = 1, nOrb + do a = 1, nOrb tmp_m(a,c,m) = eps * rho(a,c,m) enddo enddo @@ -66,19 +66,19 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER !$OMP END DO !$OMP END PARALLEL - call dgemm("N", "T", nBas*nBas, nBas*nBas, nS, 1.d0, & - tmp_m(1,1,1), nBas*nBas, rho(1,1,1), nBas*nBas, & - 0.d0, tmp(1,1,1,1), nBas*nBas) + call dgemm("N", "T", nOrb*nOrb, nOrb*nOrb, nS, 1.d0, & + tmp_m(1,1,1), nOrb*nOrb, rho(1,1,1), nOrb*nOrb, & + 0.d0, tmp(1,1,1,1), nOrb*nOrb) deallocate(tmp_m) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(a, b, aa, ab, c, d, cd, tmp_ab) & - !$OMP SHARED(nO, nBas, nR, nS, a0, lambda4, tmp, KC) + !$OMP SHARED(nO, nOrb, nR, nS, a0, lambda4, tmp, KC) !$OMP DO - do a = nO+1, nBas-nR + do a = nO+1, nOrb-nR aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO - do b = a, nBas-nR + do b = a, nOrb-nR ab = aa + b tmp_ab = lambda4 @@ -87,8 +87,8 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER endif cd = 0 - do c = nO+1, nBas-nR - do d = c, nBas-nR + do c = nO+1, nOrb-nR + do d = c, nOrb-nR cd = cd + 1 KC(ab,cd) = -tmp_ab * (tmp(a,c,b,d) + tmp(a,d,b,c)) @@ -105,12 +105,12 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER deallocate(tmp) -! do a=nO+1,nBas-nR -! do b=a,nBas-nR +! do a=nO+1,nOrb-nR +! do b=a,nOrb-nR ! ab = ab + 1 ! cd = 0 -! do c=nO+1,nBas-nR -! do d=c,nBas-nR +! do c=nO+1,nOrb-nR +! do d=c,nOrb-nR ! cd = cd + 1 ! ! chi = 0d0 @@ -125,7 +125,7 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER ! OpenMP implementation ! --- --- --- ! -! a0 = nBas - nR - nO +! a0 = nOrb - nR - nO ! lambda4 = 4.d0 * lambda ! eta2 = eta * eta ! @@ -141,11 +141,11 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER ! ! !$OMP PARALLEL DEFAULT(NONE) & ! !$OMP PRIVATE(a, b, aa, ab, c, d, cd, m, tmp_ab) & -! !$OMP SHARED(nO, nBas, nR, nS, a0, lambda4, Om_tmp, rho, KC) +! !$OMP SHARED(nO, nOrb, nR, nS, a0, lambda4, Om_tmp, rho, KC) ! !$OMP DO -! do a = nO+1, nBas-nR +! do a = nO+1, nOrb-nR ! aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO -! do b = a, nBas-nR +! do b = a, nOrb-nR ! ab = aa + b ! ! tmp_ab = lambda4 @@ -154,8 +154,8 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER ! endif ! ! cd = 0 -! do c = nO+1, nBas-nR -! do d = c, nBas-nR +! do c = nO+1, nOrb-nR +! do d = c, nOrb-nR ! cd = cd + 1 ! ! KC(ab,cd) = 0d0 @@ -184,12 +184,12 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER ! --- --- --- ! ! ab = 0 -! do a=nO+1,nBas-nR -! do b=a,nBas-nR +! do a=nO+1,nOrb-nR +! do b=a,nOrb-nR ! ab = ab + 1 ! cd = 0 -! do c=nO+1,nBas-nR -! do d=c,nBas-nR +! do c=nO+1,nOrb-nR +! do d=c,nOrb-nR ! cd = cd + 1 ! ! chi = 0d0 @@ -215,21 +215,21 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER if(ispin == 2) then - a0 = nBas - nR - nO - 1 + a0 = nOrb - nR - nO - 1 lambda4 = 4.d0 * lambda eta2 = eta * eta - allocate(tmp_m(nBas,nBas,nS)) - allocate(tmp(nBas,nBas,nBas,nBas)) + allocate(tmp_m(nOrb,nOrb,nS)) + allocate(tmp(nOrb,nOrb,nOrb,nOrb)) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(m, c, a, eps) & - !$OMP SHARED(nS, nBas, eta2, Om, rho, tmp_m) + !$OMP SHARED(nS, nOrb, eta2, Om, rho, tmp_m) !$OMP DO do m = 1, nS eps = Om(m) / (Om(m)*Om(m) + eta2) - do c = 1, nBas - do a = 1, nBas + do c = 1, nOrb + do a = 1, nOrb tmp_m(a,c,m) = eps * rho(a,c,m) enddo enddo @@ -237,24 +237,24 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER !$OMP END DO !$OMP END PARALLEL - call dgemm("N", "T", nBas*nBas, nBas*nBas, nS, 1.d0, & - tmp_m(1,1,1), nBas*nBas, rho(1,1,1), nBas*nBas, & - 0.d0, tmp(1,1,1,1), nBas*nBas) + call dgemm("N", "T", nOrb*nOrb, nOrb*nOrb, nS, 1.d0, & + tmp_m(1,1,1), nOrb*nOrb, rho(1,1,1), nOrb*nOrb, & + 0.d0, tmp(1,1,1,1), nOrb*nOrb) deallocate(tmp_m) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(a, b, aa, ab, c, d, cd) & - !$OMP SHARED(nO, nBas, nR, nS, a0, lambda4, tmp, KC) + !$OMP SHARED(nO, nOrb, nR, nS, a0, lambda4, tmp, KC) !$OMP DO - do a = nO+1, nBas-nR + do a = nO+1, nOrb-nR aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO - 1 - do b = a+1, nBas-nR + do b = a+1, nOrb-nR ab = aa + b cd = 0 - do c = nO+1, nBas-nR - do d = c+1, nBas-nR + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR cd = cd + 1 KC(ab,cd) = lambda4 * (-tmp(a,c,b,d) + tmp(a,d,b,c)) @@ -272,7 +272,7 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER ! OpenMP implementation ! --- --- --- ! -! a0 = nBas - nR - nO - 1 +! a0 = nOrb - nR - nO - 1 ! lambda4 = 4.d0 * lambda ! eta2 = eta * eta ! @@ -288,16 +288,16 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER ! ! !$OMP PARALLEL DEFAULT(NONE) & ! !$OMP PRIVATE(a, b, aa, ab, c, d, cd, m) & -! !$OMP SHARED(nO, nBas, nR, nS, a0, lambda4, Om_tmp, rho, KC) +! !$OMP SHARED(nO, nOrb, nR, nS, a0, lambda4, Om_tmp, rho, KC) ! !$OMP DO -! do a = nO+1, nBas-nR +! do a = nO+1, nOrb-nR ! aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO - 1 -! do b = a+1, nBas-nR +! do b = a+1, nOrb-nR ! ab = aa + b ! ! cd = 0 -! do c = nO+1, nBas-nR -! do d = c+1, nBas-nR +! do c = nO+1, nOrb-nR +! do d = c+1, nOrb-nR ! cd = cd + 1 ! ! KC(ab,cd) = 0d0 @@ -323,12 +323,12 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER ! --- --- --- ! ! ab = 0 -! do a=nO+1,nBas-nR -! do b=a+1,nBas-nR +! do a=nO+1,nOrb-nR +! do b=a+1,nOrb-nR ! ab = ab + 1 ! cd = 0 -! do c=nO+1,nBas-nR -! do d=c+1,nBas-nR +! do c=nO+1,nOrb-nR +! do d=c+1,nOrb-nR ! cd = cd + 1 ! ! chi = 0d0 diff --git a/src/GW/RGW_ppBSE_static_kernel_D.f90 b/src/GW/RGW_ppBSE_static_kernel_D.f90 index 16c73a5..5ff0bc0 100644 --- a/src/GW/RGW_ppBSE_static_kernel_D.f90 +++ b/src/GW/RGW_ppBSE_static_kernel_D.f90 @@ -1,4 +1,4 @@ -subroutine RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI,Om,rho,KD) +subroutine RGW_ppBSE_static_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,lambda,ERI,Om,rho,KD) ! Compute the OOOO block of the static screening W for the pp-BSE @@ -8,7 +8,7 @@ subroutine RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ER ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -17,9 +17,9 @@ subroutine RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ER integer,intent(in) :: nOO double precision,intent(in) :: eta double precision,intent(in) :: lambda - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: rho(nOrb,nOrb,nS) ! Local variables diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index d5b4146..498e508 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -44,7 +44,6 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) double precision, external :: trace_matrix - N = nOO + nVV allocate(M(N,N), Z(N,N), Om(N)) diff --git a/src/LR/ppLR_GW_davidson.f90 b/src/LR/ppLR_GW_davidson.f90 new file mode 100644 index 0000000..a86aab9 --- /dev/null +++ b/src/LR/ppLR_GW_davidson.f90 @@ -0,0 +1,402 @@ + + +! --- + +subroutine ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e, eF, n_states_diag, & + ERI, eta, rho, Om, U, W) + + + implicit none + + integer, intent(in) :: ispin + integer, intent(in) :: n_states_diag + integer, intent(in) :: nOO, nVV, nOrb, nC, nO, nR, nS + double precision, intent(in) :: lambda, eF, eta + double precision, intent(in) :: e(nOrb) + double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision, intent(in) :: rho(nOrb,nOrb,nS), Om(nS) + double precision, intent(in) :: U(nOO+nVV,n_states_diag) + double precision, intent(out) :: W(nOO+nVV,n_states_diag) + + integer :: i, j, ij, k, l, kl + integer :: a, b, c, d, ab, cd + integer :: m + integer :: state + double precision :: mat_tmp, chi, eps + + double precision, external :: Kronecker_delta + + if(ispin .eq. 1) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = a, nOrb-nR + ab = ab + 1 + + do state = 1, n_states_diag + + W(ab,state) = 0.d0 + + cd = 0 + do c = nO+1, nOrb-nR + do d = c, nOrb-nR + cd = cd + 1 + + mat_tmp = (e(a) + e(b) - eF) * Kronecker_delta(a, c) * Kronecker_delta(b, d) & + + lambda * (ERI(a,b,c,d) + ERI(a,b,d,c)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(c, d))) + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(a,c,m) * rho(b,d,m) * Om(m) / eps & + - rho(a,d,m) * rho(b,c,m) * Om(m) / eps + enddo + mat_tmp = mat_tmp + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(c, d))) + + W(ab,state) = W(ab,state) + mat_tmp * U(cd,state) + enddo + enddo + + ij = nVV + do i = nC+1, nO + do j = i, nO + ij = ij + 1 + + mat_tmp = lambda * (ERI(a,b,i,j) + ERI(a,b,j,i)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(i, j))) + + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,a,m) * rho(j,b,m) * Om(m) / eps & + - rho(i,b,m) * rho(a,j,m) * Om(m) / eps + enddo + mat_tmp = mat_tmp + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(i, j))) + + W(ab,state) = W(ab,state) - mat_tmp * U(ij,state) + enddo + enddo + + enddo ! state + enddo ! b + enddo ! a + + ! --- + + ij = nVV + do i = nC+1, nO + do j = i, nO + ij = ij + 1 + + do state = 1, n_states_diag + + W(ij,state) = 0.d0 + + ab = 0 + do a = nO+1, nOrb-nR + do b = a, nOrb-nR + ab = ab + 1 + + mat_tmp = lambda * (ERI(a,b,i,j) + ERI(a,b,j,i)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(i, j))) + + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,a,m) * rho(j,b,m) * Om(m) / eps & + - rho(i,b,m) * rho(a,j,m) * Om(m) / eps + enddo + mat_tmp = mat_tmp + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(i, j))) + + W(ij,state) = W(ij,state) + mat_tmp * U(ab,state) + enddo + enddo + + kl = nVV + do k = nC+1, nO + do l = k, nO + kl = kl + 1 + + mat_tmp = - (e(i) + e(j) - eF) * Kronecker_delta(i, k) * Kronecker_delta(j, l) & + + lambda * (ERI(i,j,k,l) + ERI(i,j,l,k)) / dsqrt( (1.d0 + Kronecker_delta(i, j)) & + * (1.d0 + Kronecker_delta(k, l))) + + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,k,m) * rho(j,l,m) * Om(m) / eps & + - rho(i,l,m) * rho(j,k,m) * Om(m) / eps + enddo + mat_tmp = mat_tmp + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(i, j)) & + * (1.d0 + Kronecker_delta(k, l))) + + W(ij,state) = W(ij,state) - mat_tmp * U(kl,state) + enddo + enddo + + enddo ! state + enddo ! j + enddo ! i + + + elseif(ispin .eq. 2) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = a+1, nOrb-nR + ab = ab + 1 + + do state = 1, n_states_diag + + W(ab,state) = 0.d0 + + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 + + mat_tmp = (e(a) + e(b) - eF) * Kronecker_delta(a, c) * Kronecker_delta(b, d) & + + lambda * (ERI(a,b,c,d) - ERI(a,b,d,c)) + + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(a,c,m) * rho(b,d,m) * Om(m) / eps & + + rho(a,d,m) * rho(b,c,m) * Om(m) / eps + enddo + mat_tmp = mat_tmp + 4.d0 * lambda * chi + + W(ab,state) = W(ab,state) + mat_tmp * U(cd,state) + enddo + enddo + + ij = nVV + do i = nC+1, nO + do j = i+1, nO + ij = ij + 1 + + mat_tmp = lambda * (ERI(a,b,i,j) - ERI(a,b,j,i)) + + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,a,m) * rho(j,b,m) * Om(m) / eps & + + rho(i,b,m) * rho(a,j,m) * Om(m) / eps + end do + mat_tmp = mat_tmp + 4.d0 * lambda * chi + + W(ab,state) = W(ab,state) - mat_tmp * U(ij,state) + enddo + enddo + + enddo ! state + enddo ! b + enddo ! a + + ! --- + + ij = nVV + do i = nC+1, nO + do j = i+1, nO + ij = ij + 1 + + do state = 1, n_states_diag + + W(ij,state) = 0.d0 + + ab = 0 + do a = nO+1, nOrb-nR + do b = a+1, nOrb-nR + ab = ab + 1 + + mat_tmp = lambda * (ERI(a,b,i,j) - ERI(a,b,j,i)) + + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,a,m) * rho(j,b,m) * Om(m) / eps & + + rho(i,b,m) * rho(a,j,m) * Om(m) / eps + end do + mat_tmp = mat_tmp + 4.d0 * lambda * chi + + W(ij,state) = W(ij,state) + mat_tmp * U(ab,state) + enddo + enddo + + kl = nVV + do k = nC+1, nO + do l = k+1, nO + kl = kl + 1 + + mat_tmp = - (e(i) + e(j) - eF) * Kronecker_delta(i, k) * Kronecker_delta(j, l) & + + lambda * (ERI(i,j,k,l) - ERI(i,j,l,k)) + + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,k,m) * rho(j,l,m) * Om(m) / eps & + + rho(i,l,m) * rho(j,k,m) * Om(m) / eps + enddo + mat_tmp = mat_tmp + 4.d0 * lambda * chi + + W(ij,state) = W(ij,state) - mat_tmp * U(kl,state) + enddo + enddo + + enddo ! state + enddo ! j + enddo ! i + + else + + print*, ' Error in ppLR_GW_HR_calc' + print*, ' ispin is not supported' + print*, ' ispin = ', ispin + stop + + endif + + return +end + +! --- + +subroutine ppLR_GW_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e, eF, ERI, eta, rho, Om, H_diag) + + implicit none + + integer, intent(in) :: ispin + integer, intent(in) :: nOO, nVV, nOrb, nC, nO, nR, nS + double precision, intent(in) :: lambda, eF, eta + double precision, intent(in) :: e(nOrb) + double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision, intent(in) :: rho(nOrb,nOrb,nS), Om(nS) + double precision, intent(out) :: H_diag(nOO+nVV) + + integer :: i, j, ij, k, l, kl + integer :: a, b, c, d, ab, cd + integer :: m + double precision :: chi, eps + + double precision, external :: Kronecker_delta + + + if(ispin .eq. 1) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = a, nOrb-nR + ab = ab + 1 + cd = 0 + do c = nO+1, nOrb-nR + do d = c, nOrb-nR + cd = cd + 1 + if(a .ne. c) cycle + if(b .ne. d) cycle + H_diag(ab) = e(a) + e(b) - eF & + + lambda * (ERI(a,b,c,d) + ERI(a,b,d,c)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(c, d))) + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(a,c,m) * rho(b,d,m) * Om(m) / eps & + - rho(a,d,m) * rho(b,c,m) * Om(m) / eps + end do + H_diag(ab) = H_diag(ab) + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(c, d))) + enddo + enddo + enddo ! b + enddo ! a + + ij = nVV + do i = nC+1, nO + do j = i, nO + ij = ij + 1 + kl = 0 + do k = nC+1, nO + do l = k, nO + kl = kl + 1 + if(i .ne. k) cycle + if(j .ne. l) cycle + H_diag(ij) = e(i) + e(j) - eF & + - lambda * (ERI(i,j,k,l) + ERI(i,j,l,k)) / dsqrt( (1.d0 + Kronecker_delta(i, j)) & + * (1.d0 + Kronecker_delta(k, l))) + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,k,m) * rho(j,l,m) * Om(m) / eps & + - rho(i,l,m) * rho(j,k,m) * Om(m) / eps + enddo + H_diag(ij) = H_diag(ij) - 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(i, j)) & + * (1.d0 + Kronecker_delta(k, l))) + enddo + enddo + enddo ! j + enddo ! i + + elseif(ispin .eq. 2) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = a+1, nOrb-nR + ab = ab + 1 + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 + if(a .ne. c) cycle + if(b .ne. d) cycle + H_diag(ab) = e(a) + e(b) - eF + lambda * (ERI(a,b,c,d) - ERI(a,b,d,c)) + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(a,c,m) * rho(b,d,m) * Om(m) / eps & + + rho(a,d,m) * rho(b,c,m) * Om(m) / eps + enddo + H_diag(ab) = H_diag(ab) + 4.d0 * lambda * chi + enddo + enddo + enddo ! b + enddo ! a + + ij = nVV + do i = nC+1, nO + do j = i+1, nO + ij = ij + 1 + kl = 0 + do k = nC+1, nO + do l = k+1, nO + kl = kl + 1 + if(i .ne. k) cycle + if(j .ne. l) cycle + H_diag(ij) = e(i) + e(j) - eF - lambda * (ERI(i,j,k,l) - ERI(i,j,l,k)) + chi = 0.d0 + do m = 1, nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,k,m) * rho(j,l,m) * Om(m) / eps & + + rho(i,l,m) * rho(j,k,m) * Om(m) / eps + end do + H_diag(ij) = H_diag(ij) - 4.d0 * lambda * chi + enddo + enddo + enddo ! j + enddo ! i + + else + + print*, ' Error in ppLR_GW_H_diag' + print*, ' ispin is not supported' + print*, ' ispin = ', ispin + stop + + endif + + return +end + +! --- + diff --git a/src/LR/ppLR_RPA_davidson.f90 b/src/LR/ppLR_RPA_davidson.f90 index bb001c9..fffea65 100644 --- a/src/LR/ppLR_RPA_davidson.f90 +++ b/src/LR/ppLR_RPA_davidson.f90 @@ -76,15 +76,15 @@ subroutine ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_ W(ij,state) = 0.d0 - cd = 0 - do c = nO+1, nOrb-nR - do d = c, nOrb-nR - cd = cd + 1 + ab = 0 + do a = nO+1, nOrb-nR + do b = a, nOrb-nR + ab = ab + 1 - mat_tmp = lambda * (ERI(c,d,i,j) + ERI(c,d,j,i)) / dsqrt( (1.d0 + Kronecker_delta(c, d)) & + mat_tmp = lambda * (ERI(a,b,i,j) + ERI(a,b,j,i)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) & * (1.d0 + Kronecker_delta(i, j))) - W(ij,state) = W(ij,state) + mat_tmp * U(cd,state) + W(ij,state) = W(ij,state) + mat_tmp * U(ab,state) enddo enddo @@ -154,14 +154,14 @@ subroutine ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_ W(ij,state) = 0.d0 - cd = 0 - do c = nO+1, nOrb-nR - do d = c+1, nOrb-nR - cd = cd + 1 + ab = 0 + do a = nO+1, nOrb-nR + do b = a+1, nOrb-nR + ab = ab + 1 - mat_tmp = lambda * (ERI(c,d,i,j) - ERI(c,d,j,i)) + mat_tmp = lambda * (ERI(a,b,i,j) - ERI(a,b,j,i)) - W(ij,state) = W(ij,state) + mat_tmp * U(cd,state) + W(ij,state) = W(ij,state) + mat_tmp * U(ab,state) enddo enddo @@ -230,14 +230,14 @@ subroutine ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_ W(ij,state) = 0.d0 - cd = 0 - do c = nO+1, nOrb-nR - do d = nO+1, nOrb-nR - cd = cd + 1 + ab = 0 + do a = nO+1, nOrb-nR + do b = nO+1, nOrb-nR + ab = ab + 1 - mat_tmp = lambda * ERI(c,d,i,j) + mat_tmp = lambda * ERI(a,b,i,j) - W(ij,state) = W(ij,state) + mat_tmp * U(cd,state) + W(ij,state) = W(ij,state) + mat_tmp * U(ab,state) enddo enddo @@ -259,6 +259,7 @@ subroutine ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_ else + print*, ' Error in ppLR_RPA_HR_calc' print*, ' ispin is not supported' print*, ' ispin = ', ispin stop @@ -445,6 +446,7 @@ subroutine ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI else + print*, ' Error in ppLR_RPA_H_diag' print*, ' ispin is not supported' print*, ' ispin = ', ispin stop diff --git a/src/LR/ppLR_davidson.f90 b/src/LR/ppLR_davidson.f90 index b8cdb10..ac50fce 100644 --- a/src/LR/ppLR_davidson.f90 +++ b/src/LR/ppLR_davidson.f90 @@ -1,7 +1,10 @@ ! --- -subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, Om, R, n_states, n_states_diag, kernel) +subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, & + supp_data_int, supp_data_int_size, & + supp_data_dbl, supp_data_dbl_size, & + Om, R, n_states, n_states_diag, kernel) ! ! Extract the low n_states @@ -22,17 +25,22 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, integer, intent(in) :: nC, nO, nR, nOrb, nOO, nVV integer, intent(in) :: n_states ! nb of physical states integer, intent(in) :: n_states_diag ! nb of states used to get n_states + integer, intent(in) :: supp_data_int_size + integer, intent(in) :: supp_data_dbl_size character(len=*), intent(in) :: kernel double precision, intent(in) :: lambda, eF double precision, intent(in) :: e(nOrb) double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + integer, intent(in) :: supp_data_int(supp_data_int_size) + double precision, intent(in) :: supp_data_dbl(supp_data_dbl_size) double precision, intent(out) :: Om(n_states) double precision, intent(out) :: R(nOO+nVV,n_states_diag) integer :: N, M integer :: iter, itermax, itertot integer :: shift1, shift2 - integer :: i, j, ij, k, l, kl, a, b, ab, c, d, cd + integer :: i, j, k, l, ab + integer :: p, q, mm, i_data, nS integer :: i_omax(n_states) logical :: converged character*(16384) :: write_buffer @@ -40,6 +48,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, double precision :: lambda_tmp double precision :: to_print(2,n_states) double precision :: mem + double precision :: eta + character(len=len(kernel)) :: kernel_name double precision, allocatable :: H_diag(:) double precision, allocatable :: W(:,:) double precision, allocatable :: U(:,:) @@ -47,6 +57,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, double precision, allocatable :: residual_norm(:) double precision, allocatable :: overlap(:) double precision, allocatable :: S_check(:,:) + double precision, allocatable :: rho_tmp(:,:,:), Om_tmp(:) double precision, external :: u_dot_u @@ -56,6 +67,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, itermax = 8 M = n_states_diag * itermax + call lower_case(trim(kernel), kernel_name) + if(M .ge. N) then print*, 'N = ', N print*, 'M = ', M @@ -71,7 +84,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, write(*,'(A40, I12)') 'Number of states = ', n_states write(*,'(A40, I12)') 'Number of states in diagonalization = ', n_states_diag write(*,'(A40, I12)') 'Number of basis functions = ', N - + write(*,'(A40, A12)') 'Kernel: ', kernel_name @@ -82,14 +95,55 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, allocate(overlap(n_states_diag)) allocate(residual_norm(n_states_diag)) - mem = 8.d0 * dble(nOrb + nOrb**4 + N*n_states) / 1d6 - write(*,'(A40, F12.4)') 'I/O mem (MB) = ', mem + mem = 8.d0 * dble(nOrb + nOrb**4 + N*n_states) & + + 8.d0 * dble(supp_data_dbl_size) + 4.d0 * dble(supp_data_int_size) - mem = 8.d0 * dble(N + N*M + N*M + M*M + M*M + M + n_states_diag + n_states_diag) / 1d6 - write(*,'(A40, F12.4)') 'tmp mem (MB) = ', mem + write(*,'(A40, F12.4)') 'I/O mem (MB) = ', mem / (1024.d0*1024.d0) + + mem = 8.d0 * dble(N + N*M + N*M + M*M + M*M + M + n_states_diag + n_states_diag) + + write(*,'(A40, F12.4)') 'tmp mem (MB) = ', mem / (1024.d0*1024.d0) - call ppLR_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag, kernel) + if(kernel_name .eq. "rpa") then + + call ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, & + ERI(1,1,1,1), H_diag(1)) + + elseif(kernel_name .eq. "gw") then + + nS = supp_data_int(1) + + allocate(rho_tmp(nOrb,nOrb,nS)) + allocate(Om_tmp(nS)) + + eta = supp_data_dbl(1) + i_data = 1 + do mm = 1, nS + do q = 1, nOrb + do p = 1, nOrb + i_data = i_data + 1 + rho_tmp(p,q,mm) = supp_data_dbl(i_data) + enddo + enddo + enddo + do mm = 1, nS + i_data = i_data + 1 + Om_tmp(mm) = supp_data_dbl(i_data) + enddo + + call ppLR_GW_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, & + ERI(1,1,1,1), eta, rho_tmp(1,1,1), Om_tmp(1), H_diag(1)) + + !! TODO + !elseif(kernel_name .eq. "gf2") then + + else + + print*, ' kernel not supported', kernel + stop + + endif !print*, "H_diag:" !do ab = 1, N @@ -187,8 +241,22 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, !enddo !deallocate(S_check) - call ppLR_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, & - ERI(1,1,1,1), U(1,shift1+1), W(1,shift1+1), kernel) + if(kernel_name .eq. "rpa") then + + call ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, n_states_diag, & + ERI(1,1,1,1), & + U(1,shift1+1), W(1,shift1+1)) + + elseif(kernel_name .eq. "gw") then + + call ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, n_states_diag, & + ERI(1,1,1,1), eta, rho_tmp(1,1,1), Om_tmp(1), & + U(1,shift1+1), W(1,shift1+1)) + + !! TODO + !elseif(kernel_name .eq. "gf2") then + + endif else @@ -353,92 +421,9 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, deallocate(overlap) deallocate(residual_norm) - return -end - -! --- - -subroutine ppLR_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W, kernel) - - implicit none - - integer, intent(in) :: ispin - integer, intent(in) :: n_states_diag - integer, intent(in) :: nOO, nVV, nOrb, nC, nO, nR - character(len=*), intent(in) :: kernel - double precision, intent(in) :: lambda, eF - double precision, intent(in) :: e(nOrb) - double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) - double precision, intent(in) :: U(nOO+nVV,n_states_diag) - double precision, intent(out) :: W(nOO+nVV,n_states_diag) - - character(len=len(kernel)) :: kernel_name - - call lower_case(trim(kernel), kernel_name) - - if(kernel_name .eq. "rpa") then - - call ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W) - - !! TODO - !elseif(kernel_name .eq. "gw") then - - ! call ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W) - - !! TODO - !elseif(kernel_name .eq. "gf2") then - - ! call ppLR_GF2_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W) - - else - - print*, ' Error in routine ppLR_HR_calc' - print*, ' kernel not supported', kernel - stop - - endif - - return -end - -! --- - -subroutine ppLR_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag, kernel) - - implicit none - - integer, intent(in) :: ispin - integer, intent(in) :: nOO, nVV, nOrb, nC, nO, nR - character(len=*), intent(in) :: kernel - double precision, intent(in) :: lambda, eF - double precision, intent(in) :: e(nOrb) - double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) - double precision, intent(out) :: H_diag(nOO+nVV) - - character(len=len(kernel)) :: kernel_name - - call lower_case(trim(kernel), kernel_name) - - if(kernel_name .eq. "rpa") then - - call ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag) - - !! TODO - !elseif(kernel_name .eq. "gw") then - - ! call ppLR_GW_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag) - - !! TODO - !elseif(kernel_name .eq. "gf2") then - - ! call ppLR_GF2_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag) - - else - - print*, ' Error in routine ppLR_H_diag' - print*, ' kernel not supported', kernel - stop - + if(kernel_name .eq. "gw") then + deallocate(rho_tmp) + deallocate(Om_tmp) endif return diff --git a/src/utils/non_sym_diag.f90 b/src/utils/non_sym_diag.f90 index 0bb19a6..82d368d 100644 --- a/src/utils/non_sym_diag.f90 +++ b/src/utils/non_sym_diag.f90 @@ -465,7 +465,7 @@ subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0) accu_nd = accu_nd + dabs(S(j,k)) enddo enddo - if(accu_nd .gt. 1d-12) then + if(accu_nd .gt. 1d-7) then print*, ' accu_nd =', accu_nd print*, ' your strategy for degenerates orbitals failed !' print*, m, 'deg on', i