diff --git a/input/options b/input/options index 7be6358..23ac893 100644 --- a/input/options +++ b/input/options @@ -4,8 +4,8 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 -# spin: singlet triplet spin_conserved spinf_flip TDA - T T T F F +# spin: singlet triplet spin_conserved spin_flip TDA + T T T T F # 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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 6f62482..d65bfcc 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -775,7 +775,7 @@ program QuAcK call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, & spin_conserved,spin_flip,linGW,eta_GW,nBas,nC,nO,nV,nR,nS, & - ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,PHF,cHF,eHF,eG0W0) + ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,PHF,cHF,eHF,eG0W0) else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index 07de370..453f326 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -1,6 +1,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, & spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS, & - ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,PHF,cHF,eHF,eGW) + ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,PHF,cHF,eHF,eGW) ! Perform unrestricted G0W0 calculation @@ -40,6 +40,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev 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 @@ -111,7 +112,7 @@ 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,1d0, & - eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) + eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) if(print_W) call print_excitation('RPA@UHF',5,nS_sc,Omega_sc) @@ -164,7 +165,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! 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,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -182,7 +183,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev if(BSE) then 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,eHF,eGW,EcRPA,EcBSE) + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & + eHF,eGW,EcRPA,EcBSE) ! if(exchange_kernel) then ! diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 index ef85215..38b36ee 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 @@ -1,5 +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,eW,eGW,EcRPA,EcBSE) + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & + eW,eGW,EcRPA,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -28,6 +29,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, 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 @@ -71,16 +73,19 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, ! 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,1d0, & - eW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA_sc,EcRPA(ispin),OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc) + eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcRPA(ispin), & + OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc) - 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) + 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) ! 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,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA_sc,EcBSE(ispin),OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcBSE(ispin), & + OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc) diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 index 82a551b..e3f5813 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 @@ -1,4 +1,6 @@ -subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A_lr) +subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,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 @@ -20,6 +22,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt 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(nSt) double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 index 476cde8..250978d 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 @@ -1,4 +1,6 @@ -subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B_lr) +subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & + Omega,rho,B_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response @@ -20,6 +22,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt 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(nSt) double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) diff --git a/src/QuAcK/unrestricted_linear_response.f90 b/src/QuAcK/unrestricted_linear_response.f90 index b8d8d52..29a6dfe 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,lambda, & - e,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY) + e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho,EcRPA,Omega,XpY,XmY) ! Compute linear response for unrestricted formalism @@ -27,6 +27,7 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, 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 @@ -53,20 +54,24 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ! Build A and B matrices - call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,A) + call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A) if(BSE) & - call unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A) + call unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A) ! Tamm-Dancoff approximation B = 0d0 if(.not. TDA) then - call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,B) + call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B) if(BSE) & - call unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B) + call unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,B) end if diff --git a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 index 75c63b0..ded4c78 100644 --- a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 +++ b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - e,ERI_aaaa,ERI_aabb,ERI_bbbb,A_lr) + e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_lr) ! Compute linear response @@ -23,6 +23,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa 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 @@ -46,7 +47,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa if(ispin == 1) then - ! alpha-alpha block + ! aaaa block ia = 0 do i=nC(1)+1,nO(1) @@ -65,7 +66,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end do end do - ! alpha-beta block + ! aabb block ia = 0 do i=nC(1)+1,nO(1) @@ -83,7 +84,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end do end do - ! beta-alpha block + ! bbaa block ia = 0 do i=nC(2)+1,nO(2) @@ -101,7 +102,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end do end do - ! beta-beta block + ! bbbb block ia = 0 do i=nC(2)+1,nO(2) @@ -128,7 +129,45 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa if(ispin == 2) then - print*,'spin-flip transition NYI' + A_lr(:,:) = 0d0 + + ! abab block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + + lambda*ERI_abab(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a) + + end do + end do + end do + end do + + ! baba block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + + lambda*ERI_abab(b,i,j,a) - (1d0 - delta_dRPA)*lambda*ERI_abab(b,i,a,j) + + end do + end do + end do + end do end if diff --git a/src/QuAcK/unrestricted_linear_response_B_matrix.f90 b/src/QuAcK/unrestricted_linear_response_B_matrix.f90 index 400a498..5ac615d 100644 --- a/src/QuAcK/unrestricted_linear_response_B_matrix.f90 +++ b/src/QuAcK/unrestricted_linear_response_B_matrix.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,B_lr) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B_lr) ! Compute linear response @@ -22,6 +22,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa 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 @@ -40,12 +41,12 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa if(dRPA) delta_dRPA = 1d0 !----------------------------------------------- -! Build A matrix for spin-conserving transitions +! Build B matrix for spin-conserving transitions !----------------------------------------------- if(ispin == 1) then - ! alpha-alpha block + ! aaaa block ia = 0 do i=nC(1)+1,nO(1) @@ -63,7 +64,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end do end do - ! alpha-beta block + ! aabb block ia = 0 do i=nC(1)+1,nO(1) @@ -81,7 +82,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end do end do - ! beta-alpha block + ! bbaa block ia = 0 do i=nC(2)+1,nO(2) @@ -99,7 +100,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end do end do - ! beta-beta block + ! bbbb block ia = 0 do i=nC(2)+1,nO(2) @@ -120,12 +121,48 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end if !----------------------------------------------- -! Build A matrix for spin-flip transitions +! Build B matrix for spin-flip transitions !----------------------------------------------- if(ispin == 2) then - print*,'spin-flip transition NYI' + B_lr(:,:) = 0d0 + + ! abab block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + B_lr(ia,jb) = lambda*ERI_abab(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_abab(i,j,b,a) + + end do + end do + end do + end do + + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + B_lr(nSa+ia,nSa+jb) = lambda*ERI_abab(j,i,b,a) - (1d0 - delta_dRPA)*lambda*ERI_abab(j,i,a,b) + + end do + end do + end do + end do end if