A free energy optimizer and method of implementation thereof in molecular dynamics simulations

Abstract

The application discloses a free energy optimizer and an implementation method thereof in a molecular dynamics simulation, and belongs to the technical field of molecular dynamics simulation. The free energy optimizer comprises an MD integrator or a Langevin integrator. Based on a replica group obtained by sampling an initial molecular conformation, a force term of the integrator is replaced or an additional bias force is added to the integrator, so that a negative free energy gradient is obtained, and the integrator is used to promote optimization of the replica group along the negative free energy gradient. The optimizer and the optimization method significantly improve the optimization performance of a molecular conformation, and realize the optimization effect of the molecular conformation in a high-dimensional CV or a full-atom coordinate space.

Description

Technical Field

[0001] This invention relates to the fields of computational chemistry, molecular simulation, molecular dynamics software engineering, and artificial intelligence-assisted scientific computing, specifically to a molecular conformation optimizer based on free energy and its implementation method in a molecular dynamics integrator. Background Technology

[0002] In the fields of computational chemistry, molecular simulation, molecular dynamics software engineering, and AI-assisted scientific computing, molecular dynamics (MD) is traditionally defined as a sampling simulation problem. Its basic paradigm is: using the potential energy gradient... Driven by the Newton equations of motion or Langevin stochastic differential equations, one (or more) continuous trajectories are advanced in atomic coordinate space; under the long-time limit and ergodic assumptions, the statistical distribution of these trajectories converges to the target Boltzmann distribution. In this paradigm, the MD / Langevin integrator acts as a "time-progressive sampler with balanced distribution"; the force acting on the atom corresponds to the negative gradient of the physical force field with respect to the potential energy function.

[0003] This invention proposes a shift from the traditional paradigm—addressing the conformational sampling task through free energy optimization, and reformulating the sampling problem on the function space of the replica set / empirical distribution as the objective of minimizing variational free energy (VFE):

[0004] The deterministic update direction of the replica group derived from this VFE target is: in , For CV projection The second term of the Jacobian equation is obtained by pulling back the logarithmic density gradient in the CV space from the replica set through kernelization to the atomic space. It should be noted that in this… Under the symbolic conventions, The dimensionless FP update direction is the force received by the MD integrator in this invention. Mathematically satisfied (This holds true only under the aforementioned notation conventions); The difference in scaling factors can be absorbed into the dynamic parameters of the MD integrator, such as step size, mass, and friction coefficient. The corresponding copy group update directions have equivalent optimization semantics under the CV subspace projection. The variational property of VFE means that it has a unique global minimum in the distribution space, i.e., the target Boltzmann distribution. The force used for single-step advancement is composed of physical force and a nonlocal force term induced by the current distribution of the entire copy set. This framework transforms the conformation sampling problem into an optimization problem with a global optimum, but its optimization path (especially in atomic coordinate space) is not trivial: variational objective While it possesses a unique global minimum at the replica group level, saddle points can occur during optimization, leading to known phenomena such as mode dropping. Summary of the Invention

[0005] To overcome the shortcomings of existing technologies, this invention designs a free energy optimizer and its implementation method in molecular dynamics simulations.

[0006] This invention provides a free energy optimizer, including an MD integrator or a Langevin integrator. Based on a replica group obtained by sampling an initial molecular conformation, the integrator's force terms are replaced or an external bias force is added to obtain a negative free energy gradient, and the integrator is then used to optimize the replica group along the negative free energy gradient.

[0007] Specifically, the negative free energy gradient is obtained by linear superposition of the potential energy driving term and the entropy driving term; wherein, the potential energy driving term is constructed based on the potential energy gradient of the molecular conformation, and the entropy driving term is constructed based on the entropy contribution corresponding to the replica group distribution.

[0008] Furthermore, the force term is replaced with an action force. That is, the negative free energy gradient: , in, For atomic coordinates, , Let be the potential energy function. reverse temperature , For the projection function of collective variables, It is obtained by nonparametric gradient estimation of the replica set, which is constructed based on the kernel trick.

[0009] Furthermore, the force Includes physical force and entropy The physical force The entropy force .

[0010] Furthermore, the free energy optimizer also includes a physical force calculation module, which calculates the physical forces by calling a machine learning potential energy model, a custom force field, or a classical molecular force field. .

[0011] Specifically, the machine learning potential energy model is an interatomic potential energy function constructed based on neural networks or graph neural networks, equivariant graph neural networks, kernel machine learning methods, etc.

[0012] Specifically, the custom force field or classical molecular force field is an interatomic potential energy function constructed based on neural networks or graph neural networks, equivariant graph neural networks, kernel machine learning methods, etc.

[0013] Furthermore, the free energy optimizer also includes an entropy force calculation module, which calculates the entropy force by calling the distributed gradient estimator with the replica set as the overall input. .

[0014] Specifically, the distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

[0015] Furthermore, the free energy optimizer also includes an adaptation module to adjust the force... As the input force required for each step of the MD integrator or Langevin integrator, the copy is updated according to the update format.

[0016] Specifically, the MD integrator or Langevin integrator is Langevin, BAOAB, Leimkuhler-Matthews, Brownian (overdamped Langevin), or Nosé-Hoover Langevin.

[0017] Furthermore, the free energy optimizer also includes a scheduling module, which maintains N replicas (N greater than or equal to 1) at the system level, calculates the entropy force based on the overall replica group, and schedules the execution of physical force calculation, negative free energy gradient synthesis, and integrator update for each replica.

[0018] Furthermore, the integrator is an MD integrator or a Langevin integrator, configured in a non-intrusive, pluggable manner.

[0019] Furthermore, entropy force terms are coupled at any or more points in the integrator's force calculation, external force, bias force, callback, plug-in, or input interface. .

[0020] Furthermore, the integrator's velocity / position update, constraint subroutine, neighbor table update, temperature controller, and one or more parallel communication functions are reused.

[0021] Furthermore, the force calculation of a single replica is extended to be called in parallel by N replicas, where N is greater than 1.

[0022] Furthermore, the entropy force This corresponds to the metadynamics-type CV bias parameterization.

[0023] Furthermore, the integrator O-step noise temperature satisfies And in Under the assumption of CV space boundary self-consistency and the standard construction conditions of kernel and pullback, the steady-state distribution of the integrator corresponds to the target behavior of metadynamics-type CV bias sampling under the parameterization.

[0024] Furthermore, the entropy force is calculated by calling the distributed gradient estimator with the replica set as the overall input. The distributed gradient estimator is constructed based on the nonparametric distribution estimation of the replica group.

[0025] Specifically, the distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

[0026] Furthermore, the entropy force Extending along the CV dimension, the kernel trick replaces the CV grid, enabling bias forces to be implemented in high-dimensional CV subspaces or all-atom coordinate spaces by combining kernel approximations / sparse approximations.

[0027] Furthermore, the entropy force By generalizing the bias expression, the kernelized density estimate is a nonparametric expression that does not force the parameterization of the sum of Gaussian hills and allows for more complex distribution geometry.

[0028] Furthermore, the entropy force Following the generalization of replica origin, bias information originates from the instantaneous estimate of the current distribution of the replica group, and the increase in the number of replicas directly translates into an improvement in the accuracy of bias estimation.

[0029] Furthermore, the physical temperature in the MD is included within the Gaussian noise amplitude of the Langevin integrator's O-step. Replace with an independently adjustable optimizer temperature hyperparameter. .

[0030] Furthermore, the present invention also discloses an implementation method in molecular dynamics simulation, comprising the following steps: based on the replica group obtained by sampling the initial molecular conformation, the force term of the integrator is replaced or an external bias force is added to obtain a negative free energy gradient, and the integrator is made to advance the optimization of the replica group along the negative free energy gradient.

[0031] Specifically, the negative free energy gradient is obtained by linear superposition of the potential energy driving term and the entropy driving term; wherein, the potential energy driving term is constructed based on the potential energy gradient of the molecular conformation, and the entropy driving term is constructed based on the entropy contribution corresponding to the replica group distribution.

[0032] Furthermore, the force term of the integrator is replaced with the action force. That is, the negative free energy gradient:

[0033] in, For atomic coordinates, , Let be the potential energy function. reverse temperature , For the projection function of collective variables, It is obtained by nonparametric gradient estimation of the replica set, which is constructed based on the kernel trick.

[0034] Furthermore, the force Includes physical force and entropy The physical force is The entropy force is .

[0035] Furthermore, the physical forces are calculated by invoking machine learning potential energy models, custom force fields, or classical molecular force fields. .

[0036] Furthermore, the machine learning potential energy model is an interatomic potential energy function constructed based on neural networks or graph neural networks, equivariant graph neural networks, kernel machine learning methods, etc.

[0037] Furthermore, the custom force field or classical molecular force field is a molecular force field based on a parameterized potential energy function.

[0038] Furthermore, the entropy force is calculated by calling the distributed gradient estimator with the replica set as the overall input. .

[0039] Specifically, the distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

[0040] Furthermore, the force As the input force required for each step of the MD integrator or Langevin integrator, the copy is updated according to the update format.

[0041] Furthermore, the MD integrator or Langevin integrator is Langevin, BAOAB, Leimkuhler-Matthews, Brownian (overdamped Langevin), or Nosé-Hoover Langevin.

[0042] Furthermore, at the system level, N replicas (N greater than or equal to 1) are maintained. The entropy force is calculated based on the overall replica group, and physical force calculation, negative free energy gradient synthesis, and integrator update are scheduled and executed for each replica.

[0043] Furthermore, the integrator is an MD integrator or a Langevin integrator, configured in a non-intrusive, pluggable manner.

[0044] Furthermore, entropy force terms are coupled at any or more points in the integrator's force calculation, external force, bias force, callback, plug-in, or input interface. .

[0045] Furthermore, the integrator's velocity / position update, constraint subroutine, neighbor table update, temperature controller, and one or more parallel communication functions are reused.

[0046] Furthermore, the force calculation of a single replica is extended to parallel calls of N replicas, where N is greater than 1.

[0047] Furthermore, the entropy force This corresponds to the metadynamics-type CV bias parameterization.

[0048] Furthermore, the integrator O-step noise temperature satisfies And in Under the assumption of CV space boundary self-consistency and the standard construction conditions of kernel and pullback, the steady-state distribution of the integrator corresponds to the target behavior of metadynamics-type CV bias sampling under the parameterization.

[0049] Furthermore, the entropy force The traditional Metadynamics method for accumulating single-trajectory time accumulation of Gaussian hills on the CV grid is replaced by nonparametric kernel density estimation and instantaneous estimation of the current distribution of the replica group.

[0050] Furthermore, the entropy force is calculated by calling the distributed gradient estimator with the replica set as the overall input. The distributed gradient estimator is constructed based on the nonparametric distribution estimation of the replica group.

[0051] Specifically, the distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

[0052] Furthermore, the entropy force By generalizing the bias expression, the kernelized density estimate is a nonparametric expression that does not force the parameterization of the sum of Gaussian hills and allows for more complex distribution geometry.

[0053] Furthermore, the entropy force Following the generalization of replica origin, bias information originates from the instantaneous estimate of the current distribution of the replica group, and the increase in the number of replicas directly translates into an improvement in the accuracy of bias estimation.

[0054] Furthermore, the physical temperature in the MD is included within the Gaussian noise amplitude of the Langevin integrator's O-step. Replace with an independently adjustable optimizer temperature hyperparameter. .

[0055] Furthermore, this invention also discloses a method for implementing a free energy-based molecular conformation optimizer in an MD integrator or a Langevin integrator, comprising the following steps: Step 1: Given the current position of the integrator ,momentum With force Progress to the predetermined value ,in, This is the integration time step; Step 2: Apply force to the integrator. ,in, For atomic coordinates, , Let be the potential energy function. reverse temperature , For the projection function of collective variables, It is obtained by nonparametric gradient estimation of the replica set, which is constructed based on the kernel trick; Step 3: The integrator performs optimization as a free energy optimizer.

[0056] Specifically, the negative free energy gradient corresponding to the force is obtained by the linear superposition of the potential energy driving term and the entropy driving term; wherein, the potential energy driving term is constructed based on the potential energy gradient of the molecular conformation, and the entropy driving term is constructed based on the entropy contribution corresponding to the replica group distribution.

[0057] Furthermore, the force Includes physical force and entropy The physical force is The entropy force is Furthermore, the physical forces are calculated by invoking machine learning potential energy models, custom force fields, or classical molecular force fields. .

[0058] Furthermore, the machine learning potential energy model is an interatomic potential energy function constructed based on neural networks or graph neural networks, equivariant graph neural networks, kernel machine learning methods, etc.

[0059] Furthermore, the custom force field or classical molecular force field is a molecular force field based on a parameterized potential energy function.

[0060] Furthermore, the entropy force is calculated by calling the distributed gradient estimator with the replica set as the overall input. .

[0061] Specifically, the distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

[0062] Furthermore, the force As the input force required for each step of the MD integrator or Langevin integrator, the copy is updated according to the update format.

[0063] Furthermore, the MD integrator or Langevin integrator is Langevin, BAOAB, Leimkuhler-Matthews, Brownian (overdamped Langevin), or Nosé-Hoover Langevin.

[0064] Furthermore, at the system level, N replicas (N greater than or equal to 1) are maintained. The entropy force is calculated based on the overall replica group, and physical force calculation, negative free energy gradient synthesis, and integrator update are scheduled and executed for each replica.

[0065] Furthermore, the integrator is an MD integrator or a Langevin integrator, configured in a non-intrusive, pluggable manner.

[0066] Furthermore, entropy force terms are coupled at any or more points in the integrator's force calculation, external force, bias force, callback, plug-in, or input interface. .

[0067] Furthermore, the integrator's velocity / position update, constraint subroutine, neighbor table update, temperature controller, and one or more parallel communication functions are reused.

[0068] Furthermore, the force calculation of a single replica is extended to parallel calls of N replicas, where N is greater than 1.

[0069] Furthermore, the entropy force This corresponds to the metadynamics-type CV bias parameterization.

[0070] Furthermore, the integrator O-step noise temperature satisfies And in Under the assumption of CV space boundary self-consistency and the standard construction conditions of kernel and pullback, the steady-state distribution of the integrator corresponds to the target behavior of metadynamics-type CV bias sampling under the parameterization.

[0071] Furthermore, the entropy force is calculated by calling the distributed gradient estimator with the replica set as the overall input. The distributed gradient estimator is constructed based on the nonparametric distribution estimation of the replica group.

[0072] Specifically, the distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

[0073] Furthermore, the entropy force Extending along the CV dimension, the kernel trick replaces the CV grid, enabling bias forces to be implemented in high-dimensional CV subspaces or all-atom coordinate spaces by combining kernel approximations / sparse approximations.

[0074] Furthermore, the entropy force By generalizing the bias expression, the kernelized density estimate is a nonparametric expression that does not force the parameterization of the sum of Gaussian hills and allows for more complex distribution geometry.

[0075] Furthermore, the entropy force Following the generalization of replica origin, bias information originates from the instantaneous estimate of the current distribution of the replica group, and the increase in the number of replicas directly translates into an improvement in the accuracy of bias estimation.

[0076] Furthermore, the physical temperature in the MD is included within the Gaussian noise amplitude of the Langevin integrator's O-step. Replace with an independently adjustable optimizer temperature hyperparameter. .

[0077] Furthermore, the present invention also discloses a system for implementing the free energy optimizer, characterized in that: the system includes a processor, a memory, and a computer program stored in the memory, wherein when the computer program is executed by the processor, the system performs the above-described method.

[0078] In addition, the present invention discloses a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, is used to implement the above-described method.

[0079] Compared with existing technologies, the free energy optimizer and its implementation method in molecular dynamics simulation provided by this invention significantly improve the optimization performance and spatial applicability of systems or methods using the optimizer. Attached Figure Description

[0080] Figure 1 A schematic diagram comparing the entropy force provided for the invention with the Gaussian hills bias in the low-dimensional CV space of Metadynamics.

[0081] Figure 2 A structural diagram of the free energy optimizer device based on the MD / Langevin integrator provided for the invention.

[0082] Figure 3 A schematic diagram of a sampling process for an embodiment of molecular conformation sampling provided for the invention.

[0083] Figure 4 A comparison diagram of the initial vs. convergent distribution of the replica group in the Ala2 vacuum system under one embodiment provided for the invention.

[0084] Figure 5 A comparison diagram of the free energy planes of the present invention and Metadynamics in an Ala2 explicit aqueous solution system provided for the invention. Detailed Implementation

[0085] To enable those skilled in the art to better understand the technical solution of the present invention, the present invention will be further described in detail below with reference to specific embodiments.

[0086] Technical Specifications • Replica: One of the molecular conformation samples that participates in the sampling; N replicas are maintained in the same simulation.

[0087] • Physical force : The negative gradient of the potential energy function / force field with respect to the atomic coordinates.

[0088] • Entropy / distribution-gradient force (in , ; or the all-atomic coordinate space version ): Estimated by the kernelized density of the replica set Derived non-physical force terms. Abbreviations appearing below. This refers to the Jacobian pullback form mentioned above, which is mathematically equivalent to the previous expression.

[0089] • Negative free-energy gradient / free-energy force In this invention, the essential meaning of "force" corresponds to the actual force received by the MD integrator; to avoid confusion with free energy... positive gradient For the sake of clarity, this manual will consistently refer to the negative free energy gradient or free energy driving force. This is distinct from the concept of physical force in traditional MD.

[0090] The free energy is defined using the classical definition in physicochemical ensemble, that is, a state function characterizing the equilibrium state of a system under a given thermodynamic ensemble; in the NVT (isothermal and isochoric) ensemble, it is the Helmholtz free energy F = U - TS, in the NPT (isothermal and isobaric) ensemble, it is the Gibbs free energy G = H - TS, and so on, where U is internal energy, H is enthalpy, S is entropy, and T is temperature. Its physical meaning is consistent across ensembles: under the corresponding isothermal conditions, the thermodynamic equilibrium state of the system corresponds to the minimum value of this free energy. In the molecular conformational context of this invention, the conformational free energy is defined as F(ξ) = -β. -1ln p(ξ), where ξ = ξ(R) is the projection of atomic coordinates R onto the collective variable (CV) subspace, p(ξ) is the equilibrium marginal distribution of the conformation on this CV (corresponding to the Boltzmann marginal distribution under NVT, the Boltzmann-Gibbs marginal distribution under NPT, and so on), and β = 1 / (k_B T) is the inverse temperature; this definition is consistent with the general definition of "conformational free energy" or "free energy surface (FES)" in statistical mechanics and molecular simulation. The "free energy optimizer" described in this invention advances the update of the replica group in the atomic coordinate space along the negative free energy gradient -∇F, causing the distribution to converge to the low-energy basin of F, that is, the thermodynamic equilibrium distribution of the system.

[0091] • MD integrator: A numerical scheme that converts forces into position / momentum updates; typical implementations include Velocity-Verlet, Leapfrog, Langevin BAOAB, etc.

[0092] • Free-energy optimizer: In this invention, the MD / Langevin integrator uses the negative free energy gradient (free energy driving force). As the driving force, the update process of the replica group in the atomic coordinate space is advanced, so that the marginal distribution of the replica group in the CV subspace approaches the distribution corresponding to the target free energy surface.

[0093] • Physical temperature Define the target Boltzmann distribution and the physical target free energy surface. (in The marginal projection of the Boltzmann distribution onto the CV space, and the estimated distribution in the current iteration. Physical temperature parameters that are conceptually independent of each other; It has been explicitly stated in the free energy target Therefore, the present invention is described in The middle already implies .

[0094] • Optional optimizer temperature In an optional embodiment of the present invention, the temperature parameter of the O-step Gaussian noise term of the Langevin integrator can be named as... and physical temperature Separate settings (in speed variable) Below, the O-step noise amplitude is If momentum variable If expressed as such, then the corresponding noise amplitude is... ). By default That is, it adopts the temperature control form of traditional physical MD; this default setting is also within the protection scope of this invention. Whether or not Separation, None of these factors affect the establishment of this invention. --As long as the actual content of the force call point has been replaced with The core innovation of this invention has been implemented.

[0095] • Operating mode coverage: The operating modes of the free energy optimizer described in this invention are covered by... and Relative value continuous coverage: The corresponding non-parametric metadynamics mode: the steady-state distribution of the integrator is That is, traditional molecular simulation form + by The entropy bias potential is represented by the sampling marginal distribution of the replica group in the CV space. The sampling distribution after flattening the physical distribution corresponds to the target behavior of a type of metadynamics-type CV bias sampling under the parameterization and self-consistency assumptions. This corresponds to the SGLD mode: free energy stochastic gradient descent, where replica groups explore the free energy surface in a stochastically controlled manner and cluster in its basin.

[0096] Regarding the problems described in the background section, if general first-order optimizers such as SGD, AdaGrad, and Adam, which are based on the cumulative momentum of coordinate components, are directly used to advance the update of the replica swarm in the atomic coordinate space, two structural problems exist: 1. Lack of exploration capabilities to overcome optimization bottlenecks: First-order GD / SGD is prone to stagnation near the saddle point of the VFE target, which can lead to mode dropping. 2. Inability to directly handle atomic systems: AdaGrad / Adam et al. accumulate first / second momentum for each coordinate component separately. These accumulated momentum are not covariant under the overall translation and rotation of the molecule, which does not satisfy the SE(3) covariance required for the whole atomic system; nor does it provide a docking mechanism with MD constraint algorithms (such as bond length and bond angle constraints). In multi-copy simulations, random translation / rotation between copies makes the optimizer momentum based on coordinate components inconsistent in meaning between different copies, requiring additional coordinate system alignment processing (such as Kabsch alignment) to ensure statistical significance.

[0097] Furthermore, as a time-progressive sampler with a balanced distribution, the effectiveness of MD relies on the validity of the ergodicity assumption. However, for real-world molecular systems with significant energy barriers, long lifecycles, and rare event structures, the physical simulation time required to achieve ergodicity can often reach the order of milliseconds (ms-s), far exceeding the achievable simulation time on current GPU / HPC clusters (typically in the order of µs-ns). This is the fundamental timescale bottleneck limiting MD's direct use as a conformational sampling method.

[0098] Furthermore, Metadynamics is a representative augmented sampling method proposed to overcome the aforementioned timescale bottleneck—it accelerates trajectory exploration in the CV space by accumulating a series of Gaussian hills (history-dependent Gaussian bias potential) in the CV space as bias potential. However, this method has the following structural limitations: 1. Limited CV dimensionality: Gaussian hills accumulate on the CV grid, which typically only supports 1-3 dimensions; grid bias is not feasible in high-dimensional CV spaces. 2. Reliance on single-track time accumulation: bias must be accumulated gradually along the time history of a single track. Bias information cannot be directly shared between replicas, and parallel expansion is not possible. 3. Fixed parameterization: The bias function is pre-fixed as a sum of Gaussian functions, making it difficult to express more complex distribution geometry.

[0099] To address the aforementioned problems, this invention provides a free energy optimizer, including an MD integrator or a Langevin integrator. Based on a replica group obtained by sampling the initial molecular conformation, the integrator's force terms are replaced or an external bias force is added to obtain a negative free energy gradient, and the integrator is then used to optimize the replica group along the negative free energy gradient.

[0100] The negative free energy gradient is obtained by linear superposition of the potential energy driving term and the entropy driving term; wherein the potential energy driving term is constructed based on the potential energy gradient of the molecular conformation, and the entropy driving term is constructed based on the entropy contribution corresponding to the replica group distribution.

[0101] In this context, the MD / Langevin integrator possesses an advantage over first-order gradient descent (GD / SGD) in its ability to explore and overcome optimization bottlenecks. In machine learning, Stochastic Gradient Langevin Dynamics (SGLD), as a stochastic optimizer, can provide theoretical support for escape saddle points under appropriate smoothness, noise intensity, and step size conditions. This exploratory capability is relevant to variational free energy (VFE) objectives. This is particularly crucial in the context of optimization—the variational property of VFE means that it has a unique global minimum in the distribution space (target Boltzmann distribution), but its optimization path has a saddle point structure, which causes the mode dropping problem reported in known methods such as Boltzmann Generator; MD, as a physical instance of LD, can correspond to the exploration capability of SGLD under appropriate conditions, thus providing a principle support for alleviating the optimization bottleneck and mode dropping of FP Transformer.

[0102] This invention adjusts the force acting on atoms in the traditional MD / Langevin integrator to a complete negative free energy gradient (or free energy driving force):

[0103] (or the equivalent form directly constructed in atomic space) ),in It is obtained by nonparametric gradient estimation of the replica set, which is constructed based on the kernel trick.

[0104] Under this scheme: • The core of this invention lies in the fact that the force actually received by the integrator includes --The force that was originally physical The force formed by individual components is expanded into a synthesis of "physical force + entropy force". = + ; • In a typical implementation, existing integrator modules such as the position / momentum update process derived from the Newton equation and the random noise / friction term derived from the Langevin equation can be reused as modular implementations.

[0105] Thus, the MD / Langevin integrator is changed from a time-progressing sampler with a balanced distribution to an optimizer that advances the replica group along the negative free energy gradient.

[0106] The above settings do not depend on any specific MD integral format (e.g., Langevin BAOAB / Leimkuhler-Matthews are all applicable), nor do they depend on any specific FP Transformer class sampler. The construction details.

[0107] Specifically, the MD / Langevin integrator covers Langevin-style stochastic dynamics schemes such as Langevin, BAOAB, Leimkuhler-Matthews, Brownian (overdamped Langevin), and Nosé-Hoover Langevin; each of the above integrators contains a random perturbation term, which serves as a sampler and is set as a free energy optimizer accordingly.

[0108] In particular, the present invention emphasizes that under this configuration, the force Appearing in The free energy target has been explicitly established. The embedded parameters correspond to the physical temperature. The temperature parameter of the Langevin integrator's own O-step noise term (regardless of the name). Or named independently ) belongs to the internal hyperparameters of the integrator side, regardless of whether the integrator ultimately behaves in SGLD mode ( Stochastic gradient descent of free energy, with replica groups exploring the free energy surface and clustering in their basins in a stochastically controlled manner) or non-parametric metadynamics mode ( , copy group The distribution, i.e., traditional molecular simulation + entropy bias, corresponds to a type of metadynamics-type CV bias sampling under this parameterization and self-consistency assumption.

[0109] With the above settings, this implementation method has the following advantages: (i) Langevin Dynamics as a paradigm of stochastic optimizers with exploratory capabilities: In the fields of statistical learning and deep learning, Stochastic Gradient Langevin Dynamics (SGLD) and its underlying variants, as a class of stochastic optimizers, can, under appropriate smoothness, noise intensity, and step size conditions, allow the intrinsic Gaussian noise term to enable the optimization trajectory to break through the saddle points that are easily trapped in first-order gradient descent (GD / SGD) with a non-zero probability. This property is theoretically supported. In the context of FP Transformer free energy optimization, the VFE objective... While its variational property guarantees a unique global minimum in the distribution space (i.e., the target Boltzmann distribution), its optimization path exhibits saddle point structures, leading to the mode dropping phenomenon reported in known methods such as the Boltzmann Generator. MD, as a physical instance of LD, can correspond to the exploration capabilities of SGLD under appropriate smoothness, noise intensity, and step size conditions, providing fundamental support for alleviating the optimization bottleneck and mode dropping problem of FP Transformer.

[0110] (ii) MD as a physical instance of Langevin Dynamics: Molecular dynamics integrators (especially high-precision Langevin splitting schemes such as Langevin, BAOAB, and Leimkuhler-Matthews) are essentially numerical discretization schemes of Langevin SDE. Therefore, under the aforementioned conditions, the MD integrator can correspond to the exploratory advantages of the SGLD style and be used as an FP Transformer optimizer. A sampler with stochasticity can serve as an effective non-convex optimizer; the free energy optimizer described in this invention is implemented using an MD / Langevin sampler with stochastic perturbation terms.

[0111] (iii) Applicable to general optimization problems with subspace distribution optimization + fast and slow degree of freedom separation structure, not limited to all-atomic and molecular systems: The role redefinition described in this invention also holds in a broader optimization context -- as long as a low-dimensional subspace (i.e., CV space) exists. The distribution optimization objective is on the original space. The existence of a separation between fast and slow degrees of freedom (CV being the slow variable and non-CV being the fast variable) allows the integrator described in this invention to serve as an effective optimizer for this objective. Atomic and molecular systems (CV = main chain dihedral angle; fast variables = bond vibrations, angular bending, solvent motion) are only one typical application scenario of this invention, but this embodiment is not limited to this.

[0112] Specifically, the free energy optimizer of the present invention further includes: The physics force calculation module calculates the physical forces by calling machine learning potential energy models, custom force fields, or classical molecular force fields. Among them, machine learning potential energy models are interatomic potential energy functions constructed based on neural networks or graph neural networks, equivariant graph neural networks, kernel machine learning methods, etc., such as: ANI, MACE, NequIP, Allegro, DeepMD, TorchMD-Net, etc.; custom force fields or classical molecular force fields are molecular force fields based on parameterized potential energy functions, such as: AMBER, ff14SB, CHARMM, GAFF, OPLS, etc.

[0113] Entropy force calculation module: Calculates entropy force by calling the kernelized distribution gradient estimator with the replica set as input. The kernelized distribution gradient estimator takes the replica set as the overall input and employs the kernel trick to construct the kernelized density estimate and calculate the entropy force. .

[0114] Negative free energy gradient synthesis module: and Synthesized in accordance with the manner described in Sovereignty 1 .

[0115] Integrator adapter module: converts the applied force As the input force required for each step of the MD integrator or Langevin integrator, the copy update is advanced according to the update format. The MD integrator or Langevin integrator is Langevin, BAOAB, Leimkuhler-Matthews, Brownian (overdamped Langevin), or Nosé-Hoover Langevin, etc.

[0116] Multi-replica scheduling module: Maintains N replicas at the system level, where N is a positive integer greater than or equal to 1. Calculates entropy force based on the overall replica group and schedules physical force calculation, negative free energy gradient synthesis, and integrator update for each replica separately.

[0117] In another embodiment of the present invention, the free energy optimizer described herein can be directly adapted non-intrusively and plug-inly based on existing MD integrators. In a typical implementation, entropy force terms can be coupled at any one or more of the force calculation, external force, bias force, callback, plugin, or integrator input interface, and modular implementations such as velocity / position updates, constraint (algorithm) subroutines, neighbor tables, temperature controllers, and parallel communication can be reused. The specific selection of coupling points, modification of data structures and parallel communication methods, fine-tuning of constraint call order, and docking methods for temperature control implementation are adapted by the implementer according to the target MD integrator. This plug-in adaptation enables the free energy optimizer of this invention to: (i) be constructed in the case of CV, kernel and Jacobian pullback according to SE(3) invariants / covariants -- inheriting SE(3)-equivariance -- position / momentum updates are written on physical coordinates and physical momentum, covariant under molecular translation and rotation, so that it does not need to rely on an additional global alignment mechanism in multi-copy parallel sampling (different from the general optimizer AdaGrad / Adam that accumulates momentum according to coordinate components); and (ii) be able to handle high-frequency degree of freedom constraints such as bond length, bond angle, and rigid water through the constraint projection / constraint process of the existing MD integrator -- the force replacement and position / velocity constraint projection process of this invention are orthogonal to each other, and the constraint processing task can be undertaken by the constraint subroutine of the existing MD integrator.

[0118] The free energy optimizer can be implemented using a non-intrusive plug-in approach based on typical MD integrators such as Langevin and BAOAB. As an example, typical implementation elements are as follows (implementers can adapt them according to the actual interface and data structure of the target MD integrator): • Entropy-force coupling interface: In a typical implementation, entropy-force terms can be coupled at any or more of the force calculation, external force, bias force, callback, plugin, or integrator input interface. This makes the force actually received by the integrator in each step correspond to... The specific coupling point is selected by the implementer based on the target integrator adaptation; • Reuse of integrator main loop and system modules: Directly reuse existing functional modules of the MD integrator such as velocity update, position update, constraint subroutine, neighbor table update, temperature control, and parallel communication; specific modification methods, data structure adjustments, parallel communication modifications, constraint call order, temperature control docking methods, etc. can be adapted by the implementer according to the target MD integrator. • Multi-replica scheduling: Extends the force calculation of a single replica to N replicas in parallel. In MD integrators such as GROMACS / AMBER that already support multi-replica simulation (such as replica-exchange), their multi-replica infrastructure can be directly reused.

[0119] Inherited Attribute (A): SE(3)-equivariance (conditional representation) -- In the case that CV, kernel and Jacobian pullback are constructed according to SE(3) invariants / covariants, the position / momentum update of the MD integrator is written on the physical coordinates and physical momentum, and is covariant to the overall molecular translation and rotation; negative free energy gradient / free energy driving force Depend on (SE(3) covariance) and (The covariance of SE(3) is maintained by combining the SE(3) invariant / covariant CV and the corresponding invariant / covariant kernel construction) together, thus maintaining the SE(3) covariance as a whole. Therefore, the free energy optimizer described in this invention does not need to rely on additional global alignment mechanisms (such as Kabsch alignment or additional COM / global rotation constraints) when multiple copies are executed in parallel, and can maintain physical symmetry and the statistical meaning of the optimizer; in contrast, general optimizers such as AdaGrad / Adam that accumulate momentum according to coordinate components are not covariant under the overall rotation of atoms. The specific validity of SE(3) covariance depends on the construction method of CV / kernel / pullback; when the implementer adopts the SE(3) invariant / covariant CV and kernel construction (typical construction methods include dihedral angle, internal distance, graph kernel reduced by symmetry, etc.), the inherited property (A) described in this invention is valid.

[0120] Inheritance Attribute (B): Constraint Algorithm Compatibility (Conditional Description) -- The force replacement of the free energy optimizer described in this invention is orthogonal to the constraint projection process of the existing MD integrator's position / velocity. Therefore, the force replacement can be processed by the constraint algorithm in the existing MD integrator (any constraint algorithm implementation, such as SHAKE, RATTLE, SETTLE, LINCS, etc.) through its constraint projection / constraint process. Specific constraint call order, whether to insert constraint projection before or after force calculation, and other implementation details can be adapted by the implementer according to the target MD integrator. The specific form of the constraint algorithm, the selection of constraint objects (which degrees of freedom are constrained), and whether the implementer introduces constraints are not limitations of this invention.

[0121] The two inherited properties described above enable the free energy optimizer of this invention to be applicable to (fully) atomic complex systems—something that general first-order optimizers (AdaGrad, Adam, etc.) cannot directly achieve.

[0122] In another embodiment of the present invention, the entropy force based on the kernelized distribution gradient of the replica group described herein... exist Under the default parameterization, it is equivalent to a non-parametric metadynamics implementation—that is, replacing the time accumulation of Gaussian hills on the CV grid with a non-parametric density estimation using kernel tricks; and extending it to a form that is not limited by traditional 1-3 dimensional CV grids (grid CV), and can be implemented in high-dimensional CV or all-atom coordinate space by combining kernel approximation / sparse approximation. It breaks through the three limitations of traditional metadynamics 1-3 dimensional CV, parameterized Gaussian, and single trajectory dependence, and can be accelerated by copy parallelism.

[0123] The entropy force described in this invention There is a parameterized correspondence between this invention and a class of metadynamics-type CV bias sampling methods, and this invention constitutes a nonparametric / high-dimensional / replica parallel generalization of metadynamics along three independent axes (CV dimension, bias expression, and replica parallelism). Specifically: when the Langevin integrator O-step noise temperature satisfies... (That is, using the default temperature control of traditional physical MD), and in Under the assumptions of marginal self-consistency in CV and the standard construction of the kernel and pullback, the steady-state distribution of the integrator corresponds to the target behavior of a type of metadynamics-type CV bias sampling under the parameterization, namely the non-parametric metadynamics mode described in this invention; this is demonstrated in three parts below. It should be noted that traditional metadynamics itself has many variants (standard metadynamics, well-tempered metadynamics, multiple-walker metadynamics, bias-exchange metadynamics, etc.), and the target distribution form and bias accumulation strategy of each variant have subtle differences; the correspondence described in this invention refers to the target distribution corresponding to the aforementioned type of metadynamics-type CV bias sampling under the parameterization and self-consistency assumptions, and does not claim to be identical to all metadynamics variants.

[0124] (i) The present invention corresponds to a parameterization of a type of metadynamics-type CV bias sampling (non-parametric metadynamics mode): negative free energy gradient Rewritten as effective potential Gradient identification ,but This corresponds to the standard form of "physical force + CV bias force" used in a type of metadynamics-type CV bias sampling.

[0125] when At that time, the steady-state distribution of the Langevin equation is:

[0126] Furthermore, under the assumption of marginal self-consistency of CV (i.e., assuming...) (This holds true when the replica group distribution reaches a self-consistent fixed point), thus the sampling marginal distribution of the replica group in the CV space satisfies:

[0127] in This represents the marginal projection of the Boltzmann distribution into the CV space. This flattened sampling distribution corresponds to the target behavior of a class of metadynamics-type CV bias sampling. Note: This correspondence is established based on... CV marginal self-consistency Kernel estimation achieves the assumptions of fixed points; when the above assumptions are not strictly true, the correspondence is approximate and does not affect the working mode definition of the non-parametric metadynamics implementation in engineering described in this document.

[0128] (ii) This invention extends the nonparametric / high-dimensional / replica parallelism of Metadynamics: the entropy force described in this invention Furthermore, instead of the traditional Metadynamics method of accumulating single-trajectory temporal accumulation of Gaussian hills on the CV grid, the nonparametric kernel density estimation and the instantaneous estimation of the current distribution of the replica set are used. The entropy force is calculated by calling the distribution gradient estimator with the replica set as the overall input. The distributed gradient estimator is constructed based on the nonparametric distribution estimation of the replica group. It employs the kernel trick to construct a kernelized density estimate, thus generalizing Metadynamics along three independent axes: • Extension along the CV dimension: The kernel trick replaces the CV grid, allowing bias forces to be implemented in high-dimensional CV subspaces or all-atomic coordinate spaces by combining kernel approximations / sparse approximations, without being limited by traditional 1-3 dimensional gridCV, breaking through the engineering limitations of traditional Metadynamics 1-3 dimensional CV; • Generalization along the bias expression form: Kernelized density estimation is a nonparametric expression that does not force the parametric form of the Gaussianhills sum, and can accommodate more complex distribution geometry; • Generalization along the source of replicas: bias information comes from the instantaneous estimate of the current distribution of the replica group, and the increase in the number of replicas directly translates into an improvement in the accuracy of bias estimation; traditional metadynamics must accumulate hills along a single trajectory over time, and the number of replicas and the rate of bias accumulation are not directly coupled.

[0129] (iii) The protected axes claimed by this invention: The core technical contribution of this invention, which distinguishes it from traditional metadynamics, lies in the generalization along three independent axes: non-parametric / high-dimensional CV / replica parallelism, which is not limited by traditional 1-3D grid CV. The scope of protection claimed by this invention corresponds to the non-parametric / high-dimensional / replica parallelism generalization itself, and the correspondence between this generalization and a type of metadynamics-type CV bias sampling under the parameterization and self-consistency assumptions. This constitutes the technical contribution of the non-parametric metadynamics implementation described in this invention in a patent sense.

[0130] In another embodiment of the present invention, the physical temperature, which is the default value used in traditional physical MD, is incorporated into the Gaussian noise amplitude of the Langevin integrator O-step. Replace with an independently adjustable optimizer temperature hyperparameter. Taking the BAOAB integrator as an example, the O-step can be rewritten as:

[0131] The default (Using traditional physical MD temperature control), at this point, the steady state of BAOAB is: That is, the traditional molecular simulation form + entropy bias potential, which corresponds to a type of metadynamics-type CV bias sampling (i.e., non-parametric metadynamics mode) under the aforementioned parameterization and self-consistency assumptions; if Then regardless Greater than or less than The integrators all exhibit free energy (stochastic) gradient descent—that is, the SGLD mode; its intrinsic Gaussian noise term, with controlled randomness, assists the replica group in exploring the free energy surface and eventually converging in its basin. It should be clarified here that the SGLD mode belongs to free energy stochastic gradient descent, not deterministic descent—its noise temperature... The non-zero values ​​correspond to non-zero random perturbations, giving the optimized trajectory the ability to explore saddle points. Its technical effect is as follows: • Parameterization switchability: Implementers can switch through a single hyperparameter Continuous adjustments, in non-parametric metadynamics mode ( ) and SGLD mode ( The integrator can be switched between these modes by adjusting a single temperature hyperparameter, and in principle, no overall changes to the integrator structure, force field, or topology configuration are required. • Application Extension: In application scenarios requiring rapid aggregation of replica groups within a free energy basin (such as conformational representation selection, structure refinement, and free energy minimization), a lower value can be selected. Accelerate convergence; in application scenarios requiring sampling according to a non-parametric metadynamics target distribution (i.e., corresponding to the traditional augmented sampling target distribution), maintain the default settings. That's all.

[0132] Accordingly, this invention also discloses a method for implementing the above-mentioned free energy optimizer in molecular dynamics simulations, see the appendix to the specification. Figure 2Based on the replica group obtained by sampling the initial molecular conformation, the force term of the integrator is replaced or an external bias force is added to obtain the negative free energy gradient, and the integrator is made to advance the optimization of the replica group along the negative free energy gradient.

[0133] Specifically, further limitations of the method can be found in the embodiments shown above.

[0134] In another embodiment of the present invention, a molecular conformation optimizer based on free energy and its implementation method in an MD integrator are disclosed, comprising the following steps: Step 1: Force interface diameter of MD integrator Any MD / Langevin integrator can be set at the numerical implementation level to: given the current position ,momentum With force Progressing according to the predetermined numerical plan to .in, This represents the integration time step.

[0135] The differences between different MD integrators (Langevin, BAOAB, etc.) lie only in how they extract... Advance to Numerical schemes; they share the same force interface. .

[0136] The core approach of this invention is to ensure that the force actually received by the integrator (expressed as a force term) corresponds to the complete negative free energy gradient (free energy driving force). The numerical schemes, constraint processes, and temperature control of existing MD integrators can be reused as modular implementations.

[0137] Step 2: Replace the content of the force In traditional molecular dynamics (MD), the integrator experiences physical forces:

[0138] In this invention, the force received by the integrator is replaced by the complete negative free energy gradient (free energy driving force):

[0139] in: • First item Calculated from existing MD force fields (consistent with physical force calculations in traditional MD); • Second item The entropy term introduced in this invention is estimated by kernel density from the current distribution of the replica set. To obtain, that is, The gradient is obtained from the nonparametric distribution gradient estimation of the replica set, which is constructed based on the kernel trick; when constructed in the CV space, it is obtained by... Jacobian is pulled back to atomic coordinate space—that is... .

[0140] exist Steady-state conditions reached after iterative convergence along the replica group Below, in atomic coordinate space That is, the marginal distribution of the replica group in the CV subspace approaches the distribution corresponding to the target free energy surface.

[0141] Step 3: Implementation of the MD / Langevin integrator as a free energy optimizer This step uses the Langevin BAOAB integrator as a representative implementation to illustrate the specific update method of the free energy optimizer described in this invention. It should be clarified first that, as an optimizer implementation, the integrator described in this invention must include a Langevin friction / noise term to maintain SGLD-style random perturbations; purely deterministic Newton integrals without random perturbation terms (such as Velocity-Verlet-like formats uncoupled with Langevin terms) lose their SGLD-style exploration capability and are not a recommended implementation of the free energy optimizer described in this invention.

[0142] 3a. Langevin BAOAB Integrator BAOAB is a high-precision Langevin SDE numerical discretization scheme proposed by Leimkuhler-Matthews (temperature parameters correspond to physical temperatures). = 1 / β):

[0143]

[0144]

[0145]

[0146]

[0147] The improved version of this invention (the force is replaced by a negative free energy gradient / free energy driving force) O-step temperature parameters remain unchanged in the original MD implementation, i.e. ):

[0148]

[0149]

[0150]

[0151]

[0152] Under this default implementation, the integrator's steady state is ,Right now The non-parametric metadynamics mode (corresponding to the traditional molecular simulation form + sampling of entropy bias potential energy) is used.

[0153] Optional implementation: Without changing the above BAOAB algebraic structure, the O-step noise temperature parameter can be changed from... Rename to Independent Hyperparameters and allow :

[0154] when When (default setting), the integrator is equivalent to the above non-parametric metadynamics form; when In this mode, the integrator exhibits free energy (stochastic) gradient descent (SGLD mode), where its intrinsic Gaussian noise term, with controlled randomness, assists the replica group in exploring the free energy surface and eventually converging in the free energy surface basin. This optional implementation further expands the operating mode coverage of the free energy optimizer described in this invention.

[0155] 3b. Other optional Langevin numerical discretization schemes Besides BAOAB, the following Langevin SDE numerical discretization scheme, while retaining the Langevin friction / noise term, can also be used as an optional implementation of the free energy optimizer described in this invention: • Velocity-Verlet + Langevin friction / noise term; • Leapfrog + Langevin friction / noise term; • Other variants of the Leimkuhler-Matthews split (OBABO, ABOBA, etc.); • Brownian dynamics / Smoluchowski integrator (a discretization scheme without inertia under the overdamped limit).

[0156] Under any of the above integration formats, the embodiments of the present invention can, in a typical implementation, couple the entropy force term at the force calculation, external force, bias force, callback, plugin, or integrator input interface. This can be achieved by implementing specific modifications, such as data structure adjustments, parallel communication modifications, constraint call order, and temperature control docking methods, which can be adapted by the implementer according to the target MD integrator.

[0157] 3c. Inheritance of constraint algorithm compatibility (conditionalization) The force replacement of the free energy optimizer described in this invention is orthogonal to the constraint projection process of position / velocity. Therefore, the force replacement can be processed by the constraint algorithm in the existing MD integrator (either SHAKE, RATTLE, SETTLE, LINCS, etc.) through its constraint projection / constraint process. The specific constraint call order, constraint insertion position, constraint algorithm implementation details, etc., are adapted by the implementer according to the target MD integrator. This constraint compatibility is on par with the SE(3)-equivariance (also conditionalized in the CV / kernel / pullback construction method). The internal implementation details of the specific constraint algorithm, the selection of constraint objects (which degrees of freedom are constrained), and whether the implementer introduces constraints are not limitations of this invention.

[0158] 3D MD Integrator Layer Plug-in Adaptation At the MD integrator level, the plug-in adaptation described in this invention is implemented by inserting an entropy force component at the force calculation return point.

[0159] The negative free energy gradient corresponding to the force is obtained by linear superposition of the potential energy driving term and the entropy driving term; wherein the potential energy driving term is constructed based on the potential energy gradient of the molecular conformation, and the entropy driving term is constructed based on the entropy contribution corresponding to the replica group distribution.

[0160] Corresponding to the above method, the present invention also proposes a computing system, electronic device, or computer-readable storage medium for implementing the free energy optimizer. The system includes one or more processors and a memory storing program instructions that, when executed by the processor, cause the processor to implement the above method.

[0161] The scenarios in which the solution is applicable include, but are not limited to: • Implement FP Transformer-class replica swarm sampling by directly adjusting existing MD integrators such as AMBER, GROMACS, OpenMM, NAMD, and LAMMPS using a non-intrusive plug-in method; • Deploy the free energy optimizer as an independent module in a self-developed MD integrator or a differentiable MD framework (such as JAX-MD, TorchMD, etc.); • By utilizing the SE(3)-equivariance property of the MD integrator, physical symmetry is maintained in multi-copy parallel sampling, thus eliminating the need for additional global coordinate system alignment processing; • As a high-dimensional, scalable, non-parametric, and replicable parallel version of traditional low-dimensional CV space history-dependent Gaussianhills / bias sampling methods (which can be implemented internally by grid or equivalent data structures), such as Metadynamics.

[0162] Two specific embodiments of the technical solution described above are as follows: Example 1: Alanine dipeptide vacuum system Subspace projection path sampling process using the alanine dipeptide (Ala2) system as an example. See appendix. Figure 3The present invention provides a schematic diagram of the sampling process in the alanine dipeptide (Ala2) system. The left side shows N initial Ala2 copies; the middle "× N cycle" box displays two parallel processing paths: the upper path sequentially performs subspace projection $\xi: R \to (\varphi, \psi) \in \mathbb{T}^2$, periodic KDE density estimation $\hat p(\theta)$, and CV spatial distribution gradient $\nabla_\theta \log \hat p$; the lower path calculates the potential energy $U(R)$ in atomic coordinate space and the potential energy gradient $-\beta \nabla_R U$. The two paths are coupled at the Jacobian pullback node $J^T$, and through linear superposition, form a unified update direction $\varphi(R) = -\beta \nabla_R U - J^T \nabla_\theta \log \hat p$, which is then processed by Langevin. The free energy optimizer advances the copy along this update direction; after iterative convergence, the marginal distribution of $(\varphi, \psi)$ on the right is obtained, which approaches the distribution corresponding to the target free energy surface (mainly low-energy basins correspond to the α-helix and β-sheet conformation regions of Ala2); see appendix. Figure 4 The Langevin free energy optimizer of this invention has a set of replicas before and after execution in the Ramachandran subspace. The distribution comparison diagram is shown in the image. The left panel, "Initial Samples," displays the initial distribution of the copies obtained from a uniform prior random initialization, scattered across the torus. Global domain; the "Sampling Results" panel on the right shows the replica distribution obtained after repeatedly running the Langevin BAOAB free energy optimizer described in this invention until convergence. It can be seen that the replicas are clustered in the standard low-energy conformation basin regions of the Ala2 system, such as α-helix and β-sheet. The background grayscale contour lines are the reference free energy surface.

[0163] This embodiment uses the Alanine dipeptide (Ala2) system as an example to demonstrate the complete sampling process when the subspace projection function ξ simultaneously possesses both non-injectivity (mapping from a high-dimensional atomic space to a low-dimensional CV space) and global non-bijectivity (periodic toroidal mapping), serving as an example to support the implementation of the subspace projection function. Ala2 is a standard benchmark system for studying main-chain conformation sampling in computational chemistry, and its main-chain dihedral angles (φ, ψ) are classic low-dimensional conformational representations. A schematic diagram of the sampling process in this embodiment can be found in the appendix to the specification. Figure 3 The sampling experiment verification results are shown in the appendix. Figure 4 .

[0164] A. Define the potential energy function and initialize the copy: Select the atomic potential energy function of the Alanine dipeptide (22 atoms) system. U(R) ,in R ∈ ( n = 22 ( ) represents all-atomic Cartesian coordinates. The potential energy function is given by the classical all-atomic force field (e.g., AMBER ff14SB). The number of conformational samples is set to be... N = 512 Sampling temperature T = 300 K (corresponding to the reverse temperature) Randomly initialize from a Gaussian distribution. N The atomic coordinates of each copy are used as the initial copy set. The distribution of the initial copies in the (φ, ψ) subspace is shown in the appendix. Figure 4 .

[0165] B. Extract collective variables: atomic coordinates for each copy The dihedral angle (φ, ψ) of the main chain is extracted as a collective variable through a subspace projection function: These dihedral angle pairs constitute a two-dimensional periodic CV space (torus). ).

[0166] The ξ satisfies the following: for any R, the output is unique (φ, ψ), making it a single-valued mapping; and ξ is not a bijective mapping—neither injective (3n = 66 dimensions reduced to 2 dimensions, multiple atomic conformations can correspond to the same (φ, ψ), nor globally invertible (…). Values ​​are taken from the torus. (It is periodic and there is no global inverse mapping) -- therefore ξ is the subspace projection function described in this invention.

[0167] C. Constructing kernel functions and distributed gradients in periodic CV space: Defining periodic distance and periodic Gaussian kernel functions. Based on replica set projection in (φ, ψ) space Construction of nucleation density estimation And calculate the log density gradient in the CV space. In this embodiment, the bandwidth employs the median technique: the median of the periodic distances between all replica pairs in the CV space is taken, and after appropriate scaling, it is used as a fixed bandwidth parameter. h .

[0168] D. Calculate the potential energy gradient in atomic space: Calculate the potential energy function value in the atomic coordinate space for each copy. and its gradient .

[0169] E. Projection Mapping Coupling and Update Direction Construction: Computing the Jacobian Matrix of the Subspace Projection Function The distributed gradient in the CV space is coupled back to the atomic coordinate space through Jacobian pullback, and linearly superimposed with the potential energy gradient term in the atomic coordinate space to form a unified update direction for each replica in the atomic coordinate space.

[0170] F. Optimizer Update: A copy update is performed using Langevin dynamics as the free energy optimizer. The "force" acting on the atoms represents the complete update direction described above. φ (Including the Jacobian pullback of the atomic space potential energy gradient term and the CV space distribution gradient term), the Langevin integrator acts as a free energy optimizer along... φ Proceed to the next stage of the dungeon.

[0171] G. Iterative Convergence: Repeat steps B through F until convergence. The convergence criterion is that the marginal distribution of the replica set in the (φ, ψ) space tends to stabilize—the change in the (φ, ψ) 2D histogram between adjacent iteration windows is below a preset threshold, or the corresponding free energy surface estimate tends to stabilize. After convergence, the marginal distribution of the replica set in the (φ, ψ) space approaches the distribution corresponding to the target free energy surface, while the atomic coordinates maintain a reasonable geometric configuration under potential energy constraints, ultimately obtaining a set that satisfies the above conditions. N A set of Ala2 conformation samples. The distribution of the copies in the (φ, ψ) subspace after convergence is shown in the appendix. Figure 4 It can be seen that the replicas are concentrated in several low-energy basin regions corresponding to the target free energy surface (corresponding to the main dihedral conformation regions such as α-helix and β-sheet in the Ala2 system), and are associated with the attached Figure 3 The process shown is consistent with the schematic diagram.

[0172] Following the system and process defined above, this section focuses on explaining the implementation of the Langevin integrator in the role of the free energy optimizer: A. System Definition • System: Alanine dipeptide (22 atoms), all-atomic Cartesian coordinates , ; • Force field: Classical all-atomic force field (such as ff14SB used in this example); • Collective variable: Main chain dihedral This constitutes a two-dimensional periodic toroidal space. ; • Number of copies Sampling temperature (correspond ); • Obtain an initial set of replicas by random initialization from the prior distribution.

[0173] B. Implementation Process Following the subspace projection path defined in steps B to F above—extracting the main chain dihedral angle as a collective variable, constructing a kernel function and distribution gradient in the periodic CV space (bandwidth using the median trick), calculating the potential gradient in atomic space, and coupling the CV space distribution gradient back to atomic space via Jacobian pullback to form a unified update direction.

[0174] In step F above, Langevin dynamics is used as the free energy optimizer to perform the copy update—this is precisely the implementation of the MD / Langevin integrator as the FP Transformer free energy optimizer described in this invention. The implementation details of this step are further supplemented as follows: • Force replacement: The "force" received by the integrator at each step determines the complete update direction. The corresponding negative free energy gradient (free energy driving force) rather than traditional physical forces Single item; • Plug-in adaptation: In this embodiment, the entropy force component is injected at any coupling point such as force calculation / external force / bias force / callback / plugin / integrator input interface under a typical implementation; existing Langevin integrator modules such as position / velocity update, O-step noise term, constraint subroutine, and neighbor table can be reused as modular implementations, and the specific modification method is adapted accordingly to the target MD integrator; • O-step noise temperature: Use from Gradually reduce (annealing) to This corresponds to the SGLD mode (stochastic gradient descent of free energy) of simulated annealing.

[0175] C. Numerical verification results (see) Figure 4 ) Figure 4 The numerical verification results of this embodiment are presented. The Initial Samples panel on the left shows the randomly initialized copies. The initial distribution in the subspace; the Sampling Results panel on the right shows the replica distribution after repeatedly executing the above process until convergence. It can be seen that the replicas are clustered in the standard low-energy conformation basin region of the Ala2 system. -helix、 -sheet, etc.), consistent with the basin structure of the target free energy surface.

[0176] Example 2: Alanine dipeptide explicit aqueous solution system compared with Metadynamics This embodiment replicates the procedure of Example 1 on an Alanine dipeptide explicit aqueous solution system and compares it with the Metadynamics reference results in the same system.

[0177] A. System Definition • System: Alanine dipeptide dissolved in an explicit water molecule box; • Force field: AMBER ff19SB (solute) + TIP3P (solvent); • Collective variable: Dihedral angle of solute backbone .

[0178] B. Implementation Process The implementation process is consistent with Example 1: CV extraction only acts on the dihedral angles of the solute backbone, and the solvent degree of freedom, as a non-CV degree of freedom, evolves naturally under the constraint of the physical force field; force replacement (in (Provided by the aforementioned kernelized distributed gradient estimator via Jacobian pullback) In a typical implementation, it can be coupled at any or more points in force calculation, external force, bias force, callback, plugin, or integrator input interface; existing Langevin integrator modules such as main loop, constraint subroutines, and neighbor tables can be reused as modular implementations, with specific modifications adapted to the target MD integrator.

[0179] C. Numerical verification results (see) Figure 5 ) • (a) Panel: Ramachandran free energy surface obtained by the Langevin free energy optimizer of this invention (simulation duration 70 ps); • (b) Panel: Metadynamics reference results in the same system (simulation duration 50 ns); • Both panels share the free energy colorbar; the conformational basin region and the transitional region are visually highly consistent, and the location and relative depth of the conformational basin are well restored.

[0180] D. Effect Comparison • Reduced simulation time: Under conditions where the accuracy of the conformational basin and transition state is comparable, the Langevin free energy optimizer described in this invention requires a single-copy propagation time / reporting physical simulation time of 70 ps (multiple copies executed in parallel, number of copies N = 512) to achieve the target sampling accuracy. Traditional Metadynamics requires a single-trajectory physical simulation time of 50 ns under the same system, a difference of over 100 × 100% under the aforementioned statistical caliber. This difference stems from the parallel provision of bias information by the instantaneous distribution of the copy group, avoiding the serial dependency of Metadynamics in accumulating bias along a single trajectory according to time history—corresponding to the aforementioned parallel generalization relationship of copies.

[0181] Notes: (i) The 70 ps mentioned above refers to the physical simulation propagation time of a single replica in the implementation of this invention, and is not the length of the aggregated replica trajectory (the aggregated length is N × 70 ps, ​​which needs to be converted according to the number of replicas N); (ii) The above data does not correspond to the actual wall-clock time consumption -- the wall-clock time consumption is related to the specific hardware configuration, force field implementation, parallel communication strategy, number of replicas N, etc., and needs to be converted according to the actual hardware conditions. It is not within the scope of comparison in this embodiment.

[0182] • Consistency of plug-in adaptation: The modification method of the existing MD integrator in this embodiment is consistent with that in embodiment 1 -- in a typical implementation, entropy force components can be injected at any or multiple coupling points such as force calculation / external force / bias force / callback / plugin / integrator input interface; the main loop, constraint subroutine, long-range electrostatics, neighbor table and other modules of the existing MD integrator can be reused as modular implementations, and the specific adaptation is made according to the target MD integrator.

[0183] • Advantages of SE(3)-equivariance in solvent-containing systems: In high-dimensional atomic spaces containing explicit solvents, the SE(3) covariance of the Langevin integrator allows the optimizer described in this invention to maintain physical symmetry without Kabsch alignment or global degree-of-freedom constraints - an advantage that is not easily obtained in conventional optimizers that accumulate momentum by coordinate components (AdaGrad / Adam).

[0184] The present invention has the following advantages over the prior art: 1. Non-intrusive plug-in reuse of the MD integrator ecosystem: Enables FP Transformer class replica group sampling to be directly modified and deployed in a plug-in manner on any existing MD tool such as AMBER, GROMACS, OpenMM, NAMD, LAMMPS, etc., directly reusing a large number of engineering capabilities such as force fields, constraints, neighbor tables, GPU / HPC optimization, and inheriting the parallel communication, I / O and other infrastructure of existing MD software.

[0185] 2. Mathematical equivalence: After replacing the force with the complete negative free energy gradient (free energy driving force), the replica group in... The marginal distribution at the self-consistent fixed point corresponds to the target free energy surface distribution and is equivalent to the FP Transformer principal curve; its The steady-state behavior under the default parameterization corresponds to a target distribution form of a type of metadynamics-type CV bias sampling, namely the non-parametric metadynamics mode, under the parameterization and self-consistency assumptions.

[0186] 3. SE(3)-equivariant optimizer characteristics (conditionalization): In the case of CV, kernel, and Jacobian pullback constructed according to SE(3) invariants / covariants, the MD / Langevin integrator is covariant for the overall translation and rotation of atoms. Multi-copy simulation does not need to rely on additional global alignment mechanisms (such as Kabsch alignment, COM, or global rotation constraints). Compared with general optimizers such as AdaGrad / Adam that accumulate momentum according to coordinate components, it has a structural advantage - enabling the optimizer described in this invention to directly handle (full) atomic complex systems.

[0187] 4. High-dimensional non-parametric generalization of Metadynamics-like methods: The entropy force described in this invention constitutes a non-parametric generalization of traditional low-dimensional CV space history-dependent Gaussian hills / bias augmented sampling methods (which can be implemented internally by a grid or equivalent data structure) along three independent axes: CV dimension, bias expression, and replication parallelism. When the kernel density estimation or entropy force construction adopts an expression equivalent to the CV space-time cumulative Gaussian hills, the entropy force of this invention can correspond to / can be realized as a metadynamics-type bias force.

[0188] 5. Integrator Independence: The role redefinition described in this invention is independent of the specific Langevin numerical discretization scheme and can be applied to existing / future Langevin-like integral formats such as Langevin BAOAB (recommended), Leimkuhler-Matthews, and Brownian / Smoluchowski integrators.

[0189] 6. Inheritance of Constraint Algorithm Compatibility (Conditionalization): The force replacement of the free energy optimizer described in this invention and the constraint projection process of position / velocity are orthogonal to each other. It can be processed by the constraint algorithm provided by existing MD software (any constraint algorithm such as SHAKE, RATTLE, SETTLE, LINCS, etc.) through its constraint projection / constraint process; this property is parallel to the (conditionalization) inheritance attribute of SE(3)-equivariance. The specific form of the constraint algorithm, the selection of constraint objects, the order of constraint invocation, and the specific adaptation method are not limited by this invention.

[0190] 7. Optional dual-temperature parameterization extension: Implementers can optionally extend the Langevin integrator O-step noise temperature parameter from... Rename to Independent Hyperparameters Thus, by adjusting a single hyperparameter, it is possible to achieve non-parametric metadynamics mode ( , corresponding to steady state Under the parameterization and self-consistency assumptions, this corresponds to a class of metadynamics-type CVbias sampling and SGLD modes ( The switching between free energy stochastic gradient descent and free energy stochastic gradient descent can be achieved by adjusting a single temperature hyperparameter, thus expanding the coverage of application scenarios.

[0191] 8. Relief of Optimization Bottlenecks Caused by SGLD-Style Random Perturbations: The Langevin integrator described in this invention, as an optimizer, under appropriate smoothness, noise intensity, and step size conditions, has an intrinsic Gaussian noise term that corresponds to the SGLD-style exploration capability, thereby improving the free energy target of the FP Transformer. The optimization path provides a fundamental support for overcoming saddle points. Although VFE has a unique global minimum (i.e., the target Boltzmann distribution), its optimization path has a saddle point structure, which leads to the mode dropping phenomenon reported in known methods such as the Boltzmann Generator. The Langevin-type optimizer described in this invention provides a fundamental means to alleviate this bottleneck by leveraging its random perturbation capability under the given conditions.

[0192] Glossary

[0193] The above are merely preferred embodiments of the present invention. It should be noted that the above preferred embodiments should not be considered as limitations on the present invention, and the scope of protection of the present invention should be determined by the scope defined in the claims. For those skilled in the art, several improvements and modifications can be made without departing from the spirit and scope of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A free energy optimizer, comprising an MD integrator or a Langevin integrator, characterized in that: Based on the replica group obtained by sampling the initial molecular conformation, the force term of the integrator is replaced or an external bias force is added to obtain the negative free energy gradient, and the integrator is made to advance the optimization of the replica group along the negative free energy gradient.

2. The free energy optimizer according to claim 1, characterized in that, The negative free energy gradient is obtained by linear superposition of the potential energy driving term and the entropy driving term; wherein, the potential energy driving term is constructed based on the potential energy gradient of the molecular conformation, and the entropy driving term is constructed based on the entropy contribution corresponding to the replica group distribution.

3. The free energy optimizer according to claim 2, characterized in that, Replace the force term with the action force. That is, the negative free energy gradient: in, For atomic coordinates, , Let be the potential energy function. reverse temperature , For the projection function of collective variables, It is obtained by nonparametric gradient estimation of the replica set, which is constructed based on the kernel trick.

4. The free energy optimizer according to claim 3, characterized in that, The force Includes physical force and entropy The physical force The entropy force .

5. The free energy optimizer according to claim 4, characterized in that, It also includes a physical force calculation module, which calculates the physical forces by calling machine learning potential energy models, custom force fields, or classical molecular force fields. .

6. The free energy optimizer according to claim 5, characterized in that, The machine learning potential energy model is an interatomic potential energy function constructed based on neural networks or graph neural networks, equivariant graph neural networks, kernel machine learning methods, etc.

7. The free energy optimizer according to claim 5, characterized in that, The custom force field or classical molecular force field is a molecular force field based on a parameterized potential energy function.

8. The free energy optimizer according to claim 4, characterized in that, It also includes an entropy force calculation module, which calculates the entropy force by calling the distributed gradient estimator with the replica set as the overall input. .

9. The free energy optimizer according to claim 8, characterized in that, The distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

10. The free energy optimizer according to claim 3, characterized in that, It also includes an adapter module to apply force. As the input force required for each step of the MD integrator or Langevin integrator, the copy is updated according to the update format.

11. The free energy optimizer according to claim 10, characterized in that, The MD integrator or Langevin integrator is Langevin, BAOAB, Leimkuhler-Matthews, Brownian (overdamped Langevin), or Nosé-Hoover Langevin.

12. The free energy optimizer according to claim 3, characterized in that, It also includes a scheduling module, which maintains N replicas (N greater than or equal to 1) at the system level, calculates the entropy force based on the overall replica group, and schedules the execution of physical force calculation, negative free energy gradient synthesis and integrator update for each replica.

13. The free energy optimizer according to claim 1, characterized in that, The integrator is an MD integrator or a Langevin integrator, configured in a non-intrusive, pluggable manner.

14. The free energy optimizer according to claim 13, characterized in that, Coupled entropy force terms at any or more of the integrator's force calculation, external force, bias force, callback, plug-in, or input interface. .

15. The free energy optimizer according to claim 14, characterized in that, The integrator's velocity / position update, constraint subroutine, neighbor table update, temperature controller, and parallel communication can be reused in one or more places.

16. The free energy optimizer according to claim 13, characterized in that, The force calculation of a single replica is extended to be called in parallel by N replicas, where N is greater than 1.

17. The free energy optimizer according to claim 4, characterized in that, The entropy force This corresponds to the metadynamics-type CV bias parameterization.

18. The free energy optimizer according to claim 17, characterized in that, The integrator O-step noise temperature satisfies And in Under the assumption of CV space boundary self-consistency and the standard construction conditions of kernel and pullback, the steady-state distribution of the integrator corresponds to the target behavior of metadynamics-type CV bias sampling under the parameterization.

19. The free energy optimizer according to claim 18, characterized in that, The entropy is calculated by calling the distribution gradient estimator with the replica set as the overall input. The distributed gradient estimator is constructed based on the nonparametric distribution estimation of the replica group.

20. The free energy optimizer according to claim 19, characterized in that, The distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

21. The free energy optimizer according to claim 19, characterized in that, The entropy force Extending along the CV dimension, the kernel trick replaces the CV grid, enabling bias forces to be implemented in high-dimensional CV subspaces or all-atom coordinate spaces by combining kernel approximations / sparse approximations.

22. The free energy optimizer according to claim 19, characterized in that, The entropy force By generalizing the bias expression, the kernelized density estimate is a nonparametric expression that does not force the parameterization of the sum of Gaussian hills and allows for more complex distribution geometry.

23. The free energy optimizer according to claim 19, characterized in that, The entropy force Following the generalization of replica origin, bias information originates from the instantaneous estimate of the current distribution of the replica group, and the increase in the number of replicas directly translates into an improvement in the accuracy of bias estimation.

24. The free energy optimizer according to claim 3, characterized in that, In the Gaussian noise amplitude of the Langevin integrator O-step, the physical temperature in MD is... Replace with an independently adjustable optimizer temperature hyperparameter. .

25. A method for implementation in molecular dynamics simulations, characterized in that: Based on the replica group obtained by sampling the initial molecular conformation, the force terms of the integrator are replaced or an external bias force is added to obtain the negative free energy gradient, and the integrator is then optimized along the negative free energy gradient.

26. The implementation method according to claim 25, characterized in that, The negative free energy gradient is obtained by linear superposition of the potential energy driving term and the entropy driving term; wherein, the potential energy driving term is constructed based on the potential energy gradient of the molecular conformation, and the entropy driving term is constructed based on the entropy contribution corresponding to the replica group distribution.

27. The implementation method according to claim 25, characterized in that, Replace the force term of the integrator with the applied force. That is, the negative free energy gradient: in, For atomic coordinates, , Let be the potential energy function. reverse temperature , For the projection function of collective variables, It is obtained by nonparametric gradient estimation of the replica set, which is constructed based on the kernel trick.

28. The implementation method according to claim 27, characterized in that, The force Includes physical force and entropy The physical force is The entropy force is .

29. The implementation method according to claim 28, characterized in that, The physical force is calculated by calling a machine learning potential energy model, a custom force field, or a classical molecular force field. .

30. The implementation method according to claim 29, characterized in that, The machine learning potential energy model is an interatomic potential energy function constructed based on neural networks or graph neural networks, equivariant graph neural networks, kernel machine learning methods, etc.

31. The implementation method according to claim 29, characterized in that, The custom force field or classical molecular force field is a molecular force field based on a parameterized potential energy function.

32. The implementation method according to claim 28, characterized in that, The entropy is calculated by calling the distribution gradient estimator with the replica set as the overall input. .

33. The implementation method according to claim 32, characterized in that, The distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

34. The implementation method according to claim 27, characterized in that, Force As the input force required for each step of the MD integrator or Langevin integrator, the copy is updated according to the update format.

35. The implementation method according to claim 34, characterized in that, The MD integrator or Langevin integrator is Langevin, BAOAB, Leimkuhler-Matthews, Brownian (overdamped Langevin), or Nosé-Hoover Langevin.

36. The implementation method according to claim 27, characterized in that, At the system level, maintain N replicas (N greater than or equal to 1), calculate the entropy force based on the replica group as a whole, and schedule and execute physical force calculation, negative free energy gradient synthesis and integrator update for each replica separately.

37. The implementation method according to claim 25, characterized in that, The integrator is an MD integrator or a Langevin integrator, configured in a non-intrusive, pluggable manner.

38. The implementation method according to claim 37, characterized in that, Coupled entropy force terms at any or more of the integrator's force calculation, external force, bias force, callback, plug-in, or input interface. .

39. The implementation method according to claim 38, characterized in that, The integrator's velocity / position update, constraint subroutine, neighbor table update, temperature controller, and parallel communication can be reused in one or more places.

40. The implementation method according to claim 37, characterized in that, The force calculation of a single copy is extended to N copies being called in parallel, where N is greater than 1.

41. The implementation method according to claim 25, characterized in that, The entropy force This corresponds to the metadynamics-type CV bias parameterization.

42. The implementation method according to claim 41, characterized in that, The integrator O-step noise temperature satisfies And in Under the assumption of CV space boundary self-consistency and the standard construction conditions of kernel and pullback, the steady-state distribution of the integrator corresponds to the target behavior of metadynamics-type CV bias sampling under the parameterization.

43. The implementation method according to claim 42, characterized in that, The entropy force The traditional Metadynamics method for accumulating single-trajectory time accumulation of Gaussian hills on the CV grid is replaced by nonparametric kernel density estimation and instantaneous estimation of the current distribution of the replica group.

44. The implementation method according to claim 43, characterized in that, The entropy is calculated by calling the distribution gradient estimator with the replica set as the overall input. The distributed gradient estimator is constructed based on the nonparametric distribution estimation of the replica group.

45. The implementation method according to claim 44, characterized in that, The distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

46. ​​The implementation method according to claim 43, characterized in that, The entropy force By generalizing the bias expression, the kernelized density estimate is a nonparametric expression that does not force the parameterization of the sum of Gaussian hills and allows for more complex distribution geometry.

47. The implementation method according to claim 43, characterized in that, The entropy force Following the generalization of replica origin, bias information originates from the instantaneous estimate of the current distribution of the replica group, and the increase in the number of replicas directly translates into an improvement in the accuracy of bias estimation.

48. The implementation method according to claim 27, characterized in that, In the Gaussian noise amplitude of the Langevin integrator O-step, the physical temperature in MD is... Replace with an independently adjustable optimizer temperature hyperparameter. .

49. A method for implementing a free energy-based molecular conformation optimizer in an MD integrator or a Langevin integrator, characterized in that, Includes the following steps: Step 1: Given the current position of the integrator ,momentum With force Progress to the predetermined value ,in, This is the integration time step; Step 2: Apply force to the integrator. ,in, For atomic coordinates, , Let be the potential energy function. reverse temperature , For the projection function of collective variables, It is obtained by nonparametric gradient estimation of the replica set, which is constructed based on the kernel trick; Step 3: The integrator performs optimization as a free energy optimizer.

50. The implementation method according to claim 49, characterized in that, The negative free energy gradient corresponding to the force is obtained by the linear superposition of the potential energy driving term and the entropy driving term; wherein, the potential energy driving term is constructed based on the potential energy gradient of the molecular conformation, and the entropy driving term is constructed based on the entropy contribution corresponding to the replica group distribution.

51. The implementation method according to claim 49, characterized in that, The force Includes physical force and entropy The physical force is The entropy force is .

52. The implementation method according to claim 51, characterized in that, The physical force is calculated by calling a machine learning potential energy model, a custom force field, or a classical molecular force field. .

53. The implementation method according to claim 52, characterized in that, The machine learning potential energy model is an interatomic potential energy function constructed based on neural networks or graph neural networks, equivariant graph neural networks, kernel machine learning methods, etc.

54. The implementation method according to claim 52, characterized in that, The custom force field or classical molecular force field is a molecular force field based on a parameterized potential energy function.

55. The implementation method according to claim 49, characterized in that, The entropy is calculated by calling the distribution gradient estimator with the replica set as the overall input. .

56. The implementation method according to claim 55, characterized in that, The distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

57. The implementation method according to claim 49, characterized in that, Force As the input force required for each step of the MD integrator or Langevin integrator, the copy is updated according to the update format.

58. The implementation method according to claim 57, characterized in that, The MD integrator or Langevin integrator is Langevin, BAOAB, Leimkuhler-Matthews, Brownian (overdamped Langevin), or Nosé-Hoover Langevin.

59. The implementation method according to claim 49, characterized in that, At the system level, maintain N replicas (N greater than or equal to 1), calculate the entropy force based on the replica group as a whole, and schedule and execute physical force calculation, negative free energy gradient synthesis and integrator update for each replica separately.

60. The implementation method according to claim 49, characterized in that, The integrator is an MD integrator or a Langevin integrator, configured in a non-intrusive, pluggable manner.

61. The implementation method according to claim 60, characterized in that, Coupled entropy force terms at any or more of the integrator's force calculation, external force, bias force, callback, plug-in, or input interface. .

62. The implementation method according to claim 61, characterized in that, The integrator's velocity / position update, constraint subroutine, neighbor table update, temperature controller, and parallel communication can be reused in one or more places.

63. The implementation method according to claim 60, characterized in that, The force calculation of a single copy is extended to N copies being called in parallel, where N is greater than 1.

64. The implementation method according to claim 49, characterized in that, The entropy force This corresponds to the metadynamics-type CV bias parameterization.

65. The implementation method according to claim 64, characterized in that, The integrator O-step noise temperature satisfies And in Under the assumption of CV space boundary self-consistency and the standard construction conditions of kernel and pullback, the steady-state distribution of the integrator corresponds to the target behavior of metadynamics-type CV bias sampling under the parameterization.

66. The implementation method according to claim 65, characterized in that, The entropy is calculated by calling the distribution gradient estimator with the replica set as the overall input. The distributed gradient estimator is constructed based on the nonparametric distribution estimation of the replica group.

67. The implementation method according to claim 66, characterized in that, The distributed gradient estimator employs the kernel trick to construct a kernelized density estimate.

68. The implementation method according to claim 66, characterized in that, The entropy force Extending along the CV dimension, the kernel trick replaces the CV grid, enabling bias forces to be implemented in high-dimensional CV subspaces or all-atom coordinate spaces by combining kernel approximations / sparse approximations.

69. The implementation method according to claim 66, characterized in that, The entropy force By generalizing the bias expression, the kernelized density estimate is a nonparametric expression that does not force the parameterization of the sum of Gaussian hills and allows for more complex distribution geometry.

70. The implementation method according to claim 66, characterized in that, The entropy force Following the generalization of replica origin, bias information originates from the instantaneous estimate of the current distribution of the replica group, and the increase in the number of replicas directly translates into an improvement in the accuracy of bias estimation.

71. The implementation method according to claim 49, characterized in that, In the Gaussian noise amplitude of the Langevin integrator O-step, the physical temperature in MD is... Replace with an independently adjustable optimizer temperature hyperparameter. .

72. A system for implementing the free energy optimizer, characterized in that: The system includes a processor, a memory, and a computer program stored in the memory, which, when executed by the processor, causes the system to perform the methods of claims 25-71.

73. A computer-readable storage medium having a computer program stored thereon, the computer program being executed by a processor to implement the method of claims 25-71.

Images