{"id":"77654d2e-9d14-4aee-af5f-881d0831e889","arxiv_id":"2506.23272","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A machine-learned water model combining a CCSD(T)-quality monomer neural network, flexible distributed charges, and cluster-fitted Lennard-Jones terms reproduces many bulk liquid properties in multi-nanosecond simulations of up to 8000 water molecules.","lead":"The authors build a water model that combines a neural network for the water molecule's internal energy, a charge model for electrostatics, and fitted van der Waals terms, then test it in simulations with thousands of molecules. The value is a generic recipe for building cheap but reasonably accurate force fields for liquids, demonstrated on water.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The workflow's central claim is validated only if the fitted pairwise LJ terms actually transfer beyond the cluster set; the weak spot is that bulk validation observables [rho, dHvap] were selection targets, so a decisive transfer test is still missing.","rationale":"The reader's verdict is CONDITIONAL with confidence HIGH, and the weakest assumption is the pairwise-additive representation. I agree with that identification. The paper is honestly written, with clear acknowledgments of the missing density maximum, over-structuring, and neglected many-body terms. The monomer PES is genuinely strong, and the cluster-hexamer benchmarks are useful diagnostics. The single most load-bearing concern is not that the pair approximation is wrong in principle—effective pair potentials can work well—but that the paper's most impressive bulk numbers (rho and dHvap) were selection targets, so they cannot validate transferability. The true predictions (density maximum, dielectric constant, compressibility, RDF shape, 2150 cm-1 band) are systematically off in ways the authors attribute to many-body physics. That attribution is plausible but not tested: the paper does not directly compare model interaction energies against 3-body benchmarks for bulk-like geometries. Therefore the central claim, that this generic workflow yields accurate condensed-phase models for systems without experimental data, remains conditional on a transferability test. The concrete test I propose—re-selecting parameters without experimental targets—would directly settle whether the cluster-fit LJ terms carry the physics or whether the experimental selection does the work. This matches the reader's weakest assumption and supports the CONDITIONAL verdict without moving it.","tokens_in":24177,"tokens_out":1911,"duration_ms":19830,"concrete_test":"Re-train the LJ fits for M-DFT and M-CC exactly as in Section 3.2 but select parameters using only the RMSE against cluster interaction energies, not rho/dHvap. Then run the same NpT 2000-molecule MD protocol and compare rho, dHvap, gOO(r), kappa, epsilon, and D0 against experiment. If the non-selected parameter set gives errors in rho/dHvap within, say, 2-3% of experiment, the bulk agreement is a genuine transfer consequence of the cluster fits. If the non-selected set drifts far from experiment, the workflow's success depends on experimental selection, and the claim of a generic ab initio-driven route to other liquids is weakened.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The central claim is that a generic cluster-fitting workflow yields a water model with accurate gas-phase spectra and qualitatively good bulk properties. The load-bearing assumption is that pairwise LJ plus kMDCM electrostatics absorbs the missing many-body physics. The paper's own Section 3.4 attributes the missing density maximum to neglect of 3- and 4-body interactions, and Section 3.3 shows hexamer interaction energies are underbound by 2-4 (M-DFT) and 4-8 (M-CC) kcal/mol relative to CCSD(T)/CBS. Those are real warning signs, but they are not conclusive: an effective pair potential does not have to reproduce cluster energetics term-by-term to yield correct bulk properties. What is conclusive is that the two bulk properties used to select among the 300 degenerate LJ fits (Table 3, Section 3.2) are the same two properties that anchor the claim of 'good bulk properties'. The density maximum, dielectric constant, and compressibility are genuine predictions and they deviate. The over-structuring in g(r) and the missing 2150 cm-1 band are also genuine predictions, with plausible physical origins (overstructured TIP5-like electrostatics, neglect of charge transfer). The decisive gap is therefore transferability: no bulk observable that was not already used in parameter selection is shown to be quantitative. The workflow is plausible and the monomer PES is independently strong (MAE 0.007 kcal/mol, VPT2 MAE 3.3 cm-1), but the paper's headline generalization to 'other liquids without experimental observables' rests on this untested transferability.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript proposes a generic, modular workflow for building condensed-phase energy functions: a neural-network monomer PES trained to CCSD(T)-F12 data, kMDCM conformationally dependent electrostatics, and Lennard-Jones 12-6 parameters fitted to cluster interaction energies from either DFT (M-DFT) or CCSD(T)-F12 (M-CC) reference data. For water as the test case, the authors select one of 300 near-degenerate LJ parameter sets by matching the experimental bulk density and heat of vaporization, then evaluate a wide range of cluster, thermodynamic, structural, dielectric, transport, and vibrational properties from MD simulations of 2000-8000 water molecules. The monomer PES is very accurate (energy MAE 0.007 kcal/mol, VPT2 frequency MAE 3.3 cm^-1), and several out-of-sample bulk properties (e.g., hydration free energy, self-diffusion, reorientation time) are in reasonable agreement with experiment, while others (density maximum, dielectric constant, compressibility, third-shell gOO structure) show systematic deviations that the authors attribute mainly to neglected many-body interactions.","tokens_in":24505,"tokens_out":5871,"duration_ms":63646,"significance":"If the workflow is ultimately validated, it would be a useful contribution because it combines high-level monomer energies, a flexible electrostatic model, and a cluster-based parametrization into a computationally efficient force field suitable for large-scale simulations. The paper ships its data and analysis code, and the monomer PES validation is a clear strength: the NN reproduces CCSD(T)-F12 energies and forces closely, conserves energy in NVE runs, and yields VPT2 frequencies within 3.3 cm^-1 of experiment. The explicit comparison of two reference levels (DFT and CCSD(T)) and the systematic evaluation of many observables are also commendable. However, the central validation is weakened because the two headline bulk properties, density and heat of vaporization, are used as selection targets rather than as predictions, and the genuinely out-of-sample bulk properties deviate in ways that the paper itself traces to missing many-body physics. The transferability claim for other liquids is therefore not yet established.","major_comments":[{"comment":"The agreement on density and heat of vaporization is a selection outcome, not an independent prediction. Table 3 explicitly labels these two quantities as 'training reference data,' and Section 3.2 states that the loss L = 0.5 MAE ρ + MAE ΔH was used to select one model from the 300 fitted LJ parameter sets. The statement in Section 3.4 that 'all these results correspond to performance on a test set of observables' is therefore misleading for the 300 K ρ and ΔHvap entries; only quantities not appearing in the selection loss (ρ(T), ΔHvap(T), κ, ε, D0, τ2, gOO(r)) are genuine predictions. Among those, the density maximum is absent, the dielectric constant is 63 (M-DFT) and 72 (M-CC) versus 78.2, the compressibility is 35.2/38.5 versus 45.8, and the third solvation shell of gOO is over-structured. The paper should either re-frame ρ and ΔHvap as calibration targets or supply a bulk observable that was not used in parameter selection and is reproduced quantitatively.","section":"Sections 3.2 and 3.4, Table 3"},{"comment":"The hexamer benchmarks show total interaction energies underbound by 2-4 kcal/mol (M-DFT) and 4-8 kcal/mol (M-CC) relative to CCSD(T)/CBS, even though the relative energies between isomers are better. The text describes this as 'rather successful' in capturing two-body contributions, but because the LJ parameters are fitted to cluster interaction energies and are intended to absorb all non-electrostatic, non-intramolecular physics, the magnitude of the underbinding is a warning that the pair term may not be capturing the effective many-body interactions present in bulk water. This is directly relevant to the transferability claim: if the pair-fit compensates for missing three- and four-body terms in clusters, nothing guarantees the compensation transfers to bulk densities. The authors should discuss this compensation explicitly and, ideally, test transferability on a second liquid or on a bulk property far outside the fitting set.","section":"Section 3.3"},{"comment":"The model selection uses only two observables (ρ and ΔHvap) to choose among 300 candidate LJ parameter sets, and Section 3.4 notes that only one of these 300 models was subsequently assessed. Because all 300 sets fit the cluster interaction energies with RMSE better than 1 kcal/mol, the out-of-sample properties reported in Tables 3 and 4 could depend substantially on which near-degenerate fit is selected. The authors should report the spread of predicted bulk properties across the candidate parameter sets, or at least a sensitivity analysis, to establish that the chosen set is representative rather than a fortuitous choice.","section":"Section 3.2, Section 3.4"}],"minor_comments":[{"comment":"There is a unit inconsistency: Section 3.4 reports 'ρcalc = 0.999 kg/m3' while the rest of the paper uses g/cm3 or g/ml, and Figure 6F labels ΔHvap in 'g/cm3', which should be kcal/mol.","section":"Section 3.4, Figure 6F"},{"comment":"The sign convention for the Lennard-Jones ε values should be stated explicitly; the tabulated values are negative, whereas the usual convention for a well depth is positive, and this may confuse readers who are not familiar with the CHARMM convention.","section":"Table S1"},{"comment":"The sentence 'This leads to a total of 36000 parameter combinations' follows from '180 sets of structures' times 200 optimizations, but the preceding text says 'Three repeats of this procedure were carried out, resulting in 180 sets of structures'; please clarify how 180 arises from the 200 random initial values and three repeats.","section":"Section 2.2"},{"comment":"The caption lists panel labels A1, A2, A3, B, and C1 but the text refers to 'Panels A1-3' without defining the layout of these subpanels; a short explanation of the Schlegel projection and its axes would improve readability.","section":"Figure 2 caption"}],"recommendation":"major_revision","confidential_remarks":"The manuscript fits the journal's scope and contains a strong monomer-PES component and a clearly described workflow. The main concern is circularity in the validation: ρ and ΔHvap are selection targets, so the paper's headline agreement on these quantities cannot support the transferability claim. The paper would be substantially improved by reframing those observables as calibration and by presenting a genuinely out-of-sample bulk property as the central validation, along with an analysis of how the choice among the 300 near-degenerate LJ fits affects the reported properties."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The thing to know: this is a solid methods paper that assembles three existing pieces (minimalist NN monomer PES, kMDCM electrostatics, cluster-fitted LJ terms) into a water model and tests it in boxes up to 8000 molecules. The new element is the workflow and its validation, not any single component. The monomer PES is genuinely strong: 0.007 kcal/mol MAE on hold-out energies, VPT2 frequencies within 3.3 cm^-1 of experiment. The MD protocol is careful, including finite-size corrections for diffusion.\n\nThe main soft spot is exactly what the reader and stress-test flagged: density and heat of vaporization are selection targets (Table 3 labels them training data), so agreement on those is not independent validation. The paper is transparent about this, which I credit. The genuinely predicted quantities are mixed: missing density maximum, dielectric constant low (63 vs 78), compressibility low, over-structured g(r), and missing IR band around 2150 cm^-1. Hexamer interaction energies are underbound by 2-8 kcal/mol. The authors attribute these to neglected 3- and 4-body interactions. That is plausible, and an effective pair potential can still produce good bulk properties without reproducing cluster energetics term-by-term. But the paper does not provide a decisive transfer test: no bulk observable that was not used in selection is quantitative. So the headline claim that the workflow generalizes to other liquids without experimental observables is not yet backed.\n\nI do not think the flaws are fatal. The paper is a solid incremental contribution, with released code and data, and the limitations are stated in the text. It deserves a serious referee. The authors should be asked to reframe the abstract so selection targets are separated from predictions, and ideally to add an explicit test of many-body transfer, e.g., 3-body benchmarks or a different liquid with known experimental properties.\n\nFor my own work, I would not build on it directly, but it is a useful reference for anyone constructing ML-based force fields. I would send it to review; the issues are fixable in revision.","headline":"A competent, honest integration of prior ML/electrostatics pieces into a water model with an excellent monomer PES; the bulk validation is partly selection, so the transferability claim is untested.","tokens_in":25075,"tokens_out":1918,"would_cite":false,"duration_ms":19924,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"A modular machine-learned energy function reproduces bulk water properties in large simulations.","keywords":["machine-learned energy functions","water potential","neural network potential energy surface","kernel-based minimal distributed charges","Lennard-Jones parametrization","condensed-phase molecular dynamics","water clusters","thermodynamic model selection"],"falsifier":"Compute explicit three-body CCSD(T)-F12 interaction energies for a set of water trimers and compare them with the residual that the fitted Lennard-Jones pair term absorbs: if the three-body contribution is comparable to or larger than that residual, the pairwise-additive core of the model cannot be the true source of its bulk accuracy.","tokens_in":23909,"feed_emoji":"💧","tokens_out":9972,"duration_ms":99927,"temperature":0.7,"pith_summary":"This paper tries to establish that a modular, machine-learned energy function can serve as a practical force field for liquid simulations: a neural network represents the CCSD(T)-F12 monomer surface, a kernel-based fluctuating charge model handles electrostatics, and Lennard-Jones (12-6) terms fitted to ab initio cluster interaction energies capture the remaining van der Waals physics. Applying this workflow to water produces two models—M-DFT, with Lennard-Jones terms fitted to density-functional cluster energies, and M-CC, fitted to coupled-cluster energies—that reproduce the liquid density and heat of vaporization within about 0.01 g/ml and 0.2 kcal/mol of experiment after selection among many equivalent parameter sets. The same models give gas-phase vibrational frequencies within a few reciprocal centimeters, correctly identify the prism isomer as the most stable water hexamer (M-CC), and yield reasonable dielectric, compressibility, diffusion, and reorientation properties in multi-nanosecond simulations of 2000 to 8000 molecules. The wider point is that such a workflow could be recycled for other condensed phases where experimental reference data are scarce, although the paper also documents a clear failure mode: neither model produces the density maximum of liquid water near 277 K, which it attributes to the neglect of three- and four-body interactions.","feed_headline":"Machine-learned water model reproduces bulk liquid properties","feed_subtitle":"Neural-network monomer, fluctuating charges, and cluster-fitted van der Waals terms pass density and dynamics tests.","key_machinery":"The machinery is a modular energy decomposition: $E = E_\\mathrm{NN}(\\mathrm{monomer}) + E_\\mathrm{kMDCM} + E_\\mathrm{LJ}$, where $E_\\mathrm{NN}$ is a small feed-forward neural network mapping the three interatomic distances of a water monomer to its energy, $E_\\mathrm{kMDCM}$ is the kernel-based minimal distributed charge model whose six geometry-dependent charges capture the fluctuating electrostatics (including lone-pair-like sites), and $E_\\mathrm{LJ}$ is a 12-6 Lennard-Jones pair potential. The Lennard-Jones parameters are fitted by subtracting the neural-network monomer energies and the kMDCM electrostatics from reference cluster interaction energies, so the pair potential absorbs whatever remains. Model selection is then driven by a loss based on bulk density and heat of vaporization, which chooses one working model out of many equally good cluster fits. This same three-part decomposition is the object the paper argues is generic and reusable.","core_discovery":"The paper's core claim is that decomposing the total energy into an intramolecular neural-network term, a geometry-dependent distributed-charge electrostatic term, and a fitted Lennard-Jones pair term yields an energy function that is both accurate enough for CCSD(T)-quality water and fast enough for boxes of about $10^{4}$ molecules. For water, reference monomer energies and forces came from CCSD(T)-F12B/aug-cc-pVTZ-F12 calculations; reference interaction energies for the Lennard-Jones fit came from ωB97X-V/def2-QZVP for clusters of 2 to 60 monomers (M-DFT) and from CCSD(T)-F12 for clusters of 2 to 4 monomers (M-CC). Because cluster interaction energies alone do not determine a unique Lennard-Jones parametrization, the authors generated hundreds of fits and selected the ones agreeing best with experimental liquid density and heat of vaporization. The resulting models reproduce the bulk density and vaporization enthalpy, are within the right range for dielectric constant, compressibility, self-diffusion, and orientational relaxation, and reproduce the spectroscopy of the monomer and the structure of small clusters. The paper is explicit that the omission of three-body and higher interactions shows up as the missing density maximum and a too-weak temperature dependence of the vaporization enthalpy.","pith_inferences":["The strong correlation between fitted Lennard-Jones σ and ε values means cluster interaction energies underdetermine the pair potential; the choice that wins on density and heat of vaporization may not be the one that transfers to solutes, surfaces, or ice, so those properties are genuine tests of the workflow rather than guarantees.","Because the paper's own diagnostic points to missing three- and four-body terms, a practical extension is to fit the same workflow with explicit trimer three-body energies added to the reference data and check whether the density maximum reappears without sacrificing the already-good bulk properties.","For liquids with stronger many-body electrostatics, such as ionic or deep-eutectic systems, the pairwise-additive Lennard-Jones term would have to absorb more physics, so the density and heat-of-vaporization selection could mask a deficiency; testing temperature dependence and structure would expose this earlier than the two thermodynamic observables used here."],"forward_implications":["The workflow provides a route to CCSD(T)-level condensed-phase simulations for liquids that lack extensive experimental force-field data, because the only experimental inputs used for parameter selection are density and heat of vaporization.","The M-CC model is efficient enough for routine 10 ns molecular dynamics with 2000 to 8000 water molecules on commodity hardware, and remains stable when O-H bonds are constrained with SHAKE at a 1 fs time step.","Because dimers dominate the cluster training set, the models capture two-body CCSD(T) energetics well, and the paper identifies explicit three-body terms as the next step for fixing the missing density maximum and the weak temperature dependence of the vaporization enthalpy.","Adding more experimental observables to the model-selection loss, such as the temperature dependence of density or radial distribution functions, is expected to refine the Lennard-Jones parameters within the same workflow."],"supporting_citations":[{"why":"Supplies the cluster-interaction-energy fitting protocol used to obtain the Lennard-Jones parameters.","marker":"[38]"},{"why":"Provides the kernel-based minimal distributed charge (kMDCM) electrostatics model used for all intermolecular electrostatics.","marker":"[35]"},{"why":"Introduces the minimalist neural-network architecture used for the water monomer potential energy surface.","marker":"[30]"},{"why":"Gives the CCSD(T)-F12 method that produces the coupled-cluster reference energies for monomers and small clusters.","marker":"[46]"},{"why":"Supplies the ωB97X-V density functional used to compute reference interaction energies for the larger water clusters.","marker":"[57]"},{"why":"Provides benchmark CCSD(T)/CBS water hexamer structures and energies used to validate the fitted models.","marker":"[52]"}],"fun_headline_variants":["ML water model matches bulk properties at scale","Neural-net water energy fits bulk liquid tests","Fast ML water model passes density and dynamics","Machine-learned water potential scales to 8000 molecules","Hybrid ML water model reproduces liquid behavior"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that pairwise-additive Lennard-Jones interactions plus geometry-dependent distributed charges, without explicit three-body terms, can represent the effective interactions of bulk water.","fun_headline_variants_meta":{"raw":{"variants":["ML water model matches bulk properties at scale","Neural-net water energy fits bulk liquid tests","Fast ML water model passes density and dynamics","Machine-learned water potential scales to 8000 molecules","Hybrid ML water model reproduces liquid behavior"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000187,"raw_usage":{"total_tokens":1399,"prompt_tokens":1085,"completion_tokens":314,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":701,"completion_tokens_details":{"reasoning_tokens":243}},"tokens_in":701,"tokens_out":314,"duration_ms":3897,"temperature":1.0,"reasoning_tokens":243,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T21:47:43.272953+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute explicit three-body CCSD(T)-F12 interaction energies for a set of water trimers and compare them with the residual that the fitted Lennard-Jones pair term absorbs: if the three-body contribution is comparable to or larger than that residual, the pairwise-additive core of the model cannot be the true source of its bulk accuracy.","supporting_citations":[{"cited_title":"D.; Meuwly, M","cited_arxiv_id":null,"evidence_quote":"Supplies the cluster-interaction-energy fitting protocol used to obtain the Lennard-Jones parameters."},{"cited_title":"Kernel-Based Minimal Distributed Charges: A Conformationally Dependent ESP-Model for Molecular Simulations","cited_arxiv_id":null,"evidence_quote":"Provides the kernel-based minimal distributed charge (kMDCM) electrostatics model used for all intermolecular electrostatics."},{"cited_title":"B.; Knizia, G.; Werner, H.-J","cited_arxiv_id":null,"evidence_quote":"Gives the CCSD(T)-F12 method that produces the coupled-cluster reference energies for monomers and small clusters."},{"cited_title":"B97X-V: A 10-parameter, range-separated hybrid, generalized gradient approximation density functional with nonlocal correlation, designed by a survival-of-the-fittest strategy","cited_arxiv_id":null,"evidence_quote":"Supplies the ωB97X-V density functional used to compute reference interaction energies for the larger water clusters."},{"cited_title":"K.; Straight, S","cited_arxiv_id":null,"evidence_quote":"Provides benchmark CCSD(T)/CBS water hexamer structures and energies used to validate the fitted models."}],"review_version":1}