From 8ec7d3d170e709b89002a5a2b5e941d2179387eb Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 12 Sep 2024 20:10:02 +0200 Subject: [PATCH] new routine to fix GTpp --- src/GT/RG0T0pp.f90 | 111 +++++---------------------- src/GT/RGTpp_phACFTD.f90 | 16 ++++ src/GT/RGTpp_phBSE.f90 | 12 ++- src/GT/evRGTpp.f90 | 41 ++-------- src/GT/qsRGTpp.f90 | 101 ++++++++++-------------- src/LR/ppLR.f90 | 10 ++- src/RPA/sort_ppRPA.f90 | 5 +- src/utils/orthogonalize_matrix.f90 | 118 +++++++++++++++++++++++++++++ 8 files changed, 218 insertions(+), 196 deletions(-) create mode 100644 src/utils/orthogonalize_matrix.f90 diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90 index 474ea95..63d4dff 100644 --- a/src/GT/RG0T0pp.f90 +++ b/src/GT/RG0T0pp.f90 @@ -40,8 +40,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Local variables - logical :: print_T = .true. - double precision :: lambda + logical :: print_T = .true. + logical :: plot_self = .false. integer :: isp_T integer :: iblock @@ -81,10 +81,6 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d write(*,*)'*********************************' write(*,*) -! Initialization - - lambda = 1d0 - ! TDA for T if(TDA_T) then @@ -125,33 +121,13 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,lambda,eHF,ERI,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,lambda,eHF,ERI,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) - deallocate(Bpp,Cpp,Dpp) - !print*, 'LAPACK:' - !print*, Om2s - !print*, Om1s - !n_states = nOOs + 5 - !n_states_diag = n_states + 4 - !allocate(Om(nOOs+nVVs), R(nOOs+nVVs,n_states_diag)) - !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 + deallocate(Bpp,Cpp,Dpp) if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s) if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s) @@ -168,34 +144,13 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,lambda,eHF,ERI,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,lambda,eHF,ERI,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) deallocate(Bpp,Cpp,Dpp) - !print*, 'LAPACK:' - !print*, Om2t - !print*, Om1t - - !n_states = nOOt + 5 - !n_states_diag = n_states + 4 - !allocate(Om(nOOt+nVVt), R(nOOt+nVVt,n_states_diag)) - !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) if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-alpha)',nOOt,Om2t) @@ -224,7 +179,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d end if call RGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & - Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) + Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) !---------------------------------------------- ! Solve the quasi-particle equation @@ -249,8 +204,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d end if -! call RGTpp_plot_self_energy(nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, & -! Om1t,rho1t,Om2t,rho2t) + if(plot_self) call RGTpp_plot_self_energy(nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, & + Om1t,rho1t,Om2t,rho2t) !---------------------------------------------- ! Dump results @@ -264,9 +219,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,lambda,eGT,ERI,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,lambda,eGT,ERI,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) @@ -278,9 +233,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,lambda,eGT,ERI,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,lambda,eGT,ERI,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) @@ -295,17 +250,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d if(dophBSE) then - call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & - Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & + call RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI,dipole_int,eHF,eGT,EcBSE) - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(1) - - end if - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au' @@ -319,29 +267,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d if(doACFDT) then - write(*,*) '--------------------------------------------------------' - write(*,*) 'Adiabatic connection version of phBSE correlation energy' - write(*,*) '--------------------------------------------------------' - write(*,*) - - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, & nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & Om2t,X2t,Y2t,rho1t,rho2t,ERI,eHF,eGT,EcBSE) - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(2) - - end if - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@RHF correlation energy (singlet) = ',EcBSE(1),' au' diff --git a/src/GT/RGTpp_phACFTD.f90 b/src/GT/RGTpp_phACFTD.f90 index 2bd6c97..3f697a4 100644 --- a/src/GT/RGTpp_phACFTD.f90 +++ b/src/GT/RGTpp_phACFTD.f90 @@ -85,6 +85,22 @@ subroutine RGTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple allocate(Aph(nS,nS),Bph(nS,nS),TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS)) allocate(Ec(nAC,nspin)) +! Hello World + + write(*,*) '-------------------------------------------------------------' + write(*,*) ' Adiabatic connection version of BSE@GTpp correlation energy ' + write(*,*) '-------------------------------------------------------------' + write(*,*) + +! eXtended BSE + + if(doXBS) then + + write(*,*) '*** scaled screening version (XBS) ***' + write(*,*) + + end if + ! Antisymmetrized kernel version if(exchange_kernel) then diff --git a/src/GT/RGTpp_phBSE.f90 b/src/GT/RGTpp_phBSE.f90 index dfe1d76..b225406 100644 --- a/src/GT/RGTpp_phBSE.f90 +++ b/src/GT/RGTpp_phBSE.f90 @@ -1,5 +1,5 @@ -subroutine RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & - Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & +subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & + Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & ERI,dipole_int,eT,eGT,EcBSE) ! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel @@ -9,6 +9,7 @@ subroutine RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, ! Input variables + logical,intent(in) :: exchange_kernel logical,intent(in) :: TDA_T logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -189,4 +190,11 @@ subroutine RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) end if + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(1) + + end if + end subroutine diff --git a/src/GT/evRGTpp.f90 b/src/GT/evRGTpp.f90 index 87b6421..d6838d1 100644 --- a/src/GT/evRGTpp.f90 +++ b/src/GT/evRGTpp.f90 @@ -222,12 +222,12 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B ! DIIS extrapolation - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGT(:)-eOld(:),eGT(:)) + if(max_diis > 1) then - ! Reset DIIS if required + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGT(:)-eOld(:),eGT(:)) - if(abs(rcond) < 1d-15) n_diis = 0 + end if ! Save quasiparticles energy for next cycle @@ -246,17 +246,10 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B if(BSE) then - call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & - Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & + call RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI,dipole_int,eGT,eGT,EcBSE) - if(exchange_kernel) then - - EcRPA(1) = 0.5d0*EcRPA(1) - EcRPA(2) = 1.5d0*EcRPA(1) - - end if - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1) @@ -270,29 +263,10 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B if(doACFDT) then - write(*,*) '------------------------------------------------------' - write(*,*) 'Adiabatic connection version of BSE correlation energy' - write(*,*) '------------------------------------------------------' - write(*,*) - - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, & - nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & + nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & Om2t,X2t,Y2t,rho1t,rho2t,ERI,eGT,eGT,EcBSE) - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(2) - - end if - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1) @@ -306,7 +280,6 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B end if - ! Testing zone if(dotest) then diff --git a/src/GT/qsRGTpp.f90 b/src/GT/qsRGTpp.f90 index 648a825..6e9b9ad 100644 --- a/src/GT/qsRGTpp.f90 +++ b/src/GT/qsRGTpp.f90 @@ -1,9 +1,6 @@ - -! --- - -subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, & - dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, & - nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) +subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA, & + dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb, & + nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GT calculation @@ -34,7 +31,8 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas, nOrb + integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC,nO,nV,nR,nS double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nOrb) @@ -154,10 +152,10 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d allocate(error_diis(nBas_Sq,max_diis)) allocate(F_diis(nBas_Sq,max_diis)) - allocate(Om1s(nVVs), X1s(nVVs,nVVs), Y1s(nOOs,nVVs), rho1s(nOrb,nOrb,nVVs)) - allocate(Om2s(nOOs), X2s(nVVs,nOOs), Y2s(nOOs,nOOs), rho2s(nOrb,nOrb,nOOs)) - allocate(Om1t(nVVt), X1t(nVVt,nVVt), Y1t(nOOt,nVVt), rho1t(nOrb,nOrb,nVVt)) - allocate(Om2t(nOOt), X2t(nVVt,nOOt), Y2t(nOOt,nOOt), rho2t(nOrb,nOrb,nOOt)) + allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),rho1s(nOrb,nOrb,nVVs)) + allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs),rho2s(nOrb,nOrb,nOOs)) + allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),rho1t(nOrb,nOrb,nVVt)) + allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt),rho2t(nOrb,nOrb,nOOt)) ! Initialization @@ -193,7 +191,7 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! AO to MO transformation of two-electron integrals - call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO) ! Compute linear response @@ -202,9 +200,9 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) @@ -215,9 +213,9 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) @@ -246,14 +244,14 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO(nBas, nOrb, S, c, Sig, Sigp) + call MOtoAO(nBas,nOrb,S,c,Sig,Sigp) ! Solve the quasi-particle equation F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Sigp(:,:) if(nBas .ne. nOrb) then - call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1)) - call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1)) + call AOtoMO(nBas,nOrb,c(1,1),F(1,1),Fp(1,1)) + call MOtoAO(nBas,nOrb,S(1,1),c(1,1),Fp(1,1),F(1,1)) endif ! Compute commutator and convergence criteria @@ -272,20 +270,20 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Diagonalize Hamiltonian in AO basis if(nBas .eq. nOrb) then - Fp = matmul(transpose(X), matmul(F, X)) + Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nOrb, cp, eGT) - c = matmul(X, cp) + call diagonalize_matrix(nOrb,cp,eGT) + c = matmul(X,cp) else - Fp = matmul(transpose(c), matmul(F, c)) + Fp = matmul(transpose(c),matmul(F,c)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nOrb, cp, eGT) - c = matmul(c, cp) + call diagonalize_matrix(nOrb,cp,eGT) + c = matmul(c,cp) endif ! Compute new density matrix in the AO basis - P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) ! Save quasiparticles energy for next cycle @@ -319,9 +317,8 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Print results call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGTpp(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, & - eGT, c, Sig, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, & - EqsGT, dipole) + call print_qsRGTpp(nBas,nOrb,nO,nSCF,Conv,thresh,eHF,eGT,c,Sig,Z, & + ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) end do !------------------------------------------------------------------------ @@ -338,34 +335,27 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis) - deallocate(Om1s, X1s, Y1s, rho1s) - deallocate(Om2s, X2s, Y2s, rho2s) - deallocate(Om1t, X1t, Y1t, rho1t) - deallocate(Om2t, X2t, Y2t, rho2t) + deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis) + deallocate(Om1s,X1s,Y1s,rho1s) + deallocate(Om2s,X2s,Y2s,rho2s) + deallocate(Om1t,X1t,Y1t,rho1t) + deallocate(Om2t,X2t,Y2t,rho2t) stop end if ! Deallocate memory - deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis) + deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis) ! Perform BSE calculation if(dophBSE) then - call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & - Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & + call RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI_MO,dipole_int_MO,eGT,eGT,EcBSE) - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(2) - - end if - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy (singlet) =',EcBSE(1) @@ -379,20 +369,8 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d if(doACFDT) then - write(*,*) '------------------------------------------------------' - write(*,*) 'Adiabatic connection version of BSE correlation energy' - write(*,*) '------------------------------------------------------' - write(*,*) - - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, & - nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & + nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcBSE) write(*,*) @@ -408,7 +386,6 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d end if - ! Testing zone if(dotest) then @@ -419,9 +396,9 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d end if - deallocate(Om1s, X1s, Y1s, rho1s) - deallocate(Om2s, X2s, Y2s, rho2s) - deallocate(Om1t, X1t, Y1t, rho1t) - deallocate(Om2t, X2t, Y2t, rho2t) + deallocate(Om1s,X1s,Y1s,rho1s) + deallocate(Om2s,X2s,Y2s,rho2s) + deallocate(Om1t,X1t,Y1t,rho1t) + deallocate(Om2t,X2t,Y2t,rho2t) end subroutine diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index ff2469a..63f5a97 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -43,7 +43,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)) @@ -70,10 +69,13 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! if((nOO .eq. 0) .or. (nVV .eq. 0)) then - ! Diagonalize the p-p matrix - if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z) + ! Diagonalize the pp matrix + + if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z) + ! Split the various quantities in p-p and h-h parts - call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2) + + call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! else diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index 85bed67..5c6714a 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -17,7 +17,6 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) integer :: pq,ab,ij ! integer :: deg1,ab_start,ab_end ! integer :: deg2,ij_start,ij_end - integer :: OO,VV double precision,allocatable :: M(:,:) double precision,allocatable :: Z1(:,:) double precision,allocatable :: Z2(:,:) @@ -211,8 +210,8 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! S1 = + matmul(transpose(Z1),matmul(M,Z1)) ! S2 = - matmul(transpose(Z2),matmul(M,Z2)) - if(nVV > 0) call orthogonalization_matrix(nVV,VV,S1,O1) - if(nOO > 0) call orthogonalization_matrix(nOO,OO,S2,O2) + if(nVV > 0) call orthogonalize_matrix(1,nVV,S1,O1) + if(nOO > 0) call orthogonalize_matrix(1,nOO,S2,O2) if(nVV > 0) call dgemm ('N', 'N', nOO+nVV,nVV,nVV, 1d0, Z1, nOO+nVV, O1, nVV,0d0, tmp1, nOO+nVV) Z1 = tmp1 diff --git a/src/utils/orthogonalize_matrix.f90 b/src/utils/orthogonalize_matrix.f90 new file mode 100644 index 0000000..c02b316 --- /dev/null +++ b/src/utils/orthogonalize_matrix.f90 @@ -0,0 +1,118 @@ +subroutine orthogonalize_matrix(ortho_type,nBas,S,X) + +! Compute the orthogonalization matrix X + + implicit none + +! Input variables + + integer,intent(in) :: nBas,ortho_type + double precision,intent(in) :: S(nBas,nBas) + +! Local variables + + logical :: debug + double precision,allocatable :: UVec(:,:),Uval(:) + double precision,parameter :: thresh = 1d-6 + + integer :: i + +! Output variables + + double precision,intent(out) :: X(nBas,nBas) + + debug = .false. + +! Type of orthogonalization ortho_type +! +! 1 = Lowdin +! 2 = Canonical +! 3 = SVD +! + + allocate(Uvec(nBas,nBas),Uval(nBas)) + + if(ortho_type == 1) then + +! write(*,*) +! write(*,*) ' Lowdin orthogonalization' +! write(*,*) + + Uvec = S + call diagonalize_matrix(nBas,Uvec,Uval) + + do i=1,nBas + + if(Uval(i) < thresh) then + + write(*,*) 'Eigenvalue',i,' is very small in Lowdin orthogonalization = ',Uval(i) + + endif + + Uval(i) = 1d0/sqrt(Uval(i)) + + enddo + + call ADAt(nBas,Uvec,Uval,X) + + elseif(ortho_type == 2) then + +! write(*,*) +! write(*,*) 'Canonical orthogonalization' +! write(*,*) + + Uvec = S + call diagonalize_matrix(nBas,Uvec,Uval) + + do i=1,nBas + + if(Uval(i) > thresh) then + + Uval(i) = 1d0/sqrt(Uval(i)) + + else + + write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' + + endif + + enddo + + call AD(nBas,Uvec,Uval) + X = Uvec + + elseif(ortho_type == 3) then + +! write(*,*) +! write(*,*) ' SVD-based orthogonalization NYI' +! write(*,*) + +! Uvec = S +! call diagonalize_matrix(nBas,Uvec,Uval) + +! do i=1,nBas +! if(Uval(i) > thresh) then +! Uval(i) = 1d0/sqrt(Uval(i)) +! else +! write(*,*) 'Eigenvalue',i,'too small for canonical orthogonalization' +! endif +! enddo +! +! call AD(nBas,Uvec,Uval) +! X = Uvec + + endif + +! Print results + + if(debug) then + + write(*,'(A28)') '----------------------' + write(*,'(A28)') 'Orthogonalization matrix' + write(*,'(A28)') '----------------------' + call matout(nBas,nBas,X) + write(*,*) + + endif + +end subroutine