diff --git a/input/methods b/input/methods index adb472a..b1f2358 100644 --- a/input/methods +++ b/input/methods @@ -9,11 +9,11 @@ # CIS* CID CISD F F F # RPA* RPAx* ppRPA - F T F + F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW - F F F + T F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 1a82a15..1371bb5 100644 --- a/input/options +++ b/input/options @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - F T T F + T F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 9c5c095..4a906c9 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -58,7 +58,9 @@ program QuAcK double precision,allocatable :: Hc(:,:) double precision,allocatable :: H(:,:) double precision,allocatable :: X(:,:) - double precision,allocatable :: dipole_int(:,:,:,:) + double precision,allocatable :: dipole_int(:,:,:) + double precision,allocatable :: dipole_int_aa(:,:,:) + double precision,allocatable :: dipole_int_bb(:,:,:) double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) integer :: ixyz @@ -233,8 +235,7 @@ program QuAcK ! Memory allocation for one- and two-electron integrals allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), & - S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas), & - dipole_int(nBas,nBas,ncart,nspin),ERI_AO(nBas,nBas,nBas,nBas)) + S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas)) ! Read integrals @@ -341,11 +342,13 @@ program QuAcK ! Read and transform dipole-related integrals - call read_dipole_integrals(nBas,dipole_int) + allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart)) + + call read_dipole_integrals(nBas,dipole_int_aa) + call read_dipole_integrals(nBas,dipole_int_bb) do ixyz=1,ncart - do ispin=1,nspin - call AOtoMO_transform(nBas,cHF(:,:,ispin),dipole_int(:,:,ixyz,ispin)) - end do + call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz)) end do ! Memory allocation @@ -399,10 +402,10 @@ program QuAcK ! Read and transform dipole-related integrals - ispin = 1 + allocate(dipole_int(nBas,nBas,ncart)) call read_dipole_integrals(nBas,dipole_int) do ixyz=1,ncart - call AOtoMO_transform(nBas,cHF,dipole_int(:,:,ixyz,ispin)) + call AOtoMO_transform(nBas,cHF,dipole_int(:,:,ixyz)) end do ! 4-index transform @@ -696,7 +699,7 @@ program QuAcK if(unrestricted) then call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & - ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,eHF) + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF) else @@ -822,9 +825,9 @@ program QuAcK call cpu_time(start_G0W0) if(unrestricted) then - 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,ERI_MO_abab,dipole_int,PHF,cHF,eHF,eG0W0) + 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,ERI_MO_abab, & + dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0) else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & @@ -849,9 +852,10 @@ program QuAcK call cpu_time(start_evGW) if(unrestricted) then - call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc, & - ERHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,PHF,cHF,eHF,eG0W0) + call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & + G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc, & + EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb, & + PHF,cHF,eHF,eG0W0) else diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index dd5fa5b..53e3e41 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,ERI_abab,dipole_int,PHF,cHF,eHF,eGW) +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,ERI_abab, & + dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eGW) ! Perform unrestricted G0W0 calculation @@ -41,7 +41,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev 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) :: dipole_int(nBas,nBas,ncart,nspin) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) ! Local variables @@ -181,7 +182,7 @@ 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,ERI_abab,eHF,eGW,EcBSE) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,eHF,eGW,EcBSE) ! if(exchange_kernel) then ! diff --git a/src/QuAcK/URPAx.f90 b/src/QuAcK/URPAx.f90 index ad58315..d7cf08c 100644 --- a/src/QuAcK/URPAx.f90 +++ b/src/QuAcK/URPAx.f90 @@ -1,5 +1,5 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,e) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e) ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism @@ -28,7 +28,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n 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) :: dipole_int(nBas,nBas,ncart,nspin) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) ! Local variables @@ -78,7 +79,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPAx ',5,nS_sc,Omega_sc) - call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_sc,dipole_int,Omega_sc,XpY_sc,XmY_sc) + call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + Omega_sc,XpY_sc,XmY_sc) deallocate(Omega_sc,XpY_sc,XmY_sc) diff --git a/src/QuAcK/UdRPA.f90 b/src/QuAcK/UdRPA.f90 index 5d1e7f6..fd1aa0e 100644 --- a/src/QuAcK/UdRPA.f90 +++ b/src/QuAcK/UdRPA.f90 @@ -1,5 +1,5 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,e) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e) ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism @@ -28,7 +28,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n 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) :: dipole_int(nBas,nBas,ncart,nspin) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) ! Local variables @@ -78,7 +79,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPA ',5,nS_sc,Omega_sc) -! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + Omega_sc,XpY_sc,XmY_sc) endif diff --git a/src/QuAcK/evUGW.f90 b/src/QuAcK/evUGW.f90 index 24f8e0e..322e13e 100644 --- a/src/QuAcK/evUGW.f90 +++ b/src/QuAcK/evUGW.f90 @@ -1,6 +1,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, & - ERHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,PHF,cHF,eHF,eG0W0) + ERHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -46,7 +46,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE 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) :: dipole_int(nBas,nBas,ncart,nspin) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) ! Local variables @@ -255,7 +256,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE 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,ERI_abab,eGW,eGW,EcBSE) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,eGW,eGW,EcBSE) ! if(exchange_kernel) then diff --git a/src/QuAcK/print_transition_vectors.f90 b/src/QuAcK/print_transition_vectors.f90 index 6986eec..76ee45e 100644 --- a/src/QuAcK/print_transition_vectors.f90 +++ b/src/QuAcK/print_transition_vectors.f90 @@ -21,6 +21,7 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int, ! Local variables + logical :: debug = .false. integer :: ia,jb,i,j,a,b integer :: ixyz integer,parameter :: maxS = 10 @@ -54,21 +55,25 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int, end do f(:,:) = sqrt(2d0)*f(:,:) - write(*,*) '------------------------' - write(*,*) ' Dipole moments (X Y Z) ' - write(*,*) '------------------------' - call matout(nS,ncart,f) - write(*,*) - - do ia=1,nS - os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) - end do - - write(*,*) '----------------------' - write(*,*) ' Oscillator strengths ' - write(*,*) '----------------------' - call matout(nS,1,os) - write(*,*) + if(debug) then + + write(*,*) '------------------------' + write(*,*) ' Dipole moments (X Y Z) ' + write(*,*) '------------------------' + call matout(nS,ncart,f) + write(*,*) + + do ia=1,nS + os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) + end do + + write(*,*) '----------------------' + write(*,*) ' Oscillator strengths ' + write(*,*) '----------------------' + call matout(nS,1,os) + write(*,*) + + end if end if diff --git a/src/QuAcK/print_unrestricted_transition_vectors.f90 b/src/QuAcK/print_unrestricted_transition_vectors.f90 index 2d6e8bf..27c49ad 100644 --- a/src/QuAcK/print_unrestricted_transition_vectors.f90 +++ b/src/QuAcK/print_unrestricted_transition_vectors.f90 @@ -1,4 +1,5 @@ -subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,nSt,dipole_int,Omega,XpY,XmY) +subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb, & + Omega,XpY,XmY) ! Print transition vectors for linear response calculation @@ -14,14 +15,18 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) integer,intent(in) :: nS(nspin) + integer,intent(in) :: nSa + integer,intent(in) :: nSb integer,intent(in) :: nSt - double precision :: dipole_int(nBas,nBas,ncart,nspin) + double precision :: dipole_int_aa(nBas,nBas,ncart) + double precision :: dipole_int_bb(nBas,nBas,ncart) double precision,intent(in) :: Omega(nSt) double precision,intent(in) :: XpY(nSt,nSt) double precision,intent(in) :: XmY(nSt,nSt) ! Local variables + logical :: debug = .false. integer :: ia,jb,i,j,a,b integer :: ixyz integer :: ispin @@ -43,35 +48,47 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n f(:,:) = 0d0 if(spin_allowed) then - do ispin=1,nspin - do ia=1,nSt - do ixyz=1,ncart - jb = 0 - do j=nC(ispin)+1,nO(ispin) - do b=nO(ispin)+1,nBas-nR(ispin) - jb = jb + 1 - f(ia,ixyz) = f(ia,ixyz) + dipole_int(j,b,ixyz,ispin)*XpY(ia,jb) - end do + do ia=1,nSt + do ixyz=1,ncart + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + f(ia,ixyz) = f(ia,ixyz) + dipole_int_aa(j,b,ixyz)*XpY(ia,jb) end do end do + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + f(ia,ixyz) = f(ia,ixyz) + dipole_int_bb(j,b,ixyz)*XpY(ia,nSa+jb) + end do + end do + end do end do - write(*,*) '----------------' - write(*,*) ' Dipole moments ' - write(*,*) '----------------' - call matout(nSt,ncart,f(:,:)) - write(*,*) - - do ia=1,nSt - os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) - end do - - write(*,*) '----------------------' - write(*,*) ' Oscillator strengths ' - write(*,*) '----------------------' - call matout(nSt,1,os(:)) - write(*,*) + if(debug) then + + write(*,*) '----------------' + write(*,*) ' Dipole moments ' + write(*,*) '----------------' + call matout(nSt,ncart,f(:,:)) + write(*,*) + + do ia=1,nSt + os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) + end do + + write(*,*) '----------------------' + write(*,*) ' Oscillator strengths ' + write(*,*) '----------------------' + call matout(nSt,1,os(:)) + write(*,*) + + end if end if @@ -92,7 +109,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n do j=nC(1)+1,nO(1) do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' -> ',b,' = ',X(jb)/sqrt(2d0) + if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb)/sqrt(2d0) end do end do @@ -100,10 +117,9 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n do j=nC(1)+1,nO(1) do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' <- ',b,' = ',Y(jb)/sqrt(2d0) + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb)/sqrt(2d0) end do end do - write(*,*) ! Spin-down transitions @@ -111,7 +127,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n do j=nC(2)+1,nO(2) do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' -> ',b,' = ',X(jb)/sqrt(2d0) + if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(jb)/sqrt(2d0) end do end do @@ -119,7 +135,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n do j=nC(2)+1,nO(2) do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' <- ',b,' = ',Y(jb)/sqrt(2d0) + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(jb)/sqrt(2d0) end do end do write(*,*) diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 index 5ac6e9c..1b7566a 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,EcBSE) + dipole_int_aa,dipole_int_bb,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -30,7 +30,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, 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) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) ! Local variables @@ -96,8 +97,9 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, 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,rho_RPA,EcBSE(ispin), & OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) - call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc) + call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) !------------------------------------------------- ! Compute the dynamical screening at the BSE level