diff --git a/src/dftbp/dftbplus/initprogram.F90 b/src/dftbp/dftbplus/initprogram.F90 index 29157a44ac..a91450e81c 100644 --- a/src/dftbp/dftbplus/initprogram.F90 +++ b/src/dftbp/dftbplus/initprogram.F90 @@ -1552,8 +1552,6 @@ subroutine initProgramVariables(this, input, env) allocate(this%speciesMass(this%nType)) this%speciesMass(:) = input%slako%mass(:) case(hamiltonianTypes%xtb) - ! TODO - ! call error("xTB calculation currently not supported") allocate(this%speciesMass(this%nType)) this%speciesMass(:) = getAtomicMass(this%speciesName) end select @@ -1600,8 +1598,6 @@ subroutine initProgramVariables(this, input, env) this%cutOff%mCutOff = max(this%cutOff%mCutOff, this%repulsive%getRCutOff()) end if case(hamiltonianTypes%xtb) - ! TODO - ! call error("xTB calculation currently not supported") this%cutOff%skCutoff = this%tblite%getRCutoff() this%cutOff%mCutoff = this%cutOff%skCutoff end select @@ -2527,8 +2523,8 @@ subroutine initProgramVariables(this, input, env) call error("PP-RPA does not support ${ERR}$") end if #:endfor - #:for VAR, ERR in [("tShellResolved","shell resolved hamiltonians"),& - & ("tDampH","H damping")] + #:for VAR, ERR in [("tShellResolved", "shell resolved hamiltonians"),& + & ("tDampH", "H damping")] if (input%ctrl%${VAR}$) then call error("PP-RPA does not support ${ERR}$") end if @@ -4611,8 +4607,7 @@ subroutine initHubbardUs_(input, orb, hamiltonianType, hubbU) @:ASSERT(size(input%slako%skHubbU, dim=2) == nSpecies) hubbU(:,:) = input%slako%skHubbU(1:orb%mShell, :) case(hamiltonianTypes%xtb) - ! TODO - ! call error("xTB calculation currently not supported") + ! handled elsewhere end select if (allocated(input%ctrl%hubbU)) then diff --git a/src/dftbp/timedep/linrespcommon.F90 b/src/dftbp/timedep/linrespcommon.F90 index 8d554e5db6..76080a9467 100644 --- a/src/dftbp/timedep/linrespcommon.F90 +++ b/src/dftbp/timedep/linrespcommon.F90 @@ -405,7 +405,7 @@ end subroutine getSqrOcc !> (spin polarized case) or \Omega^{S/T}_ia,jb * v_jb (singlet/triplet) !> !> For the standard RPA, (A+B)_ias,jbt * v_jbt needs to be computed (similar for singlet/triplet) - !> (see definitions in Marc Casida, in Recent Advances in Density Functional Methods, + !> (see definitions in Mark Casida, in Recent Advances in Density Functional Methods, !> World Scientific, 1995, Part I, p. 155.) !> Note: we actually compute sqrt(n_is-n_as) (A+B)_ias,jbt sqrt(n_jt-n_bt), with the !> occupations n_is, correct for finite T. @@ -538,8 +538,7 @@ subroutine actionAplusB(spin, wij, sym, win, nocc_ud, nvir_ud, nxoo_ud, nxvv_ud, ! product charges with the v*wn product, i.e. Q * v*wn oTmp(:) = 0.0_dp - call transChrg%qMatVec(iAtomStart, ovrXev, grndEigVecs, getIA, win, & - & vTmp * sqrOccIA, oTmp) + call transChrg%qMatVec(iAtomStart, ovrXev, grndEigVecs, getIA, win, vTmp * sqrOccIA, oTmp) if (.not.spin) then !-----------spin-unpolarized systems-------------- @@ -693,7 +692,7 @@ end subroutine actionAplusB !> Multiplies the excitation supermatrix with a supervector. !> (A-B)_ias,jbt * v_jbt is computed (and similar for singlet/triplet) - !> (see definitions in Marc Casida, in Recent Advances in Density Functional Methods, + !> (see definitions in Mark Casida, in Recent Advances in Density Functional Methods, !> World Scientific, 1995, Part I, p. 155.) !> See also Dominguez JCTC 9 4901 (2013), Kranz JCTC 13 1737 (2017) for DFTB specifics. subroutine actionAminusB(spin, wij, win, nocc_ud, nvir_ud, nxoo_ud, nxvv_ud, nxov_ud, nxov_rd,& diff --git a/src/dftbp/timedep/linrespgrad.F90 b/src/dftbp/timedep/linrespgrad.F90 index 80a3895aa2..35224e49e7 100644 --- a/src/dftbp/timedep/linrespgrad.F90 +++ b/src/dftbp/timedep/linrespgrad.F90 @@ -558,11 +558,11 @@ subroutine LinRespGrad_old(this, iAtomStart, grndEigVecs, grndEigVal, sccCalc, d & transChrg, this%testArnoldi, eval, xpy, xmy, this%onSiteMatrixElements, orb,& & tHybridXc, tZVector) case (linrespSolverTypes%stratmann) - call buildAndDiagExcMatrixStratmann(this%tSpin, this%subSpaceFactorStratmann, wij(:nxov_rd),& - & sym, win, nocc_ud, nvir_ud, nxoo_ud, nxvv_ud, nxov_ud, nxov_rd, iaTrans, getIA,& - & getIJ, getAB, iAtomStart, ovrXev, grndEigVecs, filling, sqrOccIA(:nxov_rd), gammaMat,& - & species0, this%spinW, transChrg, eval, xpy, xmy, this%onSiteMatrixElements, orb,& - & tHybridXc, lrGamma, tZVector) + call buildAndDiagExcMatrixStratmann(this%tSpin, this%subSpaceFactorStratmann,& + & wij(:nxov_rd), sym, win, nocc_ud, nvir_ud, nxoo_ud, nxvv_ud, nxov_ud, nxov_rd,& + & iaTrans, getIA, getIJ, getAB, iAtomStart, ovrXev, grndEigVecs, filling,& + & sqrOccIA(:nxov_rd), gammaMat, species0, this%spinW, transChrg, eval, xpy, xmy,& + & this%onSiteMatrixElements, orb, tHybridXc, lrGamma, tZVector) end select ! Excitation oscillator strengths for resulting states @@ -838,7 +838,7 @@ end subroutine LinRespGrad_old !> w !> [B A] Y = [0 -C] Y !> - !> (see definitions in Marc Casida, in Recent Advances in Density Functional Methods, + !> (see definitions by Mark Casida, in Recent Advances in Density Functional Methods, !> World Scientific, 1995, Part I, p. 155.) !> !> The hermitian EV problem is given by \Omega F = w^2 F, with @@ -1091,7 +1091,7 @@ end subroutine buildAndDiagExcMatrixArpack !> w !> [B A] Y = [0 -C] Y !> - !> (see definitions in Marc Casida, in Recent Advances in Density Functional Methods, + !> (see definitions by Mark Casida, in Recent Advances in Density Functional Methods, !> World Scientific, 1995, Part I, p. 155.) !> !> The RPA eqs are diagonalised by the Stratmann algorithm (JCP 109 8218 (1998). @@ -4819,14 +4819,21 @@ subroutine addNadiaGradients(sym, nxov, natom, species0, iAtomStart, norb, homo, do iSpin = 1, nSpin call symm(PS(:,:,iSpin), 'R', overlap, pc(:,:,iSpin), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) call symm(SPS(:,:,iSpin), 'L', overlap, PS(:,:,iSpin), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) - call symm(DS(:,:,iSpin), 'R', overlap, deltaRho(:,:,iSpin), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) + call symm(DS(:,:,iSpin), 'R', overlap, deltaRho(:,:,iSpin), 'U', 1.0_dp, 0.0_dp, nOrb,& + & nOrb) call symm(SDS(:,:,iSpin), 'L', overlap, DS(:,:,iSpin), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) - call symm(XS(:,:,iSpin,iState), 'R', overlap, xpyas(:,:,iSpin,iState), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) - call symm(SX(:,:,iSpin,iState), 'L', overlap, xpyas(:,:,iSpin,iState), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) - call symm(SXS(:,:,iSpin,iState), 'L', overlap, XS(:,:,iSpin,iState), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) - call symm(YS(:,:,iSpin,iState), 'R', overlap, xmyas(:,:,iSpin,iState), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) - call symm(SY(:,:,iSpin,iState), 'L', overlap, xmyas(:,:,iSpin,iState), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) - call symm(SYS(:,:,iSpin,iState), 'L', overlap, YS(:,:,iSpin,iState), 'U', 1.0_dp, 0.0_dp, nOrb, nOrb) + call symm(XS(:,:,iSpin,iState), 'R', overlap, xpyas(:,:,iSpin,iState), 'U', 1.0_dp,& + & 0.0_dp, nOrb, nOrb) + call symm(SX(:,:,iSpin,iState), 'L', overlap, xpyas(:,:,iSpin,iState), 'U', 1.0_dp,& + & 0.0_dp, nOrb, nOrb) + call symm(SXS(:,:,iSpin,iState), 'L', overlap, XS(:,:,iSpin,iState), 'U', 1.0_dp, 0.0_dp,& + & nOrb, nOrb) + call symm(YS(:,:,iSpin,iState), 'R', overlap, xmyas(:,:,iSpin,iState), 'U', 1.0_dp,& + & 0.0_dp, nOrb, nOrb) + call symm(SY(:,:,iSpin,iState), 'L', overlap, xmyas(:,:,iSpin,iState), 'U', 1.0_dp,& + & 0.0_dp, nOrb, nOrb) + call symm(SYS(:,:,iSpin,iState), 'L', overlap, YS(:,:,iSpin,iState), 'U', 1.0_dp, 0.0_dp,& + & nOrb, nOrb) end do end do @@ -4945,45 +4952,72 @@ subroutine addNadiaGradients(sym, nxov, natom, species0, iAtomStart, norb, homo, do mu = indAlpha, indAlpha1 do nu = indBeta, indBeta1 tmprs = tmprs + & - & ( 2.0_dp * (PS(mu,nu,iSpin) * DS(nu,mu,iSpin) + PS(nu,mu,iSpin) * DS(mu,nu,iSpin)) +& - & SPS(mu,nu,iSpin) * deltaRho(mu,nu,iSpin) + SPS(nu,mu,iSpin) * deltaRho(nu,mu,iSpin) +& - & pc(mu,nu,iSpin) * SDS(mu,nu,iSpin) + pc(nu,mu,iSpin) * SDS(nu,mu,iSpin) ) + & ( 2.0_dp * (PS(mu,nu,iSpin) * DS(nu,mu,iSpin)& + & + PS(nu,mu,iSpin) * DS(mu,nu,iSpin))& + & + SPS(mu,nu,iSpin) * deltaRho(mu,nu,iSpin)& + & + SPS(nu,mu,iSpin) * deltaRho(nu,mu,iSpin)& + & + pc(mu,nu,iSpin) * SDS(mu,nu,iSpin)& + & + pc(nu,mu,iSpin) * SDS(nu,mu,iSpin) ) tmprs = tmprs + & - & ( xpyas(mu,nu,iSpin,1) * SXS(mu,nu,iSpin,2) + xpyas(nu,mu,iSpin,2) * SXS(nu,mu,iSpin,1) +& - & SX(mu,nu,iSpin,1) * XS(mu,nu,iSpin,2) + SX(nu,mu,iSpin,2) * XS(nu,mu,iSpin,1) ) + & ( xpyas(mu,nu,iSpin,1) * SXS(mu,nu,iSpin,2)& + & + xpyas(nu,mu,iSpin,2) * SXS(nu,mu,iSpin,1)& + & + SX(mu,nu,iSpin,1) * XS(mu,nu,iSpin,2)& + & + SX(nu,mu,iSpin,2) * XS(nu,mu,iSpin,1) ) tmprs = tmprs + & - & ( xpyas(mu,nu,iSpin,2) * SXS(mu,nu,iSpin,1) + xpyas(nu,mu,iSpin,1) * SXS(nu,mu,iSpin,2) +& - & SX(mu,nu,iSpin,2) * XS(mu,nu,iSpin,1) + SX(nu,mu,iSpin,1) * XS(nu,mu,iSpin,2) ) + & ( xpyas(mu,nu,iSpin,2) * SXS(mu,nu,iSpin,1)& + & + xpyas(nu,mu,iSpin,1) * SXS(nu,mu,iSpin,2)& + & + SX(mu,nu,iSpin,2) * XS(mu,nu,iSpin,1)& + & + SX(nu,mu,iSpin,1) * XS(nu,mu,iSpin,2) ) tmprs = tmprs + 0.5_dp * & - & ( XS(mu,nu,iSpin,1) * XS(nu,mu,iSpin,2) + XS(nu,mu,iSpin,2) * XS(mu,nu,iSpin,1) +& - & SXS(mu,nu,iSpin,1) * xpyas(nu,mu,iSpin,2) + SXS(nu,mu,iSpin,2) * xpyas(mu,nu,iSpin,1) +& - & xpyas(mu,nu,iSpin,1) * SXS(nu,mu,iSpin,2) + xpyas(nu,mu,iSpin,2) * SXS(mu,nu,iSpin,1) +& - & SX(mu,nu,iSpin,1) * SX(nu,mu,iSpin,2) + SX(nu,mu,iSpin,2) * SX(mu,nu,iSpin,1) ) + & ( XS(mu,nu,iSpin,1) * XS(nu,mu,iSpin,2)& + & + XS(nu,mu,iSpin,2) * XS(mu,nu,iSpin,1)& + & + SXS(mu,nu,iSpin,1) * xpyas(nu,mu,iSpin,2)& + & + SXS(nu,mu,iSpin,2) * xpyas(mu,nu,iSpin,1)& + & + xpyas(mu,nu,iSpin,1) * SXS(nu,mu,iSpin,2)& + & + xpyas(nu,mu,iSpin,2) * SXS(mu,nu,iSpin,1)& + & + SX(mu,nu,iSpin,1) * SX(nu,mu,iSpin,2)& + & + SX(nu,mu,iSpin,2) * SX(mu,nu,iSpin,1) ) tmprs = tmprs + 0.5_dp * & - & ( XS(mu,nu,iSpin,2) * XS(nu,mu,iSpin,1) + XS(nu,mu,iSpin,1) * XS(mu,nu,iSpin,2) +& - & SXS(mu,nu,iSpin,2) * xpyas(nu,mu,iSpin,1) + SXS(nu,mu,iSpin,1) * xpyas(mu,nu,iSpin,2) +& - & xpyas(mu,nu,iSpin,2) * SXS(nu,mu,iSpin,1) + xpyas(nu,mu,iSpin,1) * SXS(mu,nu,iSpin,2) +& - & SX(mu,nu,iSpin,2) * SX(nu,mu,iSpin,1) + SX(nu,mu,iSpin,1) * SX(mu,nu,iSpin,2) ) + & ( XS(mu,nu,iSpin,2) * XS(nu,mu,iSpin,1)& + & + XS(nu,mu,iSpin,1) * XS(mu,nu,iSpin,2)& + & + SXS(mu,nu,iSpin,2) * xpyas(nu,mu,iSpin,1)& + & + SXS(nu,mu,iSpin,1) * xpyas(mu,nu,iSpin,2)& + & + xpyas(mu,nu,iSpin,2) * SXS(nu,mu,iSpin,1)& + & + xpyas(nu,mu,iSpin,1) * SXS(mu,nu,iSpin,2)& + & + SX(mu,nu,iSpin,2) * SX(nu,mu,iSpin,1)& + & + SX(nu,mu,iSpin,1) * SX(mu,nu,iSpin,2) ) tmprs = tmprs + & - & ( xmyas(mu,nu,iSpin,1) * SYS(mu,nu,iSpin,2) + xmyas(nu,mu,iSpin,2) * SYS(nu,mu,iSpin,1) +& - & SY(mu,nu,iSpin,1) * YS(mu,nu,iSpin,2) + SY(nu,mu,iSpin,2) * YS(nu,mu,iSpin,1) ) + & ( xmyas(mu,nu,iSpin,1) * SYS(mu,nu,iSpin,2)& + & + xmyas(nu,mu,iSpin,2) * SYS(nu,mu,iSpin,1)& + & + SY(mu,nu,iSpin,1) * YS(mu,nu,iSpin,2)& + & + SY(nu,mu,iSpin,2) * YS(nu,mu,iSpin,1) ) tmprs = tmprs + & - & ( xmyas(mu,nu,iSpin,2) * SYS(mu,nu,iSpin,1) + xmyas(nu,mu,iSpin,1) * SYS(nu,mu,iSpin,2) +& - & SY(mu,nu,iSpin,2) * YS(mu,nu,iSpin,1) + SY(nu,mu,iSpin,1) * YS(nu,mu,iSpin,2) ) + & ( xmyas(mu,nu,iSpin,2) * SYS(mu,nu,iSpin,1)& + & + xmyas(nu,mu,iSpin,1) * SYS(nu,mu,iSpin,2)& + & + SY(mu,nu,iSpin,2) * YS(mu,nu,iSpin,1)& + & + SY(nu,mu,iSpin,1) * YS(nu,mu,iSpin,2) ) tmprs = tmprs - 0.5_dp * & - & ( YS(mu,nu,iSpin,1) * YS(nu,mu,iSpin,2) + YS(nu,mu,iSpin,2) * YS(mu,nu,iSpin,1) +& - & SYS(mu,nu,iSpin,1) * xmyas(nu,mu,iSpin,2) + SYS(nu,mu,iSpin,2) * xmyas(mu,nu,iSpin,1) +& - & xmyas(mu,nu,iSpin,1) * SYS(nu,mu,iSpin,2) + xmyas(nu,mu,iSpin,2) * SYS(mu,nu,iSpin,1) +& - & SY(mu,nu,iSpin,1) * SY(nu,mu,iSpin,2) + SY(nu,mu,iSpin,2) * SY(mu,nu,iSpin,1) ) + & ( YS(mu,nu,iSpin,1) * YS(nu,mu,iSpin,2)& + & + YS(nu,mu,iSpin,2) * YS(mu,nu,iSpin,1)& + & + SYS(mu,nu,iSpin,1) * xmyas(nu,mu,iSpin,2)& + & + SYS(nu,mu,iSpin,2) * xmyas(mu,nu,iSpin,1)& + & + xmyas(mu,nu,iSpin,1) * SYS(nu,mu,iSpin,2)& + & + xmyas(nu,mu,iSpin,2) * SYS(mu,nu,iSpin,1)& + & + SY(mu,nu,iSpin,1) * SY(nu,mu,iSpin,2)& + & + SY(nu,mu,iSpin,2) * SY(mu,nu,iSpin,1) ) tmprs = tmprs - 0.5_dp * & - & ( YS(mu,nu,iSpin,2) * YS(nu,mu,iSpin,1) + YS(nu,mu,iSpin,1) * YS(mu,nu,iSpin,2) +& - & SYS(mu,nu,iSpin,2) * xmyas(nu,mu,iSpin,1) + SYS(nu,mu,iSpin,1) * xmyas(mu,nu,iSpin,2) +& - & xmyas(mu,nu,iSpin,2) * SYS(nu,mu,iSpin,1) + xmyas(nu,mu,iSpin,1) * SYS(mu,nu,iSpin,2) +& - & SY(mu,nu,iSpin,2) * SY(nu,mu,iSpin,1) + SY(nu,mu,iSpin,1) * SY(mu,nu,iSpin,2) ) + & ( YS(mu,nu,iSpin,2) * YS(nu,mu,iSpin,1)& + & + YS(nu,mu,iSpin,1) * YS(mu,nu,iSpin,2)& + & + SYS(mu,nu,iSpin,2) * xmyas(nu,mu,iSpin,1)& + & + SYS(nu,mu,iSpin,1) * xmyas(mu,nu,iSpin,2)& + & + xmyas(mu,nu,iSpin,2) * SYS(nu,mu,iSpin,1)& + & + xmyas(nu,mu,iSpin,1) * SYS(mu,nu,iSpin,2)& + & + SY(mu,nu,iSpin,2) * SY(nu,mu,iSpin,1)& + & + SY(nu,mu,iSpin,1) * SY(mu,nu,iSpin,2) ) end do end do end do @@ -5033,46 +5067,51 @@ subroutine addNadiaGradients(sym, nxov, natom, species0, iAtomStart, norb, homo, if (tHybridXc) then tmprs = 0.0_dp do ka = 1, nOrb - tmprs = tmprs +& - & ( PS(mu,ka,iSpin) * deltaRho(nu,ka,iSpin) + PS(nu,ka,iSpin) * deltaRho(mu,ka,iSpin) +& - & pc(mu,ka,iSpin) * DS(nu,ka,iSpin) + pc(nu,ka,iSpin) * DS(mu,ka,iSpin) ) *& - & (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) - - tmprs = tmprs + 0.5_dp * & - & ( xpyas(mu,ka,iSpin,1) * XS(nu,ka,iSpin,2) + xpyas(ka,mu,iSpin,2) * SX(ka,nu,iSpin,1) +& - & xpyas(nu,ka,iSpin,1) * XS(mu,ka,iSpin,2) + xpyas(ka,nu,iSpin,2) * SX(ka,mu,iSpin,1) )*& - & (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) - tmprs = tmprs + 0.5_dp * & - & ( xpyas(mu,ka,iSpin,2) * XS(nu,ka,iSpin,1) + xpyas(ka,mu,iSpin,1) * SX(ka,nu,iSpin,2) +& - & xpyas(nu,ka,iSpin,2) * XS(mu,ka,iSpin,1) + xpyas(ka,nu,iSpin,1) * SX(ka,mu,iSpin,2) )*& - & (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) - - tmprs = tmprs + 0.5_dp * & - & ( xmyas(mu,ka,iSpin,1) * YS(nu,ka,iSpin,2) + xmyas(ka,mu,iSpin,2) * SY(ka,nu,iSpin,1) +& - & xmyas(nu,ka,iSpin,1) * YS(mu,ka,iSpin,2) + xmyas(ka,nu,iSpin,2) * SY(ka,mu,iSpin,1) ) *& - & (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) - tmprs = tmprs + 0.5_dp * & - & ( xmyas(mu,ka,iSpin,2) * YS(nu,ka,iSpin,1) + xmyas(ka,mu,iSpin,1) * SY(ka,nu,iSpin,2) +& - & xmyas(nu,ka,iSpin,2) * YS(mu,ka,iSpin,1) + xmyas(ka,nu,iSpin,1) * SY(ka,mu,iSpin,2) ) *& - & (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) - - tmprs = tmprs + 0.5_dp * & - & ( XS(mu,ka,iSpin,1) * xpyas(ka,nu,iSpin,2) + XS(nu,ka,iSpin,1) * xpyas(ka,mu,iSpin,2) +& - & xpyas(mu,ka,iSpin,1) * SX(ka,nu,iSpin,2) + xpyas(nu,ka,iSpin,1) * SX(ka,mu,iSpin,2)) *& - & (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) - tmprs = tmprs + 0.5_dp * & - & ( XS(mu,ka,iSpin,2) * xpyas(ka,nu,iSpin,1) + XS(nu,ka,iSpin,2) * xpyas(ka,mu,iSpin,1) +& - & xpyas(mu,ka,iSpin,2) * SX(ka,nu,iSpin,1) + xpyas(nu,ka,iSpin,2) * SX(ka,mu,iSpin,1)) *& - & (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) - - tmprs = tmprs - 0.5_dp * & - & ( YS(mu,ka,iSpin,1) * xmyas(ka,nu,iSpin,2) + YS(nu,ka,iSpin,1) * xmyas(ka,mu,iSpin,2) +& - & xmyas(mu,ka,iSpin,1) * SY(ka,nu,iSpin,2) + xmyas(nu,ka,iSpin,1) * SY(ka,mu,iSpin,2)) *& - & (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) - tmprs = tmprs - 0.5_dp * & - & ( YS(mu,ka,iSpin,2) * xmyas(ka,nu,iSpin,1) + YS(nu,ka,iSpin,2) * xmyas(ka,mu,iSpin,1) +& - & xmyas(mu,ka,iSpin,2) * SY(ka,nu,iSpin,1) + xmyas(nu,ka,iSpin,2) * SY(ka,mu,iSpin,1)) *& - & (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) + tmprs = tmprs + (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) *& + & ( PS(mu,ka,iSpin) * deltaRho(nu,ka,iSpin)& + & + PS(nu,ka,iSpin) * deltaRho(mu,ka,iSpin)& + & + pc(mu,ka,iSpin) * DS(nu,ka,iSpin)& + & + pc(nu,ka,iSpin) * DS(mu,ka,iSpin) ) + tmprs = tmprs + 0.5_dp * (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) *& + & ( xpyas(mu,ka,iSpin,1) * XS(nu,ka,iSpin,2)& + & + xpyas(ka,mu,iSpin,2) * SX(ka,nu,iSpin,1)& + & + xpyas(nu,ka,iSpin,1) * XS(mu,ka,iSpin,2)& + & + xpyas(ka,nu,iSpin,2) * SX(ka,mu,iSpin,1) ) + tmprs = tmprs + 0.5_dp * (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) *& + & ( xpyas(mu,ka,iSpin,2) * XS(nu,ka,iSpin,1)& + & + xpyas(ka,mu,iSpin,1) * SX(ka,nu,iSpin,2)& + & + xpyas(nu,ka,iSpin,2) * XS(mu,ka,iSpin,1)& + & + xpyas(ka,nu,iSpin,1) * SX(ka,mu,iSpin,2) ) + tmprs = tmprs + 0.5_dp * (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) *& + & ( xmyas(mu,ka,iSpin,1) * YS(nu,ka,iSpin,2)& + & + xmyas(ka,mu,iSpin,2) * SY(ka,nu,iSpin,1)& + & + xmyas(nu,ka,iSpin,1) * YS(mu,ka,iSpin,2)& + & + xmyas(ka,nu,iSpin,2) * SY(ka,mu,iSpin,1) ) + tmprs = tmprs + 0.5_dp * (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) *& + & ( xmyas(mu,ka,iSpin,2) * YS(nu,ka,iSpin,1)& + & + xmyas(ka,mu,iSpin,1) * SY(ka,nu,iSpin,2)& + & + xmyas(nu,ka,iSpin,2) * YS(mu,ka,iSpin,1)& + & + xmyas(ka,nu,iSpin,1) * SY(ka,mu,iSpin,2) ) + tmprs = tmprs + 0.5_dp * (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) *& + & ( XS(mu,ka,iSpin,1) * xpyas(ka,nu,iSpin,2)& + & + XS(nu,ka,iSpin,1) * xpyas(ka,mu,iSpin,2)& + & + xpyas(mu,ka,iSpin,1) * SX(ka,nu,iSpin,2)& + & + xpyas(nu,ka,iSpin,1) * SX(ka,mu,iSpin,2)) + tmprs = tmprs + 0.5_dp * (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) *& + & ( XS(mu,ka,iSpin,2) * xpyas(ka,nu,iSpin,1)& + & + XS(nu,ka,iSpin,2) * xpyas(ka,mu,iSpin,1)& + & + xpyas(mu,ka,iSpin,2) * SX(ka,nu,iSpin,1)& + & + xpyas(nu,ka,iSpin,2) * SX(ka,mu,iSpin,1)) + tmprs = tmprs - 0.5_dp * (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) *& + & ( YS(mu,ka,iSpin,1) * xmyas(ka,nu,iSpin,2)& + & + YS(nu,ka,iSpin,1) * xmyas(ka,mu,iSpin,2)& + & + xmyas(mu,ka,iSpin,1) * SY(ka,nu,iSpin,2)& + & + xmyas(nu,ka,iSpin,1) * SY(ka,mu,iSpin,2)) + tmprs = tmprs - 0.5_dp * (lrGammaOrb(mu,ka) + lrGammaOrb(nu,ka)) *& + & ( YS(mu,ka,iSpin,2) * xmyas(ka,nu,iSpin,1)& + & + YS(nu,ka,iSpin,2) * xmyas(ka,mu,iSpin,1)& + & + xmyas(mu,ka,iSpin,2) * SY(ka,nu,iSpin,1)& + & + xmyas(nu,ka,iSpin,2) * SY(ka,mu,iSpin,1)) end do ! Factor of 2 for spin-polarized calculations tmprs2 = tmprs2 + cExchange * nSpin * dSo(n,m,xyz) * tmprs @@ -5476,8 +5515,10 @@ end subroutine fixNACVPhase !> Implements the CI optimizer of Bearpark et al. Chem. Phys. Lett. 223 269 (1994) with !> modifications introduced by Harabuchi/Hatanaka/Maeda CPL X 2019 !> Previous published results [Niehaus JCP 158 054103 (2023), TCA 140 34 (2021)] were obtained - !> with a differing version that assumed orthogonal X1 and X2 vectors, which leads to poor convergence - subroutine conicalIntersectionOptimizer(derivs, excDerivs, indNACouplings, energyShift, naCouplings, excEnergies) + !> with a differing version that assumed orthogonal X1 and X2 vectors, which leads to poor + !> convergence + subroutine conicalIntersectionOptimizer(derivs, excDerivs, indNACouplings, energyShift,& + & naCouplings, excEnergies) !> Ground state gradient (overwritten) real(dp), intent(inout) :: derivs(:,:) diff --git a/src/dftbp/timedep/transcharges.F90 b/src/dftbp/timedep/transcharges.F90 index 13dc308997..fab11a9962 100644 --- a/src/dftbp/timedep/transcharges.F90 +++ b/src/dftbp/timedep/transcharges.F90 @@ -481,7 +481,8 @@ end subroutine TTransCharges_qMatVecDs !> Transition charges right produced with a vector v * Q for spin up !> and negative transition charges right produced with a vector v * Q for spin down - !> R_ias = delta_s sum_A q_A^(ias) V_A, where delta_s = 1 for spin up and delta_s = -1 for spin down + !> R_ias = delta_s sum_A q_A^(ias) V_A, where delta_s = 1 for spin up and delta_s = -1 for spin + !> down pure subroutine TTransCharges_qVecMatDs(this, iAtomStart, sTimesGrndEigVecs, grndEigVecs, getia,& & win, vector, qProduct)