!> @brief OpenQP <-> ddX PCM reaction-field bridge for the SCF energy path. !> !> @details This module is the Fortran half of the energy-only PCM seam. It !> declares the iso_c_binding interfaces to the two production C adapter !> entry points (source/solvent_ddx_adapter.c) and orchestrates the closed !> reaction-field loop inside a single SCF Fock build: !> !> D -> phi_cav -> ddX q_cav -> V_pcm -> Fock / E_pcm !> !> The C adapter returns status 2 when OpenQP was built without OQP_ENABLE_DDX, !> so a PCM-enabled run on a non-ddX build aborts here with a clear message !> rather than silently producing a vacuum result. !> !> SOURCE CONSISTENCY (ddX forward Phi vs adjoint Psi): !> * Phi (forward solve RHS): the EXACT total solute potential phi_cav at the !> cavity points, built from the full AO density (electrostatic_potential_ !> unweighted) plus the analytic nuclear term. !> * Psi (adjoint solve source): a FULL-DENSITY source. Atom-centered real !> solid-harmonic multipoles M_lm are accumulated for l = 0..PCM_PSI_LMAX !> (=8) from the full AO density by numerical quadrature over a dedicated !> source-projection molecular grid that reproduces the reference !> ddCOSMO/ddPCM density partition: per-atom (PARENT-ATOM) point !> assignment with Becke-original (3-iteration) fuzzy-cell weights and !> Treutler-Ahlrichs sqrt(R_i/R_j) atomic-size shifting over the Becke !> Bragg-Slater table (H = 0.35 A), WITH the literature outside-sphere !> leak continuation q*rsph^(2l+1)/r^(l+1) for r>rsph !> (build_full_density_multipoles / pcm_grid_update). The moments are in !> the exact ddX harmonic convention -- the real-solid-harmonic basis is !> evaluated by ddX's OWN routines (use ddx_harmonics: ylmscale, ylmbas) !> so ddX stays an external, dynamically-linked dependency and no harmonic !> code is vendored into OpenQP; only the interior/exterior leak !> bookkeeping in pcm_accumulate_leak is OpenQP's. They are then mapped !> to psi by the ddX rule psi(lm,isph)=4*pi/((2l+1) rsph^l) M_lm(isph) !> using the production cavity radii (oqp_ddx_pcm_radii). The production !> q_cav is the ddX adjoint charge from oqp_ddx_pcm_solve_psi(psi, !> phi_cav): both the forward RHS and the adjoint source are full-density !> quantities. This is recorded by "PCM diag pcm_source_mode=full_density_ !> multipoles_lmax8_exact_phi" and "PCM diag psi_source=full_density_grid_ !> multipoles_lmax8_becke3_treutler_parent_atom_leak". !> * NOTE: the per-sphere moments are partition-defined integrals !> (Becke-original/Treutler cells), the same convention the reference !> ddPCM implementations project on their per-atom Becke grids; the two !> codes agree in the fine-grid limit, with only quadrature-mesh !> differences remaining (it is NOT claimed to be bit-identical to !> PySCF's grid-projected psi). !> * The legacy l<=2 atom-centered Mulliken multipole solve (Phi AND Psi from !> the l<=2 source) is still run as a DIAGNOSTIC only, to report the !> source-vs-exact phi residual and the q_cav shift between the old l<=2 psi !> and the new full-density psi (PCM diag q_cav_*_vs_*_rms). !> !> VALIDATED SCALAR CONVENTIONS (analytic Born-ion/ddX oracle gate): !> * phi_cav sign: phi_total = sum_k Z_k/|r-R_k| + phi_elec !> * q_cav sign/scale: ddX cavity-projected adjoint charge (ddx_get_xi) used !> directly as the external-charge vector for external_charge_potential. !> * E_pcm: -0.5 * dot_product(phi_cav, q_cav), with NO additional dielectric !> factor: ddX folds the full dielectric response into its ddPCM R_eps !> operators, so -0.5*<phi_cav, q_cav> = ddx_pcm_energy = the PHYSICAL !> solvation free energy. Proven by the Born-ion oracle (point charge q !> centered in a single sphere of radius R): -0.5*<phi,q_cav> reproduces !> -(1/2)(1-1/eps)*q^2/R to machine precision at eps = 78.3553 and eps = 2. !> An extra f(eps) = (eps-1)/eps here (as in PySCF's solvent.ddpcm) would !> double-count the dielectric scaling, by -1.3% at eps=78 and -50% at eps=2. !> The single canonical runtime path and these conventions are pinned by !> tests/test_pcm_canonical_runtime_path.py. module solvent_pcm use precision, only: dp use, intrinsic :: ieee_arithmetic, only: ieee_is_finite use iso_c_binding, only: c_int, c_double, c_char, c_null_char, & c_int64_t, c_bool use types, only: information use basis_tools, only: basis_set use messages, only: show_message, with_abort use io_constants, only: iw use mathlib, only: traceprod_sym_packed use oqp_tagarray_driver, only: OQP_SM, data_has_tags, tagarray_get_data use int1, only: electrostatic_potential_unweighted, external_charge_potential, & multipole_integrals use mod_dft_molgrid, only: dft_grid_t use mod_dft_gridint, only: xc_engine_t, xc_consumer_t, xc_options_t, run_grid_aos use mod_dft_partfunc, only: PTYPE_BECKE3 use dft, only: dft_initialize ! ddX's own harmonic routines (external, dynamically-linked LGPL library): ! ylmscale builds the real-spherical-harmonic scaling factors and ylmbas ! evaluates the normalised real solid harmonics in ddX's exact convention. ! Using them directly (rather than vendoring a copy) keeps ddX external and ! guarantees the per-l normalisation/ordering matches the model ddX builds. ! Only available when OpenQP is built with ddX (OQP_ENABLE_DDX); the whole ! full-density source projection is unreachable otherwise (the C adapter ! returns status 2 and add_pcm_reaction_field aborts before it is called). #ifdef OQP_ENABLE_DDX use ddx_harmonics, only: ylmscale, ylmbas #endif implicit none private public :: add_pcm_reaction_field ! Maximum angular momentum of the full-density adjoint source Psi. The ddX ! model is built with lmax = 8 (solvent_ddx_adapter.c::build_pcm_model), so ! nbasis = (PCM_PSI_LMAX+1)^2 = 81 must match ddx_get_n_basis(). integer, parameter :: PCM_PSI_LMAX = 8 ! Maximum Lebedev cavity points per atom (must match n_lebedev used when the ! ddX model is built in build_pcm_model(); used only to size the receive ! buffer for the cavity coordinates). integer(c_int), parameter :: MAX_CAV_PER_ATOM = 302 integer, parameter :: PCM_FD_MAX_SAMPLES = 3 real(dp), parameter :: PCM_FD_STEP = 1.0e-4_dp ! Finite-difference diagnostics below derive this sign/scale. It is kept as ! an explicit constant so the Fock convention is guarded by runtime evidence. real(dp), parameter :: PCM_QCAV_TO_FOCK_SCALE = -0.5_dp ! Grid consumer for the PCM full-density-Psi production path. It integrates ! the electronic density on a dedicated Becke-original/Treutler-shifted ! molecular grid, assigns each weighted point to its PARENT atom (the atom ! whose atomic grid generated the slice, xce%currAtom), and accumulates the ! corresponding negative electronic charge into ddX-convention real-solid- ! harmonic multipoles through PCM_PSI_LMAX. Nuclear monopoles are added by ! the driver after the grid loop. This reproduces the reference ! ddCOSMO/ddPCM source projection (per-atom Becke-partitioned moments) up ! to quadrature-mesh differences. type, extends(xc_consumer_t) :: pcm_psi_grid_consumer_t integer :: lmax = 0 integer :: nbasis = 0 integer :: natom = 0 real(dp), pointer :: xyz(:,:) => null() real(dp), allocatable :: vscales(:) real(dp), allocatable :: radii(:) real(dp), allocatable :: multipoles(:,:,:) contains procedure :: parallel_start => pcm_grid_parallel_start procedure :: parallel_stop => pcm_grid_parallel_stop procedure :: update => pcm_grid_update procedure :: postUpdate => pcm_grid_post_update procedure :: clean => pcm_grid_clean end type pcm_psi_grid_consumer_t interface integer(c_int) function oqp_ddx_pcm_cavity(natom, xyz_bohr, charges, & eps, max_cav, ncav_out, cav_xyz_out, message, message_len) & bind(C, name="oqp_ddx_pcm_cavity") import :: c_int, c_double, c_char integer(c_int), value :: natom real(c_double), intent(in) :: xyz_bohr(*) real(c_double), intent(in) :: charges(*) real(c_double), value :: eps integer(c_int), value :: max_cav integer(c_int), intent(out) :: ncav_out real(c_double), intent(out) :: cav_xyz_out(*) character(kind=c_char), intent(out) :: message(*) integer(c_int), value :: message_len end function oqp_ddx_pcm_cavity integer(c_int) function oqp_ddx_pcm_solve(natom, xyz_bohr, charges, & eps, ncav, phi_cav, q_cav_out, esolv_out, message, message_len) & bind(C, name="oqp_ddx_pcm_solve") import :: c_int, c_double, c_char integer(c_int), value :: natom real(c_double), intent(in) :: xyz_bohr(*) real(c_double), intent(in) :: charges(*) real(c_double), value :: eps integer(c_int), value :: ncav real(c_double), intent(in) :: phi_cav(*) real(c_double), intent(out) :: q_cav_out(*) real(c_double), intent(out) :: esolv_out character(kind=c_char), intent(out) :: message(*) integer(c_int), value :: message_len end function oqp_ddx_pcm_solve integer(c_int) function oqp_ddx_pcm_solve_multipole_source(natom, & xyz_bohr, cavity_charges, nmultipoles, source_multipoles, eps, ncav, & phi_source_out, q_cav_out, esolv_out, message, message_len) & bind(C, name="oqp_ddx_pcm_solve_multipole_source") import :: c_int, c_double, c_char integer(c_int), value :: natom real(c_double), intent(in) :: xyz_bohr(*) real(c_double), intent(in) :: cavity_charges(*) integer(c_int), value :: nmultipoles real(c_double), intent(in) :: source_multipoles(*) real(c_double), value :: eps integer(c_int), value :: ncav real(c_double), intent(out) :: phi_source_out(*) real(c_double), intent(out) :: q_cav_out(*) real(c_double), intent(out) :: esolv_out character(kind=c_char), intent(out) :: message(*) integer(c_int), value :: message_len end function oqp_ddx_pcm_solve_multipole_source integer(c_int) function oqp_ddx_pcm_solve_multipole_source_with_phi(natom, & xyz_bohr, cavity_charges, nmultipoles, source_multipoles, eps, ncav, & phi_cav, q_cav_out, esolv_out, message, message_len) & bind(C, name="oqp_ddx_pcm_solve_multipole_source_with_phi") import :: c_int, c_double, c_char integer(c_int), value :: natom real(c_double), intent(in) :: xyz_bohr(*) real(c_double), intent(in) :: cavity_charges(*) integer(c_int), value :: nmultipoles real(c_double), intent(in) :: source_multipoles(*) real(c_double), value :: eps integer(c_int), value :: ncav real(c_double), intent(in) :: phi_cav(*) real(c_double), intent(out) :: q_cav_out(*) real(c_double), intent(out) :: esolv_out character(kind=c_char), intent(out) :: message(*) integer(c_int), value :: message_len end function oqp_ddx_pcm_solve_multipole_source_with_phi integer(c_int) function oqp_ddx_pcm_radii(natom, charges, radii_bohr_out, & message, message_len) bind(C, name="oqp_ddx_pcm_radii") import :: c_int, c_double, c_char integer(c_int), value :: natom real(c_double), intent(in) :: charges(*) real(c_double), intent(out) :: radii_bohr_out(*) character(kind=c_char), intent(out) :: message(*) integer(c_int), value :: message_len end function oqp_ddx_pcm_radii integer(c_int) function oqp_ddx_pcm_solve_psi(natom, xyz_bohr, charges, & eps, ncav, nbasis, psi, phi_cav, q_cav_out, esolv_out, message, & message_len) bind(C, name="oqp_ddx_pcm_solve_psi") import :: c_int, c_double, c_char integer(c_int), value :: natom real(c_double), intent(in) :: xyz_bohr(*) real(c_double), intent(in) :: charges(*) real(c_double), value :: eps integer(c_int), value :: ncav integer(c_int), value :: nbasis real(c_double), intent(in) :: psi(*) real(c_double), intent(in) :: phi_cav(*) real(c_double), intent(out) :: q_cav_out(*) real(c_double), intent(out) :: esolv_out character(kind=c_char), intent(out) :: message(*) integer(c_int), value :: message_len end function oqp_ddx_pcm_solve_psi end interface contains !> @brief Add the ddX PCM reaction-field operator to the Fock matrices and !> return the (provisional) PCM energy contribution. !> @param[in] basis AO basis (read-only) !> @param[in] infos run information; uses control%pcm_epsilon, atoms, natom !> @param[in] d packed AO density blocks (nbf_tri, nfocks) !> @param[in] nfocks number of spin blocks !> @param[inout] f packed AO Fock blocks (nbf_tri, nfocks); V_pcm added !> @param[out] e_pcm PCM energy contribution (provisional) subroutine add_pcm_reaction_field(basis, infos, d, nfocks, f, e_pcm) type(basis_set), intent(in) :: basis type(information), intent(inout) :: infos real(dp), intent(in) :: d(:,:) integer, intent(in) :: nfocks real(dp), intent(inout) :: f(:,:) real(dp), intent(out) :: e_pcm integer(c_int) :: natom, ncav, max_cav, nmultipoles, nbasis integer :: nbf_tri, ii, iat, icav, rc, ndelta real(dp) :: eps, f_epsilon, esolv, esolv_source, esolv_l2, phin, dx, dy, dz, r real(dp) :: half_tr_dv, q_cav_sum, q_cav_absnorm, phi_cav_sum, phi_cav_min, phi_cav_max real(dp) :: source_charge_sum, phi_source_delta_rms, phi_source_delta_max real(dp) :: q_cav_shift_rms, q_cav_full_vs_l2_rms, psi_full_norm, mult_full_norm real(dp) :: fd_fock_scale_mean, fd_fock_scale_rms, fd_fock_scale_maxerr integer :: fd_fock_samples logical :: pcm_diag real(dp), allocatable :: xyz(:,:), charges(:) real(dp), allocatable :: cav_xyz(:), cx(:), cy(:), cz(:) real(dp), allocatable :: phi_elec(:), phi_cav(:), phi_source(:), q_cav(:), q_cav_source(:) real(dp), allocatable :: q_cav_l2(:) real(dp), allocatable :: dtot(:), vpcm(:) real(dp), allocatable :: ao_pop(:), atom_pop(:), source_charges(:) real(dp), allocatable :: ao_dip(:,:), atom_dip(:,:), ao_quad(:,:), atom_quad(:,:) real(dp), allocatable :: source_multipoles(:,:) real(dp), allocatable :: multipoles_full(:,:), psi_full(:,:), radii(:) real(dp), contiguous, pointer :: smat(:) character(kind=c_char) :: cmsg(256) character(len=256) :: fmsg character(len=*), parameter :: tags_overlap(1) = (/ character(len=80) :: OQP_SM /) e_pcm = 0.0_dp natom = int(infos%mol_prop%natom, c_int) eps = infos%control%pcm_epsilon ! DIAGNOSTIC-ONLY dielectric factor f(eps) = (eps-1)/eps. It is reported in ! the diagnostics below but is NOT applied to e_pcm or the Fock operator: ! ddX folds the COMPLETE dielectric response into its ddPCM R_eps operators, ! so its pcm_energy = 0.5*<xs,psi> (= -0.5*<phi_cav,q_cav> by the adjoint ! identity) is already the physical solvation free energy. This is proven by ! the Born-ion oracle: -0.5*<phi,q_cav> = -(1/2)(1-1/eps)q^2/R to machine ! precision. PySCF's solvent.ddpcm applies an extra 0.5*f_eps*<psi,Xvec> ! scaling on top of its R_eps solve, which is why it FAILS the same Born ! oracle; it must not be imitated here. f_epsilon = (eps - 1.0_dp) / eps nbf_tri = size(d, 1) ! Convention self-diagnostics (two extra baseline ddX solves, a finite- ! difference Fock-scale probe, and a ~40-line per-cycle log block) are gated ! to high verbosity. Production runs (verbose <= 2, the default) perform a ! single full-density ddX solve and emit no PCM diag lines. pcm_diag = (infos%control%verbose >= 3) allocate(xyz(3, natom), charges(natom)) xyz(:,:) = infos%atoms%xyz(:, 1:natom) charges(:) = infos%atoms%zn(1:natom) ! ---- Phase 1: build ddX cavity, retrieve cavity-point coordinates ------ max_cav = MAX_CAV_PER_ATOM * natom allocate(cav_xyz(3 * max_cav)) ncav = 0 rc = oqp_ddx_pcm_cavity(natom, xyz, charges, eps, max_cav, ncav, & cav_xyz, cmsg, int(size(cmsg), c_int)) if (rc /= 0) then call c_message_to_fortran(cmsg, fmsg) call show_message('PCM (ddX) cavity build failed: '//trim(fmsg), with_abort) end if allocate(cx(ncav), cy(ncav), cz(ncav)) do icav = 1, ncav cx(icav) = cav_xyz(3*(icav-1) + 1) cy(icav) = cav_xyz(3*(icav-1) + 2) cz(icav) = cav_xyz(3*(icav-1) + 3) end do ! ---- Phase 2: total solute electrostatic potential at cavity points ---- ! Electronic part from the AO density on a temporary copy of the total ! density, so the SCF density blocks are not disturbed. allocate(dtot(nbf_tri), source=0.0_dp) do ii = 1, nfocks dtot(:) = dtot(:) + d(:, ii) end do allocate(phi_elec(ncav), phi_cav(ncav)) call electrostatic_potential_unweighted(basis, cx, cy, cz, dtot, phi_elec) ! phi_total = sum_k Z_k/|r-R_k| + phi_elec, where the OpenQP ! Coulomb-potential primitive returns the electronic contribution with the ! electron-charge sign already included. do icav = 1, ncav phin = 0.0_dp do iat = 1, natom dx = cx(icav) - xyz(1, iat) dy = cy(icav) - xyz(2, iat) dz = cz(icav) - xyz(3, iat) r = sqrt(dx*dx + dy*dy + dz*dz) if (r > 1.0e-12_dp) phin = phin + charges(iat) / r end do phi_cav(icav) = phin + phi_elec(icav) end do ! ---- Phase 3: QM source -> ddX q_cav ----------------------------------- ! The PRODUCTION adjoint source Psi is FULL-DENSITY (3c): atom-centered real ! solid-harmonic multipoles for l = 0..PCM_PSI_LMAX accumulated from the AO ! density by parent-atom Becke-partitioned grid quadrature, in the exact ddX harmonic ! convention, mapped to psi by the ddX rule. The l<=2 Mulliken-multipole ! source (3a/3b) is retained ONLY as a diagnostic baseline. ! Production needs only the full-density adjoint solve (3c) below; q_cav is ! its cavity-projected adjoint charge that drives the Fock matrix and e_pcm. allocate(q_cav(ncav)) if (pcm_diag) then ! (3a/3b) DIAGNOSTIC ONLY -- the l<=2 Mulliken-multipole source, retained ! as a baseline so the q_cav shift from upgrading to the full-density Psi ! can be measured. Skipped in production: it costs two extra ddX solves per ! SCF cycle and never feeds the Fock matrix or e_pcm. allocate(ao_pop(basis%nbf), atom_pop(natom), source_charges(natom), & ao_dip(3, basis%nbf), atom_dip(3, natom), & ao_quad(6, basis%nbf), atom_quad(6, natom), source=0.0_dp) call data_has_tags(infos%dat, tags_overlap, & 'solvent_pcm:add_pcm_reaction_field', with_abort) call tagarray_get_data(infos%dat, OQP_SM, smat) call mulliken_atomic_population_from_density(basis, smat, dtot, & ao_pop, atom_pop) call mulliken_atomic_multipoles_from_density(basis, dtot, atom_pop, & ao_dip, atom_dip, ao_quad, atom_quad) source_charges(:) = charges(:) - atom_pop(:) nmultipoles = 9_c_int allocate(source_multipoles(nmultipoles, natom), source=0.0_dp) call pack_ddx_l2_multipoles(source_charges, atom_dip, atom_quad, source_multipoles) allocate(phi_source(ncav), q_cav_source(ncav), q_cav_l2(ncav)) ! (3a) legacy all-multipole solve: both Phi and Psi from the l<=2 source. rc = oqp_ddx_pcm_solve_multipole_source(natom, xyz, charges, & nmultipoles, source_multipoles, eps, ncav, phi_source, q_cav_source, & esolv_source, cmsg, int(size(cmsg), c_int)) if (rc /= 0) then call c_message_to_fortran(cmsg, fmsg) call show_message('PCM (ddX) diagnostic source solve failed: '//trim(fmsg), with_abort) end if ! (3b) previous production path: exact phi_cav + l<=2 multipole Psi. rc = oqp_ddx_pcm_solve_multipole_source_with_phi(natom, xyz, charges, & nmultipoles, source_multipoles, eps, ncav, phi_cav, q_cav_l2, esolv_l2, & cmsg, int(size(cmsg), c_int)) if (rc /= 0) then call c_message_to_fortran(cmsg, fmsg) call show_message('PCM (ddX) l2-psi diagnostic solve failed: '//trim(fmsg), with_abort) end if end if ! (3c) PRODUCTION solve -- full-density Psi (l = 0..PCM_PSI_LMAX) from the AO ! density + the EXACT total cavity potential phi_cav (Phase 2). Both the ! forward RHS and the adjoint source are full-density. q_cav is the ! cavity-projected adjoint charge (ddx_get_xi) from this solve and is what ! drives the Fock matrix and e_pcm. nbasis = int((PCM_PSI_LMAX+1)**2, c_int) allocate(multipoles_full(nbasis, natom), psi_full(nbasis, natom), radii(natom)) ! Query the production ddX cavity radii FIRST: the full-density Psi must use ! each sphere's rsph for the outside-sphere "leak" continuation (the QM ! density tail beyond the small vdW sphere), exactly as PySCF's ! cache_fake_multipoles does with (r_vdw/r)^(2l+1). rc = oqp_ddx_pcm_radii(natom, charges, radii, cmsg, int(size(cmsg), c_int)) if (rc /= 0) then call c_message_to_fortran(cmsg, fmsg) call show_message('PCM (ddX) radii query failed: '//trim(fmsg), with_abort) end if call build_full_density_multipoles(basis, infos, dtot, charges, radii, xyz, & int(natom), PCM_PSI_LMAX, multipoles_full, & mult_full_norm) call multipoles_to_psi(multipoles_full, radii, PCM_PSI_LMAX, psi_full) psi_full_norm = sqrt(sum(psi_full*psi_full)) rc = oqp_ddx_pcm_solve_psi(natom, xyz, charges, eps, ncav, nbasis, & psi_full, phi_cav, q_cav, esolv, cmsg, int(size(cmsg), c_int)) if (rc /= 0) then call c_message_to_fortran(cmsg, fmsg) call show_message('PCM (ddX) full-density psi solve failed: '//trim(fmsg), with_abort) end if ! ---- Phase 4: V_pcm AO matrix, add to Fock blocks, report E_pcm -------- allocate(vpcm(nbf_tri)) call external_charge_potential(basis, vpcm, cx, cy, cz, q_cav) ! FD self-test of dE/dphi = -0.5*q_cav about the EXACT phi_cav baseline ! (verbose >= 3 only): two extra perturbed ddX re-solves per cycle, purely ! for convention validation. Skipped in production. fd_fock_scale_mean = 0.0_dp; fd_fock_scale_rms = 0.0_dp fd_fock_scale_maxerr = 0.0_dp; fd_fock_samples = 0 if (pcm_diag) then call pcm_fock_scale_fd_diagnostic(natom, xyz, charges, eps, ncav, nbasis, & psi_full, phi_cav, q_cav, & fd_fock_scale_mean, fd_fock_scale_rms, fd_fock_scale_maxerr, & fd_fock_samples) end if ! Fock reaction-field operator. The variational PCM free energy ! E_pcm = -0.5 * <phi_cav(D), q_cav(D)> ! is quadratic in D (both phi_cav and q_cav are linear in D for the ! symmetric ddPCM response), so its derivative dE/dD = -ext(q_cav): ! the explicit factor of 1/2 cancels against the two equal D-dependent terms ! (the phi-side and psi-side contributions, equal by the symmetry of the ! continuum reaction-field kernel). The full coupling is therefore ! 2*PCM_QCAV_TO_FOCK_SCALE = -1, NOT the bare -0.5 explicit-phi factor that ! pcm_fock_scale_fd_diagnostic verifies for dE/dphi. No dielectric factor is ! applied: q_cav already carries the full eps response (see f_epsilon note). vpcm(:) = 2.0_dp * PCM_QCAV_TO_FOCK_SCALE * vpcm(:) do ii = 1, nfocks f(:, ii) = f(:, ii) + vpcm(:) end do ! PCM reaction-field (solvation) energy. The apparent surface charges q_cav ! (from the full-density-Psi exact-phi ddX solve) are contracted with the ! EXACT total solute potential at the cavity points (phi_cav = nuclear + ! electronic). The -0.5 factor is the linear-response polarization factor. ! No additional dielectric factor: -0.5*<phi_cav,q_cav> equals ddX's ! pcm_energy and the physical solvation free energy (Born-ion oracle). ! ! FOCK DERIVATIVE SCOPE: the production SCF uses the full linear-dielectric ! coupling V_pcm = -external_charge_potential(q_cav), recorded as ! fock_mode=ddpcm_physical_full_variational_coupling. The finite-difference ! probe above only verifies the explicit dE/dphi relation (-0.5*q_cav); ! future analytic gradients/response work should add a dedicated dPsi/dD ! check for the grid-projected full-density source. e_pcm = PCM_QCAV_TO_FOCK_SCALE * dot_product(phi_cav, q_cav) ! ---- Diagnostic block (validation gate; does NOT affect e_pcm or Fock) -- ! Exposes, in Fortran, the quantities needed to validate the QM SCF PCM ! conventions against a reference (see tests/test_pcm_literature_benchmarks.py ! and tests/data/pcm_literature_benchmarks.json). These prints are read-only ! summaries of arrays already computed above; the energy and Fock are ! unchanged. The host-side polarization energy 0.5*Tr[D.V_pcm] is reported ! alongside the ddX esolv so the e_pcm-vs-(1/2)Tr[D.V] bookkeeping question ! can be measured rather than assumed. psi_source records the QM source now ! used consistently for ddX phi and psi. if (pcm_diag) then half_tr_dv = 0.5_dp * traceprod_sym_packed(dtot, vpcm, basis%nbf) q_cav_sum = sum(q_cav) q_cav_absnorm = sqrt(sum(q_cav*q_cav)) phi_cav_sum = sum(phi_cav) phi_cav_min = minval(phi_cav) phi_cav_max = maxval(phi_cav) source_charge_sum = sum(source_charges) ! RMS shift in the surface charge from the full-density exact-phi production ! solve relative to the legacy all-multipole (l<=2 Phi and Psi) solve. q_cav_shift_rms = sqrt(sum((q_cav - q_cav_source)**2) / real(ncav, dp)) ! RMS shift from upgrading the adjoint source from l<=2 multipole Psi to the ! full-density Psi, both with the EXACT phi_cav: the quantitative measure of ! what the full-density Psi buys over the previous production path. q_cav_full_vs_l2_rms = sqrt(sum((q_cav - q_cav_l2)**2) / real(ncav, dp)) phi_source_delta_rms = 0.0_dp phi_source_delta_max = 0.0_dp ndelta = 0 do icav = 1, ncav if (ieee_is_finite(phi_source(icav)) .and. ieee_is_finite(phi_cav(icav))) then phi_source_delta_rms = phi_source_delta_rms + & (phi_source(icav) - phi_cav(icav))**2 phi_source_delta_max = max(phi_source_delta_max, & abs(phi_source(icav) - phi_cav(icav))) ndelta = ndelta + 1 end if end do if (ndelta > 0) then phi_source_delta_rms = sqrt(phi_source_delta_rms / real(ndelta, dp)) else phi_source_delta_rms = huge(1.0_dp) phi_source_delta_max = huge(1.0_dp) end if write(iw,'(1x,"PCM diag e_pcm=",ES22.14)') e_pcm write(iw,'(1x,"PCM diag esolv_full_density_psi=",ES22.14)') esolv write(iw,'(1x,"PCM diag esolv_l2_psi_exact_phi=",ES22.14)') esolv_l2 write(iw,'(1x,"PCM diag esolv_source_multipole=",ES22.14)') esolv_source write(iw,'(1x,"PCM diag half_tr_dv=",ES22.14)') half_tr_dv write(iw,'(1x,"PCM diag q_cav_sum=",ES22.14)') q_cav_sum write(iw,'(1x,"PCM diag q_cav_absnorm=",ES22.14)') q_cav_absnorm write(iw,'(1x,"PCM diag fock_q_scale=",ES22.14)') PCM_QCAV_TO_FOCK_SCALE write(iw,'(1x,"PCM diag f_epsilon=",ES22.14)') f_epsilon write(iw,'(1x,"PCM diag fock_q_coupling=",ES22.14)') & 2.0_dp * PCM_QCAV_TO_FOCK_SCALE write(iw,'(1x,"PCM diag fd_fock_scale_mean=",ES22.14)') fd_fock_scale_mean write(iw,'(1x,"PCM diag fd_fock_scale_rms=",ES22.14)') fd_fock_scale_rms write(iw,'(1x,"PCM diag fd_fock_scale_maxerr=",ES22.14)') fd_fock_scale_maxerr write(iw,'(1x,"PCM diag fd_fock_samples=",I0)') fd_fock_samples write(iw,'(1x,"PCM diag source_charge_sum=",ES22.14)') source_charge_sum write(iw,'(1x,"PCM diag phi_source_vs_exact_rms=",ES22.14)') phi_source_delta_rms write(iw,'(1x,"PCM diag phi_source_vs_exact_max=",ES22.14)') phi_source_delta_max write(iw,'(1x,"PCM diag phi_cav_sum=",ES22.14)') phi_cav_sum write(iw,'(1x,"PCM diag phi_cav_min=",ES22.14)') phi_cav_min write(iw,'(1x,"PCM diag phi_cav_max=",ES22.14)') phi_cav_max write(iw,'(1x,"PCM diag ncav=",I0)') int(ncav) write(iw,'(1x,"PCM diag q_cav_source_vs_exact_rms=",ES22.14)') q_cav_shift_rms write(iw,'(1x,"PCM diag q_cav_full_vs_l2_rms=",ES22.14)') q_cav_full_vs_l2_rms write(iw,'(1x,"PCM diag multipoles_full_norm=",ES22.14)') mult_full_norm write(iw,'(1x,"PCM diag psi_full_norm=",ES22.14)') psi_full_norm block integer :: ldiag, mdiag, lmdiag real(dp) :: lnorm do ldiag = 0, PCM_PSI_LMAX lnorm = 0.0_dp do mdiag = -ldiag, ldiag lmdiag = ldiag*ldiag + ldiag + 1 + mdiag lnorm = lnorm + sum(multipoles_full(lmdiag,:)**2) end do write(iw,'(1x,"PCM diag mult_l",I0,"_norm=",ES22.14)') ldiag, sqrt(lnorm) end do write(iw,'(1x,"PCM diag atom_q0=",10ES16.8)') & multipoles_full(1,:) * sqrt(4.0_dp*acos(-1.0_dp)) end block write(iw,'(1x,"PCM diag pcm_source_mode=full_density_multipoles_lmax8_exact_phi")') write(iw,'(1x,"PCM diag psi_source=full_density_grid_multipoles_lmax8_becke3_treutler_parent_atom_leak")') write(iw,'(1x,"PCM diag fock_mode=ddpcm_physical_full_variational_coupling")') end if end subroutine add_pcm_reaction_field subroutine mulliken_atomic_population_from_density(basis, smat, density, ao_pop, atom_pop) type(basis_set), intent(in) :: basis real(dp), intent(in) :: smat(:), density(:) real(dp), intent(out) :: ao_pop(:), atom_pop(:) integer :: mu, nu, ish, iatom, i0, i1, idx ao_pop(:) = 0.0_dp atom_pop(:) = 0.0_dp do mu = 1, basis%nbf do nu = 1, basis%nbf idx = packed_index(mu, nu) ao_pop(mu) = ao_pop(mu) + density(idx) * smat(idx) end do end do do ish = 1, basis%nshell iatom = basis%origin(ish) i0 = basis%ao_offset(ish) i1 = basis%ao_offset(ish) + basis%naos(ish) - 1 atom_pop(iatom) = atom_pop(iatom) + sum(ao_pop(i0:i1)) end do end subroutine mulliken_atomic_population_from_density subroutine mulliken_atomic_multipoles_from_density(basis, density, atom_pop, ao_dip, atom_dip, ao_quad, atom_quad) type(basis_set), intent(in) :: basis real(dp), intent(in) :: density(:) real(dp), intent(in) :: atom_pop(:) real(dp), intent(out) :: ao_dip(:,:), atom_dip(:,:) real(dp), intent(out) :: ao_quad(:,:), atom_quad(:,:) integer :: iat, mu, nu, ish, i0, i1, idx integer, allocatable :: ao_atom(:) real(dp), allocatable :: moment_ints(:,:) atom_dip(:,:) = 0.0_dp atom_quad(:,:) = 0.0_dp allocate(ao_atom(basis%nbf), moment_ints(size(density), 9)) moment_ints(:,:) = 0.0_dp ao_atom(:) = 0 do ish = 1, basis%nshell iat = basis%origin(ish) i0 = basis%ao_offset(ish) i1 = basis%ao_offset(ish) + basis%naos(ish) - 1 ao_atom(i0:i1) = iat end do do iat = 1, size(atom_pop) ao_dip(:,:) = 0.0_dp ao_quad(:,:) = 0.0_dp moment_ints(:,:) = 0.0_dp call multipole_integrals(basis, moment_ints, basis%atoms%xyz(:, iat), 2) do mu = 1, basis%nbf if (ao_atom(mu) /= iat) cycle do nu = 1, basis%nbf idx = packed_index(mu, nu) ao_dip(1, mu) = ao_dip(1, mu) + density(idx) * moment_ints(idx, 1) ao_dip(2, mu) = ao_dip(2, mu) + density(idx) * moment_ints(idx, 2) ao_dip(3, mu) = ao_dip(3, mu) + density(idx) * moment_ints(idx, 3) ao_quad(1, mu) = ao_quad(1, mu) + density(idx) * moment_ints(idx, 4) ao_quad(2, mu) = ao_quad(2, mu) + density(idx) * moment_ints(idx, 5) ao_quad(3, mu) = ao_quad(3, mu) + density(idx) * moment_ints(idx, 6) ao_quad(4, mu) = ao_quad(4, mu) + density(idx) * moment_ints(idx, 7) ao_quad(5, mu) = ao_quad(5, mu) + density(idx) * moment_ints(idx, 8) ao_quad(6, mu) = ao_quad(6, mu) + density(idx) * moment_ints(idx, 9) end do end do atom_dip(:, iat) = -sum(ao_dip(:, :), dim=2) atom_quad(:, iat) = -sum(ao_quad(:, :), dim=2) end do end subroutine mulliken_atomic_multipoles_from_density subroutine pack_ddx_l2_multipoles(source_charges, atom_dip, atom_quad, multipoles) real(dp), intent(in) :: source_charges(:), atom_dip(:,:), atom_quad(:,:) real(dp), intent(out) :: multipoles(:,:) integer :: iat real(dp) :: sqrt4pi, sqrt4pi_over3 real(dp), parameter :: q_xy = 1.0925484305920792_dp real(dp), parameter :: q_z2 = 0.31539156525252005_dp real(dp), parameter :: q_x2y2 = 0.5462742152960396_dp sqrt4pi = sqrt(4.0_dp * acos(-1.0_dp)) sqrt4pi_over3 = sqrt(4.0_dp * acos(-1.0_dp) / 3.0_dp) multipoles(:,:) = 0.0_dp do iat = 1, size(source_charges) ! ddX real-solid-harmonic order for l<=2 is: ! 1: charge; 2: y dipole; 3: z dipole; 4: x dipole; ! 5: xy; 6: yz; 7: z^2; 8: xz; 9: x^2-y^2. multipoles(1, iat) = source_charges(iat) / sqrt4pi multipoles(2, iat) = -atom_dip(2, iat) / sqrt4pi_over3 multipoles(3, iat) = atom_dip(3, iat) / sqrt4pi_over3 multipoles(4, iat) = -atom_dip(1, iat) / sqrt4pi_over3 multipoles(5, iat) = q_xy * atom_quad(4, iat) multipoles(6, iat) = q_xy * atom_quad(6, iat) multipoles(7, iat) = q_z2 * (-atom_quad(1, iat) - atom_quad(2, iat) + & 2.0_dp * atom_quad(3, iat)) multipoles(8, iat) = q_xy * atom_quad(5, iat) multipoles(9, iat) = q_x2y2 * (atom_quad(1, iat) - atom_quad(2, iat)) end do end subroutine pack_ddx_l2_multipoles !> @brief Accumulate the full-density ddPCM source moment of one grid point !> onto its parent sphere, with the outside-sphere exterior ("leak") !> continuation. !> !> The real-solid-harmonic basis Y_lm is evaluated by ddX's OWN ylmbas !> (use ddx_harmonics), so the per-l normalisation and (l,m) ordering match !> the multipole array ddX consumes internally (multipole_psi: !> psi(lm,isph)=4*pi/((2l+1) rsph^l) M_lm). ddX therefore stays an external, !> dynamically-linked dependency and no harmonic code is copied into OpenQP. !> Only the interior/exterior bookkeeping below -- the literature leak !> continuation q*rsph^(2l+1)/r^(l+1) for points outside their sphere -- is !> OpenQP's own and is NOT part of ddX's point-to-multipole routine: !> !> r <= rsph : M_lm += q * r^l * Y_lm (interior) !> r > rsph : M_lm += q * Y_lm * rsph^(2l+1) / r^(l+1) (exterior leak) !> = (interior term) * (rsph/r)^(2l+1) !> !> @param[in] c radius vector from the charge to the sphere centre !> @param[in] src_q charge of the source point !> @param[in] p maximal degree !> @param[in] vscales scaling factors from ylmscale, dim (p+1)**2 !> @param[in] rsph radius of the sphere the moment is accumulated on !> @param[inout] dst_m multipole coefficients, dim (p+1)**2 (accumulated) subroutine pcm_accumulate_leak(c, src_q, p, vscales, rsph, dst_m) integer, intent(in) :: p real(dp), intent(in) :: c(3), src_q, vscales((p+1)**2), rsph real(dp), intent(inout) :: dst_m((p+1)**2) real(dp) :: vylm((p+1)**2), vplm((p+1)**2), vcos(p+1), vsin(p+1) real(dp) :: rho, ctheta, stheta, cphi, sphi, t, ratio2, sqrt4pi integer :: n, ind #ifdef OQP_ENABLE_DDX if (src_q == 0.0_dp) return sqrt4pi = sqrt(4.0_dp * acos(-1.0_dp)) ! ddX's own real-solid-harmonic evaluation (external library). call ylmbas(c, rho, ctheta, stheta, cphi, sphi, p, vscales, vylm, vplm, & vcos, vsin) if (rho == 0.0_dp) then ! Exactly at the sphere centre (e.g. a nuclear charge): only l=0 survives. dst_m(1) = dst_m(1) + src_q / sqrt4pi return end if if (rho <= rsph .or. rsph <= 0.0_dp) then ! Interior point: standard bare moment q*r^l*Y_lm. t = src_q do n = 0, p ind = n*n + n + 1 dst_m(ind-n:ind+n) = dst_m(ind-n:ind+n) + t*vylm(ind-n:ind+n) t = t * rho end do else ! Exterior (leak) point: q * rsph^(2l+1)/rho^(l+1) * Y_lm ! (= the interior term q*rho^l scaled by (rsph/rho)^(2l+1)). ratio2 = rsph * (rsph/rho) t = src_q * (rsph/rho) do n = 0, p ind = n*n + n + 1 dst_m(ind-n:ind+n) = dst_m(ind-n:ind+n) + t*vylm(ind-n:ind+n) t = t * ratio2 end do end if #else ! Unreachable without ddX (callers gate on the C adapter's ddX-availability ! status); abort defensively. Unused dummy args here are intentional. call show_message('PCM full-density source projection requires a ddX-enabled & &build (OQP_ENABLE_DDX)', with_abort) #endif end subroutine pcm_accumulate_leak subroutine build_full_density_multipoles(basis, infos, density_packed, charges, & radii, xyz, natom, lmax, multipoles, mult_norm) type(basis_set), intent(in) :: basis type(information), intent(inout) :: infos real(dp), intent(in) :: density_packed(:) real(dp), intent(in), target :: charges(:), xyz(:,:) real(dp), intent(in) :: radii(:) integer, intent(in) :: natom, lmax real(dp), intent(out) :: multipoles(:,:) real(dp), intent(out) :: mult_norm type(basis_set) :: basis_work type(dft_grid_t), target :: molgrid type(xc_options_t) :: xc_opts type(pcm_psi_grid_consumer_t) :: dat real(dp), allocatable, target :: density_full(:,:) integer :: nbf, i, j, idx, maxl, nang real(dp) :: sqrt4pi ! Extra ddX ylmscale outputs we do not use, sized per its interface. real(dp) :: v4pi2lp1(lmax+1), vscales_rel((lmax+1)**2) character(kind=c_char) :: saved_xcname(20) integer(c_int64_t) :: saved_partfun, saved_bfc_algo, saved_nrad, & saved_nang, saved_radtype logical(c_bool) :: saved_pruned nbf = basis%nbf if (size(multipoles,1) /= (lmax+1)**2 .or. size(multipoles,2) /= natom) then call show_message('PCM full-density Psi multipole buffer has wrong shape', with_abort) end if allocate(density_full(nbf, nbf), source=0.0_dp) do i = 1, nbf do j = 1, nbf idx = packed_index(i, j) density_full(i,j) = density_packed(idx) * basis%bfnrm(i) * basis%bfnrm(j) end do end do basis_work = basis ! HF/reference-SCF inputs can leave infos%dft%xc_functional_name unset/garbage. ! dft_set_options still inspects the C-string even when need_functional=.false., ! so follow the existing SAP-grid pattern: temporarily blank the name while ! constructing the quadrature grid, then restore it. saved_xcname = infos%dft%xc_functional_name infos%dft%xc_functional_name = c_null_char ! The PCM source-projection grid is PINNED to the reference ddCOSMO/ddPCM ! density-partition convention, independent of any user XC-grid settings: ! * Becke's ORIGINAL fuzzy-cell partition (3 softening iterations, ! JCP 88, 2547 (1988)) -- PTYPE_BECKE3, ! * Treutler-Ahlrichs atomic-size surface shifting chi = sqrt(R_i/R_j) ! (JCP 102, 346 (1995)) over the Becke Bragg-Slater table (H = 0.35 A) ! -- dft_bfc_algo = 2, ! * unpruned 240x302 atomic grids on the standard (MHL) radial map. ! Together with the parent-atom point assignment in pcm_grid_update this ! makes the per-sphere source moments converge to the SAME partitioned ! integrals as the reference ddPCM implementations (e.g. PySCF's ! ddcosmo.make_psi_vmat on its Becke-partitioned atomic grids); the grid ! mesh itself only controls quadrature accuracy, not the partition limit. saved_partfun = infos%dft%dft_partfun saved_bfc_algo = infos%dft%dft_bfc_algo saved_nrad = infos%dft%grid_rad_size saved_nang = infos%dft%grid_ang_size saved_radtype = infos%dft%rad_grid_type saved_pruned = infos%dft%grid_pruned infos%dft%dft_partfun = int(PTYPE_BECKE3, c_int64_t) infos%dft%dft_bfc_algo = 2_c_int64_t infos%dft%grid_rad_size = 240_c_int64_t infos%dft%grid_ang_size = 302_c_int64_t infos%dft%rad_grid_type = 0_c_int64_t infos%dft%grid_pruned = .false._c_bool call dft_initialize(infos, basis_work, molgrid, verbose=.false., need_functional=.false.) infos%dft%dft_partfun = saved_partfun infos%dft%dft_bfc_algo = saved_bfc_algo infos%dft%grid_rad_size = saved_nrad infos%dft%grid_ang_size = saved_nang infos%dft%rad_grid_type = saved_radtype infos%dft%grid_pruned = saved_pruned infos%dft%xc_functional_name = saved_xcname dat%lmax = lmax dat%nbasis = (lmax+1)**2 dat%natom = natom dat%xyz => xyz allocate(dat%vscales(dat%nbasis)) #ifdef OQP_ENABLE_DDX call ylmscale(lmax, dat%vscales, v4pi2lp1, vscales_rel) #else call show_message('PCM full-density source projection requires a ddX-enabled & &build (OQP_ENABLE_DDX)', with_abort) #endif allocate(dat%radii(natom)) dat%radii(:) = radii(1:natom) call dat%pe%init(infos%mpiinfo%comm, infos%mpiinfo%usempi) maxl = maxval(basis_work%am) nang = maxl + 2 xc_opts%isGGA = .false. xc_opts%needTau = .false. xc_opts%hasBeta = .false. xc_opts%isWFVecs = .false. xc_opts%numAOs = nbf xc_opts%maxPts = molgrid%maxSlicePts xc_opts%limPts = molgrid%maxNRadTimesNAng xc_opts%numAtoms = natom xc_opts%maxAngMom = nang xc_opts%nDer = 0 xc_opts%numOccAlpha = infos%mol_prop%nelec_A xc_opts%numOccBeta = infos%mol_prop%nelec_B xc_opts%wfAlpha => density_full xc_opts%molGrid => molgrid ! No density-based weight screening for the source projection: the Becke ! partition has small but nonzero tail weights on every atom's grid, and ! the per-sphere moments should integrate them like the reference does. xc_opts%dft_threshold = 0.0_dp xc_opts%ao_threshold = infos%dft%grid_ao_threshold ! Keep AO pruning enabled using the normal DFT threshold; the consumer uses ! xce%compMOs/compRho so both pruned and unpruned paths are handled by the engine. xc_opts%ao_sparsity_ratio = infos%dft%grid_ao_sparsity_ratio if (infos%dft%grid_pruned) xc_opts%ao_sparsity_ratio = 0.0_dp call run_grid_aos(xc_opts, dat, basis_work) multipoles(:,:) = dat%multipoles(:,:,1) ! Add nuclear point charges exactly at their own atom centers. For c=0 only ! the l=0 real harmonic contributes: q * Y_00 = q/sqrt(4*pi). sqrt4pi = sqrt(4.0_dp * acos(-1.0_dp)) do i = 1, natom multipoles(1,i) = multipoles(1,i) + charges(i) / sqrt4pi end do mult_norm = sqrt(sum(multipoles*multipoles)) call dat%clean() end subroutine build_full_density_multipoles subroutine multipoles_to_psi(multipoles, radii, lmax, psi) real(dp), intent(in) :: multipoles(:,:), radii(:) integer, intent(in) :: lmax real(dp), intent(out) :: psi(:,:) integer :: iat, l, m, lm real(dp) :: pi4, denom pi4 = 4.0_dp * acos(-1.0_dp) psi(:,:) = 0.0_dp do iat = 1, size(multipoles,2) do l = 0, lmax denom = real(2*l + 1, dp) * radii(iat)**l do m = -l, l lm = l*l + l + 1 + m psi(lm,iat) = pi4 * multipoles(lm,iat) / denom end do end do end do end subroutine multipoles_to_psi subroutine pcm_grid_parallel_start(self, xce, nthreads) class(pcm_psi_grid_consumer_t), target, intent(inout) :: self class(xc_engine_t), intent(in) :: xce integer, intent(in) :: nthreads if (allocated(self%multipoles)) deallocate(self%multipoles) allocate(self%multipoles(self%nbasis, self%natom, nthreads), source=0.0_dp) end subroutine pcm_grid_parallel_start subroutine pcm_grid_parallel_stop(self) class(pcm_psi_grid_consumer_t), intent(inout) :: self if (allocated(self%multipoles)) then if (size(self%multipoles,3) > 1) then self%multipoles(:,:,1) = sum(self%multipoles, dim=3) end if call self%pe%allreduce(self%multipoles(:,:,1), size(self%multipoles(:,:,1))) end if end subroutine pcm_grid_parallel_stop subroutine pcm_grid_update(self, xce, mythread) class(pcm_psi_grid_consumer_t), intent(inout) :: self class(xc_engine_t), intent(in) :: xce integer :: mythread real(dp), allocatable :: rho(:,:) real(dp) :: qpt, c(3) integer :: ipt, iown ! PARENT-ATOM partition, exactly the reference ddCOSMO/ddPCM source ! projection (PySCF ddcosmo.make_psi_vmat): every point of the current ! slice belongs to the atom whose atomic grid generated it ! (xce%currAtom); the fuzzy-cell share of the molecular density at that ! point is already carried by the Becke-original/Treutler-shifted ! partition weight inside xce%wts (see build_full_density_multipoles). ! The multipole about the owning sphere uses the outside-sphere leak ! continuation (q*rsph^(2l+1)/r^(l+1) for r>rsph), exactly as PySCF's ! cache_fake_multipoles caps points beyond the vdW sphere, which keeps ! the QM density tail from blowing up the bare interior moment q*r^l. iown = xce%currAtom if (iown < 1 .or. iown > self%natom) then call show_message('PCM psi grid consumer: slice parent atom not set', & with_abort) end if call xce%compMOs() allocate(rho(2, xce%numPts), source=0.0_dp) call xce%compRho(rho) do ipt = 1, xce%numPts qpt = -sum(rho(:,ipt)) * xce%wts(ipt) if (qpt == 0.0_dp) cycle c(:) = xce%xyzw(ipt,1:3) - self%xyz(:,iown) call pcm_accumulate_leak(c, qpt, self%lmax, self%vscales, & self%radii(iown), & self%multipoles(:,iown,mythread)) end do end subroutine pcm_grid_update subroutine pcm_grid_post_update(self, xce, mythread) class(pcm_psi_grid_consumer_t), intent(inout) :: self class(xc_engine_t), intent(in) :: xce integer :: mythread ! No per-slice postprocessing needed; update accumulates directly. end subroutine pcm_grid_post_update subroutine pcm_grid_clean(self) class(pcm_psi_grid_consumer_t), intent(inout) :: self if (allocated(self%vscales)) deallocate(self%vscales) if (allocated(self%radii)) deallocate(self%radii) if (allocated(self%multipoles)) deallocate(self%multipoles) nullify(self%xyz) end subroutine pcm_grid_clean subroutine pcm_fock_scale_fd_diagnostic(natom, xyz, charges, eps, ncav, & nbasis, psi, phi_cav, q_cav, scale_mean, scale_rms, maxerr, nsample) integer(c_int), intent(in) :: natom, ncav, nbasis real(dp), intent(in) :: xyz(:,:), charges(:), eps real(dp), intent(in) :: psi(:,:), phi_cav(:), q_cav(:) real(dp), intent(out) :: scale_mean, scale_rms, maxerr integer, intent(out) :: nsample integer :: icav, slot, worst_slot, sample_idx(PCM_FD_MAX_SAMPLES), rc real(dp) :: sample_abs(PCM_FD_MAX_SAMPLES) real(dp) :: h, eplus, eminus, fd, scale, err, min_abs real(dp), allocatable :: phi_plus(:), phi_minus(:), q_tmp(:) character(kind=c_char) :: cmsg(256) character(len=256) :: fmsg sample_idx(:) = 0 sample_abs(:) = -1.0_dp do icav = 1, ncav if (abs(q_cav(icav)) <= 1.0e-14_dp) cycle worst_slot = 1 min_abs = sample_abs(1) do slot = 2, PCM_FD_MAX_SAMPLES if (sample_abs(slot) < min_abs) then min_abs = sample_abs(slot) worst_slot = slot end if end do if (abs(q_cav(icav)) > min_abs) then sample_abs(worst_slot) = abs(q_cav(icav)) sample_idx(worst_slot) = icav end if end do scale_mean = 0.0_dp scale_rms = 0.0_dp maxerr = 0.0_dp nsample = 0 h = PCM_FD_STEP allocate(phi_plus(ncav), phi_minus(ncav), q_tmp(ncav)) do slot = 1, PCM_FD_MAX_SAMPLES icav = sample_idx(slot) if (icav <= 0) cycle phi_plus(:) = phi_cav(:) phi_minus(:) = phi_cav(:) phi_plus(icav) = phi_plus(icav) + h phi_minus(icav) = phi_minus(icav) - h rc = oqp_ddx_pcm_solve_psi(natom, xyz, charges, eps, ncav, nbasis, & psi, phi_plus, q_tmp, eplus, cmsg, int(size(cmsg), c_int)) if (rc /= 0) then call c_message_to_fortran(cmsg, fmsg) call show_message('PCM (ddX) full-psi phi+ FD solve failed: '//trim(fmsg), with_abort) end if rc = oqp_ddx_pcm_solve_psi(natom, xyz, charges, eps, ncav, nbasis, & psi, phi_minus, q_tmp, eminus, cmsg, int(size(cmsg), c_int)) if (rc /= 0) then call c_message_to_fortran(cmsg, fmsg) call show_message('PCM (ddX) full-psi phi- FD solve failed: '//trim(fmsg), with_abort) end if fd = (eplus - eminus) / (2.0_dp * h) scale = fd / q_cav(icav) err = fd - PCM_QCAV_TO_FOCK_SCALE * q_cav(icav) scale_mean = scale_mean + scale scale_rms = scale_rms + scale * scale maxerr = max(maxerr, abs(err)) nsample = nsample + 1 end do if (nsample > 0) then scale_mean = scale_mean / real(nsample, dp) scale_rms = sqrt(scale_rms / real(nsample, dp)) else scale_mean = huge(1.0_dp) scale_rms = huge(1.0_dp) maxerr = huge(1.0_dp) end if end subroutine pcm_fock_scale_fd_diagnostic pure integer function packed_index(i, j) result(idx) integer, intent(in) :: i, j if (i >= j) then idx = i * (i - 1) / 2 + j else idx = j * (j - 1) / 2 + i end if end function packed_index !> @brief Copy a NUL-terminated C character buffer into a Fortran string. subroutine c_message_to_fortran(cmsg, fmsg) character(kind=c_char), intent(in) :: cmsg(:) character(len=*), intent(out) :: fmsg integer :: i fmsg = '' do i = 1, min(size(cmsg), len(fmsg)) if (cmsg(i) == c_null_char) exit fmsg(i:i) = cmsg(i) end do end subroutine c_message_to_fortran end module solvent_pcm