{"id":"275e0aa5-c100-4c51-ba22-4c769defd8f9","arxiv_id":"2608.08357","paper_version":1,"verdict":"ACCEPT","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"The S0 method for computing chemical potentials from structure factors is generalized to neutral multicomponent mixtures using Gaussian process integration with active learning.","lead":"This paper extends a simulation-based method for computing chemical potentials from two-component to neutral multi-component liquid mixtures, using structure factors from ordinary molecular dynamics and a Gaussian process integration step. It demonstrates the approach on a molten Fe-Cu-Ni alloy and on paracetamol solubility in water-ethanol, making multicomponent thermodynamic properties more accessible from atomistic models.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Matrix OZ extrapolation Eq. (8) with fixed k_cut is the least-secure link; Fe-Cu-Ni validation does not transfer to molecular solvents at 300 K.","rationale":"The reader's weakest assumption correctly identifies the S0 extrapolation as the load-bearing link: the entire pipeline from structure factors to chemical potentials depends on the fidelity of the matrix OZ extrapolation in Eq. (8). I agree with that assessment and with the reader's moderate confidence, but I think the missing k_cut sensitivity analysis is enough to make acceptance conditional rather than unconditional. The Fe-Cu-Ni benchmark is genuinely independent support for the method in a high-temperature metallic alloy, where the low-k structure factors are smooth and the OZ form is likely reliable. However, the paracetamol application at 303 K is the paper's showcase of practical relevance, and there is no independent check that the same fixed cutoff and truncated k^2 form are unbiased for a hydrogen-bonded molecular mixture. The two paracetamol reference points from FEP/TI only fix integration constants; they do not validate the derivative matrix that the GP integrates. Since the deposited code and trajectories make the proposed refitting test straightforward, running it would settle whether the concern lands. If the test passes, the ACCEPT verdict stands; if it fails, the central claim would need qualification. Therefore I recommend CONDITIONAL acceptance pending that robustness check.","tokens_in":14011,"tokens_out":11157,"duration_ms":103292,"concrete_test":"Use the deposited MD trajectories and S0_multi code to recompute S0 from Eq. (8) with k_cut^2 = 0.001, 0.002, 0.005, 0.01, and 0.02 (in units of 4*pi^2/Angstrom^2), and also fit an extended form S^{-1}(k) = S^{-1}(0) + k^2 L + k^4 M. Propagate each fit through Eqs. (10), (15), (18) and the GP integration to obtain G_mix for Fe-Cu-Ni at x_Cu=0.25 and paracetamol solubility in pure water and pure ethanol. If the outputs move by more than the propagated statistical uncertainty, the hand-chosen cutoff is load-bearing and a cutoff criterion or robustness statement is required.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The central claim requires that the extrapolated zero-wavevector limit of the NPT partial structure factor equals the grand-canonical fluctuation matrix B_alpha_beta / sqrt(x_alpha x_beta) used in Eq. (10). The only bridge to k=0 is the fitted matrix Ornstein-Zernike form SOZ(k) = (S^{-1}(0) + k^2 L)^{-1} in Eq. (8), with the cutoff fixed at k_cut^2 = 0.005 * 4*pi^2/Angstrom^2 for both applications and no sensitivity analysis. The Fe-Cu-Ni validation against an independent CALPHAD-style model is real support, but it is a high-temperature metallic liquid whose structure factors are smooth in the fitted k-range; the paracetamol/water/ethanol system at 303 K is molecular and hydrogen-bonded, where the same truncated k^2 form could be biased. A biased S0 propagates through Eqs. (4), (10), (15), and (18) into every derivative that the GP integrates, so the reported solubilities inherit the bias. The two FEP/TI reference points in pure solvents do not constrain the derivative matrix, only the integration constants.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript generalizes the S0 method to neutral multi-component mixtures. It derives the multi-component fluctuation relation between the zero-wavevector limit of partial structure factors and chemical potential derivatives, introduces a matrix Ornstein-Zernike extrapolation, and combines these ingredients with Gaussian-process regression using gradient observations and active learning. The method is applied to compute mixing free energies in liquid Fe-Cu-Ni and paracetamol solubilities in water-ethanol. The Fe-Cu-Ni results agree with an independent CALPHAD-style thermodynamic model and with the two-component S0 method; the paracetamol solubility in pure ethanol falls within the experimental spread, while the pure-water solubility is underpredicted and is attributed to force-field error.","tokens_in":14312,"tokens_out":9867,"duration_ms":93007,"significance":"If the S0-to-chemical-potential pipeline is reliable, the method offers a practical route to multi-component chemical potentials from equilibrium NPT simulations without particle insertion or multi-stage thermodynamic integration. The thermodynamic derivation in Section II and Appendix VI A is clean and follows standard Kirkwood-Buff theory; the Gaussian-process integration with gradient data and the active-learning selection are useful methodological advances. The manuscript provides public code and data, which is a strength. There is no circularity in the central pipeline: the derivative matrix is computed from fluctuation data, and the FEP/TI references fix only integration constants. The main uncertainty is the fidelity of the zero-wavevector extrapolation, and the validation for a high-temperature metallic liquid does not by itself guarantee the same accuracy for hydrogen-bonded molecular solutions.","major_comments":[{"comment":"The entire method rests on the identification of the extrapolated S0 with the grand-canonical fluctuation matrix B, but the only bridge from finite-k NPT structure factors to k=0 is the truncated matrix Ornstein-Zernike form in Eq. (8), used with a fixed cutoff k_cut^2 = 0.005 x 4*pi^2/Angstrom^2 for both systems. No sensitivity analysis is reported for the fitting range, the initial guess for L, or the functional form (e.g., adding a k^4 term). Because a biased S0 propagates through Eqs. (10), (15), and (18) into every chemical potential and solubility, the authors should demonstrate that S0 and the final free energies/solubilities are robust to these choices, and should include the OZ fitting uncertainty in the reported error bars.","section":"Section II A, Eq. (8); Sections III and IV"},{"comment":"The RMSE comparison between CUR and random selection is evaluated against a 231-point full grid that is generated by the same S0/GP pipeline, so it measures internal consistency and sampling efficiency rather than absolute accuracy. The absolute validation is against the independent CALPHAD-style model in Fig. 4, but only for fixed x_Cu slices. The manuscript should state this distinction explicitly, and ideally report a single quantitative error metric against the thermodynamic model over the full composition range used.","section":"Section III, Fig. 3a"},{"comment":"The experimental validation for the mixed-solvent solubility curve is largely qualitative: only the pure-ethanol Form I solubility is matched within the experimental scatter, while the pure-water solubility is below the experimental range. The force-field explanation is plausible but is not quantified here. Since the mixed-solvent trend is the paper's main molecular demonstration, a quantitative comparison against the three experimental datasets, or an independent FEP/TI point at an intermediate solvent composition, would materially strengthen the claim that the method captures the solvent-composition dependence.","section":"Section IV, Fig. 6b"}],"minor_comments":[{"comment":"The NPT structure factor at fixed particle numbers has a trivial k=0 value, so the equivalence in Eq. (4) is a grand-canonical statement that relies on the thermodynamic limit and on the OZ extrapolation; a sentence stating this explicitly where Eq. (5) is introduced would help avoid confusion.","section":"Section II A, Eq. (5)"},{"comment":"The choice of representative atoms (oxygen, nitrogen, hydroxyl-adjacent carbon) rather than molecular centers of mass should be justified; the k=0 limit is the same, but the finite-k behavior and hence the OZ fit may differ.","section":"Section IV"},{"comment":"The mapping between active-learning rounds and the number of selected compositions is not stated; please report the number of points added per round, m, and the stopping criterion used in the CUR selection.","section":"Section II E and Fig. 3"},{"comment":"The GP hyperparameters alpha, theta, and sigma_g are quoted as fixed values, but it is not stated whether they were optimized, chosen by cross-validation, or set a priori; please clarify.","section":"Section IV, Eqs. (23)-(24)"},{"comment":"At the lowest paracetamol mole fraction (x_para = 5e-4), the number of paracetamol molecules in a box of roughly 10,000 solvent molecules is small; the statistical reliability of the S_para-para and S_para-solvent structure factors in this regime should be discussed.","section":"Section IV"}],"recommendation":"major_revision","confidential_remarks":"The paper is appropriate for the journal and the method is potentially impactful. My main reservation is the lack of sensitivity analysis for the Ornstein-Zernike extrapolation, which is the most load-bearing component of the pipeline; the current error bars do not include this source of uncertainty. I would be willing to accept after a major revision that provides this analysis and clarifies the validation metrics."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Quick take: this is a genuine, useful generalization of the S0 method to arbitrary multicomponent neutral mixtures, backed by two validation tests and shipped code/data. The thermodynamic derivation in Section II and the appendix is standard Kirkwood-Buff ensemble switching, and the matrix OZ fit plus GP gradient integration with CUR active learning is a sensible toolbox. The Fe-Cu-Ni mixing free energies agree well with an independent CALPHAD/Widom model, and the paracetamol solubility in pure ethanol lands inside the experimental spread. That is real evidence the method works.\n\nThe GP-with-gradients integration and the CUR-based active learning are the genuinely new bits. The C-component generalization of the S0 relation is algebraically natural, but the matrix OZ extrapolation and the GP machinery make it practical beyond binaries.\n\nSoft spots, in order of importance. The zero-wavevector extrapolation, Eq. (8), uses a hand-chosen cutoff k_cut^2 = 0.005 × 4π^2/Å^2 with no sensitivity analysis. The stress-test note is right that the Fe-Cu-Ni test is a high-temperature metal whose structure factors are smooth, while the paracetamol system is hydrogen-bonded and molecular. That said, the pure-ethanol solubility matching experiment is a nontrivial check that the extrapolation is not badly biased in a molecular solvent. Still, I'd want a referee to ask for a cutoff sweep or a comparison with the scalar OZ form to see how much the S0 values change.\n\nSecond, the GP hyperparameters for the paracetamol system are set by hand (alpha, theta, sigma_g). The paper says these were chosen based on the underlying problem, but there's no sensitivity analysis or marginal-likelihood optimization for that system. This is minor because the Fe-Cu-Ni GP optimizes its hyperparameters and the paracetamol result is consistent with experiment, but it leaves a small question about reproducibility.\n\nThird, the GP reconstruction does not explicitly enforce Gibbs-Duhem on the integrated surface. The derivative matrix U satisfies x^T U = 0 by construction at each point, and the GP uses those derivatives, so the violation would be small, but the paper doesn't quantify it.\n\nThere are no critical red flags. The central claim holds: chemical potentials in multicomponent neutral mixtures can be obtained from moderate equilibrium NPT runs. This is a paper for practitioners of free-energy calculations, Kirkwood-Buff users, and people doing solubility prediction. It deserves a serious referee and, given the code and data, should be publishable after minor revision addressing the cutoff sensitivity and GP hyperparameter choice.","headline":"Solid extension of the S0 method to multicomponent neutral mixtures, with real validation; the OZ extrapolation is the soft spot but not a fatal one.","tokens_in":14804,"tokens_out":2320,"would_cite":true,"duration_ms":20931,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"Structure factors alone give chemical potentials for neutral mixtures.","keywords":["chemical potential","structure factors","multi-component mixtures","S0 method","Gaussian process regression","active learning","solubility","Ornstein-Zernike extrapolation"],"falsifier":"Run the identical Fe-Cu-Ni or paracetamol systems in much larger simulation boxes (or with multiple $k_{cut}$ values) and check whether the extrapolated $S^0$ matrix is independent of box size and cutoff; if $S^0$ drifts, the chemical potentials inherit that drift. A sharper test is to compute chemical potentials at a handful of off-reference compositions with an independent method, such as free-energy perturbation or thermodynamic integration, and compare with the GP-integrated values.","tokens_in":13793,"feed_emoji":"🧪","tokens_out":5087,"duration_ms":42377,"temperature":0.7,"pith_summary":"This paper extends the S0 method, which reads chemical-potential derivatives off particle-number fluctuations in equilibrium simulations, from binary solutions to mixtures with any number of neutral components. The authors derive the high-dimensional version of the fluctuation-response relations, replace the scalar Ornstein–Zernike extrapolation with a matrix fit of the whole structure-factor matrix, and integrate the resulting derivative matrix with a gradient-aware Gaussian process that actively chooses which new compositions to simulate. They validate the method on a molten Fe-Cu-Ni alloy and use it to compute paracetamol solubilities in water-ethanol, matching the experimentally observed ethanol-dependent trend. If correct, the method gives a practical route to mixing free energies and solubilities that avoids particle insertion, thermodynamic integration, and real-space Kirkwood-Buff integrals.","feed_headline":"Structure factors alone give chemical potentials for neutral mixtures","feed_subtitle":"A matrix fit to k→0 structure factors yields mixing free energies and solubilities from ordinary MD simulations.","key_machinery":"The carrying object is the matrix Ornstein–Zernike extrapolation of Eq. (8): $S_{\\mathrm{OZ}}(k) = (S^{-1}(0)+k^2 L)^{-1}$, fitted simultaneously to all diagonal and off-diagonal elements of the partial structure-factor matrix. Extrapolating to $k=0$ gives $S^0$, which equals the grand-canonical fluctuation matrix $B$ scaled by $1/\\sqrt{x_\\alpha x_\\beta}$; inverting and projecting that matrix through Eqs. (10), (15), and (18) yields the chemical-potential derivative matrix $\\Gamma^{\\mathrm{ex}}$ defined on the independent composition variables. The second carrying mechanism is Gaussian process regression that consumes both the derivative rows $\\Gamma[\\alpha,:]$ and a reference chemical potential, producing global chemical-potential surfaces with uncertainty estimates, and a CUR low-rank selection of the gradient covariance matrix that drives the active-learning loop.","core_discovery":"The central claim is that the matrix of chemical-potential derivatives with respect to independent mole fractions, $\\Gamma = U M$, can be obtained directly from partial static structure factors $S_{\\alpha\\beta}(k)$ computed in ordinary NPT simulations, and that integrating $\\Gamma$ yields excess chemical potentials over the whole composition space. The key step is the generalized matrix Ornstein–Zernike fit $S_{\\mathrm{OZ}}(k) = (S^{-1}(0)+k^2 L)^{-1}$, which extracts the zero-wavevector matrix $S^0$ without the finite-size contamination that plagues real-space Kirkwood-Buff integrals. A Gaussian process trained on both function values and gradient data performs the integration, and CUR-based active learning chooses the most informative compositions to simulate next. The paper demonstrates the scheme on a ternary alloy, where mixing free energies agree with a CALPHAD-style thermodynamic model, and on paracetamol solubility, where the computed ethanol dependence matches experiment.","pith_inferences":["A systematic convergence test across box sizes and $k_{cut}$ values would directly probe the one assumption this method leans on, since the paper uses a single cutoff for the matrix Ornstein–Zernike fit.","The same derivative matrix $\\Gamma$ could yield activity coefficients, osmotic compressibilities, or Kirkwood-Buff integrals over the entire composition space at no extra simulation cost.","For mixtures with more than three components, the GP integration's cost grows with dimension; the paper notes more simulations are needed, and warped or non-stationary kernels may matter more as composition space grows."],"forward_implications":["For a neutral mixture with $C$ components, only $C-1$ independent mole fractions are needed, and the formalism handles them through a projection matrix, so the method extends beyond three components.","Mixing free energies of fully miscible non-ideal ternary alloys can be computed from modest NPT simulations, as shown by the Fe-Cu-Ni validation against a CALPHAD-style thermodynamic model.","Solubilities of molecular crystals in mixed solvents follow from the intersection of the GP-integrated solute chemical potential with the solid-state chemical potential; the paracetamol example reproduces the experimental rise and plateau with ethanol content.","The combination of GP integration and CUR-based active learning lowers the number of compositions that must be simulated compared with random sampling, at fixed accuracy of the mixing free energy.","The paper explicitly limits itself to neutral components; extending to charged species needs a charge-aware small-k treatment, which the authors say will be the subject of a follow-up."],"supporting_citations":[{"why":"Introduces the original S0 method for binary mixtures, including the zero-wavevector structure-factor relation and the scalar Ornstein–Zernike extrapolation that this paper generalizes.","marker":"[15]"},{"why":"Kirkwood–Buff theory, source of the thermodynamic relation between particle-number fluctuations and chemical-potential derivatives that the method exploits.","marker":"[12]"},{"why":"Ben-Naim's molecular theory of solutions, used for the ensemble-switching and fluctuation-response identities in the derivation.","marker":"[22]"},{"why":"Provides the CALPHAD-style thermodynamic model and virtual semigrand canonical Widom data used as the benchmark for Fe-Cu-Ni mixing free energies.","marker":"[38]"},{"why":"Supplies the reference absolute chemical potentials and solid-state chemical potentials of paracetamol polymorphs used to fix and validate the solubility calculations.","marker":"[16]"},{"why":"Provides the CUR decomposition algorithm used by the active-learning scheme to select the most informative new compositions.","marker":"[37]"},{"why":"Gaussian process regression with gradient observations, the integration method that replaces path-based numerical integration.","marker":"[27-29]"},{"why":"Callen thermodynamics, used to show that the second term in the total-derivative relation sums to zero via the Gibbs–Duhem condition.","marker":"[26]"}],"fun_headline_variants":["S0 method extracts chemical potentials from structure factors","Structure factors alone yield mixture chemical potentials","Gaussian process integration turns S0 into chemical potentials"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The zero-wavevector limit of the NPT partial structure factors, obtained from the matrix Ornstein–Zernike fit with the chosen cutoff, must equal the true grand-canonical fluctuation matrix; any bias in that extrapolation propagates into every chemical potential and solubility.","fun_headline_variants_meta":{"raw":{"variants":["S0 method extracts chemical potentials from structure factors","Structure factors alone yield mixture chemical potentials","Gaussian process integration turns S0 into chemical potentials"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000596,"raw_usage":{"total_tokens":2754,"prompt_tokens":876,"completion_tokens":1878,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":492,"completion_tokens_details":{"reasoning_tokens":1832}},"tokens_in":492,"tokens_out":1878,"duration_ms":12625,"temperature":1.0,"reasoning_tokens":1832,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T00:07:48.853665+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the identical Fe-Cu-Ni or paracetamol systems in much larger simulation boxes (or with multiple $k_{cut}$ values) and check whether the extrapolated $S^0$ matrix is independent of box size and cutoff; if $S^0$ drifts, the chemical potentials inherit that drift. A sharper test is to compute chemical potentials at a handful of off-reference compositions with an independent method, such as free-energy perturbation or thermodynamic integration, and compare with the GP-integrated values.","supporting_citations":[{"cited_title":"Cheng, J","cited_arxiv_id":null,"evidence_quote":"Introduces the original S0 method for binary mixtures, including the zero-wavevector structure-factor relation and the scalar Ornstein–Zernike extrapolation that this paper generalizes."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Kirkwood–Buff theory, source of the thermodynamic relation between particle-number fluctuations and chemical-potential derivatives that the method exploits."},{"cited_title":"Ben-Naim,Molecular Theory of Solutions(Oxford University Press, Oxford, UK, 2006)","cited_arxiv_id":null,"evidence_quote":"Ben-Naim's molecular theory of solutions, used for the ensemble-switching and fluctuation-response identities in the derivation."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the CALPHAD-style thermodynamic model and virtual semigrand canonical Widom data used as the benchmark for Fe-Cu-Ni mixing free energies."},{"cited_title":"Reinhardt, P","cited_arxiv_id":null,"evidence_quote":"Supplies the reference absolute chemical potentials and solid-state chemical potentials of paracetamol polymorphs used to fix and validate the solubility calculations."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the CUR decomposition algorithm used by the active-learning scheme to select the most informative new compositions."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Callen thermodynamics, used to show that the second term in the total-derivative relation sums to zero via the Gibbs–Duhem condition."}],"review_version":1}