Abstract Phaseless auxiliary-field quantum Monte Carlo (AFQMC) has in several cases been found to perform well on strongly correlated systems, including compounds containing transition-metal elements. Here, we benchmark the method for three iron–sulfur clusters ([2Fe-2S], [4Fe-4S], and the FeMo cofactor) using a hierarchy of systematically improved trial states derived from coupled cluster (CC) theory, including up to quadruple excitations relative to the Hartree–Fock (HF) reference, as well as multi-Slater trial states derived from unrestricted density matrix renormalization group wave functions. Our results reveal for these systems that, as the symmetry-broken trial is improved beyond the mean-field level, the phaseless AFQMC energy can become less accurate, and in some cases even less accurate than the underlying trial projected energy itself, displaying an inverted energy pattern that is only corrected once the trial fidelity is sufficiently high. For [2Fe-2S], we show that a decreased accuracy in the AFQMC energy can coincide with a simultaneous improvement in the trial state and the walker ensemble, as measured by their fidelity relative to the exact ground state, a behavior which is allowed by the nonvariational nature of the phaseless AFQMC energy estimator. We further find that this is not solely due to the use of spin-unrestricted trial states, as the inversion persists in [2Fe-2S] when we explicitly break the symmetry of the Hamiltonian by applying a fictitious local spin-Zeeman field. Instead, we find that the energy inversion is related to the choice of measurement trial in the estimator, where using a high-order CC trial state for measurements may introduce errors that are suppressed when the measurement wave function is restricted to lower excitation subspaces. In particular, measuring the energy with the mean-field reference while guiding the walkers with a CC trial removes the inverted pattern and improves the overall accuracy across the iron–sulfur clusters, with a possible exception for [4Fe-4S] where energies with a CCSDT guide may be less accurate than those guided by CCSD for some broken-symmetry families. Taken together, our findings suggest that the relatively accurate energies obtained with an HF trial state in these systems arise from favorable error cancellation, warranting significant caution about the reliability of phaseless AFQMC with such trials for strongly correlated transition-metal systems of this kind.
Abstract Cytochrome P450 2C19 (CYP2C19) genetic variation contributes to interindividual differences in drug metabolism, yet relating sequence-level variation to substrate-specific catalytic mechanisms remains challenging. Here, we combine a generative sequence model based on the maximum entropy (MaxEnt) framework for mutation scoring, published biochemical activity data, molecular dynamics (MD) simulations, and quantum mechanics/molecular mechanics adaptive steered molecular dynamics (QM/MM ASMD) calculations to examine CYP2C19 variant effects on S-mephenytoin oxidation. In the primary correlation analysis, MaxEnt scores showed nonsignificant negative rank trends with Vmax and CLint, whereas Km showed no corresponding rank association. CYP2C19*1A (WT) was used as the reference system, and CYP2C19*23 and CYP2C19*26 were selected as contrasting variants for atomistic analysis. MD simulations revealed variant-dependent changes in conformational sampling, correlated motions, and substrate positioning, but the population of geometries satisfying a predefined reactive distance-angle criterion did not reproduce the experimental activity trend. At the production pulling rate, the Seed 1 QM/MM-ASMD profiles showed a lower C12 oxygen-addition barrier for CYP2C19*23 and a higher barrier for CYP2C19*26 relative to CYP2C19*1A. Additional CYP2C19*1A and CYP2C19*23 starting configurations yielded a lower three-seed mean barrier for CYP2C19*23 at that pulling rate, although the barrier ranges overlapped and starting-configuration dependence was appreciable, and a slower-pulling control showed substantial rate sensitivity. The sampled NIH-shift profiles exhibited apparent low barriers in all three systems; however, substantial discrepancies between partition-level nonequilibrium-work estimators precluded establishing quantitatively converged CV2 barrier magnitudes. Thus, the calculations support qualitative differences among the sampled reaction profiles but do not quantitatively reproduce the magnitude of the experimental kinetic differences. The combined sequence, MD, and QM/MM workflow provides a framework for prioritizing variants and developing substrate-specific mechanistic hypotheses while making the sampling limitations explicit.
Abstract Universal machine learning interatomic potentials (MLIPs) have made density functional theory (DFT)-level molecular dynamics highly affordable. However, their rapid proliferation makes it prohibitively expensive to benchmark every new model against the potential of mean force (PMF), the bridge between microscopic atomic fluctuations and macroscopic thermodynamic and kinetic observables. Reweighting the trajectory of a single “source” MLIP onto a suite of target potentials offers a cheaper route, but it is hindered by the statistical collapse of exponential free energy perturbation when the phase-space overlap between potentials is low. We show that for system sizes of practical interest, this collapse is essentially unavoidable, and we quantify it with an information-theoretic diagnostic, the reweighting entropy, that indicates when reweighting can be trusted. To sidestep the collapse, we introduce a systematic, scalable reweighting framework spanning exponential reweighting, the Gaussian approximation, and the energy-only estimator. When applied to the intramolecular radical cyclization of hexenyl (17 atoms), all three estimators perform well and reproduce the target ωB97M-V/def2-TZVP PMF. However, when extended to Li+ transport through a water-solvated zeolite (601 atoms), only the energy-only estimator recovers the PMFs of a diverse suite of MLIPs from a single source simulation by discarding the noisy entropic contribution to the free energy correction. An independent, brute-force umbrella sampling of one of the most divergent targets (MACE-MATPES-r2SCAN) for this system confirms it: beyond the shape of the PMF, both the reaction and activation free energies are recovered within chemical accuracy. The same analysis reveals that the studied MLIPs partition into two clusters, one of plain generalized gradient approximation (GGA) models (PBE, PBE-sol) and one of meta-GGA and/or dispersion-corrected models (PBE-D3, r2SCAN, r2SCAN-D4), reflecting differences in their training data and DFT reference. The resulting protocol thus offers an affordable and reliable route to cross-model consensus on materials properties, even when exact reweighting is statistically intractable.
Abstract We combined alchemical transfer free energy perturbation (FEP) and replica exchange with solute tempering (REST) molecular dynamics to probe the energetics of the binding of a charged peptide to a protein target. For this purpose, we have selected the coreNLS peptide KKPKKE, the conserved fragment from the nuclear localization signal (NLS) sequence in viral capsid protein binding to nuclear transporter importin-α (impα). This peptide binds to impα, natively forming a specific pose similar to that resolved experimentally. The two main outcomes of our work are as follows. First, we further develop the FEP methodology. Alchemical transfer is a rarely used version of FEP, in which an unbound peptide is gradually decoupled, while the bound peptide is coupled. Because the net charge of the system remains unchanged, no charge-related corrections within FEP are necessary. Additionally, our FEP protocol utilizes a “laissez faire” soft-walled flat-bottom spherical restraint, which does not affect a ligand as long as it occupies the major NLS binding site in impα. We showed that combining FEP with REST molecular dynamics affords exhaustive sampling of the binding process and reduces the free energy error down to 0.4 kcal/mol. Second, we found that the coreNLS binds impα with moderate affinity of ΔGb0 = −4.6 kcal/mol. We analyzed the factors contributing to binding energetics. Direct interactions between the peptide and impα favor binding and are likely to be responsible for largely native binding pose of the coreNLS. These direct binding interactions are supplemented by the gains in water–water, ion–ion, and water–ion interactions. Both, direct and indirect, enthalpic gains are about equally important, but cannot compete with severe enthalpic losses stemming from the expulsion of water and ions from the binding interface. As a result, the net binding enthalpy is unfavorable. The factor, which ultimately “saves” coreNLS binding to impα, is a strong entropic gain occurring upon water and ion release from the binding interface. Based on our investigation, we classify the mechanism governing coreNLS binding to impα as entropic. Our isothermal titration calorimetry measurements confirm the estimates of ΔGb0 provided by FEP. By connecting with our prior studies and other reports, we propose that such binding mechanism may be applicable to other charged peptides reflecting complex interplay between binding enthalpy and entropy. Interestingly, our data suggest that the ligands binding diffusively or by adopting specific poses may utilize an entropic binding mechanism.
Abstract Access to accelerated computational methods at the electronic and atomistic scales is important in materials modeling, as it will allow more realistic system descriptions. Against this background, a systematic study of a range of electronic/atomistic methods applied to cellobiose monomers and dimers is presented, with the aim of identifying methods that achieve a significant acceleration of cellulose computations compared to hybrid DFT methods while retaining adequate accuracy. We benchmark such methods, namely the semiempirical quantum-mechanical tight-binding GFNn-xTB (n = 1,2) and G-xTB methods, the MACE-OMol-0 machine-learning potential, and the CHARMM36 force field, against three DFT references (ωB97M-V, B3LYP-D3, PBE-D3) for gas-phase cellobiose monomers and dimers, focusing on structures, energies, and vibrational properties. Cellobiose monomers and dimers, as presented in this work, do not represent the structural complexity of cellulose, but they provide tractable structural models for assessing the qualities of the selected computational models regarding cellobiose–cellobiose intermolecular interactions that are also present in cellulose. We find that (i) the chain–chain and sheet–sheet dimer interactions are of comparable magnitude, suggesting that both intermolecular dispersion and hydrogen bonding make important contributions to dimer stability, so accurate treatment of noncovalent interactions the calculations is essential; (ii) GFN2-xTB and G-xTB reproduce most conformations well, with a speed-up of two to three orders-of magnitude compared to the DFT calculations, but fail for puckered, strained conformations, artificially improving apparent dimer-energy agreement with DFT; (iii) MACE reaches near-DFT accuracy with a further order of magnitude speed-up, although in a few cases fails to reproduce the hybrid DFT conformations. This is especially true for the chain–chain dimer structure, which is also reflected in a significant deviation of the chain–chain energy-distance potential curve, which highlights the importance of system-dependent model validation (as was done here) and suggests a need for model fine-tuning for cellulose systems; (iv) CHARMM36 has some conformational and energetic limitations,but being 30,000× faster than G-xTB (itself orders of magnitude faster than DFT), remains practical for qualitative studies.
