From ff58cd17c6436ee49bc4396fdb676d2217044e09 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 24 Sep 2020 14:39:37 +0200 Subject: [PATCH] UCIS --- input/basis | 55 ++++---- input/methods | 4 +- input/molecule | 4 +- input/molecule.xyz | 2 +- input/options | 8 +- src/QuAcK/UCIS.f90 | 127 ++++++++++++++++++ src/QuAcK/UG0W0.f90 | 46 +++---- src/QuAcK/UHF.f90 | 4 +- src/QuAcK/URPAx.f90 | 2 + src/QuAcK/USigmaC.f90 | 6 +- src/QuAcK/dUSigmaC.f90 | 6 +- src/QuAcK/print_UG0W0.f90 | 9 +- src/QuAcK/unrestricted_Bethe_Salpeter.f90 | 47 ++++--- .../unrestricted_Bethe_Salpeter_A_matrix.f90 | 24 ++-- .../unrestricted_Bethe_Salpeter_B_matrix.f90 | 24 ++-- src/QuAcK/unrestricted_linear_response.f90 | 18 +-- 16 files changed, 251 insertions(+), 135 deletions(-) create mode 100644 src/QuAcK/UCIS.f90 diff --git a/input/basis b/input/basis index 6f3d2a9..b2b2293 100644 --- a/input/basis +++ b/input/basis @@ -1,39 +1,30 @@ -1 10 +1 6 S 8 - 1 24350.0000000 0.0005020 - 2 3650.0000000 0.0038810 - 3 829.6000000 0.0199970 - 4 234.0000000 0.0784180 - 5 75.6100000 0.2296760 - 6 26.7300000 0.4327220 - 7 9.9270000 0.3506420 - 8 1.1020000 -0.0076450 + 1 1469.0000000 0.0007660 + 2 220.5000000 0.0058920 + 3 50.2600000 0.0296710 + 4 14.2400000 0.1091800 + 5 4.5810000 0.2827890 + 6 1.5800000 0.4531230 + 7 0.5640000 0.2747740 + 8 0.0734500 0.0097510 S 8 - 1 24350.0000000 -0.0001180 - 2 3650.0000000 -0.0009150 - 3 829.6000000 -0.0047370 - 4 234.0000000 -0.0192330 - 5 75.6100000 -0.0603690 - 6 26.7300000 -0.1425080 - 7 9.9270000 -0.1777100 - 8 1.1020000 0.6058360 + 1 1469.0000000 -0.0001200 + 2 220.5000000 -0.0009230 + 3 50.2600000 -0.0046890 + 4 14.2400000 -0.0176820 + 5 4.5810000 -0.0489020 + 6 1.5800000 -0.0960090 + 7 0.5640000 -0.1363800 + 8 0.0734500 0.5751020 S 1 - 1 2.8360000 1.0000000 -S 1 - 1 0.3782000 1.0000000 + 1 0.0280500 1.0000000 P 3 - 1 54.7000000 0.0171510 - 2 12.4300000 0.1076560 - 3 3.6790000 0.3216810 + 1 1.5340000 0.0227840 + 2 0.2749000 0.1391070 + 3 0.0736200 0.5003750 P 1 - 1 1.1430000 1.0000000 -P 1 - 1 0.3300000 1.0000000 + 1 0.0240300 1.0000000 D 1 - 1 4.0140000 1.0000000 -D 1 - 1 1.0960000 1.0000000 -F 1 - 1 2.5440000 1.0000000 - + 1 0.1239000 1.0000000 diff --git a/input/methods b/input/methods index 478a02f..ba17592 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF MOM - T F F + F T F # MP2 MP3 MP2-F12 F F F # CCD CCSD CCSD(T) @@ -13,7 +13,7 @@ # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0 evGW qsGW - F F T + T F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/molecule b/input/molecule index edeba31..058d6dd 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 5 5 0 0 + 1 2 1 0 0 # Znuc x y z - Ne 0.0 0.0 0.0 + Li 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 1c70680..c9a5a65 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,3 @@ 1 - Ne 0.0000000000 0.0000000000 0.0000000000 + Li 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/options b/input/options index 6af6fdb..55538aa 100644 --- a/input/options +++ b/input/options @@ -5,14 +5,14 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: singlet triplet spin_conserved spin_flip TDA - T T T T F + T T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 - 256 0.00001 T 5 T 0.0 F F F F F + 256 0.00001 T 5 F 0.0 F F F F F # ACFDT: AC Kx XBS - F F T + T F T # BSE: BSE dBSE dTDA evDyn - T F T T + T T T T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/UCIS.f90 b/src/QuAcK/UCIS.f90 new file mode 100644 index 0000000..d27a61b --- /dev/null +++ b/src/QuAcK/UCIS.f90 @@ -0,0 +1,127 @@ +subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,eHF) + +! Perform configuration interaction single calculation` + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + 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) :: nS(nspin) + double precision,intent(in) :: eHF(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) :: ERI_abab(nBas,nBas,nBas,nBas) + +! Local variables + + logical :: dump_matrix = .false. + logical :: dump_trans = .false. + integer :: ispin + double precision :: lambda + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: A_sc(:,:) + double precision,allocatable :: Omega_sc(:) + + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: A_sf(:,:) + double precision,allocatable :: Omega_sf(:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Configuration Interaction Singles |' + write(*,*)'************************************************' + write(*,*) + +! Adiabatic connection scaling + + lambda = 1d0 + +!----------------------------! +! Spin-conserved transitions ! +!----------------------------! + + if(spin_conserved) then + + ispin = 1 + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(A_sc(nS_sc,nS_sc),Omega_sc(nS_sc)) + + call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sc) + + if(dump_matrix) then + print*,'CIS matrix (spin-conserved transitions)' + call matout(nS_sc,nS_sc,A_sc) + write(*,*) + endif + + call diagonalize_matrix(nS_sc,A_sc,Omega_sc) + call print_excitation('UCIS ',5,nS_sc,Omega_sc) + + if(dump_trans) then + print*,'Spin-conserved CIS transition vectors' + call matout(nS_sc,nS_sc,A_sc) + write(*,*) + endif + + deallocate(A_sc,Omega_sc) + + endif + +!-----------------------! +! Spin-flip transitions ! +!-----------------------! + + if(spin_flip) then + + ispin = 2 + + ! 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 + + allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf)) + + call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,lambda,eHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sf) + + if(dump_matrix) then + print*,'CIS matrix (spin-conserved transitions)' + call matout(nS_sf,nS_sf,A_sf) + write(*,*) + endif + + call diagonalize_matrix(nS_sf,A_sf,Omega_sf) + call print_excitation('UCIS ',6,nS_sf,Omega_sf) + + if(dump_trans) then + print*,'Spin-flip CIS transition vectors' + call matout(nS_sf,nS_sf,A_sf) + write(*,*) + endif + + deallocate(A_sf,Omega_sf) + + endif + +end subroutine UCIS diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index 10b2b58..5df2b99 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -47,13 +47,16 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev logical :: print_W = .true. integer :: is integer :: ispin - double precision :: EcRPA(nspin) + double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision,allocatable :: SigC(:,:) double precision,allocatable :: Z(:,:) integer :: nS_aa,nS_bb,nS_sc - double precision,allocatable :: Omega_sc(:),XpY_sc(:,:),XmY_sc(:,:),rho_sc(:,:,:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:,:) double precision,allocatable :: eGWlin(:,:) @@ -101,8 +104,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc), & - rho_sc(nBas,nBas,nS_sc,nspin),eGWlin(nBas,nspin)) + allocate(SigC(nBas,nspin),Z(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc), & + rho_RPA(nBas,nBas,nS_sc,nspin),eGWlin(nBas,nspin)) !-------------------! ! Compute screening ! @@ -113,27 +116,27 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ispin = 1 call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) + eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - if(print_W) call print_excitation('RPA@UHF',5,nS_sc,Omega_sc) + if(print_W) call print_excitation('RPA@UHF',5,nS_sc,OmRPA) !----------------------! ! Excitation densities ! !----------------------! - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,rho_sc) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) !---------------------! ! Compute self-energy ! !---------------------! - call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Omega_sc,rho_sc,SigC) + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC) !--------------------------------! ! Compute renormalization factor ! !--------------------------------! - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Omega_sc,rho_sc,Z) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) !-----------------------------------! ! Solve the quasi-particle equation ! @@ -153,31 +156,24 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Find graphical solution of the QP equation do is=1,nspin - call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),Omega_sc, & - rho_sc,eGWlin(:,is),eGW(:,is)) + call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),OmRPA, & + rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is)) end do end if +! Compute RPA correlation energy + + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + ! Dump results call print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA) -! Compute the RPA correlation energy - - call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(ispin) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA(ispin) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - ! Free memory - deallocate(Omega_sc,XpY_sc,XmY_sc,rho_sc) + deallocate(OmRPA,XpY_RPA,XmY_RPA,rho_RPA) ! Perform BSE calculation @@ -185,7 +181,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & - eHF,eGW,EcRPA,EcBSE) + eHF,eGW,EcBSE) ! if(exchange_kernel) then ! diff --git a/src/QuAcK/UHF.f90 b/src/QuAcK/UHF.f90 index 8559d03..34eb9ef 100644 --- a/src/QuAcK/UHF.f90 +++ b/src/QuAcK/UHF.f90 @@ -166,7 +166,8 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH n_diis = min(n_diis+1,max_diis) do ispin=1,nspin - call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin)) + if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin), & + err(:,:,ispin),F(:,:,ispin)) end do ! Reset DIIS if required @@ -232,7 +233,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH ! Compute final UHF energy - call matout(nBas,2,e) call print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) end subroutine UHF diff --git a/src/QuAcK/URPAx.f90 b/src/QuAcK/URPAx.f90 index 149c62d..f351ca6 100644 --- a/src/QuAcK/URPAx.f90 +++ b/src/QuAcK/URPAx.f90 @@ -78,6 +78,7 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO call print_excitation('URPAx ',5,nS_sc,Omega_sc) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + deallocate(Omega_sc,XpY_sc,XmY_sc) endif @@ -100,6 +101,7 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO call print_excitation('URPAx ',6,nS_sf,Omega_sf) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + deallocate(Omega_sf,XpY_sf,XmY_sf) endif diff --git a/src/QuAcK/USigmaC.f90 b/src/QuAcK/USigmaC.f90 index d701cd7..c5e0817 100644 --- a/src/QuAcK/USigmaC.f90 +++ b/src/QuAcK/USigmaC.f90 @@ -18,7 +18,7 @@ double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) integer,intent(in) :: nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -34,14 +34,14 @@ double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) do i=nC+1,nO do jb=1,nS eps = w - e(i) + Omega(jb) - USigmaC = uSigmaC + rho(p,i,jb,1)**2*eps/(eps**2 + eta**2) + USigmaC = uSigmaC + rho(p,i,jb)**2*eps/(eps**2 + eta**2) end do end do do a=nO+1,nBas-nR do jb=1,nS eps = w - e(a) - Omega(jb) - USigmaC = USigmaC + rho(p,a,jb,1)**2*eps/(eps**2 + eta**2) + USigmaC = USigmaC + rho(p,a,jb)**2*eps/(eps**2 + eta**2) end do end do diff --git a/src/QuAcK/dUSigmaC.f90 b/src/QuAcK/dUSigmaC.f90 index cbd0cdb..418e131 100644 --- a/src/QuAcK/dUSigmaC.f90 +++ b/src/QuAcK/dUSigmaC.f90 @@ -18,7 +18,7 @@ double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) integer,intent(in) :: nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -34,7 +34,7 @@ double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) do i=nC+1,nO do jb=1,nS eps = w - e(i) + Omega(jb) - dUSigmaC = dUSigmaC + rho(p,i,jb,1)**2*(eps/(eps**2 + eta**2))**2 + dUSigmaC = dUSigmaC + rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2 end do end do @@ -43,7 +43,7 @@ double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) do a=nO+1,nBas-nR do jb=1,nS eps = w - e(a) - Omega(jb) - dUSigmaC = dUSigmaC + rho(p,a,jb,1)**2*(eps/(eps**2 + eta**2))**2 + dUSigmaC = dUSigmaC + rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2 end do end do diff --git a/src/QuAcK/print_UG0W0.f90 b/src/QuAcK/print_UG0W0.f90 index 3cd21cd..f678a1e 100644 --- a/src/QuAcK/print_UG0W0.f90 +++ b/src/QuAcK/print_UG0W0.f90 @@ -33,9 +33,10 @@ subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA) write(*,*)' Unrestricted one-shot G0W0 calculation (eV)' write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' + write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') & + '|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|' write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') & - '|','#','|','e_HF up','e_HF dw','|','Sig_c up','Sig_c dw','|', & - 'Z up','Z dw','|','e_QP up','e_QP dw','|' + '|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|' write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' @@ -52,8 +53,8 @@ subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA) write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'RPA@HF total energy =',ENuc + EHF + EcRPA - write(*,'(2X,A30,F15.6)') 'RPA@HF correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' write(*,*) diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 index 992ee9b..5ac6e9c 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 @@ -1,6 +1,6 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & - eW,eGW,EcRPA,EcBSE) + eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -38,10 +38,12 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, integer :: isp_W integer :: nS_aa,nS_bb,nS_sc - double precision,allocatable :: OmRPA_sc(:) - double precision,allocatable :: XpY_RPA_sc(:,:) - double precision,allocatable :: XmY_RPA_sc(:,:) - double precision,allocatable :: rho_RPA_sc(:,:,:,:) + 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 :: OmBSE_sc(:) double precision,allocatable :: XpY_BSE_sc(:,:) double precision,allocatable :: XmY_BSE_sc(:,:) @@ -53,48 +55,46 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, ! Output variables - double precision,intent(out) :: EcRPA(nspin) double precision,intent(out) :: EcBSE(nspin) -!----------------------------! -! Spin-conserved excitations ! -!----------------------------! - - isp_W = 1 - ! Memory allocation nS_aa = nS(1) nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(OmRPA_sc(nS_sc),XpY_RPA_sc(nS_sc,nS_sc),XmY_RPA_sc(nS_sc,nS_sc),rho_RPA_sc(nBas,nBas,nS_sc,nspin)) + allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) + +!--------------------------! +! Spin-conserved screening ! +!--------------------------! + + isp_W = 1 + EcRPA = 0d0 ! Compute spin-conserved RPA screening call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA_sc,rho_RPA_sc,EcRPA(isp_W), & - OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc) + eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) -! call print_excitation('RPA@UG0W0',5,nS_sc,OmRPA_sc) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb, & - XpY_RPA_sc,rho_RPA_sc) + +!----------------------------! +! Spin-conserved excitations ! +!----------------------------! if(spin_conserved) then ispin = 1 - EcBSE(ispin) = 0d0 allocate(OmBSE_sc(nS_sc),XpY_BSE_sc(nS_sc,nS_sc),XmY_BSE_sc(nS_sc,nS_sc)) ! Compute spin-conserved BSE excitation energies - OmBSE_sc(:) = OmRPA_sc(:) - call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA_sc,rho_RPA_sc,EcBSE(ispin), & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), & OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc) @@ -128,7 +128,6 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, if(spin_flip) then ispin = 2 - EcBSE(ispin) = 0d0 ! Memory allocation @@ -140,7 +139,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf)) call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA_sc,rho_RPA_sc,EcBSE(ispin), & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), & OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) call print_excitation('BSE@UG0W0',6,nS_sf,OmBSE_sf) diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 index 55d9a25..1b1e415 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & +subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism @@ -17,15 +17,15 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt - integer,intent(in) :: nSsc + integer,intent(in) :: nS_sc double precision,intent(in) :: eta double precision,intent(in) :: lambda 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) :: ERI_abab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nSsc) - double precision,intent(in) :: rho(nBas,nBas,nSsc,nspin) + double precision,intent(in) :: Omega(nS_sc) + double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin) ! Local variables @@ -55,12 +55,12 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n jb = jb + 1 chi = 0d0 - do kc=1,nSsc + do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps enddo - A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 4d0*lambda*chi + A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi enddo enddo @@ -79,12 +79,12 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n jb = jb + 1 chi = 0d0 - do kc=1,nSsc + do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps enddo - A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 4d0*lambda*chi + A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 2d0*lambda*chi enddo enddo @@ -111,12 +111,12 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n jb = jb + 1 chi = 0d0 - do kc=1,nSsc + do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps enddo - A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 4d0*lambda*chi + A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 2d0*lambda*chi end do end do @@ -135,12 +135,12 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n jb = jb + 1 chi = 0d0 - do kc=1,nSsc + do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Omega(kc)/eps enddo - A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_abab(b,i,a,j) + 4d0*lambda*chi + A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_abab(b,i,a,j) + 2d0*lambda*chi end do end do diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 index e954814..2f559ee 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & +subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,B_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response @@ -17,15 +17,15 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt - integer,intent(in) :: nSsc + integer,intent(in) :: nS_sc double precision,intent(in) :: eta double precision,intent(in) :: lambda 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) :: ERI_abab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nSsc) - double precision,intent(in) :: rho(nBas,nBas,nSsc,nspin) + double precision,intent(in) :: Omega(nS_sc) + double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin) ! Local variables @@ -55,12 +55,12 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n jb = jb + 1 chi = 0d0 - do kc=1,nSsc + do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps enddo - B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 4d0*lambda*chi + B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi enddo enddo @@ -80,12 +80,12 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n jb = jb + 1 chi = 0d0 - do kc=1,nSsc + do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps enddo - B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 4d0*lambda*chi + B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 2d0*lambda*chi enddo enddo @@ -113,12 +113,12 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n jb = jb + 1 chi = 0d0 - do kc=1,nSsc + do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps enddo - B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 4d0*lambda*chi + B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 2d0*lambda*chi end do end do @@ -137,12 +137,12 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n jb = jb + 1 chi = 0d0 - do kc=1,nSsc + do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Omega(kc)/eps enddo - B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_abab(b,j,i,a) + 4d0*lambda*chi + B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_abab(b,j,i,a) + 2d0*lambda*chi end do end do diff --git a/src/QuAcK/unrestricted_linear_response.f90 b/src/QuAcK/unrestricted_linear_response.f90 index b61cb38..e816162 100644 --- a/src/QuAcK/unrestricted_linear_response.f90 +++ b/src/QuAcK/unrestricted_linear_response.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & - e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,EcRPA,Omega,XpY,XmY) +subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,Omega,XpY,XmY) ! Compute linear response for unrestricted formalism @@ -21,7 +21,7 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt - integer,intent(in) :: nSsc + 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) @@ -29,8 +29,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega_RPA(nSsc) - double precision,intent(in) :: rho_RPA(nBas,nBas,nSsc,nspin) + double precision,intent(in) :: OmRPA(nS_sc) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin) ! Local variables @@ -61,8 +61,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A) if(BSE) & - call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,A) + call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,A) ! Tamm-Dancoff approximation @@ -73,8 +73,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B) if(BSE) & - call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,B) + call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,B) end if