diff --git a/source/tdhf_mrsf_lib.F90 b/source/tdhf_mrsf_lib.F90 index 84b977c..e17c94e 100644 --- a/source/tdhf_mrsf_lib.F90 +++ b/source/tdhf_mrsf_lib.F90 @@ -335,6 +335,52 @@ subroutine mrinivec(infos,ea,eb,bvec_mo,xm,nvec) end subroutine mrinivec +!> Transform MRSF response vectors from MO to AO basis +!> +!> This subroutine performs the transformation of MRSF-TDDFT response amplitudes +!> X^(k), where k is singlet or triplet, from molecular orbital (MO) representation +!> to atomic orbital (AO) basis. The AO-basis response matrices are needed +!> for contraction with two-electron integrals and Fock matrix contributions +!> in the response equations. +!> +!> Physical context: +!> In MRSF-TDDFT, the response space is constructed from MS=+/-1 triplet references +!> to eliminate spin contamination in target singlet and triplet excited states. +!> +!> Orbital spaces in MRSF-TDDFT: +!> - C (Closed): Doubly-occupied orbitals (indices i,j,k,l) +!> - O (Open): Singly-occupied orbitals O1=HOMO-1, O2=HOMO (indices u,v,w,z) +!> - V (Virtual): Unoccupied orbitals (indices a,b,c,d) +!> +!> For each reference state k, the response amplitudes X^(k)_pq represent orbital +!> excitations between different orbital spaces. The response configurations are: +!> - Type I (OO): Open-to-Open transitions (O1<->O2) +!> - Type II (CO): Closed-to-Open transitions (C->O1, C->O2) +!> - Type III (OV): Open-to-Virtual transitions (O1->V, O2->V) +!> - Type IV (CV): Closed-to-Virtual transitions (C->V) +!> +!> The six response components in this subroutine correspond to: +!> 1. bo2v: O2(HOMO, alpha) -> V(beta) - OV block +!> 2. bo1v: O1(HOMO-1, alpha) -> V(beta) - OV block +!> 3. bco1: C(alpha) -> O1(HOMO-1, beta) - CO block +!> 4. bco2: C(alpha) -> O2(HOMO, beta) - CO block +!> 5. o21v: Mixed OV component coupling O1 and O2 with V +!> 6. co12: Mixed CO component coupling C with O1 and O2 +!> +!> Transformation scheme: +!> For each component, we transform X^(k)_pq (MO basis) to P^(k)_(mu,nu) (AO basis): +!> P^(k)_(mu,nu) = sum_pq C_(mu,p) X^(k)_pq C_(nu,q) +!> where C are MO coefficient matrices (va for alpha-spin, vb for beta-spin). +!> +!> Spin-pairing coupling between MS=+1 and MS=-1 reference states is realized +!> through specific linear combinations of MO coefficients va and vb (with proper +!> signs and 1/sqrt(2) normalization), following Slater-Condon rules. This enables +!> proper description of singlet and triplet target states from the mixed-reference formalism. +!> +!> Reference: Lee et al., J. Chem. Phys. 150, 184111 (2019), Eq. 2.11-2.18 +!> +!> \author Konstantin Komarov (constlike@gmail.com) +!> subroutine mrsfcbc(infos,va,vb,bvec,fmrsf) use messages, only: show_message, with_abort @@ -351,7 +397,7 @@ subroutine mrsfcbc(infos,va,vb,bvec,fmrsf) fmrsf real(kind=dp), allocatable, dimension(:,:) :: & - tmp, tmp1, tmp2 + tmp real(kind=dp), pointer, dimension(:,:) :: & bo2v, bo1v, bco1, bco2, ball, co12, o21v integer :: nocca, noccb, mrst, i, j, m, nbf, lr1, lr2, ok @@ -375,53 +421,252 @@ subroutine mrsfcbc(infos,va,vb,bvec,fmrsf) lr1 = nocca-1 lr2 = nocca - allocate(tmp(nbf,nbf), & - tmp1(nbf,4), & - tmp2(nbf,4), & - source=0.0_dp, stat=ok) + allocate(tmp(nbf,max(1,noccb)), source=0.0_dp, stat=ok) if( ok/=0 ) call show_message('Cannot allocate memory',with_abort) - do j = nocca+1, nbf - tmp1(:,1) = tmp1(:,1)+vb(:,j)*bvec(lr2,j) - tmp1(:,2) = tmp1(:,2)+vb(:,j)*bvec(lr1,j) - tmp1(:,3) = tmp1(:,3)+vb(:,j)*bvec(lr2,j) - tmp1(:,4) = tmp1(:,4)+vb(:,j)*bvec(lr1,j) - end do - - do i = 1, noccb - tmp2(:,1) = tmp2(:,1)+va(:,i)*bvec(i,lr1) - tmp2(:,2) = tmp2(:,2)+va(:,i)*bvec(i,lr2) - tmp2(:,3) = tmp2(:,3)+va(:,i)*bvec(i,lr1) - tmp2(:,4) = tmp2(:,4)+va(:,i)*bvec(i,lr2) - end do - - do m = 1, nbf - bo2v(:,m) = bo2v(:,m)+va(:,lr2)*tmp1(m,1) - bo1v(:,m) = bo1v(:,m)+va(:,lr1)*tmp1(m,2) - bco1(:,m) = bco1(:,m)+vb(m,lr1)*tmp2(:,1) - bco2(:,m) = bco2(:,m)+vb(m,lr2)*tmp2(:,2) - end do - - do m = 1, nbf - o21v(m,:) = o21v(m,:)+va(:,lr1)*tmp1(m,3) & - -va(:,lr2)*tmp1(m,4) - co12(m,:) = co12(m,:)+vb(m,lr2)*tmp2(:,3) & - -vb(m,lr1)*tmp2(:,4) - end do - + !----------------------------------------------------------------------- + ! Component 1: bo2v - O2(HOMO, alpha) -> V(beta) excitations (OV block) + !----------------------------------------------------------------------- + ! Physical meaning: This component represents Open-to-Virtual (OV) excitations + ! from the HOMO orbital of alpha-spin (O2, lr2 = nocca) to virtual orbitals + ! of beta-spin. In MRSF theory, this corresponds to Type III response + ! configurations. + ! + ! MO->AO transformation: P^bo2v_(mu,nu) = C^alpha_(mu,HOMO) * X_(HOMO,a) * C^beta_(nu,a) + ! where a runs over virtual beta-orbitals (nocca+1:nbf) + ! + ! Step 1: Intermediate vector tmp = sum_a C^beta_(mu,a) X_(HOMO,a) + ! tmp_mu = sum_{a in virt_beta} C^beta_(mu,a) * X_(HOMO,a) + call dgemm('n', 't', nbf, 1, nbf-nocca, & + 1.0_dp, vb(:,nocca+1), nbf, & + bvec(lr2:lr2,nocca+1), nbf, & + 0.0_dp, tmp(:,1), nbf) + + ! Step 2: Outer product to form AO-basis matrix + ! P^bo2v_(mu,nu) += C^alpha_(mu,HOMO) * tmp_nu + call dgemm('n', 't', nbf, nbf, 1, & + 1.0_dp, va(:,lr2:lr2), nbf, & + tmp(:,1:1), nbf, & + 1.0_dp, bo2v, nbf) + + !----------------------------------------------------------------------- + ! Component 2: bo1v - O1(HOMO-1, alpha) -> V(beta) excitations (OV block) + !----------------------------------------------------------------------- + ! Physical meaning: This component represents Open-to-Virtual (OV) excitations + ! from HOMO-1 orbital of alpha-spin (O1, lr1 = nocca-1) to virtual orbitals + ! of beta-spin. Together with bo2v, this forms the complete set of Type III + ! spin-flip excitations from the two singly-occupied MOs (O1 and O2) in the + ! MS=+/-1 triplet reference. + ! + ! MO->AO transformation: P^bo1v_(mu,nu) = C^alpha_(mu,HOMO-1) * X_(HOMO-1,a) * C^beta_(nu,a) + ! + ! Step 1: Intermediate vector tmp = sum_a C^beta_(mu,a) X_(HOMO-1,a) + ! tmp_mu = sum_{a in virt_beta} C^beta_(mu,a) * X_(HOMO-1,a) + call dgemm('n', 't', nbf, 1, nbf-nocca, & + 1.0_dp, vb(:,nocca+1), nbf, & + bvec(lr1:lr1,nocca+1), nbf, & + 0.0_dp, tmp(:,1), nbf) + + ! Step 2: Outer product to form AO-basis matrix + ! P^bo1v_(mu,nu) += C^alpha_(mu,HOMO-1) * tmp_nu + call dgemm('n', 't', nbf, nbf, 1, & + 1.0_dp, va(:,lr1:lr1), nbf, & + tmp(:,1:1), nbf, & + 1.0_dp, bo1v, nbf) + + !----------------------------------------------------------------------- + ! Component 3: bco1 - C(alpha) -> O1(HOMO-1, beta) excitations (CO block) + !----------------------------------------------------------------------- + ! Physical meaning: This component represents Closed-to-Open (CO) excitations + ! from doubly-occupied alpha-orbitals (C, 1:noccb) to the HOMO-1 orbital of + ! beta-spin (O1, lr1). In MRSF theory, this corresponds to Type II response + ! configurations. These are spin-flip de-excitations that complement the + ! bo1v/bo2v excitations, maintaining the symmetry of the response space. + ! + ! MO->AO transformation: P^bco1_(mu,nu) = C^alpha_(mu,i) * X_(i,HOMO-1) * C^beta_(nu,HOMO-1) + ! where i runs over doubly-occupied orbitals (1:noccb) + ! + ! Step 1: Intermediate vector tmp = sum_i C^alpha_(mu,i) X_(i,HOMO-1) + ! tmp_mu = sum_{i in occ_alpha} C^alpha_(mu,i) * X_(i,HOMO-1) + call dgemm('n', 'n', nbf, 1, noccb, & + 1.0_dp, va, nbf, & + bvec(1:noccb,lr1:lr1), nbf, & + 0.0_dp, tmp(:,1), nbf) + + ! Step 2: Outer product to form AO-basis matrix + ! P^bco1_(mu,nu) += tmp_mu * C^beta_(nu,HOMO-1) + call dgemm('n', 't', nbf, nbf, 1, & + 1.0_dp, tmp(:,1:1), nbf, & + vb(:,lr1:lr1), nbf, & + 1.0_dp, bco1, nbf) + + !----------------------------------------------------------------------- + ! Component 4: bco2 - C(alpha) -> O2(HOMO, beta) excitations (CO block) + !----------------------------------------------------------------------- + ! Physical meaning: This component represents Closed-to-Open (CO) excitations + ! from doubly-occupied alpha-orbitals (C, 1:noccb) to the HOMO orbital of + ! beta-spin (O2, lr2). In MRSF theory, this corresponds to Type II response + ! configurations. Together with bco1, this completes the set of spin-flip + ! de-excitations from doubly-occupied orbitals to the two singly-occupied MOs. + ! + ! MO->AO transformation: P^bco2_(mu,nu) = C^alpha_(mu,i) * X_(i,HOMO) * C^beta_(nu,HOMO) + ! where i runs over doubly-occupied orbitals (1:noccb) + ! + ! Step 1: Intermediate vector tmp = sum_i C^alpha_(mu,i) X_(i,HOMO) + ! tmp_mu = sum_{i in occ_alpha} C^alpha_(mu,i) * X_(i,HOMO) + call dgemm('n', 'n', nbf, 1, noccb, & + 1.0_dp, va, nbf, & + bvec(1:noccb,lr2:lr2), nbf, & + 0.0_dp, tmp(:,1), nbf) + + ! Step 2: Outer product to form AO-basis matrix + ! P^bco2_(mu,nu) += tmp_mu * C^beta_(nu,HOMO) + call dgemm('n', 't', nbf, nbf, 1, & + 1.0_dp, tmp(:,1:1), nbf, & + vb(:,lr2:lr2), nbf, & + 1.0_dp, bco2, nbf) + + !----------------------------------------------------------------------- + ! Component 5: o21v - Mixed (O1<->O2)(alpha) x V(beta) (OV block) + !----------------------------------------------------------------------- + ! Physical meaning: This is a mixed Open-to-Virtual (OV) component coupling + ! both HOMO and HOMO-1 alpha-orbitals (O1 and O2) with virtual beta-orbitals. + ! The subtraction ensures proper antisymmetry and represents coherent + ! superpositions of spin-flip excitations. This component is essential for + ! the correct description of spin-adapted states in MRSF theory, arising + ! from the coupling between different OV response configurations. + ! + ! MO->AO transformation: + ! P^o21v_(mu,nu) = sum_a [ + ! C^alpha_(mu,HOMO-1) * X_(HOMO,a) - C^alpha_(mu,HOMO) * X_(HOMO-1,a) + ! ] * C^beta_(nu,a) + ! + ! Step 1: Intermediate vector from HOMO -> virt_beta amplitudes + ! tmp_mu = sum_{a in virt_beta} C^beta_(mu,a) * X_(HOMO,a) + call dgemm('n', 't', nbf, 1, nbf-nocca, & + 1.0_dp, vb(:,nocca+1), nbf, & + bvec(lr2:lr2,nocca+1), nbf, & + 0.0_dp, tmp(:,1), nbf) + + ! Combined outer products with subtraction + ! P^o21v_(mu,nu) += C^alpha_(mu,HOMO-1) * tmp_nu (positive contribution) + call dgemm('n', 't', nbf, nbf, 1, & + 1.0_dp, tmp(:,1:1), nbf, & + va(:,lr1:lr1), nbf, & + 1.0_dp, o21v, nbf) + + ! Step 1: Intermediate vector from HOMO-1 -> virt_beta amplitudes + ! tmp_mu = sum_{a in virt_beta} C^beta_(mu,a) * X_(HOMO-1,a) + call dgemm('n', 't', nbf, 1, nbf-nocca, & + 1.0_dp, vb(:,nocca+1), nbf, & + bvec(lr1:lr1,nocca+1), nbf, & + 0.0_dp, tmp(:,1), nbf) + + ! P^o21v_(mu,nu) -= C^alpha_(mu,HOMO) * tmp_nu (negative contribution) + call dgemm('n', 't', nbf, nbf, 1, & + -1.0_dp, tmp(:,1:1), nbf, & + va(:,lr2:lr2), nbf, & + 1.0_dp, o21v, nbf) + + !----------------------------------------------------------------------- + ! Component 6: co12 - C(alpha) x Mixed (O1<->O2)(beta) (CO block) + !----------------------------------------------------------------------- + ! Physical meaning: This is a mixed Closed-to-Open (CO) component coupling + ! doubly-occupied alpha-orbitals with both HOMO and HOMO-1 beta-orbitals + ! (O1 and O2). The subtraction ensures proper antisymmetry and represents + ! coherent superpositions of spin-flip de-excitations. Together with o21v, + ! this maintains the full symmetry of the MRSF response space under orbital + ! permutations, arising from coupling between different CO configurations. + ! + ! MO->AO transformation: + ! P^co12_(mu,nu) = sum_i [ + ! C^beta_(mu,HOMO) * X_(i,HOMO-1) - C^beta_(mu,HOMO-1) * X_(i,HOMO) + ! ] * C^alpha_(nu,i) + ! + ! Step 1: Intermediate vector from occ_alpha -> HOMO-1_beta amplitudes + ! tmp_mu = sum_{i in occ_alpha} C^alpha_(mu,i) * X_(i,HOMO-1) + call dgemm('n', 'n', nbf, 1, noccb, & + 1.0_dp, va, nbf, & + bvec(1:noccb,lr1:lr1), nbf, & + 0.0_dp, tmp(:,1), nbf) + + ! Combined outer products with subtraction + ! P^co12_(mu,nu) += C^beta_(mu,HOMO) * tmp_nu (positive contribution) + call dgemm('n', 't', nbf, nbf, 1, & + 1.0_dp, vb(:,lr2:lr2), nbf, & + tmp(:,1:1), nbf, & + 1.0_dp, co12, nbf) + + ! Step 1: Intermediate vector from occ_alpha -> HOMO_beta amplitudes + ! tmp_mu = sum_{i in occ_alpha} C^alpha_(mu,i) * X_(i,HOMO) + call dgemm('n', 'n', nbf, 1, noccb, & + 1.0_dp, va, nbf, & + bvec(1:noccb,lr2:lr2), nbf, & + 0.0_dp, tmp(:,1), nbf) + + ! P^co12_(mu,nu) -= C^beta_(mu,HOMO-1) * tmp_nu (negative contribution) + call dgemm('n', 't', nbf, nbf, 1, & + -1.0_dp, vb(:,lr1:lr1), nbf, & + tmp(:,1:1), nbf, & + 1.0_dp, co12, nbf) + + !----------------------------------------------------------------------- + ! Sum the four primary components into the total response matrix + !----------------------------------------------------------------------- + ! Physical meaning: Combine bo2v, bo1v, bco1, and bco2 (the four components + ! without mixed character) into the total AO-basis response matrix ball. + ! These represent the OV and CO blocks (Type II and Type III configurations). + ! The mixed components o21v and co12 are not included here as they are + ! handled separately in the spin-dependent sections below. ball = ball + bo2v + bo1v + bco1 + bco2 + !----------------------------------------------------------------------- + ! Additional general contribution: C(alpha) x V(beta) block (CV) + !----------------------------------------------------------------------- + ! Physical meaning: Transform the general Closed-to-Virtual (CV) block of + ! response amplitudes (doubly-occupied_alpha -> virtual_beta) to AO basis. + ! This represents Type IV spin-flip excitations from the doubly-occupied + ! core orbitals. These excitations do not involve the singly-occupied + ! frontier orbitals (O1, O2) and represent the CV response configurations. + ! + ! Transformation: P^ball_(mu,nu) += sum_ia C^alpha_(mu,i) * X_(i,a) * C^beta_(nu,a) + ! where i in doubly-occupied (1:noccb), a in virtual_beta (nocca+1:nbf) + ! + ! Step 1: Intermediate tmp_(mu,i) = sum_a C^beta_(mu,a) * X_(i,a) call dgemm('n', 't', nbf, noccb, nbf-nocca, & 1.0_dp, vb(:,nocca+1), nbf, & bvec(:,nocca+1), nbf, & - 0.0_dp, tmp, nbf) + 0.0_dp, tmp(:,1:noccb), nbf) + ! Step 2: Outer product P^ball_(mu,nu) += sum_i C^alpha_(mu,i) * tmp_(nu,i) call dgemm('n', 't', nbf, nbf, noccb, & 1.0_dp, va, nbf, & - tmp, nbf, & + tmp(:,1:noccb), nbf, & 1.0_dp, ball, nbf) + !----------------------------------------------------------------------- + ! Spin-dependent corrections for (O1<->O2) x (O1<->O2) block (OO) + !----------------------------------------------------------------------- + ! Physical meaning: Add corrections for the Open-to-Open (OO) Type I response + ! configurations that depend on the target spin state (singlet mrst=1 or + ! triplet mrst=3). These involve the special elements X_(HOMO,HOMO-1), + ! X_(HOMO-1,HOMO), and X_(HOMO-1,HOMO-1) that couple the two singly-occupied + ! MOs (O1 and O2). The 1/sqrt(2) factor ensures proper normalization of + ! spin-adapted states in the MRSF formalism. + ! + ! Spin-pairing coupling between MS=+1 and MS=-1 triplet references (following + ! Slater-Condon rules) is implemented via specific linear combinations of + ! MO coefficients va and vb. The different sign patterns for singlet/triplet + ! reflect the different spin symmetries: + ! - Singlet: antisymmetric spatial wavefunction (subtraction) + ! - Triplet: symmetric spatial wavefunction (addition) if (mrst==1) then + ! Singlet state corrections (mrst=1): + ! Three terms that couple HOMO and HOMO-1 orbitals: + ! 1. X_(HOMO,HOMO-1) * C^alpha_HOMO (outer) C^beta_HOMO-1 + ! 2. X_(HOMO-1,HOMO) * C^alpha_HOMO-1 (outer) C^beta_HOMO + ! 3. X_(HOMO-1,HOMO-1) * (C^alpha_HOMO-1 (outer) C^beta_HOMO-1 - C^alpha_HOMO (outer) C^beta_HOMO) / sqrt(2) + ! The subtraction in term 3 ensures proper singlet spin coupling. do m = 1, nbf ball(:,m) = ball(:,m) & +va(:,lr2)*bvec(lr2,lr1)*vb(m,lr1) & @@ -430,6 +675,11 @@ subroutine mrsfcbc(infos,va,vb,bvec,fmrsf) *bvec(lr1,lr1)*isqrt2 end do else if (mrst==3) then + ! Triplet state corrections (mrst=3): + ! Single term for diagonal HOMO-1,HOMO-1 element: + ! X_(HOMO-1,HOMO-1) * (C^alpha_HOMO-1 (outer) C^beta_HOMO-1 + C^alpha_HOMO (outer) C^beta_HOMO) / sqrt(2) + ! The addition ensures proper triplet spin coupling. + ! Off-diagonal OO terms are zero for triplet states. do m = 1, nbf ball(:,m) = ball(:,m) & +(va(:,lr1)*vb(m,lr1)+va(:,lr2)*vb(m,lr2)) & @@ -453,120 +703,38 @@ subroutine mrsfcbc(infos,va,vb,bvec,fmrsf) return end subroutine mrsfcbc -! subroutine umrsfcbc(infos, va, vb, bvec, mrsf_density) - -! use messages, only: show_message, with_abort -! use types, only: information -! use io_constants, only: iw -! use precision, only: dp - -! implicit none - -! type(information), intent(in) :: infos -! real(kind=dp), intent(in), dimension(:,:) :: & -! va, vb, bvec -! real(kind=dp), intent(inout), target, dimension(:,:,:) :: & -! mrsf_density - -! real(kind=dp), allocatable, dimension(:,:) :: & -! tmp, tmp1, tmp2 -! real(kind=dp), pointer, dimension(:,:) :: & -! bo2va, bo2vb, bo1va, bo1vb, bco1a, bco1b, & -! bco2a, bco2b, ball, co12, o21v -! integer :: nocca, noccb, mrst, i, j, m, nbf, lr1, lr2, ok -! logical :: debug_mode -! real(kind=dp), parameter :: isqrt2 = 1.0_dp/sqrt(2.0_dp) - -! ball => mrsf_density(11,:,:) -! bo2va => mrsf_density(1,:,:) -! bo2vb => mrsf_density(2,:,:) -! bo1va => mrsf_density(3,:,:) -! bo1vb => mrsf_density(4,:,:) -! bco1a => mrsf_density(5,:,:) -! bco1b => mrsf_density(6,:,:) -! bco2a => mrsf_density(7,:,:) -! bco2b => mrsf_density(8,:,:) -! o21v => mrsf_density(9,:,:) -! co12 => mrsf_density(10,:,:) - -! nbf = infos%basis%nbf -! nocca = infos%mol_prop%nelec_A -! noccb = infos%mol_prop%nelec_B -! mrst = infos%tddft%mult -! debug_mode = infos%tddft%debug_mode - -! lr1 = nocca-1 -! lr2 = nocca -! allocate(tmp(nbf,nbf), & -! tmp1(nbf,4), & -! tmp2(nbf,4), & -! source=0.0_dp, stat=ok) -! if( ok/=0 ) call show_message('Cannot allocate memory',with_abort) - -! do j = nocca+1, nbf -! tmp1(:,1) = tmp1(:,1)+va(:,j)*bvec(nocca,j) -! tmp1(:,2) = tmp1(:,2)+vb(:,j)*bvec(nocca,j) -! tmp1(:,3) = tmp1(:,3)+va(:,j)*bvec(nocca-1,j) -! tmp1(:,4) = tmp1(:,4)+vb(:,j)*bvec(nocca-1,j) -! end do - -! do i = 1, nocca-2 -! tmp2(:,1) = tmp2(:,1)+va(:,i)*bvec(i,nocca-1) -! tmp2(:,2) = tmp2(:,2)+vb(:,i)*bvec(i,nocca-1) -! tmp2(:,3) = tmp2(:,3)+va(:,i)*bvec(i,nocca) -! tmp2(:,4) = tmp2(:,4)+vb(:,i)*bvec(i,nocca) -! end do - -! do m = 1, nbf -! ball(:,m) = ball(:,m)+tmp1(m,2)*va(:,nocca) & -! +tmp1(m,4)*va(:,nocca-1) & -! +tmp2(:,1)*vb(m,nocca-1) & -! +tmp2(:,3)*vb(m,nocca) - -! bo2va(:,m) = bo2va(:,m)+tmp1(m,1)*va(:,nocca) -! bo2vb(:,m) = bo2vb(:,m)+tmp1(m,2)*vb(:,nocca) -! bo1va(:,m) = bo1va(:,m)+tmp1(m,3)*va(:,nocca-1) -! bo1vb(:,m) = bo1vb(:,m)+tmp1(m,4)*vb(:,nocca-1) -! o21v(:,m) = o21v(:,m)+tmp1(m,2)*va(:,nocca-1) & -! -tmp1(m,4)*va(:,nocca) - -! bco1a(:,m) = bco1a(:,m)+tmp2(:,1)*va(m,nocca-1) -! bco1b(:,m) = bco1b(:,m)+tmp2(:,2)*vb(m,nocca-1) -! bco2a(:,m) = bco2a(:,m)+tmp2(:,3)*va(m,nocca) -! bco2b(:,m) = bco2b(:,m)+tmp2(:,4)*vb(m,nocca) -! co12(:,m) = co12(:,m)+tmp2(:,1)*vb(m,nocca) & -! -tmp2(:,3)*vb(m,nocca-1) -! end do - -! call dgemm('n','t', nbf, nocca-2, nbf-nocca, & -! 1.0_dp, vb(:,nocca+1), nbf, & -! bvec(:,nocca+1), nbf, & -! 0.0_dp, tmp, nbf) - -! call dgemm('n','t', nbf, nbf, nocca-2, & -! 1.0_dp, va, nbf, & -! tmp, nbf, & -! 1.0_dp, ball, nbf) - -! if (mrst==1) then -! do m = 1, nbf -! ball(:,m) = ball(:,m) & -! +va(:,noca)*bvec(noca,noca-1)*vb(m,noca-1) & -! +va(:,noca-1)*bvec(noca-1,noca)*vb(m,noca) & -! +(va(:,noca-1)*vb(m,noca-1)-va(:,noca)*vb(m,noca)) & -! *bvec(noca-1,noca-1)*isqrt2 -! end do -! elseif (mrsft) then -! do m = 1, nbf -! ball(:,m) = ball(:,m) -! +(va(:,noca-1)*vb(m,noca-1)+va(:,noca)*vb(m,noca)) -! *bvec(noca-1,noca-1)*sqrt2 -! end do -! endif - -! return -! end subroutine umrsfcbc +!> Transform MRSF Fock-like matrices from AO to MO basis +!> +!> This subroutine performs the transformation of MRSF-TDDFT Fock-like matrices +!> (or generalized density contributions) from atomic orbital (AO) basis back to +!> molecular orbital (MO) representation. +!> +!> Physical context: +!> After contracting the response amplitudes (in AO basis) with two-electron +!> integrals and other operators, we obtain Fock-like matrices P^(k)_(mu,nu) in +!> AO basis. These must be transformed back to MO basis to extract elements +!> corresponding to specific orbital transitions in the MRSF response space. +!> +!> Orbital spaces (same as in mrsfcbc): +!> - C (Closed): Doubly-occupied orbitals +!> - O (Open): Singly-occupied O1=HOMO-1, O2=HOMO +!> - V (Virtual): Unoccupied orbitals +!> +!> Transformation scheme: +!> For the general contribution: F^MO_pq = sum_(mu,nu) C^alpha_(mu,p) P^AO_(mu,nu) C^beta_(nu,q) +!> +!> Then, specific corrections are added for each of the six response components +!> corresponding to different blocks of the MRSF response space: +!> 1. Section 3: Corrections from ado1v (OV) and aco12 (CO mixed) -> C x O2 block +!> 2. Section 4: Corrections from ado2v (OV) and aco12 (CO mixed) -> C x O1 block +!> 3. Section 5: Corrections from adco2 (CO) and ao21v (OV mixed) -> O1 x V block +!> 4. Section 6: Corrections from adco1 (CO) and ao21v (OV mixed) -> O2 x V block +!> +!> Reference: Lee et al., J. Chem. Phys. 150, 184111 (2019), Eq. 2.11-2.18 +!> +!> \author Konstantin Komarov (constlike@gmail.com) +!> subroutine mrsfmntoia(infos, fmrsf, pmo, va, vb, ivec) use precision, only: dp @@ -614,6 +782,19 @@ subroutine mrsfmntoia(infos, fmrsf, pmo, va, vb, ivec) lr1 = noca-1 lr2 = noca + !----------------------------------------------------------------------- + ! Initial AO->MO transformation of general Fock-like contribution + !----------------------------------------------------------------------- + ! Physical meaning: Transform the general (summed) Fock-like matrix agdlr + ! from AO basis to MO basis. This represents the baseline contribution + ! before adding specific corrections for each response component. + ! + ! Transformation: F^MO_pq = C^alpha^T * P^AO * C^beta + ! Step 1: wrk = C^alpha^T * agdlr (transform first index) + ! Step 2: scr = wrk * C^beta (transform second index) + ! + ! Result: scr contains the general MO-basis matrix that will be corrected + ! by the specific response component contributions in sections 3-6. call dgemm('t','n',nbf,nbf,nbf, & one, va, nbf, & agdlr,nbf, & @@ -627,35 +808,166 @@ subroutine mrsfmntoia(infos, fmrsf, pmo, va, vb, ivec) ! ----- (m,n) to (i+,n) ----- wrk = scr -! 3 - call dgemv('n',nbf,nbf,one,ado1v,nbf,vb(:,lr2),1,zero,tmp,1) - call dgemv('n',nbf,nbf,one,aco12,nbf,vb(:,lr1),1,one,tmp,1) - do i = 1, noca-2 - wrk(i,lr2) = wrk(i,lr2) + dot_product(va(:,i),tmp(:)) - end do -! 4 - call dgemv('n',nbf,nbf,one,aco12,nbf,vb(:,lr2),1,zero,tmp,1) - call dgemv('n',nbf,nbf,one,ado2v,nbf,vb(:,lr1),1,-one,tmp,1) - do i = 1, noca-2 - wrk(i,lr1) = wrk(i,lr1) + dot_product(va(:,i),tmp(:)) - end do -! 5 - call dgemv('t',nbf,nbf,one,adco2,nbf,va(:,lr1),1,zero,tmp,1) - call dgemv('t',nbf,nbf,one,ao21v,nbf,va(:,lr2),1, one,tmp,1) - do j = noca+1, nbf - wrk(lr1,j) = wrk(lr1,j) + dot_product(vb(:,j),tmp(:)) - end do -! 6 - call dgemv('t',nbf,nbf,one,ao21v,nbf,va(:,lr1),1,zero,tmp,1) - call dgemv('t',nbf,nbf,one,adco1,nbf,va(:,lr2),1,-one,tmp,1) - do j = noca+1, nbf - wrk(lr2,j) = wrk(lr2,j) + dot_product(vb(:,j),tmp(:)) - end do - + !----------------------------------------------------------------------- + ! Section 3: Corrections for C(alpha) -> O2(HOMO, beta) response element + !----------------------------------------------------------------------- + ! Physical meaning: Add corrections to F^MO_(i,HOMO) from two sources: + ! 1. ado1v: O1(HOMO-1, alpha) -> V(beta) excitation component (OV block) + ! 2. aco12: C(alpha) x Mixed (O1<->O2)(beta) component (CO mixed block) + ! + ! These represent the backflow of MO to AO, accounting for coupling + ! between different excitation types (OV and CO blocks). + ! The result corrects wrk(i,lr2) for doubly-occupied + ! orbitals i in the Closed space. + ! + ! Combined transformation: + ! tmp = ado1v * C^beta_HOMO + aco12 * C^beta_HOMO-1 + ! F^MO_(i,HOMO) += C^alpha_i^T * tmp (for i=1:noca-2) + ! + ! Step 1: Contract ado1v with HOMO_beta MO coefficient + ! tmp_mu = sum_nu P^ado1v_(mu,nu) * C^beta_(nu,HOMO) + call dgemm('n','n',nbf,1,nbf, & + one, ado1v, nbf, & + vb(:,lr2:lr2), nbf, & + zero, tmp, nbf) + ! Step 2: Add contribution from aco12 with HOMO-1_beta MO coefficient + ! tmp_mu += sum_nu P^aco12_(mu,nu) * C^beta_(nu,HOMO-1) + call dgemm('n','n',nbf,1,nbf, & + one, aco12, nbf, & + vb(:,lr1:lr1), nbf, & + one, tmp, nbf) + ! Step 3: Project onto doubly-occupied alpha-orbitals + ! F^MO_(i,HOMO) += sum_mu C^alpha_(mu,i) * tmp_mu (i=1:noca-2) + call dgemm('t','n',noca-2,1,nbf, & + one, va, nbf, & + tmp, nbf, & + one, wrk(1:noca-2,lr2:lr2), noca-2) + !----------------------------------------------------------------------- + ! Section 4: Corrections for C(alpha) -> O1(HOMO-1, beta) response element + !----------------------------------------------------------------------- + ! Physical meaning: Add corrections to F^MO_(i,HOMO-1) from two sources: + ! 1. ado2v: O2(HOMO, alpha) -> V(beta) excitation component (OV block, positive) + ! 2. aco12: C(alpha) x Mixed (O1<->O2)(beta) component (CO mixed, negative) + ! + ! The subtraction in this section ensures proper antisymmetry between + ! HOMO and HOMO-1 contributions, maintaining consistency with the mixed + ! character of the aco12 and ao21v components. This antisymmetry arises from + ! the spin-pairing coupling realized via linear combinations of va and vb. + ! + ! Combined transformation: + ! tmp = ado2v * C^beta_HOMO-1 - aco12 * C^beta_HOMO + ! F^MO_(i,HOMO-1) += C^alpha_i^T * tmp (for i=1:noca-2) + ! + ! Step 1: Contract aco12 with HOMO_beta MO coefficient + ! tmp_mu = sum_nu P^aco12_(mu,nu) * C^beta_(nu,HOMO) + call dgemm('n','n',nbf,1,nbf, & + one, aco12, nbf, & + vb(:,lr2:lr2), nbf, & + zero, tmp, nbf) + ! Step 2: Add ado2v contribution and subtract aco12 contribution + ! tmp_mu = sum_nu P^ado2v_(mu,nu) * C^beta_(nu,HOMO-1) - tmp_mu + call dgemm('n','n',nbf,1,nbf, & + one, ado2v, nbf, & + vb(:,lr1:lr1), nbf, & + -one, tmp, nbf) + ! Step 3: Project onto doubly-occupied alpha-orbitals + ! F^MO_(i,HOMO-1) += sum_mu C^alpha_(mu,i) * tmp_mu (i=1:noca-2) + call dgemm('t','n',noca-2,1,nbf, & + one, va, nbf, & + tmp, nbf, & + one, wrk(1:noca-2,lr1:lr1), noca-2) + !----------------------------------------------------------------------- + ! Section 5: Corrections for O1(HOMO-1, alpha) -> V(beta) response element + !----------------------------------------------------------------------- + ! Physical meaning: Add corrections to F^MO_(HOMO-1,a) from two sources: + ! 1. adco2: C(alpha) -> O2(HOMO, beta) excitation component (CO block) + ! 2. ao21v: Mixed (O1<->O2)(alpha) x V(beta) component (OV mixed block) + ! + ! This section computes contributions to the O1(HOMO-1, alpha) -> V(beta) block + ! by contracting the transposed AO matrices with alpha-spin MO coefficients, + ! then projecting onto virtual beta-orbitals. + ! + ! Combined transformation: + ! tmp = adco2^T * C^alpha_HOMO-1 + ao21v^T * C^alpha_HOMO + ! F^MO_(HOMO-1,a) += C^beta_a^T * tmp (for a in virt_beta) + ! + ! Step 1: Contract adco2^T with HOMO-1_alpha MO coefficient + ! tmp_mu = sum_nu P^adco2_(nu,mu) * C^alpha_(nu,HOMO-1) + call dgemm('t','n',nbf,1,nbf, & + one, adco2, nbf, & + va(:,lr1:lr1), nbf, & + zero, tmp, nbf) + ! Step 2: Add contribution from ao21v^T with HOMO_alpha MO coefficient + ! tmp_mu += sum_nu P^ao21v_(nu,mu) * C^alpha_(nu,HOMO) + call dgemm('t','n',nbf,1,nbf, & + one, ao21v, nbf, & + va(:,lr2:lr2), nbf, & + one, tmp, nbf) + ! Step 3: Project onto virtual beta-orbitals + ! F^MO_(HOMO-1,a) += sum_mu C^beta_(mu,a) * tmp_mu (a=noca+1:nbf) + call dgemm('t','n',nbf-noca,1,nbf, & + one, vb(:,noca+1), nbf, & + tmp, nbf, & + one, wrk(lr1:lr1,noca+1:nbf), nbf-noca) + !----------------------------------------------------------------------- + ! Section 6: Corrections for O2(HOMO, alpha) -> V(beta) response element + !----------------------------------------------------------------------- + ! Physical meaning: Add corrections to F^MO_(HOMO,a) from two sources: + ! 1. adco1: C(alpha) -> O1(HOMO-1, beta) excitation component (CO block, positive) + ! 2. ao21v: Mixed (O1<->O2)(alpha) x V(beta) component (OV mixed, negative) + ! + ! The subtraction in this section ensures proper antisymmetry between + ! HOMO and HOMO-1 contributions, complementing Section 5 and maintaining + ! consistency with the mixed character of ao21v. This antisymmetry arises from + ! the spin-pairing coupling realized via linear combinations of va and vb. + ! + ! Combined transformation: + ! tmp = adco1^T * C^alpha_HOMO - ao21v^T * C^alpha_HOMO-1 + ! F^MO_(HOMO,a) += C^beta_a^T * tmp (for a in virt_beta) + ! + ! Step 1: Contract ao21v^T with HOMO-1_alpha MO coefficient + ! tmp_mu = sum_nu P^ao21v_(nu,mu) * C^alpha_(nu,HOMO-1) + call dgemm('t','n',nbf,1,nbf, & + one, ao21v, nbf, & + va(:,lr1:lr1), nbf, & + zero, tmp, nbf) + ! Step 2: Add adco1 contribution and subtract ao21v contribution + ! tmp_mu = sum_nu P^adco1_(nu,mu) * C^alpha_(nu,HOMO) - tmp_mu + call dgemm('t','n',nbf,1,nbf, & + one, adco1, nbf, & + va(:,lr2:lr2), nbf, & + -one, tmp, nbf) + ! Step 3: Project onto virtual beta-orbitals + ! F^MO_(HOMO,a) += sum_mu C^beta_(mu,a) * tmp_mu (a=noca+1:nbf) + call dgemm('t','n',nbf-noca,1,nbf, & + one, vb(:,noca+1), nbf, & + tmp, nbf, & + one, wrk(lr2:lr2,noca+1:nbf), nbf-noca) + + !----------------------------------------------------------------------- + ! Spin-dependent corrections for (O1,O1) diagonal element (OO block) + !----------------------------------------------------------------------- + ! Physical meaning: Apply spin-state-dependent transformations to the + ! F^MO_(HOMO-1,HOMO-1) diagonal element based on whether we are computing + ! a singlet (mrst=1) or triplet (mrst=3) excited state. + ! + ! The 1/sqrt(2) factor and the addition/subtraction of diagonal elements ensure + ! proper normalization and spin coupling for the MRSF response equations. + ! This reflects the spin-pairing coupling between MS=+1 and MS=-1 triplet + ! reference states (via linear combinations of va and vb), following Slater-Condon rules. if (mrst==1) then + ! Singlet state (mrst=1): + ! F^MO_(HOMO-1,HOMO-1) = (scr_(HOMO-1,HOMO-1) - scr_(HOMO,HOMO)) / sqrt(2) + ! The subtraction reflects the antisymmetric spin coupling in singlets. + ! The (HOMO,HOMO) element is zeroed as it's not part of the singlet response space. wrk(lr1,lr1) = (scr(lr1,lr1)-scr(lr2,lr2))*sqrt2 wrk(lr2,lr2) = 0.0_dp else if (mrst==3) then + ! Triplet state (mrst=3): + ! F^MO_(HOMO-1,HOMO-1) = (scr_(HOMO-1,HOMO-1) + scr_(HOMO,HOMO)) / sqrt(2) + ! The addition reflects the symmetric spin coupling in triplets. + ! All other OO coupling elements (HOMO,HOMO-1), (HOMO-1,HOMO), and + ! (HOMO,HOMO) are zeroed as they're outside the triplet response space. wrk(lr1,lr1) = (scr(lr1,lr1)+scr(lr2,lr2))*sqrt2 wrk(lr2,lr1) = 0.0_dp wrk(lr1,lr2) = 0.0_dp