module hf_hessian_mod implicit none character(len=*), parameter :: module_name = "hf_hessian_mod" contains !############################################################################### subroutine hf_hessian_C(c_handle) bind(C, name="hf_hessian") use c_interop, only: oqp_handle_t, oqp_handle_get_info use types, only: information type(oqp_handle_t) :: c_handle type(information), pointer :: inf inf => oqp_handle_get_info(c_handle) call hf_hessian(inf) end subroutine hf_hessian_C !############################################################################### subroutine hf_hessian(infos) ! Native OpenQP HF/DFT Hessian CPHF response prepass. ! ! This routine deliberately exercises the production Fortran CPHF/CPKS PCG ! solver for every Cartesian nuclear perturbation used by a ground-state ! analytic Hessian. It builds the closed-shell occupied-virtual RHS from ! OpenQP derivative integrals and the current OpenQP SCF density/MOs, then ! calls cphf_solve on the full 3N RHS block and stores the native Hessian ! matrix in OQP::hf_hessian for the Python frequency driver. use precision, only: dp use types, only: information use basis_tools, only: basis_set use oqp_tagarray_driver, only: tagarray_get_data, OQP_DM_A, OQP_VEC_MO_A, OQP_E_MO_A, & OQP_hf_hessian, TA_TYPE_REAL64 use mathlib, only: unpack_matrix, pack_matrix use grd1, only: der_overlap_matrix, der_kinetic_matrix, der_nucattr_matrix, hess_nn use fock_deriv_mod, only: fock_deriv_contract use scf_addons, only: fock_jk use cphf_mod, only: cphf_solve use io_constants, only: iw use messages, only: show_message, WITH_ABORT implicit none type(information), target, intent(inout) :: infos type(basis_set), pointer :: basis real(kind=dp), contiguous, pointer :: dmat_a(:), mo_a(:,:), eps(:) real(kind=dp), allocatable :: pfull(:,:), probe(:,:), gx(:,:) real(kind=dp), allocatable :: dSa(:,:,:,:), dTa(:,:,:,:), dVa(:,:,:,:) real(kind=dp), allocatable :: Sx(:,:), hx(:,:), F0x(:,:), Gd0(:,:) real(kind=dp), allocatable :: d0(:,:), d0p(:,:), gp(:,:), gfull(:,:) real(kind=dp), allocatable :: bvec(:,:), uvec(:,:), scr(:,:), col(:,:), hess_native(:,:) real(kind=dp), contiguous, pointer :: hess_store(:,:) real(kind=dp) :: hfscale integer :: nbf, nbf2, nocc, nvir, natom, ncart integer :: i, j, a, mu, nu, ia, icart, kc, cc ! Unsupported-feature guards (apply to ALL references, RHF/RKS included). ! Effective-core-potential (ECP) second derivatives ARE supported: RHF/UHF ! contract the ECP skeleton d^2 V_ECP/dR^2 analytically (add_ecphess, libecpint ! deriv order 2) plus the ECP core-derivative in the CPHF response; ROHF folds ! the ECP gradient (add_ecpder) into its semi-numerical resp_grad. ! Range-separated (CAM/LC) functionals are also supported: the 2e derivative ! integrals are erfc-attenuation capable, so grd2_hess_driver (skeleton), ! grd2_driver (fock_deriv_contract response) and fock_jk (cphf) all run the ! long-range Coulomb + short-range erfc-exchange two-pass split when ! infos%dft%cam_flag is set. ! Open-shell (UHF/ROHF) dispatch. The body below is the closed-shell ! (RHF/RKS) kernel: it reads only the alpha density/MOs (OQP_DM_A, mo_a, eps) ! and treats nocc as doubly occupied, so it must never run on an open-shell ! SCF. UHF (scftype==2) -> hf_hessian_uhf, ROHF (scftype==3) -> hf_hessian_rohf ! (both HF and DFT, finite-difference validated). if (infos%control%scftype == 2) then call hf_hessian_uhf(infos) return else if (infos%control%scftype == 3) then call hf_hessian_rohf(infos) return else if (infos%control%scftype > 3) then call show_message('Native analytic Hessian supports RHF/RKS, UHF (HF) '// & 'and ROHF (HF) references only for this scftype. Use [hess] '// & 'type=numerical.', WITH_ABORT) end if basis => infos%basis basis%atoms => infos%atoms nbf = basis%nbf nbf2 = nbf*(nbf+1)/2 nocc = infos%mol_prop%nocc nvir = nbf - nocc natom = size(basis%atoms%xyz, 2) ncart = 3*natom hfscale = 1.0_dp if (infos%control%hamilton >= 20) hfscale = infos%dft%hfscale open(unit=iw, file=infos%log_filename, position="append") write(iw,'(/,A)') 'PyOQP: Native OpenQP HF/DFT Hessian CPHF response prepass' write(iw,'(A,I6,A,I6,A,I6,A,I6)') ' nbf=', nbf, ' nocc=', nocc, ' nvir=', nvir, ' rhs=', ncart write(iw,'(A)') ' Storing native OpenQP HF/DFT analytic Hessian matrix in OQP::hf_hessian.' if (nocc <= 0 .or. nvir <= 0 .or. ncart <= 0) then write(iw,'(A)') ' Native CPHF prepass skipped: empty occupied/virtual/nuclear space.' close(iw) return end if call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a) call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo_a) call tagarray_get_data(infos%dat, OQP_E_MO_A, eps) allocate(pfull(nbf,nbf)); call unpack_matrix(dmat_a, pfull) allocate(dSa(nbf,nbf,3,natom), dTa(nbf,nbf,3,natom), dVa(nbf,nbf,3,natom)) call der_overlap_matrix(basis, dSa) call der_kinetic_matrix(basis, dTa) call der_nucattr_matrix(basis, basis%atoms%xyz, & basis%atoms%zn - basis%ecp_zn_num, dVa) ! ECP-screened point charge ! der_* matrices are returned in the UNNORMALIZED basis; bring them into the ! same normalized (bfnrm) convention as the MO coefficients / density so the ! CPHF RHS and the response contractions are correct for d/f functions ! (bfnrm /= 1). Invisible for s/p-only bases (e.g. STO-3G). block integer :: kc2, cc2, mu2, nu2 do kc2 = 1, natom do cc2 = 1, 3 do nu2 = 1, nbf do mu2 = 1, nbf dSa(mu2,nu2,cc2,kc2) = dSa(mu2,nu2,cc2,kc2)*basis%bfnrm(mu2)*basis%bfnrm(nu2) dTa(mu2,nu2,cc2,kc2) = dTa(mu2,nu2,cc2,kc2)*basis%bfnrm(mu2)*basis%bfnrm(nu2) dVa(mu2,nu2,cc2,kc2) = dVa(mu2,nu2,cc2,kc2)*basis%bfnrm(mu2)*basis%bfnrm(nu2) end do end do end do end do end block ! ECP first-derivative integrals enter the core-Hamiltonian derivative ! dHcore/dR (added into dVa, the nuclear-attraction derivative tensor), so the ! ECP contributes to the CPHF right-hand side and the orbital-relaxation ! response exactly as point-charge nuclear attraction does. libecpint returns ! these already in the OpenQP normalized convention, hence added AFTER the ! bfnrm scaling above. No-op for non-ECP bases. block use ecp_tool, only: ecp_deriv_ints real(kind=dp), allocatable :: dVecp(:,:,:,:) allocate(dVecp(nbf,nbf,3,natom)) call ecp_deriv_ints(basis, basis%atoms%xyz, dVecp) dVa = dVa + dVecp deallocate(dVecp) end block allocate(scr(nbf,nbf), col(nbf,nbf)) allocate(Sx(nbf,nbf), hx(nbf,nbf), F0x(nbf,nbf), Gd0(nbf,nbf)) allocate(probe(nbf,nbf), gx(3,natom)) allocate(d0(nbf,nbf), d0p(nbf2,1), gp(nbf2,1), gfull(nbf,nbf)) allocate(bvec(nocc*nvir,ncart), uvec(nocc*nvir,ncart), source=0.0_dp) icart = 0 do kc = 1, natom do cc = 1, 3 icart = icart + 1 call mo_transform(mo_a, dSa(:,:,cc,kc), nbf, scr, col, Sx) scr = dTa(:,:,cc,kc) + dVa(:,:,cc,kc) call mo_transform(mo_a, scr, nbf, col, F0x, hx) F0x = hx do a = 1, nvir do i = 1, nocc do mu = 1, nbf do nu = 1, nbf probe(mu,nu) = 0.5_dp*( mo_a(mu,nocc+a)*mo_a(nu,i) + mo_a(mu,i)*mo_a(nu,nocc+a) ) end do end do call fock_deriv_contract(infos, basis, pfull, probe, hfscale, gx) F0x(i,nocc+a) = hx(i,nocc+a) + 2.0_dp*gx(cc,kc) end do end do ! --- XC contribution to the CPKS right-hand side (DFT only) ------------- ! The A-matrix (cphf_apbx) includes the XC kernel fxc, so the perturbation ! RHS must carry BOTH XC pieces or the relaxed response dPx is wrong (the ! HF response, ~2x too large for DFT): ! (i) skeleton dVxc/dR (fixed orbitals, basis+grid move) -> in F0x ! (ii) fxc[d0], d0 = reorthonormalization density -> in Gd0 ! Both enter B_ai with the SAME (minus) sign as the other Fock terms, so ! they are captured together by ONE central FD of the XC Fock matrix ! (dftexcor) along the combined path: geometry R +/- h AND occupied MOs ! reorthonormalized by dmo_i = -1/2 sum_j C_j S^x_ji. dftexcor handles all ! density/spin scale factors internally, so no manual convention factors. if (infos%control%hamilton == 20) then block use dft, only: dft_initialize, dftclean, dftexcor use mod_dft_molgrid, only: dft_grid_t type(dft_grid_t) :: mgr real(dp), allocatable :: dmoR(:,:), mop(:,:), frp(:), frm(:), dVxcR(:,:), hxcR(:,:) real(dp) :: hxr, telr, tknr, exr integer :: ir, jr allocate(dmoR(nbf,nocc), mop(nbf,nbf), frp(nbf2), frm(nbf2), dVxcR(nbf,nbf), hxcR(nbf,nbf)) hxr = 1.0d-3 dmoR = 0.0_dp do ir = 1, nocc do jr = 1, nocc dmoR(:,ir) = dmoR(:,ir) - 0.5_dp*mo_a(:,jr)*Sx(jr,ir) end do end do basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) + hxr call basis%init_shell_centers() call dft_initialize(infos, basis, mgr) mop = mo_a; mop(:,1:nocc) = mo_a(:,1:nocc) + hxr*dmoR frp = 0.0_dp call dftexcor(basis, mgr, 1, frp, frp, mop, mop, nbf, nbf2, exr, telr, tknr, infos) call dftclean(infos) basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) - 2*hxr call basis%init_shell_centers() call dft_initialize(infos, basis, mgr) mop = mo_a; mop(:,1:nocc) = mo_a(:,1:nocc) - hxr*dmoR frm = 0.0_dp call dftexcor(basis, mgr, 1, frm, frm, mop, mop, nbf, nbf2, exr, telr, tknr, infos) call dftclean(infos) basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) + hxr call basis%init_shell_centers() call unpack_from_packed((frp - frm)/(2*hxr), dVxcR, nbf) call mo_transform(mo_a, dVxcR, nbf, scr, col, hxcR) do a = 1, nvir do i = 1, nocc F0x(i,nocc+a) = F0x(i,nocc+a) + hxcR(i,nocc+a) end do end do deallocate(dmoR, mop, frp, frm, dVxcR, hxcR) end block end if d0 = 0.0_dp do i = 1, nocc do j = 1, nocc do mu = 1, nbf do nu = 1, nbf d0(mu,nu) = d0(mu,nu) - 2.0_dp*Sx(i,j)*mo_a(mu,i)*mo_a(nu,j) end do end do end do end do call pack_matrix(d0, d0p(:,1)) gp = 0.0_dp call fock_jk(basis, d=d0p, f=gp, scale_exch=hfscale, infos=infos) call unpack_from_packed(gp(:,1), gfull, nbf) call mo_transform(mo_a, gfull, nbf, scr, col, Gd0) ia = 0 do a = 1, nvir do i = 1, nocc ia = ia + 1 bvec(ia,icart) = -F0x(i,nocc+a) + eps(i)*Sx(i,nocc+a) - Gd0(i,nocc+a) end do end do end do end do call cphf_solve(infos, ncart, bvec, uvec) ! ===== CPHF orbital-relaxation response ===== ! H^resp_xy = 4 Tr[F^x dm1^y] - 4 Tr[S^x (eps.dm1^y)] - 2 Tr[s1oo^x mo_e1^y] ! dm1^y_pq = sum_k dC^y_pk C_qk (one-sided) ! mo_e1^y_kl = (h^y + G[P]^y + G[dP^y])^MO_kl - 1/2 (eps_k+eps_l) s1oo^y_kl ! F^x = h^x + G[P]^x; dC^y from the validated CPHF amplitudes U^y. The first ! two terms equal Tr[dP^y F^x] and the eps-weighted overlap term; the third ! is the FULL occ-occ energy-weighted term (the off-diagonal part is what a ! diagonal dε approximation misses). 2e traces use fock_deriv_contract ! (=1/2 Tr[M G[P]^x]) and fock_jk (G[dP^y]). allocate(hess_native(ncart,ncart), source=0.0_dp) block real(dp), allocatable :: sflat(:,:,:), hflat(:,:,:) real(dp), allocatable :: dCx(:,:,:), dPx(:,:,:), Gdp(:,:,:) real(dp), allocatable :: s1oo(:,:,:), hMOoo(:,:,:), GdpMOoo(:,:,:), moe1a(:,:,:) real(dp), allocatable :: Mi(:,:), gxy(:,:), A2(:,:), tGP(:,:), hresp(:,:) real(dp), allocatable :: s1(:,:), s2(:,:), bMO(:,:), dpp(:,:), gpp(:,:), gfl(:,:) real(dp), allocatable :: cocc(:,:), tmpno(:,:) real(dp) :: a1v, a3v, t3a, dcsx integer :: x, yy, ii, jj, kk, ll, aa, ia2, mu2, nu2, ccx, kcx allocate(sflat(nbf,nbf,ncart), hflat(nbf,nbf,ncart)) do x = 1, ncart ccx = mod(x-1,3)+1; kcx = (x-1)/3+1 sflat(:,:,x) = dSa(:,:,ccx,kcx) hflat(:,:,x) = dTa(:,:,ccx,kcx) + dVa(:,:,ccx,kcx) end do allocate(cocc(nbf,nocc)); cocc = mo_a(:,1:nocc) ! occ-occ MO blocks of S^x and h^x allocate(s1oo(nocc,nocc,ncart), hMOoo(nocc,nocc,ncart), source=0.0_dp) allocate(s1(nbf,nbf), s2(nbf,nbf), bMO(nbf,nbf), tmpno(nbf,nocc)) do x = 1, ncart call dgemm('n','n',nbf,nocc,nbf,1.0_dp,sflat(:,:,x),nbf,cocc,nbf,0.0_dp,tmpno,nbf) call dgemm('t','n',nocc,nocc,nbf,1.0_dp,cocc,nbf,tmpno,nbf,0.0_dp,s1oo(:,:,x),nocc) call dgemm('n','n',nbf,nocc,nbf,1.0_dp,hflat(:,:,x),nbf,cocc,nbf,0.0_dp,tmpno,nbf) call dgemm('t','n',nocc,nocc,nbf,1.0_dp,cocc,nbf,tmpno,nbf,0.0_dp,hMOoo(:,:,x),nocc) end do ! relaxed orbital derivative dC^y, density dP^y (total), response Fock G[dP^y] allocate(dCx(nbf,nocc,ncart), dPx(nbf,nbf,ncart), Gdp(nbf,nbf,ncart), source=0.0_dp) allocate(GdpMOoo(nocc,nocc,ncart), source=0.0_dp) allocate(dpp(nbf2,1), gpp(nbf2,1), gfl(nbf,nbf)) do yy = 1, ncart ia2 = 0 do aa = 1, nvir do ii = 1, nocc ia2 = ia2 + 1 dCx(:,ii,yy) = dCx(:,ii,yy) + mo_a(:,nocc+aa)*uvec(ia2,yy) end do end do do ii = 1, nocc do jj = 1, nocc dCx(:,ii,yy) = dCx(:,ii,yy) - 0.5_dp*mo_a(:,jj)*s1oo(jj,ii,yy) end do end do do ii = 1, nocc do mu2 = 1, nbf do nu2 = 1, nbf dPx(mu2,nu2,yy) = dPx(mu2,nu2,yy) & + 2.0_dp*(dCx(mu2,ii,yy)*mo_a(nu2,ii) + mo_a(mu2,ii)*dCx(nu2,ii,yy)) end do end do end do call pack_matrix(dPx(:,:,yy), dpp(:,1)) gpp = 0.0_dp call fock_jk(basis, d=dpp, f=gpp, scale_exch=hfscale, infos=infos) call unpack_from_packed(gpp(:,1), gfl, nbf); Gdp(:,:,yy) = gfl call dgemm('n','n',nbf,nocc,nbf,1.0_dp,gfl,nbf,cocc,nbf,0.0_dp,tmpno,nbf) call dgemm('t','n',nocc,nocc,nbf,1.0_dp,cocc,nbf,tmpno,nbf,0.0_dp,GdpMOoo(:,:,yy),nocc) end do ! mo_e1 without the G[P]^y part (added via Mi trick in term3) allocate(moe1a(nocc,nocc,ncart)) do yy = 1, ncart do ll = 1, nocc do kk = 1, nocc moe1a(kk,ll,yy) = hMOoo(kk,ll,yy) + GdpMOoo(kk,ll,yy) & - 0.5_dp*(eps(kk)+eps(ll))*s1oo(kk,ll,yy) end do end do end do ! 2e traces: A2(x,y)=Tr[dP^y G[P]^x]; tGP(x,y)=Tr[M^x G[P]^y] ! with M^x = sum_kl s1oo^x_kl C_k C_l^T allocate(gxy(3,natom), A2(ncart,ncart), tGP(ncart,ncart), Mi(nbf,nbf), source=0.0_dp) do yy = 1, ncart gxy = 0.0_dp call fock_deriv_contract(infos, basis, pfull, dPx(:,:,yy), hfscale, gxy) A2(:,yy) = 2.0_dp*reshape(gxy, [ncart]) end do do x = 1, ncart call dgemm('n','n',nbf,nocc,nocc,1.0_dp,cocc,nbf,s1oo(:,:,x),nocc,0.0_dp,tmpno,nbf) call dgemm('n','t',nbf,nbf,nocc,1.0_dp,tmpno,nbf,cocc,nbf,0.0_dp,Mi,nbf) gxy = 0.0_dp call fock_deriv_contract(infos, basis, pfull, Mi, hfscale, gxy) tGP(x,:) = 2.0_dp*reshape(gxy, [ncart]) end do ! assemble response hresp(x,y) = 4Tr[F^x dm1^y]-4Tr[S^x eps.dm1^y]-2Tr[s1oo^x mo_e1^y] ! = (Tr[dP^y h^x] + A2) - 4 A3 - 2 (sum_kl s1oo^x_kl moe1a^y_kl) - 2 tGP allocate(hresp(ncart,ncart), source=0.0_dp) do x = 1, ncart do yy = 1, ncart a1v = sum(dPx(:,:,yy)*hflat(:,:,x)) a3v = 0.0_dp do ii = 1, nocc dcsx = 0.0_dp do mu2 = 1, nbf do nu2 = 1, nbf dcsx = dcsx + dCx(mu2,ii,yy)*sflat(mu2,nu2,x)*mo_a(nu2,ii) end do end do a3v = a3v + eps(ii)*dcsx end do t3a = 0.0_dp do ll = 1, nocc do kk = 1, nocc t3a = t3a + s1oo(kk,ll,x)*moe1a(kk,ll,yy) end do end do hresp(x,yy) = (a1v + A2(x,yy)) - 4.0_dp*a3v - 2.0_dp*t3a - 2.0_dp*tGP(x,yy) end do end do hess_native = 0.5_dp*(hresp + transpose(hresp)) ! --- DFT exchange-correlation second-derivative contribution ----------- ! The XC part of the Hessian is obtained by central finite differencing the ! analytic XC nuclear gradient (derexc_blk) over geometry while displacing ! the density by the analytic relaxed density derivative dP^y. This adds ! both the XC skeleton (d2Exc/dR2 at fixed density) and the XC response ! (through dP^y) in one shot, with no re-SCF. The HF-exchange fraction is ! already in the Coulomb/exchange terms above (hfscale); derexc_blk ! supplies the remaining DFT exchange-correlation functional. if (infos%control%hamilton == 20) then block use dft, only: dft_initialize, dftclean, dftexcor use mod_dft_gridint_grad, only: derexc_blk use mod_dft_molgrid, only: dft_grid_t type(dft_grid_t) :: mg real(dp), allocatable :: dap(:,:), dedp(:,:), dedm(:,:) real(dp), allocatable :: mop(:,:), frp(:), frm(:), dFxc(:,:), dFoo(:,:) real(dp), allocatable :: tmpn(:,:), dHse(:,:), dHt3(:,:) real(dp) :: hx, tele, tkin, eexc integer :: yy2, ccy, kcy, nang, x2, kk2, ll2 hx = 1.0d-3; nang = maxval(basis%am) + 2 ! XC contribution split into a skeleton+density-response term and an ! energy-weighting term, realised through the OpenQP moving-grid XC ! machinery so it stays consistent with the OpenQP numerical Hessian: ! dHse : skeleton + density-response (term1). Central FD of the analytic ! XC gradient (derexc) along the relaxed path R+lambda, P+lambda*dP. ! This is the genuine total derivative d/dR[g_XC(R,P(R))] of the ! OpenQP XC gradient, so the moving-grid weight derivatives are ! handled identically to the SCF/numerical-gradient convention. ! dHt3 : -2 Tr[s1oo^x (vxc^y+fxc[dP^y])_oo], the XC part of the ! energy-weighted (mo_e1) term, from the FD of the XC Fock ! matrix (dftexcor) along the same relaxed orbital path. allocate(dap(nbf,nbf), dedp(3,natom), dedm(3,natom)) allocate(mop(nbf,nbf), frp(nbf2), frm(nbf2), dFxc(nbf,nbf), dFoo(nocc,nocc)) allocate(tmpn(nbf,nocc), dHse(ncart,ncart), dHt3(ncart,ncart)) dHt3 = 0.0_dp ! warm-up to flush any stale grid state left by the CPHF solver call dft_initialize(infos, basis, mg); call dftclean(infos) do yy2 = 1, ncart ccy = mod(yy2-1,3)+1; kcy = (yy2-1)/3+1 basis%atoms%xyz(ccy,kcy) = basis%atoms%xyz(ccy,kcy) + hx call basis%init_shell_centers() call dft_initialize(infos, basis, mg) dap = pfull + hx*dPx(:,:,yy2); dedp = 0.0_dp ! skeleton + density response call derexc_blk(basis, mg, dap, dap, dedp, tele, tkin, nang, nbf, & infos%dft%grid_density_cutoff, .false., infos) mop = mo_a; mop(:,1:nocc) = mo_a(:,1:nocc) + hx*dCx(:,:,yy2) call dftexcor(basis, mg, 1, frp, frp, mop, mop, nbf, nbf2, eexc, tele, tkin, infos) call dftclean(infos) basis%atoms%xyz(ccy,kcy) = basis%atoms%xyz(ccy,kcy) - 2*hx call basis%init_shell_centers() call dft_initialize(infos, basis, mg) dap = pfull - hx*dPx(:,:,yy2); dedm = 0.0_dp call derexc_blk(basis, mg, dap, dap, dedm, tele, tkin, nang, nbf, & infos%dft%grid_density_cutoff, .false., infos) mop = mo_a; mop(:,1:nocc) = mo_a(:,1:nocc) - hx*dCx(:,:,yy2) call dftexcor(basis, mg, 1, frm, frm, mop, mop, nbf, nbf2, eexc, tele, tkin, infos) call dftclean(infos) basis%atoms%xyz(ccy,kcy) = basis%atoms%xyz(ccy,kcy) + hx call basis%init_shell_centers() dHse(:,yy2) = reshape((dedp - dedm)/(2*hx), [ncart]) ! term3: -2 s1oo^x (vxc^y + fxc[dP^y])_oo call unpack_from_packed((frp - frm)/(2*hx), dFxc, nbf) call dgemm('n','n',nbf,nocc,nbf,1.0_dp,dFxc,nbf,mo_a,nbf,0.0_dp,tmpn,nbf) call dgemm('t','n',nocc,nocc,nbf,1.0_dp,mo_a,nbf,tmpn,nbf,0.0_dp,dFoo,nocc) do x2 = 1, ncart do ll2 = 1, nocc do kk2 = 1, nocc dHt3(x2,yy2) = dHt3(x2,yy2) - 2.0_dp*s1oo(kk2,ll2,x2)*dFoo(kk2,ll2) end do end do end do end do hess_native = hess_native + 0.5_dp*(dHse + transpose(dHse)) & + 0.5_dp*(dHt3 + transpose(dHt3)) deallocate(dap, dedp, dedm, mop, frp, frm, dFxc, dFoo, tmpn, dHse, dHt3) end block end if deallocate(sflat, hflat, dCx, dPx, Gdp, s1oo, hMOoo, GdpMOoo, moe1a, & Mi, gxy, A2, tGP, hresp, s1, s2, bMO, dpp, gpp, gfl, cocc, tmpno) end block call hess_nn(basis%atoms, basis%ecp_zn_num, hess_native) ! --- One-electron + Pulay second-derivative skeleton (fixed density) ------ ! Mirrors the production HF gradient assembly (hf_1e_grad): the analytic ! Hessian skeleton is d/dx of [grad_ee_overlap(W) + grad_ee_kinetic(P) ! + grad_en(P)] evaluated at the fixed converged density, i.e. the ! second-derivative integral contractions hess_ee_overlap / hess_ee_kinetic ! / hess_en. This is distinct from (and additive to) the CPHF response ! term above; the 2e ERI second-derivative skeleton is added separately. block use grd1, only: eijden, hess_ee_overlap, hess_ee_kinetic, hess_en use ecp_tool, only: add_ecphess real(kind=dp), allocatable :: wlag(:), pden(:), hcc(:,:) allocate(wlag(nbf2), pden(nbf2), hcc(ncart,ncart), source=0.0_dp) call eijden(wlag, nbf, infos) ! energy-weighted (Lagrangian) density pden = dmat_a ! total density (closed-shell RHF) call hess_ee_overlap(basis, wlag, hess_native) ! overlap / Pulay call hess_ee_kinetic(basis, pden, hess_native) ! kinetic call hess_en(basis, basis%atoms%xyz, & basis%atoms%zn - basis%ecp_zn_num, pden, hess_native, hess_cc=hcc) call add_ecphess(basis, basis%atoms%xyz, pden, hess_native) ! ECP skeleton (if any) deallocate(wlag, pden, hcc) end block ! --- Two-electron (ERI) second-derivative skeleton (fixed density) -------- ! d^2/dR^2 of the analytic 2e gradient contraction at the converged density, ! i.e. sum P P d^2/dR^2 [ (ij|kl) - 1/4 c_x (ik|jl) ]. Validated against a ! finite difference of grd2_driver (see grd2_hess_selftest). Additive to the ! CPHF response and 1e skeleton above. block use grd2, only: grd2_hess_driver, grd2_compute_data_t use hf_gradient_mod, only: grd2_rhf_compute_data_t type(grd2_rhf_compute_data_t) :: gcomp gcomp = grd2_rhf_compute_data_t( da = dmat_a, hfscale = hfscale, nbf = nbf ) call gcomp%init() call gcomp%build_cart(basis) call grd2_hess_driver(infos, basis, hess_native, gcomp) call gcomp%clean() end block call infos%dat%alloc_or_die(OQP_hf_hessian, (/ ncart, ncart /), hess_store, & description='Native OpenQP HF/DFT analytic Hessian matrix') hess_store = hess_native write(iw,'(A)') 'PyOQP: Native OpenQP HF/DFT Hessian matrix stored' close(iw) deallocate(pfull, dSa, dTa, dVa, scr, col, Sx, hx, F0x, Gd0, probe, gx, & d0, d0p, gp, gfull, bvec, uvec, hess_native) end subroutine hf_hessian !############################################################################### subroutine hf_hessian_uhf(infos) ! Native open-shell (UHF) analytic HF Hessian. ! ! Mirrors the closed-shell hf_hessian response assembly per spin, summed over ! s in {alpha, beta} with single (not doubled) occupation factors. Each spin ! uses its own MO set C^s, orbital energies eps^s and density P^s; the ! two-electron couplings are open-shell (Coulomb from the total density ! P = Pa + Pb, exchange from the spin density P^s): ! ! B^s_ia = -(h^x_ia + G^{s,x}[P]_ia) + eps^s_i S^x_ia - G^s[d0]_ia , ! d0^s = -sum_ij S^x,s_ij C^s_i C^s_j^T (reorthonormalization), ! G^s[.] = J[.^a + .^b] - c_x K[.^s] (scf_addons::fock_jk), ! G^{s,x}[P] via fock_deriv_mod::fock_deriv_contract_os (Coulomb P, exch P^s). ! ! The 3N right-hand sides are solved with cphf_mod::cphf_solve_uhf, and the ! orbital-relaxation response is assembled as (per spin, summed): ! ! H^resp_xy = sum_s [ Tr[dP^s,y h^x] + Tr[dP^s,y G^{s,x}[P]] ] ! - 2 sum_s sum_i eps^s_i (dC^s,y_i . S^x . C^s_i) ! - sum_s sum_kl s1oo^s,x_kl moe1^s,y_kl ! - sum_s Tr[Mi^s,x G^{s,y}[P]] , ! ! moe1^s,y_kl = h^x_kl(MO) + G^s[dP^y]_kl(MO) - 1/2(eps^s_k+eps^s_l) s1oo^s,y_kl , ! Mi^s,x = sum_kl s1oo^s,x_kl C^s_k C^s_l^T . ! ! The fixed-density skeleton (1e total density + open-shell Lagrangian W, 2e ! via grd2_uhf_compute_data_t) and the nuclear-repulsion term are added on ! top, exactly as in hess_skel_open_selftest. HF only (the UKS f_xc response ! is not finite-difference validated). use precision, only: dp use types, only: information use basis_tools, only: basis_set use oqp_tagarray_driver, only: tagarray_get_data, OQP_DM_A, OQP_DM_B, & OQP_VEC_MO_A, OQP_VEC_MO_B, OQP_E_MO_A, OQP_E_MO_B, OQP_hf_hessian, TA_TYPE_REAL64 use mathlib, only: unpack_matrix, pack_matrix use grd1, only: der_overlap_matrix, der_kinetic_matrix, der_nucattr_matrix, hess_nn use fock_deriv_mod, only: fock_deriv_contract_os use scf_addons, only: fock_jk use cphf_mod, only: cphf_solve_uhf use io_constants, only: iw implicit none type(information), target, intent(inout) :: infos !> Per-spin work container (alpha/beta have different nocc/nvir). type :: uhf_spin_t real(dp), allocatable :: mo(:,:) ! MO coefficients (nbf,nbf) real(dp), allocatable :: eps(:) ! orbital energies (nbf) real(dp), allocatable :: p(:,:) ! spin AO density (nbf,nbf) integer :: nocc = 0, nvir = 0, loff = 0 ! occ/vir count, CPHF block offset real(dp), allocatable :: s1oo(:,:,:) ! occ-occ MO of S^x (nocc,nocc,ncart) real(dp), allocatable :: hoo(:,:,:) ! occ-occ MO of h^x real(dp), allocatable :: g2e(:,:) ! G^{s,x}[P]_ia for all coords (nocc*nvir,ncart) real(dp), allocatable :: dCx(:,:,:) ! relaxed dC (nbf,nocc,ncart) real(dp), allocatable :: dPx(:,:,:) ! relaxed spin density derivative real(dp), allocatable :: gdpoo(:,:,:) ! occ-occ MO of G^s[dP^y] real(dp), allocatable :: moe1(:,:,:) ! occ-occ energy-weighted derivative end type type(basis_set), pointer :: basis real(dp), contiguous, pointer :: dma(:), dmb(:), moa(:,:), mob(:,:), epsa(:), epsb(:) real(dp), contiguous, pointer :: hess_store(:,:) real(dp), allocatable :: ptot(:,:), dSa(:,:,:,:), dTa(:,:,:,:), dVa(:,:,:,:) real(dp), allocatable :: sflat(:,:,:), hflat(:,:,:) real(dp), allocatable :: bvec(:,:), uvec(:,:), hess_native(:,:) real(dp), allocatable :: scr(:,:), tmp(:,:), gx(:,:), probe(:,:) real(dp), allocatable :: SxMO(:,:), hxMO(:,:), d0a(:,:), d0b(:,:) real(dp), allocatable :: dpck(:,:), fpck(:,:), gfull(:,:) real(dp), allocatable :: Gd0(:,:), Mi(:,:) real(dp), allocatable :: A2(:,:), tGP(:,:), hresp(:,:) type(uhf_spin_t) :: sp(2) real(dp) :: hfscale, a1v, a3v, t3a, dcsx integer :: nbf, nbf2, natom, ncart, nocca, noccb, nvira, nvirb, la, lb, ltot integer :: s, i, j, a, ia, icart, kc, cc, x, yy, kk, ll, mu, nu basis => infos%basis basis%atoms => infos%atoms nbf = basis%nbf nbf2 = nbf*(nbf+1)/2 natom = size(basis%atoms%xyz, 2) ncart = 3*natom nocca = infos%mol_prop%nelec_A noccb = infos%mol_prop%nelec_B nvira = nbf - nocca nvirb = nbf - noccb la = nocca*nvira lb = noccb*nvirb ltot = la + lb hfscale = 1.0_dp if (infos%control%hamilton >= 20) hfscale = infos%dft%hfscale write(iw,'(/,A)') 'PyOQP: Native OpenQP open-shell (UHF) HF Hessian CPHF response prepass' write(iw,'(A,I6,A,I6,A,I6,A,I6,A,I6)') ' nbf=', nbf, ' nocca=', nocca, & ' noccb=', noccb, ' rhs=', ncart, ' ltot=', ltot write(iw,'(A)') ' Storing native OpenQP open-shell HF analytic Hessian in OQP::hf_hessian.' if (ncart <= 0 .or. (la <= 0 .and. lb <= 0)) then write(iw,'(A)') ' UHF CPHF prepass skipped: empty occupied/virtual/nuclear space.' return end if call tagarray_get_data(infos%dat, OQP_DM_A, dma) call tagarray_get_data(infos%dat, OQP_DM_B, dmb) call tagarray_get_data(infos%dat, OQP_VEC_MO_A, moa) call tagarray_get_data(infos%dat, OQP_VEC_MO_B, mob) call tagarray_get_data(infos%dat, OQP_E_MO_A, epsa) call tagarray_get_data(infos%dat, OQP_E_MO_B, epsb) ! per-spin containers sp(1)%nocc = nocca; sp(1)%nvir = nvira; sp(1)%loff = 0 sp(2)%nocc = noccb; sp(2)%nvir = nvirb; sp(2)%loff = la allocate(sp(1)%mo(nbf,nbf), sp(1)%eps(nbf), sp(1)%p(nbf,nbf)) allocate(sp(2)%mo(nbf,nbf), sp(2)%eps(nbf), sp(2)%p(nbf,nbf)) sp(1)%mo = moa; sp(1)%eps = epsa sp(2)%mo = mob; sp(2)%eps = epsb call unpack_matrix(dma, sp(1)%p) call unpack_matrix(dmb, sp(2)%p) allocate(ptot(nbf,nbf)); ptot = sp(1)%p + sp(2)%p ! derivative integrals (normalized into the bfnrm convention of the MOs) allocate(dSa(nbf,nbf,3,natom), dTa(nbf,nbf,3,natom), dVa(nbf,nbf,3,natom)) call der_overlap_matrix(basis, dSa) call der_kinetic_matrix(basis, dTa) call der_nucattr_matrix(basis, basis%atoms%xyz, & basis%atoms%zn - basis%ecp_zn_num, dVa) ! ECP-screened point charge block integer :: kc2, cc2, mu2, nu2 do kc2 = 1, natom do cc2 = 1, 3 do nu2 = 1, nbf do mu2 = 1, nbf dSa(mu2,nu2,cc2,kc2) = dSa(mu2,nu2,cc2,kc2)*basis%bfnrm(mu2)*basis%bfnrm(nu2) dTa(mu2,nu2,cc2,kc2) = dTa(mu2,nu2,cc2,kc2)*basis%bfnrm(mu2)*basis%bfnrm(nu2) dVa(mu2,nu2,cc2,kc2) = dVa(mu2,nu2,cc2,kc2)*basis%bfnrm(mu2)*basis%bfnrm(nu2) end do end do end do end do end block ! ECP first-derivative integrals -> core-Hamiltonian derivative dHcore/dR (see ! the RHF kernel for the rationale). Already in the normalized convention, so ! added after the bfnrm scaling. No-op for non-ECP bases. block use ecp_tool, only: ecp_deriv_ints real(dp), allocatable :: dVecp(:,:,:,:) allocate(dVecp(nbf,nbf,3,natom)) call ecp_deriv_ints(basis, basis%atoms%xyz, dVecp) dVa = dVa + dVecp deallocate(dVecp) end block ! flat (ncart) AO views of S^x and h^x = (T+V)^x allocate(sflat(nbf,nbf,ncart), hflat(nbf,nbf,ncart)) do x = 1, ncart cc = mod(x-1,3)+1; kc = (x-1)/3+1 sflat(:,:,x) = dSa(:,:,cc,kc) hflat(:,:,x) = dTa(:,:,cc,kc) + dVa(:,:,cc,kc) end do ! occ-occ and occ MO transforms needed by the response assembly allocate(scr(nbf,nbf), tmp(nbf,nbf), SxMO(nbf,nbf), hxMO(nbf,nbf)) do s = 1, 2 allocate(sp(s)%s1oo(sp(s)%nocc, sp(s)%nocc, ncart), source=0.0_dp) allocate(sp(s)%hoo (sp(s)%nocc, sp(s)%nocc, ncart), source=0.0_dp) do x = 1, ncart call mo_transform(sp(s)%mo, sflat(:,:,x), nbf, scr, tmp, SxMO) call mo_transform(sp(s)%mo, hflat(:,:,x), nbf, scr, tmp, hxMO) sp(s)%s1oo(:,:,x) = SxMO(1:sp(s)%nocc,1:sp(s)%nocc) sp(s)%hoo (:,:,x) = hxMO(1:sp(s)%nocc,1:sp(s)%nocc) end do end do ! ===== CPHF right-hand sides B^s (occ-vir) for all 3N perturbations ===== allocate(bvec(ltot,ncart), uvec(ltot,ncart), source=0.0_dp) allocate(probe(nbf,nbf), gx(3,natom)) allocate(d0a(nbf,nbf), d0b(nbf,nbf), gfull(nbf,nbf), Gd0(nbf,nbf)) allocate(dpck(nbf2,2), fpck(nbf2,2)) ! 2e response-Fock skeleton G^{s,x}[P]_ia for ALL 3N coordinates. The ! occ-vir probe C^s_a C^s_i^T is geometry-independent, so a single open-shell ! derivative-Fock contraction per occ-vir pair yields every Cartesian ! component at once (avoids an ncart-fold redundant grd2 sweep). do s = 1, 2 allocate(sp(s)%g2e(sp(s)%nocc*sp(s)%nvir, ncart), source=0.0_dp) do a = 1, sp(s)%nvir do i = 1, sp(s)%nocc do mu = 1, nbf do nu = 1, nbf probe(mu,nu) = 0.5_dp*( sp(s)%mo(mu,sp(s)%nocc+a)*sp(s)%mo(nu,i) & + sp(s)%mo(mu,i)*sp(s)%mo(nu,sp(s)%nocc+a) ) end do end do gx = 0.0_dp call fock_deriv_contract_os(infos, basis, ptot, sp(s)%p, probe, hfscale, gx) ia = (a-1)*sp(s)%nocc + i sp(s)%g2e(ia,:) = reshape(gx, [ncart]) end do end do end do icart = 0 do kc = 1, natom do cc = 1, 3 icart = icart + 1 ! reorthonormalization density per spin: d0^s = -sum_ij S^x,s_ij C^s_i C^s_j^T d0a = 0.0_dp; d0b = 0.0_dp do s = 1, 2 call mo_transform(sp(s)%mo, dSa(:,:,cc,kc), nbf, scr, tmp, SxMO) do i = 1, sp(s)%nocc do j = 1, sp(s)%nocc do mu = 1, nbf do nu = 1, nbf if (s == 1) then d0a(mu,nu) = d0a(mu,nu) - SxMO(i,j)*sp(s)%mo(mu,i)*sp(s)%mo(nu,j) else d0b(mu,nu) = d0b(mu,nu) - SxMO(i,j)*sp(s)%mo(mu,i)*sp(s)%mo(nu,j) end if end do end do end do end do end do call pack_matrix(d0a, dpck(:,1)) call pack_matrix(d0b, dpck(:,2)) fpck = 0.0_dp call fock_jk(basis, d=dpck, f=fpck, scale_exch=hfscale, infos=infos) do s = 1, 2 call mo_transform(sp(s)%mo, dSa(:,:,cc,kc), nbf, scr, tmp, SxMO) call mo_transform(sp(s)%mo, dTa(:,:,cc,kc)+dVa(:,:,cc,kc), nbf, scr, tmp, hxMO) call unpack_from_packed(fpck(:,s), gfull, nbf) ! G^s[d0] call mo_transform(sp(s)%mo, gfull, nbf, scr, tmp, Gd0) do a = 1, sp(s)%nvir do i = 1, sp(s)%nocc ia = (a-1)*sp(s)%nocc + i bvec(sp(s)%loff+ia,icart) = & -(hxMO(i,sp(s)%nocc+a) + sp(s)%g2e(ia,icart)) & + sp(s)%eps(i)*SxMO(i,sp(s)%nocc+a) & - Gd0(i,sp(s)%nocc+a) end do end do end do ! --- XC contribution to the CPKS right-hand side (UKS only) ------------- ! Mirror the closed-shell RKS RHS XC: one central FD of the spin XC Fock ! matrices (open-shell dftexcor) along R +/- h AND occupied MOs reorthonor- ! malized by dmoR^s = -1/2 sum_j C^s_j S^x,s_ji captures both the XC ! skeleton dVxc/dR and f_xc[d0]; subtract the vir-occ MO blocks from B ! (which carries -F0x). if (infos%control%hamilton >= 20) then block use dft, only: dft_initialize, dftclean, dftexcor use mod_dft_molgrid, only: dft_grid_t type(dft_grid_t) :: mgr real(dp), allocatable :: dmoa(:,:), dmob(:,:), mopa(:,:), mopb(:,:) real(dp), allocatable :: SxMOa(:,:), SxMOb(:,:) real(dp), allocatable :: frap(:), frbp(:), fram(:), frbm(:), dvx(:,:), hxc(:,:) real(dp) :: hxr, exr, telr, tknr integer :: ir, jr allocate(dmoa(nbf,nocca), dmob(nbf,noccb), mopa(nbf,nbf), mopb(nbf,nbf)) allocate(SxMOa(nbf,nbf), SxMOb(nbf,nbf)) allocate(frap(nbf2), frbp(nbf2), fram(nbf2), frbm(nbf2), dvx(nbf,nbf), hxc(nbf,nbf)) hxr = 1.0d-3 call mo_transform(sp(1)%mo, dSa(:,:,cc,kc), nbf, scr, tmp, SxMOa) call mo_transform(sp(2)%mo, dSa(:,:,cc,kc), nbf, scr, tmp, SxMOb) dmoa = 0.0_dp do ir = 1, nocca do jr = 1, nocca dmoa(:,ir) = dmoa(:,ir) - 0.5_dp*SxMOa(jr,ir)*sp(1)%mo(:,jr) end do end do dmob = 0.0_dp do ir = 1, noccb do jr = 1, noccb dmob(:,ir) = dmob(:,ir) - 0.5_dp*SxMOb(jr,ir)*sp(2)%mo(:,jr) end do end do basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) + hxr call basis%init_shell_centers() call dft_initialize(infos, basis, mgr) mopa = sp(1)%mo; mopa(:,1:nocca) = sp(1)%mo(:,1:nocca) + hxr*dmoa mopb = sp(2)%mo; mopb(:,1:noccb) = sp(2)%mo(:,1:noccb) + hxr*dmob frap = 0.0_dp; frbp = 0.0_dp call dftexcor(basis, mgr, int(infos%control%scftype), frap, frbp, mopa, mopb, & nbf, nbf2, exr, telr, tknr, infos) call dftclean(infos) basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) - 2*hxr call basis%init_shell_centers() call dft_initialize(infos, basis, mgr) mopa = sp(1)%mo; mopa(:,1:nocca) = sp(1)%mo(:,1:nocca) - hxr*dmoa mopb = sp(2)%mo; mopb(:,1:noccb) = sp(2)%mo(:,1:noccb) - hxr*dmob fram = 0.0_dp; frbm = 0.0_dp call dftexcor(basis, mgr, int(infos%control%scftype), fram, frbm, mopa, mopb, & nbf, nbf2, exr, telr, tknr, infos) call dftclean(infos) basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) + hxr call basis%init_shell_centers() call unpack_from_packed((frap - fram)/(2*hxr), dvx, nbf) call mo_transform(sp(1)%mo, dvx, nbf, scr, tmp, hxc) do a = 1, nvira do i = 1, nocca ia = (a-1)*nocca + i bvec(ia,icart) = bvec(ia,icart) - hxc(i,nocca+a) end do end do call unpack_from_packed((frbp - frbm)/(2*hxr), dvx, nbf) call mo_transform(sp(2)%mo, dvx, nbf, scr, tmp, hxc) do a = 1, nvirb do i = 1, noccb ia = (a-1)*noccb + i bvec(la+ia,icart) = bvec(la+ia,icart) - hxc(i,noccb+a) end do end do deallocate(dmoa, dmob, mopa, mopb, SxMOa, SxMOb, frap, frbp, fram, frbm, dvx, hxc) end block end if end do end do call cphf_solve_uhf(infos, ncart, bvec, uvec) ! ===== open-shell CPHF orbital-relaxation response ===== ! relaxed dC^s, spin density derivative dP^s do s = 1, 2 allocate(sp(s)%dCx(nbf, sp(s)%nocc, ncart), source=0.0_dp) allocate(sp(s)%dPx(nbf, nbf, ncart), source=0.0_dp) do yy = 1, ncart do a = 1, sp(s)%nvir do i = 1, sp(s)%nocc ia = (a-1)*sp(s)%nocc + i sp(s)%dCx(:,i,yy) = sp(s)%dCx(:,i,yy) & + sp(s)%mo(:,sp(s)%nocc+a)*uvec(sp(s)%loff+ia, yy) end do end do do i = 1, sp(s)%nocc do j = 1, sp(s)%nocc sp(s)%dCx(:,i,yy) = sp(s)%dCx(:,i,yy) - 0.5_dp*sp(s)%mo(:,j)*sp(s)%s1oo(j,i,yy) end do end do do i = 1, sp(s)%nocc do mu = 1, nbf do nu = 1, nbf sp(s)%dPx(mu,nu,yy) = sp(s)%dPx(mu,nu,yy) & + sp(s)%dCx(mu,i,yy)*sp(s)%mo(nu,i) + sp(s)%mo(mu,i)*sp(s)%dCx(nu,i,yy) end do end do end do end do allocate(sp(s)%gdpoo(sp(s)%nocc, sp(s)%nocc, ncart), source=0.0_dp) allocate(sp(s)%moe1 (sp(s)%nocc, sp(s)%nocc, ncart), source=0.0_dp) end do ! G^s[dP^y] (couples both spins via fock_jk) -> occ-occ MO block do yy = 1, ncart call pack_matrix(sp(1)%dPx(:,:,yy), dpck(:,1)) call pack_matrix(sp(2)%dPx(:,:,yy), dpck(:,2)) fpck = 0.0_dp call fock_jk(basis, d=dpck, f=fpck, scale_exch=hfscale, infos=infos) do s = 1, 2 call unpack_from_packed(fpck(:,s), gfull, nbf) call mo_transform(sp(s)%mo, gfull, nbf, scr, tmp, hxMO) sp(s)%gdpoo(:,:,yy) = hxMO(1:sp(s)%nocc,1:sp(s)%nocc) end do end do ! energy-weighted derivative occ-occ block (without the G[P]^y piece, which ! is folded into tGP via the Mi^x probe below) do s = 1, 2 do yy = 1, ncart do ll = 1, sp(s)%nocc do kk = 1, sp(s)%nocc sp(s)%moe1(kk,ll,yy) = sp(s)%hoo(kk,ll,yy) + sp(s)%gdpoo(kk,ll,yy) & - 0.5_dp*(sp(s)%eps(kk)+sp(s)%eps(ll))*sp(s)%s1oo(kk,ll,yy) end do end do end do end do ! 2e response traces, summed over spin: ! A2(x,y) = sum_s Tr[dP^s,y G^{s,x}[P]] ! tGP(x,y) = sum_s Tr[Mi^s,x G^{s,y}[P]], Mi^s,x = sum_kl s1oo^s,x_kl C^s_k C^s_l^T allocate(A2(ncart,ncart), tGP(ncart,ncart), Mi(nbf,nbf), source=0.0_dp) do yy = 1, ncart do s = 1, 2 gx = 0.0_dp call fock_deriv_contract_os(infos, basis, ptot, sp(s)%p, sp(s)%dPx(:,:,yy), hfscale, gx) A2(:,yy) = A2(:,yy) + reshape(gx, [ncart]) end do end do do x = 1, ncart do s = 1, 2 Mi = 0.0_dp do ll = 1, sp(s)%nocc do kk = 1, sp(s)%nocc do mu = 1, nbf do nu = 1, nbf Mi(mu,nu) = Mi(mu,nu) + sp(s)%s1oo(kk,ll,x)*sp(s)%mo(mu,kk)*sp(s)%mo(nu,ll) end do end do end do end do gx = 0.0_dp call fock_deriv_contract_os(infos, basis, ptot, sp(s)%p, Mi, hfscale, gx) tGP(x,:) = tGP(x,:) + reshape(gx, [ncart]) end do end do ! assemble H^resp_xy allocate(hresp(ncart,ncart), source=0.0_dp) do x = 1, ncart do yy = 1, ncart a1v = 0.0_dp; a3v = 0.0_dp; t3a = 0.0_dp do s = 1, 2 a1v = a1v + sum(sp(s)%dPx(:,:,yy)*hflat(:,:,x)) do i = 1, sp(s)%nocc dcsx = 0.0_dp do mu = 1, nbf do nu = 1, nbf dcsx = dcsx + sp(s)%dCx(mu,i,yy)*sflat(mu,nu,x)*sp(s)%mo(nu,i) end do end do a3v = a3v + sp(s)%eps(i)*dcsx end do do ll = 1, sp(s)%nocc do kk = 1, sp(s)%nocc t3a = t3a + sp(s)%s1oo(kk,ll,x)*sp(s)%moe1(kk,ll,yy) end do end do end do hresp(x,yy) = (a1v + A2(x,yy)) - 2.0_dp*a3v - t3a - tGP(x,yy) end do end do allocate(hess_native(ncart,ncart)) hess_native = 0.5_dp*(hresp + transpose(hresp)) ! --- DFT (UKS) exchange-correlation second-derivative contribution -------- ! Open-shell analog of the closed-shell RKS XC block: ! dHse : XC skeleton + density-response, central FD of the analytic ! open-shell XC gradient (derexc_blk) along R +/- h, P^s +/- h dP^s. ! dHt3 : -2 sum_s sum_kl s1oo^s,x_kl (vxc^s,y + fxc[dP^s,y])_kl, the XC part ! of the energy-weighted term, from the FD of the spin XC Fock ! (dftexcor) along R, C^s +/- h dC^s. ! The HF-exchange fraction is already in the Coulomb/exchange terms (hfscale). if (infos%control%hamilton >= 20) then block use dft, only: dft_initialize, dftclean, dftexcor use mod_dft_gridint_grad, only: derexc_blk use mod_dft_molgrid, only: dft_grid_t type(dft_grid_t) :: mg real(dp), allocatable :: dapa(:,:), dapb(:,:), dedp(:,:), dedm(:,:) real(dp), allocatable :: mopa(:,:), mopb(:,:), frap(:), frbp(:), fram(:), frbm(:) real(dp), allocatable :: dFxc(:,:), tmpn(:,:), dHse(:,:), dHt3(:,:), dFoo(:,:) real(dp) :: hx, exr, telr, tknr integer :: yy2, ccy, kcy, nang, x2, kk2, ll2, ss hx = 1.0d-3; nang = maxval(basis%am) + 2 allocate(dapa(nbf,nbf), dapb(nbf,nbf), dedp(3,natom), dedm(3,natom)) allocate(mopa(nbf,nbf), mopb(nbf,nbf), frap(nbf2), frbp(nbf2), fram(nbf2), frbm(nbf2)) allocate(dFxc(nbf,nbf), dHse(ncart,ncart), dHt3(ncart,ncart)) dHse = 0.0_dp; dHt3 = 0.0_dp call dft_initialize(infos, basis, mg); call dftclean(infos) ! warm-up do yy2 = 1, ncart ccy = mod(yy2-1,3)+1; kcy = (yy2-1)/3+1 basis%atoms%xyz(ccy,kcy) = basis%atoms%xyz(ccy,kcy) + hx call basis%init_shell_centers() call dft_initialize(infos, basis, mg) dapa = sp(1)%p + hx*sp(1)%dPx(:,:,yy2); dapb = sp(2)%p + hx*sp(2)%dPx(:,:,yy2) dedp = 0.0_dp call derexc_blk(basis, mg, dapa, dapb, dedp, telr, tknr, nang, nbf, & infos%dft%grid_density_cutoff, .true., infos) mopa = sp(1)%mo; mopa(:,1:nocca) = sp(1)%mo(:,1:nocca) + hx*sp(1)%dCx(:,:,yy2) mopb = sp(2)%mo; mopb(:,1:noccb) = sp(2)%mo(:,1:noccb) + hx*sp(2)%dCx(:,:,yy2) frap = 0.0_dp; frbp = 0.0_dp call dftexcor(basis, mg, int(infos%control%scftype), frap, frbp, mopa, mopb, & nbf, nbf2, exr, telr, tknr, infos) call dftclean(infos) basis%atoms%xyz(ccy,kcy) = basis%atoms%xyz(ccy,kcy) - 2*hx call basis%init_shell_centers() call dft_initialize(infos, basis, mg) dapa = sp(1)%p - hx*sp(1)%dPx(:,:,yy2); dapb = sp(2)%p - hx*sp(2)%dPx(:,:,yy2) dedm = 0.0_dp call derexc_blk(basis, mg, dapa, dapb, dedm, telr, tknr, nang, nbf, & infos%dft%grid_density_cutoff, .true., infos) mopa = sp(1)%mo; mopa(:,1:nocca) = sp(1)%mo(:,1:nocca) - hx*sp(1)%dCx(:,:,yy2) mopb = sp(2)%mo; mopb(:,1:noccb) = sp(2)%mo(:,1:noccb) - hx*sp(2)%dCx(:,:,yy2) fram = 0.0_dp; frbm = 0.0_dp call dftexcor(basis, mg, int(infos%control%scftype), fram, frbm, mopa, mopb, & nbf, nbf2, exr, telr, tknr, infos) call dftclean(infos) basis%atoms%xyz(ccy,kcy) = basis%atoms%xyz(ccy,kcy) + hx call basis%init_shell_centers() dHse(:,yy2) = reshape((dedp - dedm)/(2*hx), [ncart]) ! dHt3: -2 sum_s s1oo^s,x (dVxc^s,y + fxc[dP^s,y])_oo do ss = 1, 2 if (ss == 1) then call unpack_from_packed((frap - fram)/(2*hx), dFxc, nbf) else call unpack_from_packed((frbp - frbm)/(2*hx), dFxc, nbf) end if allocate(tmpn(nbf,sp(ss)%nocc), dFoo(sp(ss)%nocc,sp(ss)%nocc)) call dgemm('n','n', nbf, sp(ss)%nocc, nbf, 1.0_dp, dFxc, nbf, sp(ss)%mo, nbf, 0.0_dp, tmpn, nbf) call dgemm('t','n', sp(ss)%nocc, sp(ss)%nocc, nbf, 1.0_dp, sp(ss)%mo, nbf, tmpn, nbf, 0.0_dp, dFoo, sp(ss)%nocc) do x2 = 1, ncart do ll2 = 1, sp(ss)%nocc do kk2 = 1, sp(ss)%nocc dHt3(x2,yy2) = dHt3(x2,yy2) - 1.0_dp*sp(ss)%s1oo(kk2,ll2,x2)*dFoo(kk2,ll2) end do end do end do deallocate(tmpn, dFoo) end do end do hess_native = hess_native + 0.5_dp*(dHse + transpose(dHse)) & + 0.5_dp*(dHt3 + transpose(dHt3)) deallocate(dapa, dapb, dedp, dedm, mopa, mopb, frap, frbp, fram, frbm, dFxc, dHse, dHt3) end block end if ! nuclear repulsion call hess_nn(basis%atoms, basis%ecp_zn_num, hess_native) ! --- one-electron + Pulay second-derivative skeleton (fixed density) ------ block use grd1, only: eijden, hess_ee_overlap, hess_ee_kinetic, hess_en use ecp_tool, only: add_ecphess real(dp), allocatable :: wlag(:), pden(:), hcc(:,:) allocate(wlag(nbf2), pden(nbf2), hcc(ncart,ncart), source=0.0_dp) call eijden(wlag, nbf, infos) ! open-shell Lagrangian W pden = dma + dmb ! total density (Pa + Pb) call hess_ee_overlap(basis, wlag, hess_native) call hess_ee_kinetic(basis, pden, hess_native) call hess_en(basis, basis%atoms%xyz, & basis%atoms%zn - basis%ecp_zn_num, pden, hess_native, hess_cc=hcc) call add_ecphess(basis, basis%atoms%xyz, pden, hess_native) ! ECP skeleton (if any) deallocate(wlag, pden, hcc) end block ! --- two-electron (ERI) second-derivative skeleton (fixed density) -------- block use grd2, only: grd2_hess_driver, grd2_compute_data_t use hf_gradient_mod, only: grd2_uhf_compute_data_t type(grd2_uhf_compute_data_t) :: gcomp gcomp = grd2_uhf_compute_data_t( da = dma, db = dmb, hfscale = hfscale, nbf = nbf ) call gcomp%init() call gcomp%build_cart(basis) call grd2_hess_driver(infos, basis, hess_native, gcomp) call gcomp%clean() end block call infos%dat%alloc_or_die(OQP_hf_hessian, (/ ncart, ncart /), hess_store, & description='Native OpenQP open-shell (UHF) HF analytic Hessian matrix') hess_store = hess_native write(iw,'(A)') 'PyOQP: Native OpenQP open-shell (UHF) HF Hessian matrix stored' deallocate(ptot, dSa, dTa, dVa, sflat, hflat, bvec, uvec, scr, tmp, SxMO, hxMO, & probe, gx, d0a, d0b, gfull, Gd0, dpck, fpck, & A2, tGP, Mi, hresp, hess_native) end subroutine hf_hessian_uhf !############################################################################### subroutine hf_hessian_rohf(infos) ! Native open-shell (ROHF) analytic HF Hessian (HF only). ! ! ROHF uses a SINGLE MO set with a docc/socc/virt partition, so the orbital ! response is solved over the ROHF rotation space (cphf_solve_rohf) rather ! than the UHF spin blocks. The ROHF energy has the same functional form as ! UHF in terms of (Pa, Pb), so the Hessian decomposes identically into ! H = E_nn'' + skeleton(1e total density + open-shell W, 2e via grd2_uhf) ! + response(orbital relaxation), ! where the skeleton + nuclear repulsion are exactly hess_skel_open. ! ! The orbital-relaxation response is evaluated SEMI-NUMERICALLY, reusing the ! validated analytic open-shell gradient: with the relaxed orbital derivative ! dC^b (from the ROHF CPHF amplitudes) the response is the central finite ! difference, AT FIXED GEOMETRY, of the density/Lagrangian-dependent gradient ! along the orbital path C_occ +/- h dC^b: ! H^resp(:,b) = [ g(C + h dC^b) - g(C - h dC^b) ] / 2h , ! g(C') = grad_ee_overlap(W') + grad_ee_kinetic(P') + grad_en(P') ! + grad_2e(Pa', Pb') , W' = -(Pa' Fa' Pa' + Pb' Fb' Pb') , ! with Fa'/Fb' rebuilt from the perturbed densities (Hcore + fock_jk). This ! captures BOTH the relaxed-density and the energy-weighted (W) response ! through the gradient's own W build (eijden convention), so no ROHF-specific ! Lagrangian-derivative algebra is required. The CPHF right-hand side is the ! non-canonical Pulay form (orbital energies replaced by the full Fock occ-occ ! blocks), reducing to the validated UHF RHS in the canonical limit. use precision, only: dp use types, only: information use basis_tools, only: basis_set use oqp_tagarray_driver, only: tagarray_get_data, OQP_DM_A, OQP_DM_B, & OQP_VEC_MO_A, OQP_FOCK_A, OQP_FOCK_B, OQP_Hcore, OQP_hf_hessian, TA_TYPE_REAL64 use mathlib, only: unpack_matrix, pack_matrix, orthogonal_transform_sym use grd1, only: der_overlap_matrix, der_kinetic_matrix, der_nucattr_matrix, hess_nn, & grad_ee_overlap, grad_ee_kinetic, grad_en_hellman_feynman, grad_en_pulay use grd2, only: grd2_driver, grd2_compute_data_t use hf_gradient_mod, only: grd2_uhf_compute_data_t use fock_deriv_mod, only: fock_deriv_contract_os use scf_addons, only: fock_jk use cphf_mod, only: cphf_solve_rohf, rohf_pack_trial, rohf_unpack_trial use io_constants, only: iw use messages, only: show_message, WITH_ABORT implicit none type(information), target, intent(inout) :: infos type(basis_set), pointer :: basis real(dp), contiguous, pointer :: dma(:), dmb(:), mo(:,:), focka(:), fockb(:), hcore(:) real(dp), contiguous, pointer :: hess_store(:,:) real(dp), allocatable :: pa(:,:), pb(:,:), ptot(:,:) real(dp), allocatable :: dSa(:,:,:,:), dTa(:,:,:,:), dVa(:,:,:,:) real(dp), allocatable :: faMO(:,:), fbMO(:,:) real(dp), allocatable :: scr(:,:), tmp(:,:), SxMO(:,:), hxMO(:,:), probe(:,:) real(dp), allocatable :: ga2e(:,:,:), gb2e(:,:,:) real(dp), allocatable :: d0a(:,:), d0b(:,:), dpck(:,:), fpck(:,:), gfull(:,:), Gd0(:,:) real(dp), allocatable :: ba(:,:), bb(:,:), bvec(:,:), uvec(:,:) real(dp), allocatable :: xa(:,:), xb(:,:), dCa(:,:), dCb(:,:), gp(:,:), gm(:,:) real(dp), allocatable :: zneff(:), hess_native(:,:), hresp(:,:) real(dp), allocatable :: faop(:), fbop(:) integer, allocatable :: iecp_atom(:) real(dp) :: hfscale, hstep, gx(3, size(infos%atoms%xyz,2)) integer :: nbf, nbf2, natom, ncart, nocca, noccb, nvira, nvirb, offset, ltot integer :: i, j, a, icart, kc, cc, x, mu, nu, ie, nec basis => infos%basis basis%atoms => infos%atoms nbf = basis%nbf nbf2 = nbf*(nbf+1)/2 natom = size(basis%atoms%xyz, 2) ncart = 3*natom nocca = infos%mol_prop%nelec_A noccb = infos%mol_prop%nelec_B nvira = nbf - nocca nvirb = nbf - noccb offset = nocca - noccb ltot = noccb*(offset + nvira) + offset*nvira hfscale = 1.0_dp if (infos%control%hamilton >= 20) hfscale = infos%dft%hfscale hstep = 1.0d-3 write(iw,'(/,A)') 'PyOQP: Native OpenQP open-shell (ROHF) HF Hessian CPHF response prepass' write(iw,'(A,I6,A,I6,A,I6,A,I6,A,I6)') ' nbf=', nbf, ' nocca=', nocca, & ' noccb=', noccb, ' rhs=', ncart, ' rotdim=', ltot write(iw,'(A)') ' Storing native OpenQP open-shell (ROHF) HF analytic Hessian in OQP::hf_hessian.' if (ncart <= 0 .or. ltot <= 0) then write(iw,'(A)') ' ROHF CPHF prepass skipped: empty rotation/nuclear space.' return end if call tagarray_get_data(infos%dat, OQP_DM_A, dma) call tagarray_get_data(infos%dat, OQP_DM_B, dmb) call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo) call tagarray_get_data(infos%dat, OQP_FOCK_A, focka) call tagarray_get_data(infos%dat, OQP_FOCK_B, fockb) call tagarray_get_data(infos%dat, OQP_Hcore, hcore) allocate(pa(nbf,nbf), pb(nbf,nbf), ptot(nbf,nbf)) call unpack_matrix(dma, pa); call unpack_matrix(dmb, pb); ptot = pa + pb allocate(zneff(natom)); zneff = basis%atoms%zn - basis%ecp_zn_num ! Map each atom to its ECP-centre index in ecp_coord (which is sized ! 3*num_ecps, i.e. one (x,y,z) triple per ECP centre, NOT per atom). The ! semi-numerical resp_grad displaces atoms one Cartesian at a time and must ! move the matching ECP centre in lockstep; iecp_atom(kc)=0 means atom kc ! carries no ECP (its centre must not be touched). allocate(iecp_atom(natom)); iecp_atom = 0 if (basis%ecp_params%is_ecp) then nec = size(basis%ecp_params%n_expo) do ie = 1, nec do i = 1, natom if (all(abs(basis%ecp_params%ecp_coord(3*(ie-1)+1:3*ie) & - basis%atoms%xyz(:,i)) < 1.0e-6_dp)) then iecp_atom(i) = ie exit end if end do end do ! Every ECP centre must have been matched to an atom: an unmapped centre ! would stay fixed while its atom is displaced, silently corrupting the ! semi-numerical response. Abort loudly instead. if (count(iecp_atom > 0) /= nec) then call show_message('hf_hessian (ROHF): could not map every ECP centre '// & 'to an atom (coordinate mismatch > 1e-6 bohr); analytic Hessian '// & 'would be wrong - use [hess] type=numerical for this system.', WITH_ABORT) end if end if ! derivative integrals (normalized into the bfnrm/MO convention) allocate(dSa(nbf,nbf,3,natom), dTa(nbf,nbf,3,natom), dVa(nbf,nbf,3,natom)) call der_overlap_matrix(basis, dSa) call der_kinetic_matrix(basis, dTa) call der_nucattr_matrix(basis, basis%atoms%xyz, & basis%atoms%zn - basis%ecp_zn_num, dVa) ! ECP-screened point charge block integer :: kc2, cc2, mu2, nu2 do kc2 = 1, natom do cc2 = 1, 3 do nu2 = 1, nbf do mu2 = 1, nbf dSa(mu2,nu2,cc2,kc2) = dSa(mu2,nu2,cc2,kc2)*basis%bfnrm(mu2)*basis%bfnrm(nu2) dTa(mu2,nu2,cc2,kc2) = dTa(mu2,nu2,cc2,kc2)*basis%bfnrm(mu2)*basis%bfnrm(nu2) dVa(mu2,nu2,cc2,kc2) = dVa(mu2,nu2,cc2,kc2)*basis%bfnrm(mu2)*basis%bfnrm(nu2) end do end do end do end do end block ! ECP first-derivative integrals -> core-Hamiltonian derivative dHcore/dR, ! feeding the non-canonical CPHF RHS (hxMO below). The ECP skeleton + response ! is then completed by add_ecpder inside resp_grad (semi-numerical). Already ! normalized, so added after the bfnrm scaling. No-op for non-ECP bases. block use ecp_tool, only: ecp_deriv_ints real(dp), allocatable :: dVecp(:,:,:,:) allocate(dVecp(nbf,nbf,3,natom)) call ecp_deriv_ints(basis, basis%atoms%xyz, dVecp) dVa = dVa + dVecp deallocate(dVecp) end block ! occ-occ Fock blocks (MO) of the converged spin Fock matrices (non-canonical) allocate(scr(nbf,nbf), tmp(nbf,nbf), SxMO(nbf,nbf), hxMO(nbf,nbf)) allocate(faMO(nbf,nbf), fbMO(nbf,nbf)) call unpack_matrix(focka, scr); call mo_transform(mo, scr, nbf, tmp, hxMO, faMO) call unpack_matrix(fockb, scr); call mo_transform(mo, scr, nbf, tmp, hxMO, fbMO) ! 2e response-Fock skeleton G^{s,x}[P]_ai for all coordinates (per spin) allocate(ga2e(nvira,nocca,ncart), gb2e(nvirb,noccb,ncart), source=0.0_dp) allocate(probe(nbf,nbf)) do a = 1, nvira do i = 1, nocca do mu = 1, nbf do nu = 1, nbf probe(mu,nu) = 0.5_dp*( mo(mu,nocca+a)*mo(nu,i) + mo(mu,i)*mo(nu,nocca+a) ) end do end do gx = 0.0_dp call fock_deriv_contract_os(infos, basis, ptot, pa, probe, hfscale, gx) ga2e(a,i,:) = reshape(gx, [ncart]) end do end do do a = 1, nvirb do i = 1, noccb do mu = 1, nbf do nu = 1, nbf probe(mu,nu) = 0.5_dp*( mo(mu,noccb+a)*mo(nu,i) + mo(mu,i)*mo(nu,noccb+a) ) end do end do gx = 0.0_dp call fock_deriv_contract_os(infos, basis, ptot, pb, probe, hfscale, gx) gb2e(a,i,:) = reshape(gx, [ncart]) end do end do ! ===== CPHF right-hand sides (non-canonical Pulay form), packed ===== allocate(d0a(nbf,nbf), d0b(nbf,nbf), gfull(nbf,nbf), Gd0(nbf,nbf)) allocate(dpck(nbf2,2), fpck(nbf2,2)) allocate(ba(nvira,nocca), bb(nvirb,noccb)) allocate(bvec(ltot,ncart), uvec(ltot,ncart), source=0.0_dp) icart = 0 do kc = 1, natom do cc = 1, 3 icart = icart + 1 call mo_transform(mo, dSa(:,:,cc,kc), nbf, scr, tmp, SxMO) call mo_transform(mo, dTa(:,:,cc,kc)+dVa(:,:,cc,kc), nbf, scr, tmp, hxMO) ! reorthonormalization densities d0^s = -sum_ij S^x_ij C_i C_j (per spin occ) d0a = 0.0_dp; d0b = 0.0_dp do i = 1, nocca do j = 1, nocca do mu = 1, nbf do nu = 1, nbf d0a(mu,nu) = d0a(mu,nu) - SxMO(i,j)*mo(mu,i)*mo(nu,j) end do end do end do end do do i = 1, noccb do j = 1, noccb do mu = 1, nbf do nu = 1, nbf d0b(mu,nu) = d0b(mu,nu) - SxMO(i,j)*mo(mu,i)*mo(nu,j) end do end do end do end do call pack_matrix(d0a, dpck(:,1)); call pack_matrix(d0b, dpck(:,2)) fpck = 0.0_dp call fock_jk(basis, d=dpck, f=fpck, scale_exch=hfscale, infos=infos) ! Non-canonical Pulay RHS. The reorthonormalization Fock-coupling is the ! occupied-projected anticommutator of S^x and the spin Fock: ! B^s_ai = -(h^x + G2e + G[d0])_ai ! + sum_{j in occ} ( S^x_aj F^s_ji + F^s_aj S^x_ji ) . ! The first sum is the usual eps_i S^x_ai in the canonical (diagonal-Fock) ! limit; the second vanishes there (F^s_aj is a vir-occ Fock element) and ! supplies the non-canonical correction needed for the socc rotations. call unpack_from_packed(fpck(:,1), gfull, nbf) call mo_transform(mo, gfull, nbf, scr, tmp, Gd0) do i = 1, nocca do a = 1, nvira ba(a,i) = -(hxMO(i,nocca+a) + ga2e(a,i,icart) + Gd0(i,nocca+a)) & + dot_product(SxMO(nocca+a,1:nocca), faMO(1:nocca,i)) & + dot_product(faMO(nocca+a,1:nocca), SxMO(1:nocca,i)) end do end do ! beta block call unpack_from_packed(fpck(:,2), gfull, nbf) call mo_transform(mo, gfull, nbf, scr, tmp, Gd0) do i = 1, noccb do a = 1, nvirb bb(a,i) = -(hxMO(i,noccb+a) + gb2e(a,i,icart) + Gd0(i,noccb+a)) & + dot_product(SxMO(noccb+a,1:noccb), fbMO(1:noccb,i)) & + dot_product(fbMO(noccb+a,1:noccb), SxMO(1:noccb,i)) end do end do ! --- XC contribution to the CPKS right-hand side (ROKS only) ----------- ! Central FD of the spin XC Fock (open-shell dftexcor) along R +/- h AND ! occupied MOs reorthonormalized by dmoR^s = -1/2 sum_j C_j S^x_ji; the XC ! skeleton dVxc/dR + f_xc[d0], subtracted from B (which carries -F0x). if (infos%control%hamilton >= 20) then block use dft, only: dft_initialize, dftclean, dftexcor use mod_dft_molgrid, only: dft_grid_t type(dft_grid_t) :: mgr real(dp), allocatable :: dmoa(:,:), dmob(:,:), mopa(:,:), mopb(:,:) real(dp), allocatable :: frap(:), frbp(:), fram(:), frbm(:), dvx(:,:), hxc(:,:) real(dp) :: hxr, exr, telr, tknr integer :: ir, jr allocate(dmoa(nbf,nocca), dmob(nbf,noccb), mopa(nbf,nbf), mopb(nbf,nbf)) allocate(frap(nbf2), frbp(nbf2), fram(nbf2), frbm(nbf2), dvx(nbf,nbf), hxc(nbf,nbf)) hxr = 1.0d-3 dmoa = 0.0_dp do ir = 1, nocca do jr = 1, nocca dmoa(:,ir) = dmoa(:,ir) - 0.5_dp*SxMO(jr,ir)*mo(:,jr) end do end do dmob = 0.0_dp do ir = 1, noccb do jr = 1, noccb dmob(:,ir) = dmob(:,ir) - 0.5_dp*SxMO(jr,ir)*mo(:,jr) end do end do basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) + hxr call basis%init_shell_centers() call dft_initialize(infos, basis, mgr) mopa = mo; mopa(:,1:nocca) = mo(:,1:nocca) + hxr*dmoa mopb = mo; mopb(:,1:noccb) = mo(:,1:noccb) + hxr*dmob frap = 0.0_dp; frbp = 0.0_dp call dftexcor(basis, mgr, int(infos%control%scftype), frap, frbp, mopa, mopb, & nbf, nbf2, exr, telr, tknr, infos) call dftclean(infos) basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) - 2*hxr call basis%init_shell_centers() call dft_initialize(infos, basis, mgr) mopa = mo; mopa(:,1:nocca) = mo(:,1:nocca) - hxr*dmoa mopb = mo; mopb(:,1:noccb) = mo(:,1:noccb) - hxr*dmob fram = 0.0_dp; frbm = 0.0_dp call dftexcor(basis, mgr, int(infos%control%scftype), fram, frbm, mopa, mopb, & nbf, nbf2, exr, telr, tknr, infos) call dftclean(infos) basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) + hxr call basis%init_shell_centers() call unpack_from_packed((frap - fram)/(2*hxr), dvx, nbf) call mo_transform(mo, dvx, nbf, scr, tmp, hxc) do i = 1, nocca do a = 1, nvira ba(a,i) = ba(a,i) - hxc(i,nocca+a) end do end do call unpack_from_packed((frbp - frbm)/(2*hxr), dvx, nbf) call mo_transform(mo, dvx, nbf, scr, tmp, hxc) do i = 1, noccb do a = 1, nvirb bb(a,i) = bb(a,i) - hxc(i,noccb+a) end do end do deallocate(dmoa, dmob, mopa, mopb, frap, frbp, fram, frbm, dvx, hxc) end block end if call rohf_pack_trial(bvec(:,icart), ba, bb, nbf, nocca, noccb) end do end do call cphf_solve_rohf(infos, ncart, bvec, uvec) ! ===== semi-numerical orbital-relaxation response ===== ! Build the relaxed alpha/beta orbital derivatives independently, UHF-style: ! dCa_i = sum_a C^{vir_a}_a xa(a,i) - 1/2 sum_{j in docc+socc} S^x_ji C_j ! dCb_i = sum_a C^{vir_b}_a xb(a,i) - 1/2 sum_{j in docc} S^x_ji C_j ! The socc-docc rotation lives in xb (socc is beta-virtual), so it relaxes Pb ! and leaves Pa invariant (it is an alpha occ-occ rotation) -- exactly the ! ROHF physics, with no socc-docc cross term needed in dCa. allocate(xa(nvira,nocca), xb(nvirb,noccb), dCa(nbf,nocca), dCb(nbf,noccb)) allocate(gp(3,natom), gm(3,natom), hresp(ncart,ncart), source=0.0_dp) allocate(faop(nbf2), fbop(nbf2)) if (infos%control%hamilton >= 20) then ! flush stale grid state from the CPHF solver block use dft, only: dft_initialize, dftclean use mod_dft_molgrid, only: dft_grid_t type(dft_grid_t) :: mgw call dft_initialize(infos, basis, mgw); call dftclean(infos) end block end if do x = 1, ncart cc = mod(x-1,3)+1; kc = (x-1)/3+1 call rohf_unpack_trial(uvec(:,x), xa, xb, nbf, nocca, noccb) call mo_transform(mo, dSa(:,:,cc,kc), nbf, scr, tmp, SxMO) dCa = 0.0_dp do i = 1, nocca do a = 1, nvira dCa(:,i) = dCa(:,i) + mo(:,nocca+a)*xa(a,i) end do do j = 1, nocca dCa(:,i) = dCa(:,i) - 0.5_dp*SxMO(j,i)*mo(:,j) end do end do dCb = 0.0_dp do i = 1, noccb do a = 1, nvirb dCb(:,i) = dCb(:,i) + mo(:,noccb+a)*xb(a,i) end do do j = 1, noccb dCb(:,i) = dCb(:,i) - 0.5_dp*SxMO(j,i)*mo(:,j) end do end do call resp_grad( 1.0_dp, gp) call resp_grad(-1.0_dp, gm) hresp(:,x) = reshape((gp - gm)/(2.0_dp*hstep), [ncart]) end do ! The central difference of the ELECTRONIC gradient over geometry AND the ! relaxed orbital path already contains the full electronic Hessian (skeleton ! + orbital-relaxation response); only the (orbital-independent) nuclear ! repulsion second derivative is added analytically. allocate(hess_native(ncart,ncart)) hess_native = 0.5_dp*(hresp + transpose(hresp)) call hess_nn(basis%atoms, basis%ecp_zn_num, hess_native) call infos%dat%alloc_or_die(OQP_hf_hessian, (/ ncart, ncart /), hess_store, & description='Native OpenQP open-shell (ROHF) HF analytic Hessian matrix') hess_store = hess_native write(iw,'(A)') 'PyOQP: Native OpenQP open-shell (ROHF) HF Hessian matrix stored' deallocate(pa, pb, ptot, dSa, dTa, dVa, faMO, fbMO, scr, tmp, & SxMO, hxMO, probe, ga2e, gb2e, d0a, d0b, dpck, fpck, gfull, Gd0, & ba, bb, bvec, uvec, xa, xb, dCa, dCb, gp, gm, hresp, zneff, & hess_native, faop, fbop, iecp_atom) contains !> Electronic gradient (1e + 2e + Pulay-W; NO nuclear repulsion) at the !> geometry displaced by sgn*hstep in coordinate (cc,kc) AND the alpha-occ !> MOs displaced by sgn*hstep*dC (host-associated cc,kc,dC,hstep). Central !> differencing over sgn therefore captures the electronic skeleton AND the !> orbital-relaxation response together: the one-electron Hamiltonian, all !> gradient integrals, the densities Pa'/Pb' and the energy-weighted density !> W' = -(Pa' Fa' Pa' + Pb' Fb' Pb') (eijden convention, Fock rebuilt as !> Hcore' + fock_jk) are all evaluated at the displaced point, so no !> ROHF-specific Lagrangian-derivative algebra is needed. subroutine resp_grad(sgn, gout) use int1, only: omp_hst real(dp), intent(in) :: sgn real(dp), intent(out) :: gout(:,:) real(dp), allocatable :: cocc(:,:), pap(:,:), pbp(:,:) real(dp), allocatable, target :: paP_tri(:), pbP_tri(:) real(dp), allocatable :: ptP_tri(:), wlag(:), ta(:), hc(:), sm(:), tm(:) real(dp) :: tol integer :: ii, ij type(grd2_uhf_compute_data_t) :: gc allocate(cocc(nbf,nocca), pap(nbf,nbf), pbp(nbf,nbf)) allocate(paP_tri(nbf2), pbP_tri(nbf2), ptP_tri(nbf2), wlag(nbf2), ta(nbf2)) allocate(hc(nbf2), sm(nbf2), tm(nbf2)) ! displace geometry and rebuild the one-electron Hamiltonian there. The ECP ! center (ecp_coord) is a separate array from atoms%xyz, so it must be moved ! in lockstep or the displaced add_ecpint/add_ecpder would see the basis and ! the ECP at mismatched centers (catastrophic for the ECP atom). basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) + sgn*hstep if (iecp_atom(kc) > 0) & basis%ecp_params%ecp_coord(3*(iecp_atom(kc)-1)+cc) = & basis%ecp_params%ecp_coord(3*(iecp_atom(kc)-1)+cc) + sgn*hstep call basis%init_shell_centers() tol = log(10.0d0)*20.0_dp call omp_hst(basis, basis%atoms%xyz, basis%atoms%zn - basis%ecp_zn_num, & hc, sm, tm, logtol=tol, comm=infos%mpiinfo%comm, usempi=infos%mpiinfo%usempi) ! NB: the ECP one-electron potential is deliberately NOT added to hc here. ! The full ECP gradient (operator + basis-centre/Pulay derivatives) is the ! analytic add_ecpder below; folding the ECP into the spin Fock used to build ! the energy-weighted density W' would double-count its Pulay contribution ! (verified: doing so gives ~1.5e-2 vs the numerical Hessian, omitting it ! gives ~2e-5). ! relaxed orbitals -> perturbed densities and spin Fock matrices cocc(:,1:nocca) = mo(:,1:nocca) + sgn*hstep*dCa call dgemm('n','t', nbf, nbf, nocca, 1.0_dp, cocc, nbf, cocc, nbf, 0.0_dp, pap, nbf) cocc(:,1:noccb) = mo(:,1:noccb) + sgn*hstep*dCb call dgemm('n','t', nbf, nbf, noccb, 1.0_dp, cocc, nbf, cocc, nbf, 0.0_dp, pbp, nbf) call pack_matrix(pap, paP_tri); call pack_matrix(pbp, pbP_tri) ptP_tri = paP_tri + pbP_tri dpck(:,1) = paP_tri; dpck(:,2) = pbP_tri fpck = 0.0_dp call fock_jk(basis, d=dpck, f=fpck, scale_exch=hfscale, infos=infos) faop = hc + fpck(:,1); fbop = hc + fpck(:,2) gout = 0.0_dp ! DFT (ROKS): add the XC potential to the spin Focks (so W' is the full KS ! energy-weighted density) and the explicit open-shell XC gradient to gout. ! Both are evaluated at the displaced geometry with the relaxed orbitals, so ! the geometry+orbital FD gives the full KS Hessian (Pulay/W XC + explicit XC) ! with no separate analytic XC term. if (infos%control%hamilton >= 20) then block use dft, only: dft_initialize, dftclean, dftexcor use mod_dft_gridint_grad, only: derexc_blk use mod_dft_molgrid, only: dft_grid_t type(dft_grid_t) :: mg real(dp), allocatable :: mopa(:,:), mopb(:,:), fra(:), frb(:), dedft(:,:) real(dp) :: exr, telr, tknr integer :: nang allocate(mopa(nbf,nbf), mopb(nbf,nbf), fra(nbf2), frb(nbf2), dedft(3,natom)) nang = maxval(basis%am) + 2 call dft_initialize(infos, basis, mg) mopa = mo; mopa(:,1:nocca) = mo(:,1:nocca) + sgn*hstep*dCa mopb = mo; mopb(:,1:noccb) = mo(:,1:noccb) + sgn*hstep*dCb fra = 0.0_dp; frb = 0.0_dp call dftexcor(basis, mg, int(infos%control%scftype), fra, frb, mopa, mopb, & nbf, nbf2, exr, telr, tknr, infos) faop = faop + fra; fbop = fbop + frb dedft = 0.0_dp call derexc_blk(basis, mg, pap, pbp, dedft, telr, tknr, nang, nbf, & infos%dft%grid_density_cutoff, .true., infos) call dftclean(infos) gout = gout + dedft deallocate(mopa, mopb, fra, frb, dedft) end block end if call orthogonal_transform_sym(nbf, nbf, faop, pap, nbf, ta) call orthogonal_transform_sym(nbf, nbf, fbop, pbp, nbf, wlag) wlag = -wlag - ta ij = 0 do ii = 1, nbf ij = ij + ii wlag(ij) = 0.5_dp*wlag(ij) end do call grad_ee_overlap(basis, wlag, gout) call grad_ee_kinetic(basis, ptP_tri, gout) call grad_en_hellman_feynman(basis, basis%atoms%xyz, zneff, ptP_tri, gout) call grad_en_pulay(basis, basis%atoms%xyz, zneff, ptP_tri, gout) ! ECP gradient at the displaced geometry/density: central FD over the ! geometry+orbital path then yields BOTH the ECP skeleton second derivative ! and the ECP orbital-relaxation response. No-op for non-ECP bases. block use ecp_tool, only: add_ecpder call add_ecpder(basis, basis%atoms%xyz, ptP_tri, gout) end block gc = grd2_uhf_compute_data_t( da = paP_tri, db = pbP_tri, hfscale = hfscale, nbf = nbf ) call gc%init() call gc%build_cart(basis) call grd2_driver(infos, basis, gout, gc) call gc%clean() ! restore geometry (and the ECP center moved above) basis%atoms%xyz(cc,kc) = basis%atoms%xyz(cc,kc) - sgn*hstep if (iecp_atom(kc) > 0) & basis%ecp_params%ecp_coord(3*(iecp_atom(kc)-1)+cc) = & basis%ecp_params%ecp_coord(3*(iecp_atom(kc)-1)+cc) - sgn*hstep call basis%init_shell_centers() deallocate(cocc, pap, pbp, paP_tri, pbP_tri, ptP_tri, wlag, ta, hc, sm, tm) end subroutine resp_grad end subroutine hf_hessian_rohf !############################################################################### subroutine mo_transform(c_mo, a_ao, n, s1, s2, b_mo) use precision, only: dp real(kind=dp), intent(in) :: c_mo(:,:), a_ao(:,:) integer, intent(in) :: n real(kind=dp), intent(inout) :: s1(:,:), s2(:,:), b_mo(:,:) call dgemm('t','n', n, n, n, 1.0_dp, c_mo, n, a_ao, n, 0.0_dp, s1, n) call dgemm('n','n', n, n, n, 1.0_dp, s1, n, c_mo, n, 0.0_dp, b_mo, n) end subroutine mo_transform !############################################################################### subroutine unpack_from_packed(gpk, gfu, n) use precision, only: dp real(kind=dp), intent(in) :: gpk(:) real(kind=dp), intent(inout) :: gfu(:,:) integer, intent(in) :: n integer :: ii, jj, ij ij = 0 do ii = 1, n do jj = 1, ii ij = ij + 1 gfu(ii,jj) = gpk(ij); gfu(jj,ii) = gpk(ij) end do end do end subroutine unpack_from_packed end module hf_hessian_mod