{"id":"51823419-81f1-40f1-bd37-3561e63d1ebe","arxiv_id":"2412.15370","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"MLACS is a production Python package that iteratively trains linear MLIP surrogates with active learning and MBAR reweighting to sample the DFT canonical ensemble at 50 to 100 times lower DFT cost.","lead":"This paper describes MLACS, a Python package that accelerates ab initio molecular dynamics by training a machine learning potential on the fly and reweighting configurations so that ensemble averages match DFT. It reports that the method reaches near-DFT accuracy (about 1 meV/atom) while reducing the number of DFT single-point calculations by one to two orders of magnitude.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"MeV-level free-energy claims rely on an in-sample cumulant correction; out-of-sample validation is needed.","rationale":"The reader's weakest-assumption field identifies expressiveness of the linear MLIP family, which is a real limitation with no general error bound. I partially agree: expressiveness is fundamental, but the paper provides external AIMD comparisons for Mg and Au phonons that give some empirical support. The more acute and checkable flaw is that the free-energy correction ΔF_int→AI is computed in-sample: the same weighted configurations that trained the MLIP are used to estimate the residual correction to the ab initio free energy. This is a systematic bias, not merely a missing error bound, and it directly supports the headline meV-level free-energy/phase-diagram claims. The reader's rationale does mention 'reported free energy corrections are computed in-sample,' so there is overlap, but it is not the reader's stated weakest assumption. I do not change the CONDITIONAL verdict: the software and the canonical-sampling benchmarks are valuable, but the free-energy accuracy claim needs an out-of-sample test before it can be taken at face value.","tokens_in":38377,"tokens_out":10837,"duration_ms":100734,"concrete_test":"Run a held-out validation for the fcc-Au thermodynamic points in Fig. 11. After the MLACS loop converges, generate 50–100 independent configurations by continuing MLMD with the final surrogate potential (or by drawing from the final MBAR weights), compute DFT energies for these configurations, and evaluate Eq. (36) on this out-of-sample set. Compare with the in-sample ΔF_int→AI values and report the kurtosis of ΔV to test the second-order truncation. If out-of-sample |ΔF_int→AI| exceeds 1 meV/atom at any point where Fig. 11 claims ≤1 meV/atom, the free-energy accuracy claim must be revised or replaced by an explicit out-of-sample correction procedure.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central free-energy claim (abstract, §5.3, and Fig. 11) rests on the correction ΔF_int→AI computed from Eq. (36), a second-order cumulant expansion of ΔV = V − eV over the MBAR-reweighted database. These are exactly the configurations and DFT energies used to fit the MLIP at every self-consistent step (Alg. 1, Eq. (17)). The residuals are therefore in-sample: active learning selects configurations to minimize weighted least-squares errors, so the empirical mean and variance of ΔV understate the true discrepancy between the surrogate and DFT distributions. The near-zero values in Fig. 11 (10⁻⁴ meV/atom below 2000 K, ~1 meV/atom at extremes) are evidence that the surrogate fits its own training set, not that the free-energy correction is converged. Moreover, Eq. (36) is a Gaussian approximation to the exact free-energy perturbation identity (33); the paper reports no test of Gaussianity of ΔV, particularly for gold at 10,000 K and 1 TPa where anharmonicity is large. Since the paper claims 'ab initio accuracy (≤1 meV/atom)' for free energies and phase diagrams, and that claim is supported almost entirely by this in-sample correction, the meV-level free-energy results are not established by the data shown.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents MLACS, a Python package that accelerates ab initio canonical sampling by replacing AIMD with a self-consistent variational loop. A linear MLIP (SNAP, MTP, or ACE) is iteratively trained on actively selected DFT energies, forces, and stresses, with MBAR reweighting used to account for the evolving surrogate distribution and to compute canonical averages. The package also implements free-energy calculations via nonequilibrium thermodynamic integration from an Einstein crystal or Uhlenbeck-Ford reference, followed by a cumulant-expansion correction for the residual difference between the MLIP and the DFT target. Applications include canonical sampling of Mg, phonon spectra and a phase diagram for Au, geometry optimization of AuCu alloys, NEB-based vacancy diffusion in Ag, liquid water at 400 K, and UO2 phonons. The central claims are that the self-consistent MLIP sampling reproduces AIMD distributions at roughly a 50-fold reduction in DFT single-point calculations, and that free energies are obtained with ab initio accuracy of about 1 meV/atom.","tokens_in":38664,"tokens_out":7392,"duration_ms":72097,"significance":"The empirical part of the paper is strong and practically valuable. The Mg benchmark compares MLACS directly with AIMD and shows agreement within statistical uncertainty; the Au phonon comparison at two thermodynamic states demonstrates the claimed acceleration; and the AuCu excess-enthalpy results agree with direct DFT optimizations. The open-source package, tutorials, and tests are a further strength. However, the theoretical derivation has two uncontrolled steps, and the meV-level free-energy claim rests on an in-sample cumulant correction rather than on an externally validated error estimate. The paper is likely to be useful to the computational materials community once the free-energy validation and the derivation gaps are addressed.","major_comments":[{"comment":"The derivation of Eq. (7) is incomplete as written. The condition <V>_{\\tilde V} = <\\tilde V>_{\\tilde V} is not a consequence of the preceding gradient calculation; it is imposed 'without loss of generality' because adding a constant to \\tilde V leaves the canonical distribution unchanged. This invariance argument requires that the linear parameterization in Eq. (2) actually contains an adjustable constant (or an equivalent projection), and it must be reconciled with the fact that the fixed point of Eq. (7) is otherwise not constrained to satisfy the offset condition. Please state the constant-offset assumption explicitly and show that Eq. (7) follows from the full stationarity condition of the KLD under that assumption.","section":"Section 1.1, Eq. (7)"},{"comment":"The gradient of the Fisher divergence is truncated by dropping the second term on the right-hand side of Eq. (11) with the justification that it is quadratic in the force (or stress) error and therefore small near the fixed point. This is an uncontrolled approximation for a linear MLIP whose force and stress errors do not vanish at the fixed point, and the magnitude of the dropped term is nowhere quantified or tested. Since Eq. (12) is the actual fitting rule used in Algorithm 1 for force and stress data, the authors should either retain the term, bound it, or provide a numerical check, for example by comparing fits with and without the term on a representative system.","section":"Section 1.1, Eqs. (11)-(12)"},{"comment":"The free-energy correction DeltaF_int_to_AI is computed from a second-order cumulant expansion using the same DFT data that were actively selected for fitting the MLIP. The data are therefore in-sample: the active-learning fit minimizes the weighted energy residuals, so the empirical mean and variance of DeltaV = V - \\tilde V understate the true surrogate-to-DFT discrepancy. The near-zero values in Fig. 11 are not independent evidence of 1 meV/atom free-energy accuracy, and no test of Gaussianity of DeltaV is provided for the extreme Au thermodynamic states. The meV-level claim should be supported by out-of-sample evaluation, for example hold-out configurations or independent thermodynamic-integration calculations at a few state points, and by a check of the cumulant truncation.","section":"Section 2.1, Eqs. (33)-(36) and Fig. 11"},{"comment":"The statement that the method 'prove[s] that a self-consistent active learning strategy using a MLIP enables to sample the BO surface as accurately as using AIMD simulations' overstates what the derivation establishes. The variational calculation shows that the KLD-optimal distribution within the chosen linear family is a fixed point of Eq. (7), but it does not by itself bound the residual KLD or the error of canonical averages. The accuracy claims rest on the empirical benchmarks, which are convincing for the tested systems, but the theoretical wording should be qualified accordingly.","section":"Section 1.1 and abstract"}],"minor_comments":[{"comment":"Reporting RMSE and MAE values to sub-meV and sub-meV/A precision without statistical error bars or the number of samples is potentially misleading; adding error estimates or at least the database size for each weighting policy would help.","section":"Section 5.1, Table 3"},{"comment":"The main object is described as 'Mlas' in one place; this should read 'Mlacs'. The same section also contains the phrase 'Funtions', which should be 'Functions'.","section":"Section 4.4"},{"comment":"The caption contains the typo 'pannel' in 'Left pannel' and 'Right panel'; please correct it.","section":"Figure 10 caption"},{"comment":"The sentence 'the free energy of the anharmonic surrogate system Fint (which is exacly the same quantity as eF0)' contains the typo 'exacly' and should read 'exactly'.","section":"Section 2.1"}],"recommendation":"major_revision","confidential_remarks":"This is a software and methods paper for CPC. The empirical demonstrations are the strongest part of the manuscript and justify publication after revision, but the current free-energy evidence is in-sample and the derivation of the fitting equations has uncontrolled approximations. The theoretical gaps are repairable, so I would not reject; however, the meV-level free-energy claim should not be accepted on the present Fig. 11 alone."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Andreas,\n\nI've read the MLACS paper. The package is the real contribution here: a cleanly orchestrated Python implementation of the variational inference method from their 2022 PRB, with MBAR reweighting, NETI, NEB/PAFI, delta learning, and a sensible object model for plugging in different descriptors and DFT codes. The applications section carries the paper. The Mg benchmark is a fair direct comparison to AIMD: MLACS converges to the same mean energy with an order of magnitude fewer DFT steps and even reproduces the bimodal energy distribution that raw MLIP MD misses. The Au phonon comparisons at 1000 K and 8000 K are also convincing external validation of the canonical sampling. AuCu, silver vacancy, water, and UO2 are useful demonstrations of the package's versatility, even if the comparisons are less rigorous.\n\nThe soft spot is the free energy story. The correction in Eq. (36) is evaluated on the same configurations used to train the surrogate at every step. The near-zero values in the right panel of Fig. 11 are essentially a measure of the training residuals, not of the discrepancy between the surrogate and DFT distributions. For the gold phase diagram, no independent check is offered — no AIMD-based TI at even one thermodynamic point — so the '≤1 meV/atom' free energy accuracy claim is not established by the data shown. The paper also brushes past two technical points in the derivation: the constant-offset condition that lets them drop the energy mean in Eq. (7), and the term discarded from the force/stress gradient in Eq. (11). Both are probable, but they are not justified in the text.\n\nAll that said, the central sampling claim is credible because the phonon benchmarks are external. The free energy issue is a limitation, not a fatal flaw. I'd accept this for peer review with the request that the authors either provide out-of-sample free energy validation for a couple of points or tone down the abstract and conclusion to match what is actually demonstrated. The package is useful, the writing is clear, and the benchmarks are informative. Worth a serious referee.\n\nBest,","headline":"A useful MLACS software release with strong AIMD sampling benchmarks, but the meV-level free energy claims rest on an in-sample correction that needs independent validation.","tokens_in":39254,"tokens_out":3239,"would_cite":true,"duration_ms":32971,"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":"The MLACS package claims a self-consistent MLIP loop samples the Born-Oppenheimer surface as accurately as AIMD, with a ~50x reduction in DFT cost and energies within ~1 meV/atom.","keywords":["machine learning interatomic potential","canonical sampling","variational inference","Kullback-Leibler divergence","MBAR reweighting","free energy","active learning","anharmonicity"],"falsifier":"Run MLACS on a strongly anharmonic or liquid system with a known multimodal canonical distribution, for example liquid water or a high-temperature bcc metal, and compare the reweighted energy histogram, pair distribution functions, and free energy against a long AIMD trajectory and an independent thermodynamic-integration reference. If the histogram misses a second peak, the pair distribution functions deviate beyond the AIMD statistical error, or the free energy drifts by more than about 1 meV/atom, the expressiveness assumption is falsified; a cheaper diagnostic is the second-order cumulant correction $\\Delta F_{\\mathrm{int}\\to\\mathrm{AI}}$, which should stay below 1 meV/atom on a converged loop.","tokens_in":38150,"feed_emoji":"⚛️","tokens_out":9913,"duration_ms":73587,"temperature":0.7,"pith_summary":"The paper presents MLACS, a production software package that moves almost all of the sampling in an ab initio molecular dynamics calculation onto a machine-learned surrogate potential while keeping the ab initio system as the simulated object. Its central claim is that a self-consistent active-learning loop over the parameters of a linear machine-learning interatomic potential, with every stored configuration reweighted by MBAR, converges to the DFT canonical distribution closely enough that energies come within about 1 meV/atom of AIMD while using one to two orders of magnitude fewer DFT single-point calculations. This matters because finite-temperature properties of anharmonic solids, liquids, alloys, and f-electron materials are currently out of reach for direct AIMD over the long trajectories needed for convergence. The theory connects the sampling problem to variational inference: minimizing the Kullback-Leibler divergence between surrogate and DFT distributions is the Gibbs-Bogoliubov inequality, and the MLIP update reduces to a weighted least-squares fit on energies, forces, and stresses.","feed_headline":"MLACS samples the DFT ensemble with 50x fewer DFT calls","feed_subtitle":"A machine-learning surrogate plus MBAR reweighting keeps energies within 1 meV/atom of AIMD at a fraction of the cost.","key_machinery":"The load-bearing object is the self-consistent parameter update for a linear surrogate potential $\\tilde V_\\gamma(R)=\\sum_k \\tilde D_k(R)\\gamma_k$. Minimizing the Kullback-Leibler divergence between the surrogate and DFT distributions gives the fixed-point equation $\\gamma = \\langle \\tilde D^T \\tilde D \\rangle^{-1} \\langle \\tilde D^T V \\rangle$; using forces and stresses through the Fisher divergence gives the analogous equation with $\\nabla_\\eta \\tilde D$ and $\\nabla_\\eta V$. The averages are taken with respect to the surrogate distribution itself, which is why the solve is circular and is iterated as an active-learning loop. The second load-bearing component is MBAR reweighting, which assigns weights so that configurations generated under earlier surrogate potentials remain valid samples for the current one; this reuse is what makes the small database sufficient and keeps the fit from being dominated by off-target configurations.","core_discovery":"On the paper's own terms, the central discovery is that the surrogate distribution should approximate the DFT canonical distribution at a single thermodynamic point rather than the Born-Oppenheimer surface globally, and that this approximation problem is exactly a variational-inference problem. The authors show that minimizing the Kullback-Leibler divergence $D_{\\mathrm{KL}}(q_{\\gamma}\\Vert p)$ is equivalent to minimizing the Gibbs-Bogoliubov free energy, so the optimal surrogate parameters solve the self-consistent least-squares system $\\gamma = \\langle \\tilde D^T \\tilde D \\rangle^{-1} \\langle \\tilde D^T V \\rangle$, with a corresponding force- and stress-based equation obtained from the Fisher divergence. Iterating this solve with short surrogate molecular dynamics runs, adding one DFT-evaluated configuration per cycle, and reweighting the full database by MBAR yields a weighted configuration set that reproduces the DFT canonical ensemble. The demonstrated consequence is that phonon spectra, free energies, and phase boundaries are obtained with near-DFT accuracy, about 1 meV/atom, at one to two orders of magnitude lower DFT cost, including a factor of about fifty for bcc-gold phonon frequencies.","pith_inferences":["Because MBAR reweighting also works when applied a posteriori to an existing configuration database, the same machinery could mine previously computed DFT trajectories for new thermodynamic conditions without additional ab initio runs; the paper demonstrates the reweighting but not this large-scale reuse.","A nonlinear machine-learning potential could replace the linear model inside the same self-consistent loop, but the closed-form least-squares update would be lost and the optimization would no longer be a single matrix solve; this is a natural testbed for whether the expressiveness limitation identified here is the real bottleneck.","The current theory treats nuclei classically; combining the variational loop with path-integral sampling would extend the same KL-minimization idea to quantum nuclear fluctuations, which are relevant for light elements and hydrogen-rich materials.","The second-order cumulant free-energy difference could serve as an on-the-fly convergence criterion, since it measures exactly how far the surrogate distribution is from the target; the package now relies on user-chosen property thresholds instead."],"forward_implications":["Finite-temperature phonon spectra of anharmonic and dynamically unstable crystals can be computed at AIMD accuracy with roughly fifty times fewer DFT steps, as demonstrated for bcc gold.","Free-energy differences and phase diagrams over hundreds of thermodynamic points become feasible at about 1 meV/atom accuracy, because only the cheap surrogate is driven through nonequilibrium thermodynamic integration.","Geometry relaxation and minimum-energy-path searches on large cells converge several times faster than direct DFT optimization and reproduce DFT excess enthalpies, volumes, and migration barriers.","Liquid-state sampling works when delta-learning corrections are added, with the water example reproducing AIMD pair distribution functions over a much longer trajectory.","The same self-consistent loop converts approximate MLIP databases into reweighted ab initio-quality ensembles, so the package can serve as a data-generation engine for building larger training sets."],"supporting_citations":[{"why":"Defines the original MLACS method and its variational-inference derivation, which this paper generalizes into a production package.","marker":"[1]"},{"why":"Provides MBAR, the statistically optimal reweighting that assigns weights to configurations generated by different surrogate potentials.","marker":"[42]"},{"why":"Supplies the Fisher-divergence formulation used to fit the surrogate from forces and stresses.","marker":"[37]"},{"why":"Defines the spectral-neighbor-analysis descriptor family used in most of the paper's demonstrations.","marker":"[32]"},{"why":"Establishes the self-consistent harmonic approximation whose Gibbs-Bogoliubov variational structure MLACS generalizes to arbitrary linear MLIPs.","marker":"[13]"},{"why":"Supplies the nonequilibrium thermodynamic-integration protocol used to carry the surrogate free energy from a reference model to the ab initio system.","marker":"[72]"},{"why":"Underpins the cumulant expansion used to correct the surrogate free energy to the ab initio value.","marker":"[44]"}],"fun_headline_variants":["MLACS: MLIP surrogate matches DFT at 50x fewer calls","MLACS variational sampling: 1 meV accuracy, 50x speedup","MLACS: DFT-free canonical sampling matches AIMD accuracy","MLACS: reweighted MLIP sampling gets 1 meV/atom at 50x faster"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the chosen family of linear machine-learning potentials, built from the selected descriptors, is expressive enough to closely approximate the true quantum-mechanical distribution at the thermodynamic point of interest; if those descriptors cannot represent the relevant energy landscape, the self-consistent loop converges to a biased ensemble no matter how well it converges.","fun_headline_variants_meta":{"raw":{"variants":["MLACS: MLIP surrogate matches DFT at 50x fewer calls","MLACS variational sampling: 1 meV accuracy, 50x speedup","MLACS: DFT-free canonical sampling matches AIMD accuracy","MLACS: reweighted MLIP sampling gets 1 meV/atom at 50x faster"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001105,"raw_usage":{"total_tokens":4634,"prompt_tokens":1001,"completion_tokens":3633,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":617,"completion_tokens_details":{"reasoning_tokens":3548}},"tokens_in":617,"tokens_out":3633,"duration_ms":24964,"temperature":1.0,"reasoning_tokens":3548,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-11T11:29:36.018550+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run MLACS on a strongly anharmonic or liquid system with a known multimodal canonical distribution, for example liquid water or a high-temperature bcc metal, and compare the reweighted energy histogram, pair distribution functions, and free energy against a long AIMD trajectory and an independent thermodynamic-integration reference. If the histogram misses a second peak, the pair distribution functions deviate beyond the AIMD statistical error, or the free energy drifts by more than about 1 meV/atom, the expressiveness assumption is falsified; a cheaper diagnostic is the second-order cumulant correction $\\Delta F_{\\mathrm{int}\\to\\mathrm{AI}}$, which should stay below 1 meV/atom on a converged loop.","supporting_citations":[{"cited_title":"Castellano, F","cited_arxiv_id":null,"evidence_quote":"Defines the original MLACS method and its variational-inference derivation, which this paper generalizes into a production package."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Establishes the self-consistent harmonic approximation whose Gibbs-Bogoliubov variational structure MLACS generalizes to arbitrary linear MLIPs."}],"review_version":1}