pith. sign in

arxiv: 2605.25081 · v1 · pith:GI7II7UGnew · submitted 2026-05-24 · ⚛️ physics.chem-ph · cond-mat.dis-nn· cond-mat.soft· physics.comp-ph

Separable Force Matching of PBE0 Hybrid-Functional Reference Forces for Path-Integral Simulations of Liquid Water

Jan Kessler , Thomas Spura , Kristof Karhan , Thomas D. K\"uhne This is my paper

Pith reviewed 2026-06-29 23:41 UTC · model grok-4.3

classification ⚛️ physics.chem-ph cond-mat.dis-nncond-mat.softphysics.comp-ph
keywords force matchingPBE0path-integral simulationsflexible water modelBuckingham potentialLennard-Jonesradial distribution functionsliquid water
0
0 comments X

The pith

Separable nonlinear least-squares fitting to PBE0 forces produces a stable flexible four-site water model for path-integral simulations.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper develops a force-matching protocol that takes PBE0 hybrid-functional reference forces from liquid water configurations and fits them to a classical flexible four-site potential. It reformulates the parametrization as a separable nonlinear least-squares problem so that all linear amplitudes are eliminated analytically for each trial set of nonlinear parameters. This reduces the dimension of the search, limits compensation between parameters, and permits a controlled comparison between Lennard-Jones and Buckingham oxygen-oxygen repulsion terms. The resulting models reproduce the target force distributions and give radial distribution functions that align with first-principles calculations and neutron-scattering data while remaining numerically stable under path-integral sampling.

Core claim

The projected optimization reproduces target force distributions and yields radial distribution functions close to first-principles and neutron-scattering benchmarks. The Buckingham representation gives a more flexible short-range repulsive wall than a purely Lennard-Jones form, and the final flexible model remains stable in path-integral simulations.

What carries the argument

Separable nonlinear least-squares optimization that analytically eliminates all linear force-field amplitudes for every trial set of nonlinear shape parameters, applied to PBE0 reference forces from CP2K.

If this is right

  • The Buckingham short-range repulsion term supplies greater flexibility than a Lennard-Jones form for the oxygen-oxygen wall.
  • Radial distribution functions from the fitted model stay close to both first-principles and experimental benchmarks.
  • The final flexible four-site model executes stable path-integral molecular dynamics runs without numerical instability.
  • The workflow supplies a transparent route from hybrid-functional forces to simulation-ready water potentials.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The separable formulation could be applied to force-matching for other molecular liquids using the same hybrid-functional references.
  • Direct head-to-head comparison of repulsion functional forms becomes feasible once linear-parameter compensation is removed.
  • Lower effective search dimension may reduce the chance of unphysical parameter sets in future water-model derivations.

Load-bearing premise

PBE0 reference forces computed for liquid water configurations are sufficiently accurate and representative that a classical flexible four-site model fitted to them will reproduce the structural and dynamical behavior of real liquid water under path-integral sampling.

What would settle it

If the fitted model's radial distribution functions deviate substantially from independent PBE0 molecular dynamics trajectories or from neutron-scattering measurements, the accuracy claim would be falsified.

Figures

Figures reproduced from arXiv: 2605.25081 by Jan Kessler, Kristof Karhan, Thomas D. K\"uhne, Thomas Spura.

Figure 1
Figure 1. Figure 1: Benchmark of the force-matched water model against the PBE-based reference using the original and modified optimization strategies. The projected optimization eliminates the linear coefficients at each nonlinear step and thereby separates the choice of functional form from the determination of the optimal amplitudes. where the prime excludes self terms when α = β. For orthorhombic simulation cells, efficie… view at source ↗
Figure 2
Figure 2. Figure 2: Site–site radial distribution functions for the PBE validation case. The comparison tests whether the modified force-matching procedure preserves the structural quality of the first-principles reference while retaining the efficiency of a flexible analytic water potential [PITH_FULL_IMAGE:figures/full_fig_p010_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Comparison of Lennard-Jones and Buckingham descriptions of the O–O short-range interaction for the hybrid-functional parametrization. The Buckingham form separates the exponential repulsive wall from the dispersion-like attractive contribution and is therefore more naturally compatible with the projected linear-parameter optimization. broadening expected when quantum nuclear fluctuations are included. Thes… view at source ↗
Figure 4
Figure 4. Figure 4: Radial distribution functions of the final force-matched water potential. The figure illustrates the structural quality retained after separating the linear and nonlinear parameters in the optimization. between intra- and intermolecular nuclear quantum effects. The fitted Buckingham interaction therefore provides a physically motivated route to improve the short-range structure without abandoning the simpl… view at source ↗
Figure 5
Figure 5. Figure 5: Path-integral RDFs for the final fitted potential. The flexible model remains stable under quantum-nuclear sampling and retains the intramolecular broadening required for q-TIP4P/F-like simula￾tions. 5 Conclusions We have reformulated the parametrization of first-principles-based flexible water potentials as a separable nonlinear least-squares problem and combined this optimization strategy with PBE0 hybri… view at source ↗
read the original abstract

Force-matched water models provide a practical route from first-principles reference data to long classical and path-integral molecular simulations. Previous flexible four-site potentials in the spirit of q-TIP4P/F showed that fitting analytic models to density-functional-theory forces can reproduce key structural features of liquid water while retaining the efficiency required for quantum-nuclear sampling. Here we introduce two refinements aimed at making this strategy more accurate and more reproducible for molecular simulation. First, the production fit is based on PBE0 hybrid-functional reference forces and therefore includes the Hartree--Fock exact-exchange contribution in the electronic-structure target. Second, the parametrization is formulated as a separable nonlinear least-squares problem in which all linear force-field amplitudes are eliminated analytically for every trial set of nonlinear shape parameters. The resulting force-matching protocol lowers the dimension of the nonlinear search, reduces compensation between heterogeneous parameters, and enables a controlled comparison of Lennard-Jones and Buckingham oxygen--oxygen repulsion terms. Applied to CP2K reference forces for liquid water, the projected optimization reproduces target force distributions and yields radial distribution functions close to first-principles and neutron-scattering benchmarks. The Buckingham representation gives a more flexible short-range repulsive wall than a purely Lennard-Jones form, and the final flexible model remains stable in path-integral simulations. The method provides a transparent workflow for deriving simulation-ready water potentials from accurate hybrid density-functional reference forces.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

3 major / 0 minor

Summary. The manuscript introduces a separable nonlinear least-squares force-matching protocol to parametrize a flexible four-site water model against PBE0 hybrid-functional forces computed in CP2K for liquid water. Linear amplitudes are eliminated analytically for each trial set of nonlinear shape parameters, enabling a controlled comparison of Lennard-Jones versus Buckingham oxygen-oxygen repulsion. The resulting models are reported to reproduce target force distributions, yield radial distribution functions close to both the PBE0 reference and neutron-scattering benchmarks, and remain stable under path-integral sampling.

Significance. If the quantitative claims hold, the work supplies a reproducible, lower-dimensional workflow for deriving classical potentials directly from hybrid-DFT reference forces that are suitable for quantum-nuclear simulations. The separable formulation is a clear methodological advance over earlier ad-hoc fitting of q-TIP4P/F-style models, and the explicit LJ-versus-Buckingham comparison is useful for the community.

major comments (3)
  1. [Abstract] Abstract and results section: the central claims that the projected optimization 'reproduces target force distributions' and 'yields radial distribution functions close to first-principles and neutron-scattering benchmarks' are presented without any quantitative error metrics (RMSE, MAE, Kolmogorov-Smirnov statistics), error bars, or measures of agreement. This absence makes it impossible to judge the degree of reproduction or improvement relative to prior models.
  2. [Methods] Methods and results: no description is given of the training/test split, the number of independent configurations, or any cross-validation procedure used in the force-matching optimization. Without this information the generalization of the fitted parameters beyond the reference set cannot be assessed.
  3. [Results] Results: the assertion that the Buckingham form provides a 'more flexible short-range repulsive wall' and that the final model 'remains stable in path-integral simulations' is stated without supporting quantitative diagnostics (e.g., force-component histograms, energy conservation plots, or diffusion constants) that would allow independent verification.

Simulated Author's Rebuttal

3 responses · 0 unresolved

We thank the referee for the constructive comments on our manuscript. We address each major point below and have revised the manuscript to incorporate additional quantitative information and clarifications where the original presentation was incomplete.

read point-by-point responses

  1. Referee: [Abstract] Abstract and results section: the central claims that the projected optimization 'reproduces target force distributions' and 'yields radial distribution functions close to first-principles and neutron-scattering benchmarks' are presented without any quantitative error metrics (RMSE, MAE, Kolmogorov-Smirnov statistics), error bars, or measures of agreement. This absence makes it impossible to judge the degree of reproduction or improvement relative to prior models.

    Authors: We agree that explicit quantitative metrics strengthen the claims. The revised manuscript now reports RMSE values for the force components relative to the PBE0 reference, MAE and Kolmogorov-Smirnov statistics for the oxygen-oxygen, oxygen-hydrogen, and hydrogen-hydrogen radial distribution functions against both the reference data and neutron-scattering benchmarks, together with statistical error bars obtained from block averaging of independent trajectories. revision: yes

  2. Referee: [Methods] Methods and results: no description is given of the training/test split, the number of independent configurations, or any cross-validation procedure used in the force-matching optimization. Without this information the generalization of the fitted parameters beyond the reference set cannot be assessed.

    Authors: The revised Methods section now specifies the size and provenance of the reference dataset (configurations drawn from the CP2K PBE0 trajectory), describes the selection of a held-out test set of configurations not used during optimization, and reports force-matching errors on that test set to demonstrate generalization. Formal cross-validation was omitted because of the high cost of the hybrid-functional reference calculations; this limitation is now stated explicitly. revision: yes

  3. Referee: [Results] Results: the assertion that the Buckingham form provides a 'more flexible short-range repulsive wall' and that the final model 'remains stable in path-integral simulations' is stated without supporting quantitative diagnostics (e.g., force-component histograms, energy conservation plots, or diffusion constants) that would allow independent verification.

    Authors: We accept that additional diagnostics improve verifiability. The revised Results section includes force-component histograms contrasting the Lennard-Jones and Buckingham short-range terms, total-energy conservation traces from the path-integral trajectories, and self-diffusion constants for the final model. These quantities directly support the statements on repulsive-wall flexibility and long-term stability. revision: yes

Circularity Check

0 steps flagged

Fitting to independent PBE0 forces; minor background self-citation not load-bearing

full rationale

The paper presents a force-matching protocol that optimizes a flexible four-site water model against external CP2K PBE0 hybrid-functional forces on liquid water configurations. The separable nonlinear least-squares formulation analytically eliminates linear amplitudes for each trial set of nonlinear parameters, but this is a standard numerical technique that does not render the fitted outputs equivalent to the reference data by definition. Target force distributions are reproduced by construction of the fit, yet the paper additionally validates derived radial distribution functions against independent neutron-scattering data. No equation reduces the reported structural or dynamical observables to quantities defined solely by the fitted parameters. Background reference to prior q-TIP4P/F-style models constitutes at most a minor self-citation that is not invoked to establish uniqueness theorems or to forbid alternative functional forms. The procedure therefore remains self-contained against external benchmarks, consistent with a low circularity score.

Axiom & Free-Parameter Ledger

1 free parameters · 0 axioms · 0 invented entities

The central claim rests on the validity of force-matching to PBE0 forces computed in CP2K and on the assumption that the chosen four-site functional form plus separability can capture the essential physics without additional many-body terms.

free parameters (1)
  • nonlinear shape parameters
    The nonlinear parameters that define the functional form (e.g., exponents or widths in the repulsion and dispersion terms) remain free and are optimized numerically; linear amplitudes are eliminated analytically.

pith-pipeline@v0.9.1-grok · 5806 in / 1199 out tokens · 26706 ms · 2026-06-29T23:41:25.632169+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

34 extracted references · 27 canonical work pages

  1. [1]

    J. L. F. Abascal, C. Vega, A general purpose model for the condensed phases of water: TIP4P/2005,J. Chem. Phys.123, 234505 (2005).https://doi.org/10.1063/1.2121687

    work page doi:10.1063/1.2121687 2005

  2. [2]

    C. Vega, J. L. F. Abascal, Simulating water with rigid non-polarizable models: a general perspective,Phys. Chem. Chem. Phys.13, 19663–19688 (2011). https://doi.org/10.1039/ C1CP22168J

    2011

  3. [3]

    D. Marx, M. Parrinello, Ab initio path-integral molecular dynamics,Z. Phys. B Condens. Matter95, 143–144 (1994).https://doi.org/10.1007/BF01312185

    work page doi:10.1007/bf01312185 1994

  4. [4]

    D. Marx, M. Parrinello, Ab initio path integral molecular dynamics: basic ideas,J. Chem. Phys.104, 4077–4082 (1996).https://doi.org/10.1063/1.471221

    work page doi:10.1063/1.471221 1996

  5. [5]

    M. E. Tuckerman, D. Marx, M. L. Klein, M. Parrinello, Efficient and general algorithms for path integral Car–Parrinello molecular dynamics,J. Chem. Phys.104, 5579–5588 (1996). https://doi.org/10.1063/1.471771

    work page doi:10.1063/1.471771 1996

  6. [6]

    B. Chen, I. Ivanov, M. L. Klein, M. Parrinello, Hydrogen bonding in water,Phys. Rev. Lett. 91, 215503 (2003).https://doi.org/10.1103/PhysRevLett.91.215503

    work page doi:10.1103/physrevlett.91.215503 2003

  7. [7]

    D. Marx, M. E. Tuckerman, J. Hutter, M. Parrinello, The nature of the hydrated excess proton in water,Nature397, 601–604 (1999).https://doi.org/10.1038/17579

    work page doi:10.1038/17579 1999

  8. [8]

    Spura, C

    T. Spura, C. John, S. Habershon, T. D. Kühne, Nuclear quantum effects in liquid water from path-integral simulations using an ab initio force-matching approach,Mol. Phys.113, 808–822 (2015).https://doi.org/10.1080/00268976.2014.981231

    work page doi:10.1080/00268976.2014.981231 2015

  9. [9]

    Kessler, H

    J. Kessler, H. Elgabarty, T. Spura, K. Karhan, P. Partovi-Azar, A. A. Hassanali, T. D. Kühne, Structure and dynamics of the instantaneous water/vapor interface revisited by path-integral and ab initio molecular dynamics simulations,J. Phys. Chem. B119, 10079–10086 (2015). https://doi.org/10.1021/acs.jpcb.5b04185 13

    work page doi:10.1021/acs.jpcb.5b04185 2015

  10. [10]

    On-the-fly

    T. Spura, H. Elgabarty, T. D. Kühne, “On-the-fly” coupled cluster path-integral molecular dynamics: impact of nuclear quantum effects on the protonated water dimer,Phys. Chem. Chem. Phys.17, 14355–14359 (2015).https://doi.org/10.1039/C4CP05192K

    work page doi:10.1039/c4cp05192k 2015

  11. [11]

    Ercolessi, J

    F. Ercolessi, J. B. Adams, Interatomic potentials from first-principles calculations: the force-matching method,Europhys. Lett.26, 583–588 (1994). https://doi.org/10.1209/ 0295-5075/26/8/005

    1994

  12. [12]

    Köster, T

    A. Köster, T. Spura, G. Rutkai, J. Kessler, H. Wiebeler, J. Vrabec, T. D. Kühne, Assessing the accuracy of improved force-matched water models derived from ab initio molecular dynamics simulations,J. Comput. Chem.37, 1828–1838 (2016).https://doi.org/10.1002/ jcc.24398

    2016

  13. [13]

    W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L. Klein, Comparison of simple potential functions for simulating liquid water,J. Chem. Phys.79, 926–935 (1983). https://doi.org/10.1063/1.445869

    work page doi:10.1063/1.445869 1983

  14. [14]

    J. Tao, J. P. Perdew, V. N. Staroverov, G. E. Scuseria, Climbing the density functional ladder: nonempirical meta-generalized gradient approximation designed for molecules and solids, Phys. Rev. Lett.91, 146401 (2003).https://doi.org/10.1103/PhysRevLett.91.146401

    work page doi:10.1103/physrevlett.91.146401 2003

  15. [15]

    S.Grimme, J.Antony, S.Ehrlich, H.Krieg, Aconsistentandaccurateabinitioparametrization of density functional dispersion correction (DFT-D) for the 94 elements H–Pu,J. Chem. Phys.132, 154104 (2010).https://doi.org/10.1063/1.3382344

    work page doi:10.1063/1.3382344 2010

  16. [16]

    J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett.77, 3865–3868 (1996).https://doi.org/10.1103/PhysRevLett.77.3865

    work page doi:10.1103/physrevlett.77.3865 1996

  17. [17]

    Adamo, V

    C. Adamo, V. Barone, Toward reliable density functional methods without adjustable pa- rameters: the PBE0 model,J. Chem. Phys.110, 6158–6170 (1999).https://doi.org/10. 1063/1.478522

    1999

  18. [18]

    Todorova, A

    T. Todorova, A. P. Seitsonen, J. Hutter, I.-F. W. Kuo, C. J. Mundy, Molecular dynamics simulation of liquid water: hybrid density functionals,J. Phys. Chem. B110, 3685–3691 (2006).https://doi.org/10.1021/jp055127v

    work page doi:10.1021/jp055127v 2006

  19. [19]

    T. D. Kühneet al., CP2K: An electronic structure and molecular dynamics software package – Quickstep: efficient and accurate electronic structure calculations,J. Chem. Phys.152, 194103 (2020).https://doi.org/10.1063/5.0007045

    work page doi:10.1063/5.0007045 2020

  20. [20]

    Iannuzziet al., The CP2K program package made simple,J

    M. Iannuzziet al., The CP2K program package made simple,J. Phys. Chem. B130, 1237–1310 (2026).https://doi.org/10.1021/acs.jpcb.5c05851

    work page doi:10.1021/acs.jpcb.5c05851 2026

  21. [21]

    Guidon, J

    M. Guidon, J. Hutter, J. VandeVondele, Auxiliary density matrix methods for Hartree–Fock exchange calculations,J. Chem. Theory Comput.6, 2348–2364 (2010). https://doi.org/ 10.1021/ct1002225

    work page doi:10.1021/ct1002225 2010

  22. [22]

    G. H. Golub, V. Pereyra, The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate,SIAM J. Numer. Anal.10, 413–432 (1973). https: //doi.org/10.1137/0710036

    work page doi:10.1137/0710036 1973

  23. [23]

    https://doi.org/10.1088/0266-5611/ 19/2/201

    G.H.Golub, V.Pereyra, Separablenonlinearleastsquares: thevariableprojectionmethodand its applications,Inverse Probl.19, R1–R26 (2003). https://doi.org/10.1088/0266-5611/ 19/2/201

    work page doi:10.1088/0266-5611/ 2003

  24. [24]

    D. P. O’Leary, B. W. Rust, Variable projection for nonlinear least squares problems,Comput. Optim. Appl.54, 579–593 (2013).https://doi.org/10.1007/s10589-012-9492-9 14

    work page doi:10.1007/s10589-012-9492-9 2013

  25. [25]

    Habershon, T

    S. Habershon, T. E. Markland, D. E. Manolopoulos, Competing quantum effects in the dynamics of a flexible water model,J. Chem. Phys.131, 024501 (2009). https://doi.org/ 10.1063/1.3167790

    work page doi:10.1063/1.3167790 2009

  26. [26]

    T. E. Markland, D. E. Manolopoulos, An efficient ring polymer contraction scheme for imaginary time path integral simulations,J. Chem. Phys.129, 024105 (2008). https://doi. org/10.1063/1.2953308

    work page doi:10.1063/1.2953308 2008

  27. [27]

    T. E. Markland, D. E. Manolopoulos, A refined ring polymer contraction scheme for systems with electrostatic interactions,Chem. Phys. Lett.464, 256–261 (2008).https://doi.org/ 10.1016/j.cplett.2008.09.019

    work page doi:10.1016/j.cplett.2008.09.019 2008

  28. [28]

    Available at:https://github.com/DCM-UPB/RPMD_Mainz

    DCM-UPB, RPMD_Mainz: ring-polymer molecular dynamics code for flexible water models, GitHub repository. Available at:https://github.com/DCM-UPB/RPMD_Mainz

  29. [29]

    Lippert, J

    G. Lippert, J. Hutter, M. Parrinello, A hybrid Gaussian and plane wave density functional scheme,Mol. Phys.92, 477–488 (1997).https://doi.org/10.1080/002689797170220

    work page doi:10.1080/002689797170220 1997

  30. [30]

    R. P. Feynman, A. R. Hibbs,Quantum Mechanics and Path Integrals, McGraw–Hill, New York (1965)

    1965

  31. [31]

    M. E. Tuckerman, B. J. Berne, G. J. Martyna, M. L. Klein, Efficient molecular dynamics and hybrid Monte Carlo algorithms for path integrals,J. Chem. Phys.99, 2796–2808 (1993). https://doi.org/10.1063/1.465188

    work page doi:10.1063/1.465188 1993

  32. [32]

    Ceriotti, M

    M. Ceriotti, M. Parrinello, T. E. Markland, D. E. Manolopoulos, Efficient stochastic thermostatting of path integral molecular dynamics,J. Chem. Phys.133, 124104 (2010). https://doi.org/10.1063/1.3489925

    work page doi:10.1063/1.3489925 2010

  33. [33]

    K. A. F. Röhrig, T. D. Kühne, Optimal calculation of the pair correlation function for an orthorhombic system,Phys. Rev. E87, 045301 (2013). https://doi.org/10.1103/ PhysRevE.87.045301

    2013

  34. [34]

    A. K. Soper, The radial distribution functions of water as derived from radiation total scattering experiments: is there anything we can say for sure?,ISRN Phys. Chem.2013, 279463 (2013).https://doi.org/10.1155/2013/279463 15

    work page doi:10.1155/2013/279463 2013