diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index f4d0258..108cd8d 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -53,7 +53,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons double precision :: EcAC(nspin) double precision,allocatable :: SigC(:,:) double precision,allocatable :: Z(:,:) - integer :: nS_aa,nS_bb,nS_sc + integer :: nSa,nSb,nSt + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) @@ -91,12 +94,12 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons ! Memory allocation - nS_aa = nS(1) - nS_bb = nS(2) - nS_sc = nS_aa + nS_bb + nSa = nS(1) + nSb = nS(2) + nSt = nSa + nSb allocate(SigC(nBas,nspin),Z(nBas,nspin),eGWlin(nBas,nspin),eGW(nBas,nspin), & - Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc),rho(nBas,nBas,nS_sc,nspin)) + Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt),rho(nBas,nBas,nSt,nspin)) !-------------------! ! Compute screening ! @@ -106,16 +109,18 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons ispin = 1 - call phULR(ispin,dRPA,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + + call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) - if(print_W) call print_excitation_energies('phRPA@UHF',5,nS_sc,Om) + if(print_W) call print_excitation_energies('phRPA@UHF',5,nSt,Om) !----------------------! ! Excitation densities ! !----------------------! - call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) + call UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) !------------------------------------------------! ! Compute self-energy and renormalization factor ! @@ -123,11 +128,11 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons if(regularize) then do is=1,nspin - call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is),Om,rho(:,:,:,is)) + call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eHF(:,is),Om,rho(:,:,:,is)) end do end if - call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Om,rho,SigC,Z,EcGM) + call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eHF,Om,rho,SigC,Z,EcGM) !-----------------------------------! ! Solve the quasi-particle equation ! @@ -155,7 +160,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons if(is==1) write(*,*)' Spin-up orbitals ' if(is==2) write(*,*)' Spin-down orbitals ' - call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is), & + call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nSt,eHF(:,is), & Om,rho(:,:,:,is),eGWlin(:,is),eHF(:,is),eGW(:,is),Z(:,is)) end do @@ -163,8 +168,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons ! Compute RPA correlation energy - call phULR(ispin,dRPA,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + + call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) ! Dump results diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index 6e8b14d..7f0d408 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -46,6 +46,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Local variables + logical :: dRPA logical :: linear_mixing integer :: is integer :: ispin @@ -63,8 +64,11 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, double precision,allocatable :: eGW(:,:) double precision,allocatable :: eOld(:,:) double precision,allocatable :: Z(:,:) - integer :: nS_aa,nS_bb,nS_sc + integer :: nSa,nSb,nSt double precision,allocatable :: SigC(:,:) + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) @@ -92,6 +96,11 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, write(*,*) end if +! Initialization + + EcRPA = 0d0 + dRPA = .true. + ! Linear mixing linear_mixing = .false. @@ -99,13 +108,13 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Memory allocation - nS_aa = nS(1) - nS_bb = nS(2) - nS_sc = nS_aa + nS_bb + nSa = nS(1) + nSb = nS(2) + nSt = nSa + nSb allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin), & - Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc),rho(nBas,nBas,nS_sc,nspin), & - error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin)) + Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt), & + rho(nBas,nBas,nSt,nspin),error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin)) ! Initialization @@ -126,16 +135,18 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, do while(Conv > thresh .and. nSCF <= maxSCF) - ! Compute screening + ! Compute screening - call phULR(ispin,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + + call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) !----------------------! ! Excitation densities ! !----------------------! - call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) + call UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) !------------------------------------------------! ! Compute self-energy and renormalization factor ! @@ -143,11 +154,11 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, if(regularize) then do is=1,nspin - call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eGW(:,is),Om,rho(:,:,:,is)) + call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eGW(:,is),Om,rho(:,:,:,is)) end do end if - call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,Om,rho,SigC,Z,EcGM) + call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM) !-----------------------------------! ! Solve the quasi-particle equation ! @@ -171,7 +182,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, if(is==1) write(*,*)' Spin-up orbitals ' if(is==2) write(*,*)' Spin-down orbitals ' - call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is), & + call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nSt,eHF(:,is), & Om,rho(:,:,:,is),eOld(:,is),eOld(:,is),eGW(:,is),Z(:,is)) end do diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index ee91ec2..8f3ca4b 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -63,7 +63,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, integer :: ixyz integer :: is integer :: n_diis - integer :: nS_aa,nS_bb,nS_sc + integer :: nSa,nSb,nSt double precision :: dipole(ncart) double precision :: ET(nspin) @@ -80,6 +80,9 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, double precision,external :: trace_matrix double precision,allocatable :: err_diis(:,:,:) double precision,allocatable :: F_diis(:,:,:) + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) @@ -131,14 +134,14 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Memory allocation - nS_aa = nS(1) - nS_bb = nS(2) - nS_sc = nS_aa + nS_bb + nSa = nS(1) + nSb = nS(2) + nSt = nSa + nSb - allocate(eGW(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), & - Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), & - Z(nBas,nspin),Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc),rho(nBas,nBas,nS_sc,nspin), & - err(nBas,nBas,nspin),err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) + allocate(Aph(nSt,nSt),Bph(nSt,nSt),eGW(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin), & + F(nBas,nBas,nspin),Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin), & + SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),Z(nBas,nspin),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt), & + rho(nBas,nBas,nSt,nspin),err(nBas,nBas,nspin),err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) ! Initialization @@ -198,14 +201,16 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Compute linear response - call phULR(ispin,dRPA,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + + call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) !----------------------! ! Excitation densities ! !----------------------! - call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) + call UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) !------------------------------------------------! ! Compute self-energy and renormalization factor ! @@ -213,11 +218,11 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, if(regularize) then do is=1,nspin - call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eGW(:,is),Om,rho(:,:,:,is)) + call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eGW(:,is),Om,rho(:,:,:,is)) end do end if - call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,Om,rho,SigC,Z,EcGM) + call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM) ! Make correlation self-energy Hermitian and transform it back to AO basis diff --git a/src/LR/phULR.f90 b/src/LR/phULR.f90 index 9a0018a..dc5a76a 100644 --- a/src/LR/phULR.f90 +++ b/src/LR/phULR.f90 @@ -1,5 +1,4 @@ -subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,Om,XpY,XmY) +subroutine phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) ! Compute linear response for unrestricted formalism @@ -8,33 +7,16 @@ subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, ! Input variables - integer,intent(in) :: ispin - logical,intent(in) :: dRPA logical,intent(in) :: TDA - logical,intent(in) :: BSE - integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt - integer,intent(in) :: nS_sc - double precision,intent(in) :: lambda - double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) - - double precision,intent(in) :: OmRPA(nS_sc) - double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin) + double precision,intent(in) :: Aph(nSt,nSt) + double precision,intent(in) :: Bph(nSt,nSt) ! Local variables double precision,external :: trace_matrix - double precision,allocatable :: Aph(:,:) - double precision,allocatable :: Bph(:,:) double precision,allocatable :: ApB(:,:) double precision,allocatable :: AmB(:,:) double precision,allocatable :: AmBSq(:,:) @@ -50,21 +32,18 @@ subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, ! Memory allocation - allocate(Aph(nSt,nSt),Bph(nSt,nSt),ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt)) + allocate(ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt)) ! Build A and B matrices - call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) - - if(BSE) & - call UGW_phBSE_static_kernel_A(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Aph) +! if(BSE) & +! call UGW_phBSE_static_kernel_A(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & +! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,Aph) ! Tamm-Dancoff approximation if(TDA) then - Bph(:,:) = 0d0 XpY(:,:) = Aph(:,:) call diagonalize_matrix(nSt,XpY,Om) XpY(:,:) = transpose(XpY(:,:)) @@ -72,11 +51,9 @@ subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, else - call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) - - if(BSE) & - call UGW_phBSE_static_kernel_B(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Bph) +! if(BSE) & +! call UGW_phBSE_static_kernel_B(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & +! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,Bph) ! Build A + B and A - B matrices @@ -90,10 +67,6 @@ subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, if(minval(Om) < 0d0) & call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') - ! do ia=1,nSt - ! if(Om(ia) < 0d0) Om(ia) = 0d0 - ! end do - call ADAt(nSt,AmB,1d0*sqrt(Om),AmBSq) call ADAt(nSt,AmB,1d0/sqrt(Om),AmBIv) @@ -104,10 +77,6 @@ subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, if(minval(Om) < 0d0) & call print_warning('You may have instabilities in linear response: negative excitations!!') - ! do ia=1,nSt - ! if(Om(ia) < 0d0) Om(ia) = 0d0 - ! end do - Om = sqrt(Om) XpY = matmul(transpose(Z),AmBSq) diff --git a/src/RPA/phRPA.f90 b/src/RPA/phRPA.f90 index 7eeb51f..14cb769 100644 --- a/src/RPA/phRPA.f90 +++ b/src/RPA/phRPA.f90 @@ -40,9 +40,9 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS ! Hello world write(*,*) - write(*,*)'***********************************************' - write(*,*)'| Random-phase approximation calculation |' - write(*,*)'***********************************************' + write(*,*)'*********************************' + write(*,*)'* Restricted ph-RPA Calculation *' + write(*,*)'*********************************' write(*,*) ! TDA diff --git a/src/RPA/phUACFDT.f90 b/src/RPA/phUACFDT.f90 index 2174029..8ff3f60 100644 --- a/src/RPA/phUACFDT.f90 +++ b/src/RPA/phUACFDT.f90 @@ -39,10 +39,12 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin double precision,allocatable :: Ec(:,:) double precision :: EcRPA - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:,:) + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:,:) integer :: nS_aa,nS_bb,nS_sc double precision,allocatable :: Om_sc(:) @@ -90,11 +92,16 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) nS_sf = nS_ab + nS_ba - allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) + allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc),Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc),rho(nBas,nBas,nS_sc,nspin)) - call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call phULR_A(ispin,.true.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,.true.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + + call phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,Om,XpY,XmY) + + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) + + deallocate(Aph,Bph) ! Spin-conserved manifold @@ -102,7 +109,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin ispin = 1 - allocate(Om_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc),Om_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) write(*,*) '------------------------' write(*,*) 'Spin-conserved manifold ' @@ -119,14 +126,14 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin if(doXBS) then - call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) +! call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & +! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) end if - call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sc,XpY_sc,XmY_sc) +! call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & +! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcAC(ispin),Om_sc,XpY_sc,XmY_sc) call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, & ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,Ec(iAC,ispin)) @@ -144,7 +151,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin write(*,*) '-----------------------------------------------------------------------------------' write(*,*) - deallocate(Om_sc,XpY_sc,XmY_sc) + deallocate(Aph,Bph,Om_sc,XpY_sc,XmY_sc) end if @@ -156,7 +163,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin ! Memory allocation - allocate(Om_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + allocate(Aph(nS_sf,nS_sf),Bph(nS_sf,nS_sf),Om_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) write(*,*) '--------------------' write(*,*) ' Spin-flip manifold ' @@ -173,14 +180,14 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin if(doXBS) then - call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) +! call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & +! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) end if - call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sf,XpY_sf,XmY_sf) +! call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & +! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcAC(ispin),Om_sf,XpY_sf,XmY_sf) call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, & ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,Ec(iAC,ispin)) @@ -198,7 +205,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin write(*,*) '-----------------------------------------------------------------------------------' write(*,*) - deallocate(Om_sf,XpY_sf,XmY_sf) + deallocate(Aph,Bph,Om_sf,XpY_sf,XmY_sf) end if diff --git a/src/RPA/phURPA.f90 b/src/RPA/phURPA.f90 index a5588af..356f12f 100644 --- a/src/RPA/phURPA.f90 +++ b/src/RPA/phURPA.f90 @@ -33,27 +33,26 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n ! Local variables + logical :: dRPA integer :: ispin - integer :: nS_aa,nS_bb,nS_sc - double precision,allocatable :: Om_sc(:) - double precision,allocatable :: XpY_sc(:,:) - double precision,allocatable :: XmY_sc(:,:) + integer :: nSa,nSb,nSt + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) - integer :: nS_ab,nS_ba,nS_sf - double precision,allocatable :: Om_sf(:) - double precision,allocatable :: XpY_sf(:,:) - double precision,allocatable :: XmY_sf(:,:) - double precision :: rho_sc,rho_sf + double precision :: rho double precision :: EcRPA(nspin) ! Hello world write(*,*) - write(*,*)'**************************************************************' - write(*,*)'| Unrestricted direct random phase approximation calculation |' - write(*,*)'**************************************************************' + write(*,*)'***********************************' + write(*,*)'* Unrestricted ph-RPA Calculation *' + write(*,*)'***********************************' write(*,*) ! TDA @@ -65,8 +64,8 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n ! Initialization + dRPA = .true. EcRPA(:) = 0d0 - EcRPA(:) = 0d0 ! Spin-conserved transitions @@ -76,19 +75,20 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n ! Memory allocation - nS_aa = nS(1) - nS_bb = nS(2) - nS_sc = nS_aa + nS_bb + nSa = nS(1) + nSb = nS(2) + nSt = nSa + nSb - allocate(Om_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + allocate(Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt)) - call phULR(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Om_sc,rho_sc,EcRPA(ispin),Om_sc,XpY_sc,XmY_sc) - call print_excitation_energies('phRPA@UHF',5,nS_sc,Om_sc) - call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & - c,S,Om_sc,XpY_sc,XmY_sc) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) - deallocate(Om_sc,XpY_sc,XmY_sc) + call phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) + call print_excitation_energies('phRPA@UHF',5,nSt,Om) + call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY) + + deallocate(Aph,Bph,Om,XpY,XmY) endif @@ -100,19 +100,20 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n ! Memory allocation - nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) - nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) - nS_sf = nS_ab + nS_ba + nSa = (nO(1) - nC(1))*(nV(2) - nR(2)) + nSb = (nO(2) - nC(2))*(nV(1) - nR(1)) + nSt = nSa + nSb - allocate(Om_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + allocate(Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt)) - call phULR(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Om_sf,rho_sf,EcRPA(ispin),Om_sf,XpY_sf,XmY_sf) - call print_excitation_energies('phRPA@UHF',6,nS_sf,Om_sf) - call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & - c,S,Om_sf,XpY_sf,XmY_sf) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) - deallocate(Om_sf,XpY_sf,XmY_sf) + call phULR(TDA,nSa,nSa,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) + call print_excitation_energies('phRPA@UHF',6,nSt,Om) + call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY) + + deallocate(Aph,Bph,Om,XpY,XmY) endif diff --git a/src/RPA/phURPAx.f90 b/src/RPA/phURPAx.f90 index aa2cb9b..0e4c568 100644 --- a/src/RPA/phURPAx.f90 +++ b/src/RPA/phURPAx.f90 @@ -33,27 +33,25 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC, ! Local variables + logical :: dRPA integer :: ispin - integer :: nS_aa,nS_bb,nS_sc - double precision,allocatable :: Omega_sc(:) - double precision,allocatable :: XpY_sc(:,:) - double precision,allocatable :: XmY_sc(:,:) + integer :: nSa,nSb,nSt + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) - integer :: nS_ab,nS_ba,nS_sf - double precision,allocatable :: Omega_sf(:) - double precision,allocatable :: XpY_sf(:,:) - double precision,allocatable :: XmY_sf(:,:) - - double precision :: rho_sc,rho_sf + double precision :: rho double precision :: EcRPA(nspin) ! Hello world write(*,*) - write(*,*)'*********************************************************************' - write(*,*)'| Unrestricted random phase approximation calculation with exchange |' - write(*,*)'*********************************************************************' + write(*,*)'************************************' + write(*,*)'* Unrestricted ph-RPAx Calculation *' + write(*,*)'************************************' write(*,*) ! TDA @@ -66,30 +64,31 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC, ! Initialization + dRPA = .false. EcRPA(:) = 0d0 - EcRPA(:) = 0d0 ! Spin-conserved transitions - if(spin_conserved) then + if(spin_conserved) then ispin = 1 ! Memory allocation - nS_aa = nS(1) - nS_bb = nS(2) - nS_sc = nS_aa + nS_bb + nSa = nS(1) + nSb = nS(2) + nSt = nSa + nSb - allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + allocate(Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt)) - call phULR(ispin,.false.,TDA,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) - call print_excitation_energies('phRPAx@UHF',5,nS_sc,Omega_sc) - call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & - c,S,Omega_sc,XpY_sc,XmY_sc) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) - deallocate(Omega_sc,XpY_sc,XmY_sc) + call phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) + call print_excitation_energies('phRPA@UHF',5,nSt,Om) + call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY) + + deallocate(Aph,Bph,Om,XpY,XmY) endif @@ -101,22 +100,30 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC, ! Memory allocation - nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) - nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) - nS_sf = nS_ab + nS_ba + nSa = (nO(1) - nC(1))*(nV(2) - nR(2)) + nSb = (nO(2) - nC(2))*(nV(1) - nR(1)) + nSt = nSa + nSb - allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + allocate(Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt)) - call phULR(ispin,.false.,TDA,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) - call print_excitation_energies('phRPAx@UHF',6,nS_sf,Omega_sf) - call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & - c,S,Omega_sf,XpY_sf,XmY_sf) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) - deallocate(Omega_sf,XpY_sf,XmY_sf) + call phULR(TDA,nSa,nSa,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) + call print_excitation_energies('phRPA@UHF',6,nSt,Om) + call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY) + + deallocate(Aph,Bph,Om,XpY,XmY) endif + if(exchange_kernel) then + + EcRPA(1) = 0.5d0*EcRPA(1) + EcRPA(2) = 1.5d0*EcRPA(2) + + end if + if(exchange_kernel) then EcRPA(1) = 0.5d0*EcRPA(1)