Biomolecular aqueous solutions represent a crucial class of condensed matter and serve as the fundamental basis of living organisms. Investigating their structural and chemical properties holds immense scientific and practical significance. Molecular dynamics （MD） simulations are a powerful tool for studying biomolecular systems， where accuracy critically depends on the precision of the molecular force field. Traditional force fields employ fixed atomic charges， neglecting polarization effects and charge transfer. Over recent decades， significant efforts have been devoted to developing polarizable force fields. The atom-bond electronegativity equalization method （ABEEM） polarizable force field effectively captures molecular polarization and charge transfer. This review outlines the ABEEM methodology， with emphasis on the ABEEM-7P water model and the application of the ABEEM polarizable force field is provided. The work has carried out using both the ABEEM-7P water model and the ABEEM polarizable molecular force field for molecular dynamics simulations of biomolecule aqueous systems，and the results have been given by representative case studies.
Contents
1 Introduction
2 Theoretical methods overview
2.1 Born-oppenheimer approximation for solving the time-independent Schrödinger equation
2.2 Molecular mechanics （MM）
2.3 Common force fields
2.4 ABEEM polarizable force field
2.5 Molecular dynamic simulation
3 Molecular dynamics simulation in biomolecular aqueous solution system
3.1 Application of ABEEM-7P water model and polarization force field
3.2 Molecular dynamics simulation example of ABEEM PF biomolecule aqueous solution system
Water molecules, as the most common solvent, interact with numerous other compounds and occupy a central position in both the study of molecular aqueous solutions and the development of molecular force fields. Whether in quantum mechanical simulations or classical molecular dynamics simulations, water systems remain one of the primary research focuses. Water is the major intracellular solvent (accounting for 70%–80% of cellular volume), and biomolecules such as proteins and nucleic acids can maintain their native conformations and functional activities—such as enzymatic catalysis and DNA replication—only in aqueous solutions. Non-aqueous environments (e.g., vacuum) can lead to protein misfolding or loss of function. Therefore, water serves as an essential solvent for simulating the dynamic behavior of biomolecules under physiological conditions, and the selection of an appropriate water model can significantly enhance the accuracy of such simulations.
Molecular dynamics simulations are indispensable tools for studying the properties of biological macromolecules such as proteins and nucleic acids［1-3］. However, their effectiveness depends on the accuracy of the molecular force fields that describe the interactions among all atoms. Therefore, it is essential to develop new force fields for investigating the structure and dynamics of biomolecular aqueous systems. Based on the fundamental approximations of molecular mechanics, as long as a potential energy function can adequately represent the various interactions within a system, the motions of its constituent atoms and molecules can be described using classical mechanics derived from that potential energy function. Such a potential energy function, along with its associated parameters, is referred to as a molecular force field. Early conventional molecular force fields (FFs), in order to simulate biomolecular systems with moderate computational resources, typically employed fixed atomic charges as parameters. This approach, however, has inherent limitations［4-9］, because it neglects polarization and charge transfer between charged atoms—two phenomena that are crucial for the recognition and stability of biological systems［5, 9］. In contrast, polarizable force fields allow for dynamic changes in atomic charges, thereby better capturing polarization effects. Over the past several decades, numerous polarizable force fields have been developed［10］, including AMBER［11-13］, AMOEBA［14-15］, CHARMM Drude［6］, CHARMM floating charges［16-20］, SIBFA［21］, GEM［21］, and ABEEM［22-31］. Consequently, the development of polarizable force fields for molecular dynamics simulations of biomolecular systems is highly necessary.
This review summarizes the relevant work of the ABEEM Polarizable Force Field (ABEEM PFF) in molecular dynamics simulations of biomolecular aqueous solutions. It introduces the ABEEM-7P water model and the applications of ABEEM PFF in biomolecular systems. Through several molecular dynamics simulation examples, it demonstrates that ABEEM PFF can accurately predict dynamic properties such as structure, charge distribution, and energy of biomolecules in aqueous solutions. Finally, it provides a prospect for the future development of ABEEM PFF.
2 Theoretical Methods
2.1 Solution of the Time-Independent Schrödinger Equation Using the Born–Oppenheimer Approximation
Molecules are often systems consisting of several atomic nuclei and multiple electrons, whose wave functions satisfy the time-independent Schrödinger equation. This equation serves as the fundamental equation governing the electronic structure of molecules. Based on the concept that electrons, which move rapidly, can immediately adjust themselves in response to even a slight change in molecular geometry, Born and Oppenheimer simplified the time-independent Schrödinger equation for molecular systems by separating nuclear motion from electronic motion. This approach is known as the fixed-nuclei approximation or the Born–Oppenheimer approximation. Under this approximation, the wave function of a molecular system can be expressed as the product of the electronic motion wave function and the nuclear motion wave function . Correspondingly, the nonrelativistic time-independent Schrödinger equation for the molecular system can be decomposed into two equations, which are the most important equations in the theory of electronic structure in quantum chemistry. These are the equation for electronic motion (Equation 1)and the equation for nuclear motion (Equation 2).
in the equation: represents the coordinates of the Ath nucleus; denotes the coordinates of all electrons; stands for the coordinates of all nuclei; represents the electronic wave function at a given nuclear coordinate ; denotes the nuclear wave function, and E represents the total energy of the system. Obviously, under a fixed set of nuclear coordinates , the total electronic energy of the system is determined by Equation (1) and also appears in Equation (2), which is referred to as the potential energy surface. The potential energy surface is indispensable when discussing molecular structures, chemical reactions, and related topics.
2.2 Molecular Mechanics (MM)
MM is a computational method based on classical mechanics. This method uses empirical or semi-empirical parameters to calculate the structure and energy of molecules, commonly referred to as the force field (FF) approach. Unlike quantum mechanical methods, it does not require solving the wave function of the system. Instead, it relies on the Born–Oppenheimer approximation to separate the motions of electrons and nuclei, considering only the motion of atoms and treating the system’s energy as a function of nuclear positions. The potential energy surface described by molecular force fields generally includes terms for bond stretching vibrations, bond angle bending vibrations, dihedral torsional potentials, non-planar dihedral potentials, as well as electrostatic and van der Waals interactions. These force field parameters are typically obtained through fitting to experimental data or quantum chemical calculations. By employing force field methods, one can compute various properties of large molecular systems, including energy, structure, and thermodynamic properties. This approach provides a powerful tool for studying molecular properties and finds wide-ranging applications, particularly in exploring the behavior of large molecular systems.
2.3 Common Force Fields
In the 1970s, computer simulations began to be widely applied in the study of organic and biomolecular systems. Molecular dynamics and Monte Carlo statistical mechanics were employed to simulate systems such as liquid water, aqueous solutions containing simple molecules and ions, organic liquids, and small proteins in vacuum. The primary focus of these studies was on the accurate description of intra- and intermolecular energies. Given the size of the systems and the available computational resources, classical force fields were typically used for simulations. These force fields form the foundation of molecular mechanics research for organic molecules, and differences among them mainly arise from variations in the force field models employed, including whether hydrogen bonding interactions or polarization effects are considered. Furthermore, various force fields also differ in their mathematical functional forms and parameter settings, contributing to their diversity. Based on how electrostatic interactions are handled, force fields can be broadly classified into two main types: fixed-charge force fields and polarizable force fields. A key characteristic of fixed-charge force fields is that the atomic charges within the molecular system remain constant throughout the course of a dynamical simulation. This means that even if the molecular structure undergoes changes, the charge distribution on the atoms remains unchanged. Due to its inability to accurately account for charge polarization and redistribution within molecules, this model has certain limitations when simulating molecular behavior. When dealing with interactions between highly charged ions, interactions between ions and π-electron systems, as well as ion conductivity and permeability［14］, the use of polarizable force fields provides a more accurate description of these phenomena［32］.
2.3.1 Fixed Charge Force Field
The concept of the fixed-charge force field (FFF) can be traced back to the late 1960s and early 1970s. Its core idea is to treat atomic charges within a molecule as fixed values, which are primarily derived from quantum mechanical (QM) calculations and experimental data. The AMBER (Assisted Model Building with Energy Refinement) force field［33-34］, as one of the earliest fixed-charge force field models, has been mainly used to simulate protein and nucleic acid systems. Subsequently, other force field models such as CHARMM (Chemistry at HARvard Macromolecular Mechanics)［35-36］, OPLS (Optimized Potentials for Liquid Simulations)［37-39］, and GROMOS (GROningen MOlecular Simulation)［40-41］ emerged in succession.
The AMBER force field was developed by Kollman et al.［33］. Its advantages lie in its ability to rapidly and accurately compute the properties of biomacromolecular systems, such as those containing proteins, nucleic acids, and polysaccharides, while requiring a relatively small number of force field parameters. In the potential energy function framework of the AMBER force field, the total potential energy can be decomposed into bonded and nonbonded interactions, as expressed in Equation (3).
the CHARMM force field, developed by the Karplus group［35］, can yield molecular structures and energetic properties that are nearly consistent with experimental results when applied to organic small molecules and macromolecular systems in solution. The force field parameters are derived from experimental data and *ab initio* calculations. The potential function of the CHARMM force field is shown in Equation (4), where each term represents the bond-length potential, bond-angle potential, Urey-Bradley interaction, dihedral-angle potential, out-of-plane interaction potential, dispersion interaction, van der Waals interaction, and Coulombic interaction, respectively.
to better simulate the thermodynamic properties of liquids, Jorgensen et al. ［37］developed the OPLS force field. This force field includes both all-atom and united-atom versions, primarily designed for simulating solvent-containing systems. The bond-related parameters are derived from the AMBER force field, while the non-bonded parameters are calibrated based on experimental data.
Currently, in addition to force fields such as AMBER, CHARMM, and OPLS, the GROMOS force field is also one of the most important fixed-charge molecular force fields used for biomolecular simulations. This force field is a united-atom molecular force field that has been carefully designed to enhance computational efficiency. It is suitable for simulating alkanes, proteins, and nucleic acids in the condensed phase, with the potential energy function expressed as Equation (5)..
computational simplicity is the most prominent advantage of fixed-charge force fields, making them suitable for large-scale molecular simulations and widely used in molecular dynamics simulations. For example, in drug design research, the most significant contribution of fixed-charge force fields lies in their ability to predict binding sites between drug molecules and proteins, thereby enabling the screening of potential drug candidates. In studies of protein folding and enzyme catalysis, fixed-charge force fields can provide important computational results at the molecular level. Furthermore, fixed-charge force fields are also valuable for investigating interfacial phenomena in materials chemistry, as well as phase transition behaviors and thermodynamic properties of nanomaterials.
In the process of continuous improvement and development, these force field models have gradually demonstrated their unique stylistic characteristics and applications. Although fixed-charge force fields have shown good applicability in many areas, they still have certain shortcomings and limitations. This is because the atomic charges in fixed-charge force fields remain constant and do not account for charge polarization effects. Consequently, when describing intermolecular interactions under certain conditions, such models may lack accuracy. For example, the neglect of polarization effects results in no difference between the central molecular charge distribution calculated in the gas phase and that obtained from solvent calculations, leading to incorrect conformational predictions. When molecules are subjected to strong electric fields or highly polar environments, the simulation results may deviate significantly from the actual situation. To overcome these limitations and more accurately describe intermolecular interactions, more complex polarizable force fields have been developed.
2.3.2 Polarization Force Field
Polarizable force fields (PFFs) are capable of explicitly accounting for polarization effects that arise in real-world environments during simulations. Polarizable force fields primarily represent polarization effects through two models: The first model incorporates induced dipoles (multipole moments)［15］or Drude oscillators (also known as the charge-on-spring model or shell model)［6, 42-43］to reflect the redistribution of charge within each atom; the other is a floating charge model based on charge transfer between atoms［17, 25］. In the induced dipole model, each polarizable site generates an induced dipole moment μindunder the influence of an external electric field E, where μind = αE, and αrepresents the polarizability［15］. These induced dipoles also contribute to the total electric field; therefore, polarization is not additive and is typically determined through self-consistent iterative methods. In the Drude oscillator model［6, 42-43］, a harmonic spring connects a Drude particle carrying a partial charge to the atom. The displacement of the Drude particle gives rise to an induced dipole moment. Although extended Lagrangian methods have been employed in MD simulations to approximate the exact solution, the position of the Drude particle still requires self-consistent iteration to ensure it remains in its ground state. The floating charge model is based on the principle of electronegativity equilibration, whereby atomic charges are determined by balancing the electronegativity or chemical potential at each site. Both the induced dipole and Drude oscillator models allow for the generation of a dipole at each site. Given that the displacement of the Drude oscillator is usually very small, it is often approximated as a point dipole. The floating charge model permits charge transfer between atoms. However, simple floating charge models do not account for out-of-plane polarization or the polarization of single-atom ions. These effects can be represented in the floating charge model by introducing additional virtual sites around each atom［28］.
Currently, the polarizable force fields applied to biological systems include AMBER［11-12］,AMOEBA［14-15］,CHARMM Drude［6］,CHARMM-FQ［16-17］,SIBFA (Sum of interaction between fragments ab initio) ［21］,GEM (Gaussian electrostatic model) ［21］, and ABEEM PFF ［22-23］, among others.
AMBER ff02pol［12］ is one of the earliest polarizable force fields used for proteins and nucleic acids, employing point charges and a simple undamped induced dipole model. Subsequently, ff12pol introduced a model incorporating the Thole damping function, thereby improving the accuracy of intermolecular interaction energies［11］.
The AMOEBA polarizable force field uses atomic induced dipoles to simulate polarization and employs multipole moments up to the quadrupole moment to represent permanent electrostatic interactions. The AMOEBA force field has been applied to simulate water, ions, organic molecules, and protein systems［44-49］. Mu et al.［50］used the AMOEBA force field to capture the σ-hole effect. Zhang et al.［51］developed an AMOEBA force field suitable for DNA and RNA and extensively validated it through molecular dynamics simulations lasting up to 35 μs. By simulating the crystal structures of DNA double strands, RNA double strands, and hairpin loops, their results were consistent with those determined by NMR, with a root-mean-square deviation of less than 2.0 Å (1 Å = 0.1 nm).
The CHARMM Drude force field employs the Drude oscillator model to describe polarization and off-center charges, thereby capturing the anisotropy of charge distribution［6］. The Drude force field covers most biomolecules, including small organic molecules, proteins/polypeptides, DNA［52］, and lipids［53］. Lin et al.［54］introduced displaced charges into aliphatic and aromatic halogenated molecules, enabling halogen atoms to exhibit anisotropic behavior. Furthermore, they improved the Drude force field by incorporating van der Waals parameters into the Drude particles.
SIBFA［21］is an ab initio polarizable force field composed of the sum of electrostatic multipole, short-range repulsion, polarization, charge transfer, and dispersion contributions. The values of each term are determined by reproducing the corresponding quantum mechanical terms. This force field was initially developed to handle Zn2+ proteins［55］and was subsequently extended to halogen compounds［56］and nucleic acids［57］.
2.4 ABEEM Polarizable Molecular Force Field
Our research group, based on DFT and the Electronegativity equalization method (EEM), proposed and developed the Atom–Bond Electronegativity Equalization Method model［58-59］, subsequently introducing the ABEEM Polarizable Force Field (ABEEM PFF). The distinctive feature of ABEEM PFF lies in its partitioning of the molecular single-electron density into atomic, σ-bond, π-bond, and lone-pair electron sites. This multi-site model captures the anisotropy of charge distribution around atoms. ABEEM PFF determines the charge distribution of the target molecule by solving geometry-dependent linear equations, thereby avoiding self-consistent iterative procedures and enabling rapid computation of charge distributions. Consequently, in ABEEM PFF, atomic charges fluctuate with changes in the surrounding environment during simulations, effectively capturing polarization effects and charge-transfer phenomena.
2.4.1 Charge Distribution of Molecules and Molecular Systems
The charge distribution of molecules and molecular systems—that is, the charges distributed on atoms or regions—is a highly significant physical quantity. It is indispensable for understanding both intra- and intermolecular interactions, constructing molecular force fields, and performing molecular dynamics simulations. Atomic charges, as crucial parameters for characterizing electron density distributions, play an important role in theoretical chemistry research. Different computational methods possess distinct physical foundations and application scopes; for example, wavefunction-based methods can reflect quantum mechanical characteristics, whereas molecular potential energy simulation methods are more suitable for building molecular force fields. Researchers continuously optimize these molecular force fields to more accurately describe the charge distribution properties within molecules and conduct molecular dynamics simulations, thereby providing a reliable theoretical basis for predicting molecular properties and designing new molecules.
The methods for calculating atomic charges are mainly divided into two major categories—sharp boundary and fuzzy boundary—depending on the way in which molecular regions are delineated, as illustrated in Figure 1. The sharp-boundary model employs quantum chemical topology (QCT) based on real-space scalar functions to partition the molecule entirely into atomic regions. By integrating the electron density over each atomic region, one can obtain the number of electrons associated with that atom; subtracting this value from the nuclear charge yields the atomic charge within the molecule. This category includes QTAIM charges［60-61］, PAEM charges［62-63］, and KSpot charges［64］, among others. In contrast, the fuzzy-boundary model primarily accounts for the overlap of electron densities between atoms and can be further classified as follows: 1. Methods based on wavefunctions, such as Mulliken charges［65-70］, NPA charges［71-73］, Hirshfeld charges［74-80］, MMPA charges［81］, Löwdin charges［82-85］, and AOIM charges［86］. 2. Atomic charges computed based on the principle of electronegativity equalization, including EEM charges［87-88］, QEq charges［89］, and ABEEM charges［59, 90］. 3. Atomic charges derived from fitting electrostatic potential, such as MK charges［91-93］, CHELPG charges［94］, ESP charges［95］, and RESP charges［96］. 4. Experimental atomic charges, such as APT charges obtained from atomic polarizability tensors［97-102］, IR charges inferred from infrared spectroscopy data［100, 103］, and charges derived from X-ray diffraction studies. Additionally, the dipole moments of small molecules can also be utilized to calculate atomic charges［104］. Other methods yielding atomic charges include MMFF94 charges, among others.
ABEEM (Atom-bond electronegativity equalization method) ［23, 105］ is an atom-bond electronegativity equalization method used to calculate the charge distribution on atoms, bonds, and lone pairs within a molecule. Based on the principle of electronegativity equalization, this method determines the charges of individual atoms, bonds, and lone pairs by solving a set of linear equations. ABEEM specifically accounts for hydrogen-bonding interaction regions, where the interactions between lone pairs and hydrogen atoms are described using hydrogen-bond fitting functions. The charge interactions in these regions are crucial for understanding hydrogen bonding among water molecules. The model includes several adjustable parameters, such as the positions of bond charges and lone pairs, as well as the interaction parameters for lone pairs that form hydrogen bonds with hydrogen atoms. These parameters are determined by fitting experimental data or high-accuracy quantum chemical calculation results.
Figure 2 is the ABEEM site map of the thymine molecule, where yellow elliptical spheres represent lone-pair (lp) electron sites located at the covalent radii of double-bonded atoms, with their angles determined by the hybridization orbital types. Light-blue elliptical spheres denote π-bond sites perpendicular to the σ-bond plane, situated at the covalent radii of double-bonded atoms, above and below them. Green elliptical spheres indicate σ-bond sites, positioned at the ratio of covalent radii between two bonded atoms. Red, light-gray, and white spheres respectively represent oxygen, carbon, and hydrogen atoms. From Figure 2, we can see that thymine (T) has a total of 50 charge sites, including 15 atomic sites on the nuclei, 15 σ-bond sites, 16 π-bond sites, and 4 lp-electron sites. The ABEEM model can clearly illustrate the positions of various sites within the molecule; as multiple sites change with environmental conditions, charges fluctuate among all sites, thereby manifesting polarization effects and charge-transfer effects.
In the ABEEM model, based on the effective electronegativity in DFT, =, then the effective electronegativity of each site can be expressed as Equation (6)..
according to the principle of electronegativity equilibrium, the effective electronegativity at each site in a molecule is equal to the effective electronegativity of the entire molecule, as expressed in Equation (7)..
meanwhile, the sum of charges at all sites equals the total charge of the molecule or system, as shown in Equation (8).
in the equation: Rij represents the distance between sites i and j; qi and qj are the point charges at sites i and j, respectively. When the absolute value of the charge at site i exceeds 1.0e—for example, for transition metal atoms in metalloproteins and phosphorus atoms in nucleic acids—it is necessary to set a valence charge qi* so that the absolute value of the relative point charge at site i, |qi′|=|qi-qi *|, is less than 1.0e［23］. Therefore, once the valence electronegativity , valence hardness (), and electrostatic coefficient (kij) for each site in the system have been determined, the charges of all sites under different chemical environments can be obtained by solving a system of linear equations.
2.4.3 ABEEM Polarizable Molecular Force Field Potential Energy Function
When the ABEEM method is combined with a molecular force field, the resulting ABEEM polarizable force field (ABEEM PFF) has its potential energy function expressed as Equation (9)..
in the equation: and represent, respectively, the potential energy of bond-length stretching vibrations and the potential energy of bond-angle bending vibrations. The bond-length stretching vibration potential energy of the water molecule is expressed using the Morse potential function, whereas other bond-length stretching vibration potentials are described by a harmonic oscillator model. and represent, respectively, the torsional potential energy associated with dihedral angles and the nonplanar bending potential term. EvdWand Eelecdenote, respectively, the van der Waals and electrostatic interaction potential energies, whose explicit expressions are given in Equations (10)through (13).
in the equation: De represents the potential well at the minimum potential energy; is a parameter related to the force constant kb. Its specific expression is: , where kr is the force constant for each bond. kr and kθ represent the force constants for bond-length stretching vibrations and bond-angle bending vibrations, respectively; r and θ denote the actual values of bond length and bond angle, respectively; req and θeq indicate the equilibrium bond length and equilibrium bond angle under the equilibrium geometry; ν1, ν2, and ν3 represent the expansion force constants of the dihedral torsional potential term, where ϕ denotes the actual dihedral angle value; ν denotes the force constant of the nonplanar bending potential term; is the actual value of the nonplanar dihedral angle.
Non-bonded interaction terms, including the Lennard-Jones potential term (Equation (14)) and the Coulombic term (Equation (15)).
Equation (14) In this equation, the potential energy parameters of atoms a and b are equal to the geometric mean of the two atomic potential energy parameters, which are represented by Equation (16) and (17), respectively.
in the equation: a and b are in a 1-2 or 1-3 relationship; the coefficient fab is 0.0; for interactions of the 1-4 type, fab is 0.5; whereas for all other intramolecular and intermolecular interactions, fab is 1.0.
In the non-bonded potential terms, in addition to the Lennard-Jones potential term, there is also the Coulomb potential term. In Equation (15), qi and qj represent the point charges within the molecule; Rij is the distance between the point charges qi and qj; and kij is an overall correction factor of 0.57. When i and j are capable of forming a hydrogen bond, we employ the hydrogen-bond fitting function kij = kHB, where kHB can be expressed as Equation (18):
in the equation: k′, ri,j, C, and D are parameters, while Rlp,H represents the distance between a lone pair of electrons and a hydrogen atom. During explicit water molecular dynamics simulations, hydrogen bonds exist among water molecules. The function kHB is used to describe the functional relationship of these hydrogen bonds, more effectively characterizing the electrostatic interactions among water molecules.
2.5 Molecular Dynamics Simulation
Molecular dynamics simulation (MD) is currently the most widely used computational method for studying large and complex systems. Since 1970, with the rapid advancement of molecular mechanics, numerous force fields have been systematically developed for biochemical molecular systems, polymers, metals, and nonmetallic materials, greatly enhancing both the capability and accuracy of calculating the structural, thermodynamic, and spectroscopic properties of complex systems. Molecular dynamics simulation is a computational approach that has been developed specifically to apply these force fields together with Newton’s laws of motion. Its advantages include high precision, the simultaneous acquisition of both dynamical and thermodynamic statistical data for a system, and broad applicability across various systems and property types.
Consider a system of motion containing N molecules. The energy of the system is the sum of the kinetic energies of the molecules in the system and the total potential energy of the system. Its total potential energy is a function of the positions of all atoms within the molecules, expressed as . Typically, the potential energy can be divided into two parts: the van der Waals (vdW) interactions between atoms both within and between molecules, and the internal potential energy (int). This is expressed in Equation (19).
the van der Waals interaction can generally be approximated as the summation of van der Waals interactions between each pair of atoms, while the intramolecular potential energy is the sum of potential energies corresponding to various internal coordinates (such as bond stretching, bond angle bending, etc.).
The force acting on any atom i in the system is given by the gradient of the potential energy, as shown in Equation (20). The acceleration of atom i is given in Equation (21). By integrating with respect to time, one can predict the velocity and position of atom i after a time interval t, as expressed in Equation (22).
Equation (22) where: and denote the position and velocity of the particle, respectively; the superscript "0" indicates the initial values of each physical quantity.
The fundamental principle of molecular dynamics simulations involves first calculating the potential energy of the system based on the positions of all molecules within it. Next, the forces and accelerations acting on each atom are computed. By setting t in the equations equal to , one can determine the positions and velocities of all molecules after a time increment of . The above steps are then repeated: the potential energy is recalculated from the new atomic positions, the forces and accelerations are determined once more, and the subsequent positions and velocities of the molecules are predicted after another time increment of . Through this iterative process, data on the positions, velocities, and accelerations of the molecules at various time points can be obtained. Typically, the molecular positions at each time step are referred to as the trajectory.
3 Molecular Dynamics Simulation of Biomolecular Aqueous Solutions
Molecular dynamics simulation is an important tool for analyzing the behavior of biomolecules in aqueous solutions. The ABEEM-7P water model and polarizable force field have been successfully applied to biomolecular systems such as proteins and nucleic acids. This chapter will provide a detailed introduction to their applications, as well as the work conducted using the ABEEM PFF in molecular dynamics simulations of biomolecular aqueous solutions.
3.1 The ABEEM-7P Water Model and the Application of the ABEEM Polarizable Molecular Force Field
3.1.1 ABEEM-7P water molecule model
Common simple water molecule models［106-114］ From the early ST2 and ST4 models, they gradually evolved into the Simple Point Charge (SPC) force field and the Transferable Intermolecular Potential (TIPnP) function. These models all employ a rigid water molecule structure, characterized by fixed bond lengths and bond angles. The potential energy surface between water molecules is described using van der Waals interactions, as well as three-, four-, or five-point electrostatic interactions. As shown in Figure 3, these include commonly used three-, four-, and five-point simple water molecule force field models.
These models all utilize a rigid water molecule structure and fixed charges as parameters, making it impossible to describe changes in electrostatic interactions arising from alterations in the electric field surrounding atoms (or molecules). In contrast, we have developed a new computational model for water molecular systems: the Transferable Intermolecular Potential Seven Points Approach Including Fluctuating Charge and Flexible Model (ABEEM-7P)［115］. The ABEEM-7P model employs additional sites to calculate electrostatic interactions between water molecules; Figure 4illustrates the monomeric water structure of the ABEEM-7P model. A single water molecule is divided into seven sites, with the three atomic nuclei strictly positioned at the geometric coordinates of the hydrogen and oxygen atoms. The oxygen atom adopts sp3hybridization to form a tetrahedral electron distribution. The two σ-bond charge density maxima are located along the O—H bonds, their precise positions determined by the Pauling covalent radii ratio (0.66 Å for oxygen and 0.28 Å for hydrogen). The bond angle between the two σ-bonds is 104.5°, a result of the spatial repulsion exerted by the oxygen atom’s lone pairs on the bonding orbitals. Furthermore, two lone pair (lp) sites are symmetrically distributed around the oxygen nucleus at a solid angle of 109.47°, at a distance of 0.74 Å from the oxygen nucleus center (approximately 40% of the oxygen atom’s van der Waals radius). This spatial configuration strictly adheres to the principle of minimizing repulsive energy among the lpsites.
In the ABEEM-7P model, the Morse potential function is used to describe the bond-length stretching vibrations. Due to hydrogen-bond formation in liquid water systems, the H—O bond exists in different forms, such as OH- and H3O+. The Morse potential function can describe bond stretching vibrations over a wider range, from the equilibrium bond length to bond dissociation, which is extremely important for a detailed analysis of the vibrational modes of water molecular systems.
The equilibrium bond length req of the O—H bond and the equilibrium bond angle θeq of H—O—H in the ABEEM-7P model are set to the experimental values of 0.9572 Å and 104.52°, respectively. The force field parameters of the ABEEM-7P model are listed in Table 1, where the dissociation energy D of the oxygen–hydrogen bond is 529.6 kcal/mol, and the angular force constant kθ is 34.05 kcal/mol·(°)2. α is a parameter related to the force constant of the bond length (, with kb being the bond force constant). In addition, the force field parameters for the TIP4P-FQ［116］, TIP5P［107］, and POL5［117］ models are also provided.
The monomer water charge distribution and dipole moment calculated using the ABEEM-7P water molecular model. Compared with other water molecular force fields, the ABEEM-7P water molecular model fully considers conformational changes and provides the charge distribution of all atoms, bonds, and lone-pair electrons. Meanwhile, the dipole moment obtained from the ABEEM-7P model, 1.8551, is in very good agreement with the experimental value of 1.855［115］. This excellent agreement is closely related to the reasonable partial charges at each site; therefore, an accurate description of electrostatic interactions can precisely characterize the dipole moment of water molecules.
3.1.2 Water Clusters
Water molecular clusters are one of the popular research topics in numerous theoretical and experimental studies. An important objective of studying water molecular clusters is to provide a quantitative description of aqueous solution behavior at the molecular level. Advanced research on water molecular clusters is expected to offer powerful tools for deeply revealing fundamental issues in fields such as chemistry, life sciences, and information science. The significant theoretical implications and broad application prospects of studying the structure and function of water molecular clusters have made this area one of the active frontiers of scientific research.
From a structural perspective, water molecular clusters can be classified into two types: two-dimensional structures (such as cyclic structures) and three-dimensional structures (such as cage-like, prismatic, or cubic structures). It should be further noted that there are four types of monomers in water molecular clusters: one proton donor and one acceptor ("da" type), one donor and two acceptors ("daa" type), two donors and one acceptor ("dda" type), and two donors and two acceptors ("ddaa" type).
Figure 5 presents the various low-energy structures of (H2O) n (n=2–10). Among the five low-energy structures of heptameric water, structure 7A, which resembles a cubic arrangement with one corner missing, is the lowest-energy configuration, followed by structure 7B. These findings are consistent with those reported in ab initio calculations［119-120］, and experimental spectroscopic data［121］also confirm the existence of both 7A and 7B. It can be observed that 7A and 7B can be derived from octameric water with S4 symmetry by removing either a “dda” or a “daa” type monomeric water molecule; alternatively, they can be obtained by inserting a “da” type water molecule onto or beneath the triangular base of the hexameric prism structure. Notably, all structures shown in Figure 5(except for structure D2d) possess energetically equivalent isomers that can be generated through reflection［120］. These isomers are non-superimposable on the original structures, indicating chiral symmetry. Structure 7A has 10 hydrogen bonds, with a binding energy of −230.25 kJ/mol. Structure 7B also features 10 hydrogen bonds, but its binding energy is −227.65 kJ/mol, slightly higher than that of 7A by 2.60 kJ/mol. This difference is very close to the approximately 1.90 kJ/mol discrepancy reported by Kim et al.［119］, who employed the MP2/6-311++G**//HF/6-311++G** method with a 50% basis set superposition error correction (BSSE), as well as the MP2/6-311++G**//B3LYP/6-311++G** approach incorporating 50% BSSE. Furthermore, it is similar to the results obtained by Sathyamurthy et al.［122］using the HF/6-31G(d,p)//HF/6-31G(d,p) level of theory. The double-ring structures R4 and R5 can both be regarded as arising from the hexameric book-like arrangement by the insertion of an additional water molecule; they exhibit comparable binding energies and the same number of hydrogen bonds.
This model predicted the dipole moments and structural parameters of various low-energy structures in heptameric water. For the dipole moments, the computational results obtained from this model show excellent agreement with those reported by Kim et al.［119］ at the MP2/TZ2P++//B3LYP/6-311++G** level, as well as with our own calculations performed at the HF/6-31G(d,p)//HF/6-31G(d,p) level. In addition, the ABEEM-7P model has been applied to the study of larger water clusters (H2O)n (n=20, 24, 28, 32, 35). The results indicate that for these larger water clusters, the ABEEM polarizable molecular force field model can also yield results consistent with ab initio calculations and other semi-empirical methods in terms of cluster geometry and binding energies. These findings will contribute to a deeper understanding of medium-sized, large, and even larger water clusters, provide important insights into the properties of aqueous solutions, and offer reliable evidence for further validating the rationality, accuracy of parameters, and transferability of the ABEEM polarizable molecular force field model.
3.1.3 Peptides and Proteins
Research on the conformational distribution of biomolecules is crucial for understanding their functions. The accurate search for stable conformations remains a challenging task. In 2015, Xu et al.［123］employed the ABEEM PFF to explore the stable conformations of Val-dipeptide and Val-tripeptide, identifying six and thirty-four stable conformations, respectively. These results were consistent with those obtained using the first-principles B3LYP/6-311++G(d,p) method. In contrast, the OPLS/AA-FCFF and AMBER-FCFF 99sb force fields could only locate fewer or even fewer stable conformations.
In 2011, Zhao et al. ［105］ compared the solvation free energies of 18 amino acid side chains (excluding GLY and PRO) calculated using a combined ABEEM–GB approach with results obtained from other methods ［124］. Statistical analysis of the deviations between these computational results and experimental data revealed that our computed solvation free energies exhibited the smallest deviation from experiment, at 0.17 kcal/mol. This finding indicates that the ABEEM method, when coupled with the GB model, is highly suitable for calculating solvation free energies.
To further test the correction coefficients used in the ABEEM model for calculating charges in a continuous solvent, we selected 20 dipeptides with a total of 120 stable conformations. Under the PCM continuum solvent model, the solvation free energies of these 120 conformations were computed using the MP2/AUG-CC-PVDZ and B3LYP/AUG-CC-PVDZ quantum mechanical methods. We then analyzed the linear relationships between the solvation free energies of the dipeptides calculated by ABEEM/GB and those obtained from the two quantum mechanical methods. Specifically, the linear correlation coefficient between ABEEM/GB and MP2/AUG-CC-PVDZ was 0.9951, with a slope of 1.0625; while that between ABEEM/GB and B3LYP/AUG-CC-PVDZ was 0.9945, with a slope of 0.9783. It is evident that our method’s ability to predict the solvation free energies of dipeptides is in excellent agreement with both ab initio and DFT approaches, thereby providing strong support for subsequent simulations of biological systems.
In 2011, Chen et al. ［125］combined the ABEEM polarizable molecular force field with the LIE method to develop the ABEEM/MM LIE method for estimating absolute protein–ligand binding free energies. The parameters of this method were calibrated using a training set consisting of 24 HIV-1 protease–inhibitor complexes, resulting in both a two-parameter model and a three-parameter model that computed the interaction energies and solvation energies for these 24 receptor–ligand complexes. The results indicated that the three-parameter model exhibited higher predictive accuracy than the two-parameter model, with RMSDs and mean absolute deviations of 0.77 and 0.70 kcal/mol, as well as 0.55 and 0.56 kcal/mol, respectively. To avoid overfitting, leave-one-out cross-validation was performed on the models. The cross-validation results showed that the three-parameter model yielded RMSDs of 0.82 and 0.78 kcal/mol for and , while the two-parameter model had RMSDs of 0.91 and 0.93 kcal/mol. The multiple correlation coefficients were 0.67 and 0.87 for the two-parameter and three-parameter models, respectively. Furthermore, the third parameter = 8.67 kcal/mol fell within the experimentally observed range of 7 to 11 kcal/mol ［126–127］. The linear correlations between the free energies predicted by the two-parameter and three-parameter models for the 24 complexes and the experimental data were 0.91 and 0.93, respectively, further demonstrating the superior predictive accuracy of the three-parameter model.
To further validate the model's predictive capability, the model was applied to seven different HIV-1 protease–inhibitor complexes, and the results were also satisfactory. The RMSD between the binding energies calculated by the three-parameter model and the experimental values was 1.24 kcal/mol. At present, this method exhibits good generalizability for relative binding free energy prediction; however, the parameters need to be recalibrated according to different targets. The radius of gyration and hydrogen-bond distance can be used as auxiliary indicators for evaluating binding energy.
In 2014, Yang et al.［23］developed ABEEM parameters for all transition metal atoms in metalloproteins. By introducing reference charges to account for the potentially large charges of transition metals, they calibrated the valence state electronegativity and valence state hardness parameters using linear regression and the method of least squares. Using 700 biomolecules containing transition metal atoms as model systems, they optimized the charge-related parameters and validated them by comparison with ab initio results at the HF/STO-3G level. Taking 25 metalloproteins with molybdenum in different oxidation states as examples, they established a linear relationship between the ABEEM charges of the model molecules and those calculated by the ab initio method (HF/STO-3G). All slope values were close to 1.0, and the correlation coefficient R exceeded 0.97 in every case, indicating an excellent linear correlation between the two charge sets. Furthermore, using five metalloproteins with PDB IDs 3bdja, 1ti4, 1vlfj, 1g8ja, and 1g8ka as examples, they compared the structures optimized by the ABEEM polarizable molecular force field with those obtained using the MP2/6-311G ab initio method. The results showed that the two sets of structures were highly similar, with root-mean-square deviations (RMSDs) in bond lengths of 0.0396, 0.2158, 0.2008, 0.1309, and 0.0238 Å, respectively, and RMSDs in bond angles of 7.93°, 8.91°, 9.25°, 6.96°, and 6.06°, respectively. These findings demonstrate that the ABEEM polarizable molecular force field can be effectively used to optimize the geometries of metalloproteins, thereby laying a foundation for subsequent molecular dynamics simulations of metalloproteins.
3.1.4 Nucleic Acids
DMP is an important model molecule for the phosphodiester bond in nucleic acid molecules. In 2009, Wang et al.［128］constructed a DMP model containing 35 sites using the ABEEM method and investigated the properties of gas-phase DMP as well as DMP–(H₂O)ₙ(n=1–3) clusters. Compared with the QM method, the ABEEM method yielded RMSDs of 0.027 Å, 2.4°, and 1.1° for bond lengths, bond angles, and dihedral angles, respectively, when optimizing the DMP structure. The hydrogen-bonding structures and interaction energies of DMP–(H₂O)ₙ(n=1–3) clusters were analyzed. Specifically, the RMSDs for the hydrogen-bond angle and hydrogen-bond length of the DMP–H2O dimer optimized by the ABEEM and QM methods were 1.2° and 0.020 Å, respectively. Furthermore, the ABEEM method predicted the interaction energies of 11 DMP–(H₂O)ₙ (n=1–3) clusters, with deviations from MP2/6-311++G** calculations not exceeding 1.94 kcal/mol, indicating that the ABEEM PFF can accurately describe the hydrogen-bonding interactions in DMP–H2O complexes.
Through the parameterization of the ABEEM PFF potential function, molecular dynamics (MD) simulations were conducted for an aqueous solution of DMP anions in an explicit water environment. It was found that a distinct first hydration shell formed around the non-esterified oxygen atoms of the DMP anion, with an average coordination number of 5.2［128］. Notably, we observed that two hydrogen atoms of a single water molecule could simultaneously form bifurcated hydrogen bonds with two non-esterified oxygen atoms of the DMP anion. This represents the most significant difference between the ABEEM PFF and AMBER ff99［13］force fields in MD simulations, and it also demonstrates that the floating charge force field can provide richer information compared to traditional fixed charge force fields.
The biological functions of nucleic acids (such as DNA) are closely related to their conformation. In 2009, Wang et al.［129］obtained the dihedral torsion potential energy parameters for the ABEEM PFF sugar–phosphate backbone and glycosidic bond, and enhanced the description of hydrogen-bonding interactions through a hydrogen-bond fitting function. The ABEEM PFF can effectively describe the structure, interaction energies, and dipole moments of nucleic acid base molecules. For methylated bases (adenine A, guanine G, cytosine C, and thymine T), the root-mean-square deviation (RMSD) between the optimized non-hydrogen atomic bond lengths obtained using the ABEEM method and crystallographic data is 0.005–0.006 Å, which is comparable to the results from the CHARMM27 force field［20］. The RMSD for hydrogen-related bond angles ranges from 0.4° to 3.2°; in particular, guanine exhibits a relatively large deviation (3.2°) due to the planarization assumption of the amino group, whereas other bases show smaller deviations (e.g., thymine only 0.4°). This indicates that the ABEEM PFF can accurately represent the structural characteristics of nucleic acid base molecules. Furthermore, the calculated dipole moments of methylated bases using the ABEEM PFF are in good agreement with experimental values and high-precision theoretical methods［130–132］, with an RMSD of 0.27D relative to MP2 theory, outperforming the CHARMM27 force field［61］, which has an RMSD of 0.84D. Meanwhile, the ABEEM PFF can accurately reproduce the potential energy surfaces of key dihedral angles in nucleic acid systems［129］. For example, the dihedral angle γ exhibits three minima at 43.1°, 191.0°, and 298.4° in the C3'-endo conformation, which are essentially consistent with the values of 44.1°, 191.1°, and 298.3° obtained from ab initio calculations. Analysis of the potential energy surface for the glycosidic bond dihedral angle χ reveals that the anti conformation is more stable, in line with the crystal structures.
Through the study of 38 structures and their energies of Uracil-(H2O)n (n = 1~7) clusters[31], combined with high-accuracy quantum chemical calculations (MP2/6-311++G**), the hydrogen-bonding interactions between water molecules and uracil were investigated. By comparing the interaction energies obtained from the ABEEM PFF method and those from *ab initio* calculations, as shown in Figure 6, the root-mean-square deviation (RMSD) of all energy values is 0.96 kcal/mol, with a linear correlation coefficient of 0.997, indicating that the ABEEM PFF model can effectively simulate the interactions between the base and water.
Furthermore, in 2014, Liu et al. ［133］carried out a study on the charge distribution, energy, and structure of the 8-oxo-adenine (8-oxo-A) base pair, using quantum mechanical calculations at the B3LYP/6-31+G* level as a benchmark. Compared with the quantum mechanical results, the linear correlation coefficients for charge distribution, dipole moment, and interaction energy obtained from ABEEM PFF calculations reached 0.99, 0.96, and 0.96, respectively, while the root-mean-square deviations (RMSDs) for bond lengths and bond angles were only 0.005 Å and 0.5°, respectively. Building upon these findings, in 2021 they further investigated the molecular mechanism by which the repair enzyme hOGG1 recognizes oxidized bases (8-oxo-G) through hydrogen bonding and π–π stacking interactions ［30］,with particular emphasis on elucidating the dynamic mediating role of water molecules during the recognition process. Specifically, Gln315 in hOGG1 forms hydrogen bonds with 8-oxo-G via the sequential exchange of two water molecules (A and B), while Phe319 and Cys253 engage in bottom-side π–π stacking on either side of the planar base, thereby forming a sandwich-like structure. The charge polarization effect provides an important signal to drive the water molecule exchange and maintain the enzyme’s recognition of the oxidized base. Together, these two mediating water molecules pull Gln315 toward the recognition of 8-oxo-G through hydrogen-bonding interactions. The charge polarization effect of solvent water molecule C (with an absolute charge increase of approximately 0.2 e) facilitates the exchange between A and B, thereby sustaining the enzyme’s recognition of the oxidized base. By comparing the results of quantum mechanical (QM) calculations with those obtained from simulations using different force fields, the accuracy of ABEEM PFF in terms of geometric structure, charge distribution, and hydrogen-bonding energies was validated. When the geometric parameters (bond lengths, bond angles, and dihedral angles) of the optimized model molecules were compared across ABEEM PFF, OPLS/AA, and AMBER OL15 force fields ［134］and QM (using the M06-2X/6-311++G** method), ABEEM PFF exhibited the smallest RMSDs (0.017 Å for bond lengths, 2.20° for bond angles, and 2.26° for dihedral angles), indicating its superior fit to the geometric structure. In contrast, AMBER OL15 showed a relatively large deviation in dihedral angle simulation (with an RMSD of 9.50°), leading to distortions in both the methyl conformation and the planar structure. The average absolute deviation (AAD) between the dipole moments calculated by ABEEM PFF and those obtained from QM was 1.45 D,which was significantly lower than the 3.70 D and 2.95 D observed for OPLS/AA and AMBER OL15, respectively (representing deviations 2–3 times greater) ［134］. Moreover, the AAD and RMSD of hydrogen-bonding energies calculated by ABEEM PFF, as compared to those obtained using the MP2/aug-cc-pVDZ method, were 2.61 kcal/mol and 3.30 kcal/mol, respectively—both markedly better than the corresponding values for OPLS/AA and AMBER OL15 ［134］. For OPLS/AA, the AAD and RMSD were 10.35 kcal/mol and 15.65 kcal/mol, respectively, whereas for AMBER OL15, they were 13.50 kcal/mol and 25.52 kcal/mol. These results demonstrate that ABEEM PFF can accurately describe molecular geometry, electrostatic interactions within hydrogen-bonding regions, and charge distribution.
3.2 An Example of Molecular Dynamics Simulation of Biomolecular Aqueous Solutions Using ABEEM PFF
3.2.1 Molecular Dynamics Simulations of GA88 and GB88
In 2013, Wang et al. ［135］selected two molecules, GA88 (PDB code: 2JWS) and GB88 (PDB code: 2JWU), and performed molecular dynamics simulations in aqueous solution at 295 K and in vacuum using the ABEEM PFF and OPLS-AA fixed-charge force field methods. The root-mean-square deviations between the structures of the two proteins obtained from the two methods and the experimental structures were compared, and the radii of gyration, hydrogen bonding, and charge distributions were analyzed. The results showed that the RMSDs of non-hydrogen atoms in GA88 and GB88 obtained with ABEEM PFF were smaller than those obtained with the OPLS-AA fixed-charge force field; the radius of gyration simulations well captured the protein’s “electrostatic compaction” phenomenon. The RMSDs between the average hydrogen bond lengths in the NMR structures and those simulated by ABEEM PFF for proteins GA88 and GB88 were 0.012 nm and 0.008 nm, respectively, while the corresponding values obtained with the OPLS-AA force field were 0.012 nm and 0.011 nm, which were essentially consistent. The ABEEM model not only considers atomic-region charges but also incorporates charges associated with chemical bonds and lone-pair electron regions, allowing the charges of each region to be recalculated promptly as the system environment changes. In contrast, when simulations are conducted using the OPLS-AA fixed-charge force field, the charges of all atoms remain unchanged regardless of the environmental conditions. Therefore, ABEEM PFF is better able to simulate protein structures and capture the polarization effects between proteins and surrounding water molecules.
3.2.2 Molecular Dynamics Simulation of Boric Acid-β-Lactamase
Boric acid, as a β-lactamase inhibitor, holds potential in the treatment of bioinfections and tumors. In 2020, Lu et al.［136］developed the ABEEM PFF force field tailored for boron-containing systems, thereby filling the gap in existing force fields that lacked boron-specific parameters. The ABEEM PFF is advantageous in providing a reasonable charge distribution, charge transfer, and specialized treatment of hydrogen-bonding interactions. By selecting five classes comprising 47 boric acid–β-lactamase complexes as model molecules, the relevant parameters of the ABEEM PFF were determined and validated. Comparison with quantum mechanical (QM) calculations revealed that the linear correlation coefficient between the charges computed by the two methods exceeded 0.97; the average absolute deviation (AAD) of the dipole moments was 0.97D, and the root-mean-square deviation (RMSD) was only 1.13D, indicating that the charge distribution calculated by the ABEEM PFF is reasonable. A comparison of the optimized geometries obtained from both methods confirmed the rationality of the Morse parameters, bond-angle parameters, and dihedral-angle parameters within the ABEEM PFF. The binding energies of ten boric acid–β-lactamase model complexes and ten boric acid–β-lactamase–water complexes were investigated, and the results obtained using the ABEEM PFF exhibited good correlation with those derived from the MP2 method. Building upon these findings, molecular dynamics (MD) simulations were conducted on five large-scale boric acid–β-lactamase complexes (BZB-1C3BJ, CTB-1FSWJ, NPB-1KDSJ, SM3-1MY8J, and 0NB-4E3IJ) over a simulation time of 4.0 ns. The structures generated from the MD simulations were compared with their initial structures as deposited in the Protein Data Bank (PDB). The root-mean-square deviations (RMSDs) of bond lengths were 0.041, 0.039, 0.041, 0.030, and 0.044 Å, respectively, while the RMSDs of bond angles were 3.89°, 3.03°, 3.99°, 3.73°, and 3.66°, respectively. These results demonstrate that the ABEEM PFF can accurately simulate stable structures of boric acid–β-lactamase complexes. However, due to the relatively low energy barriers associated with dihedral torsions, these complexes are highly susceptible to environmental influences, leading to significant deviations from their initial conformations. This observation underscores the ABEEM PFF’s capability to effectively account for polarization effects arising in realistic environments.
3.2.3 Molecular Dynamics Simulations of Magnesium-Ion Proteins and Magnesium-Ion Aqueous Solutions
Magnesium ions are important biological metal ions that participate in many vital physiological processes. In 2025, Zhang［137］et al. developed the ABEEM PFF bonding model (ABEEM-BM) and non-bonding model (ABEEM-NBM) for magnesium ions to simulate their behavior in magnesium–protein complexes and aqueous solutions. Forty magnesium–protein fragments were extracted from PDB entries, retaining the magnesium ion, ligand amino acid side chains, and water molecules; constrained optimization was then employed to obtain stable structures. Using the ABEEM-BM, ABEEM-NBM, OPLS-AA［38］, AMBER99［138］, and CHARMM22［35］force fields, 10-ns MD simulations were performed on these 40 magnesium–protein fragments in aqueous solution. Compared with the PDB structures, the RMSDs (bond lengths/angles) obtained from the MD simulations showed that the average bond length RMSDs for ABEEM-BM, ABEEM-NBM, OPLS-AA, AMBER99, and CHARMM22 were 0.088, 2.490, 2.697, 2.831, and 2.594 Å, respectively. The corresponding average angle RMSDs were 5.99°, 24.38°, 24.96°, 25.86°, and 24.01°. The bond length and bond angle RMSDs of ABEEM-BM were significantly lower than those of ABEEM-NBM, OPLS-AA, AMBER99, and CHARMM22 force fields, indicating that ABEEM-BM can yield a stable conformation of magnesium–protein molecules.
Using ABEEM-BM and ABEEM-NBM to perform MD simulations of magnesium ion aqueous solutions, we investigated the dynamic properties of magnesium ions in aqueous solution. The analysis of the radial distribution function (RDF) yielded the following results: the first peak of the Mg—O RDF in ABEEM-BM appeared at 2.050 Å, which is consistent with the experimental value (2.090 Å)［139］and the AIMD result (2.100 Å)［140］; whereas in ABEEM-NBM, it appeared at 1.850 Å. The integrated areas of the first Mg—O peaks obtained from both methods were 6, indicating that the coordination number of magnesium ions in aqueous solution is 6, which matches the experimental and AIMD results. The water exchange rate constants for ABEEM-BM and ABEEM-NBM were 4.7×105and 5.0×106 s-1, respectively. ABEEM-BM was closer to the experimental value (5.3×105~6.7×105 s⁻¹), with an error margin smaller than that of microMg (TIP3P)［141］. The average residence time (τ) of water molecules within the first hydration shell of magnesium ions was experimentally determined to be 1.0×10-6 s［142］. Compared with ABEEM-NBM (2.0×10-7 s), the τcalculated using ABEEM-BM (2.1×10-6 s) was closer to the experimental data. These results demonstrate that ABEEM-BM, while maintaining computational efficiency, has improved the accuracy of kinetic simulations of magnesium ion aqueous solutions. Figure 7illustrates how the ABEEM model dynamically captures the charge transfer process between magnesium ions and ligands; the partial charge range of magnesium ions was 1.07–1.24|e| (ABEEM-BM), which is in excellent agreement with the QM results (1.06–1.42|e|). This work lays a foundation for future studies utilizing the ABEEM PFF to investigate the structure and properties of hydrated magnesium proteins.
3.2.4 Heme-containing enzymes catalyze the hydroxyamination of styrene to produce enantiomers
In the design of heme-containing enzymes that regulate catalytic reactivity and selectivity, steric effects and electrostatic effects play crucial roles, while polarizable force fields can effectively describe the changes in electrostatic effects arising from structural alterations within this regulatory mechanism［143］.
In 2024, Huang et al.［144］focused on elucidating the critical functional role of designed cytochrome c(P450) in catalyzing the hydrolysis of carbon-centered Fe–porphyrin complexes. By fitting charge parameters of ABEEM PFF associated with the ammonia hydroxylation reaction using model molecules and employing the ModSeminario［145］method combined with a scanning approach, bond-length stretching and bond-angle bending parameters were determined. The interaction between Fe and its axial ligand was incorporated into a Morse potential function for description. Conformational searches were conducted to identify the structure of the reactive intermediate. Mutated residues imposed constraints on the reactive intermediate via van der Waals interactions and hydrogen bonding involving residues 100S and 101P, thereby influencing its structural orientation. Simultaneously, hydrogen-bonding interactions between the reactive intermediate and water molecules collectively facilitated the formation of the Rproduct. Furthermore, QM/MM-ABEEM PFF and FEP calculations revealed that, in the 3Rts structure, a stable sandwich-like arrangement could form between the reaction complex and neighboring amino acid residues. Compared with the crowded packing observed in the 3Sts configuration, the 3Rts structure exhibited greater stability (Figure 8). The free-energy barrier for the Renantiomer was approximately 10.0 kcal/mol lower than that for the Senantiomer, rendering the formation of the Renantiomer more thermodynamically favorable［144］. These findings provided mechanistic insights into designing cytochrome cwith superior stereoselectivity for the styrene-ammonia hydroxylation reaction.
3.2.5 DNA Variation and Repair Mechanisms
The repair enzymes within living organisms have a restorative effect on damaged DNA; their essence lies in the enzymatic catalytic reaction between the enzyme and the damaged DNA. Elucidating the microscopic mechanism of this reaction will provide assistance for the treatment of related diseases. In 2020, Liu et al.［146］performed molecular dynamics simulations under aqueous conditions using ABEEM PFF to investigate more than twenty reaction pathways associated with the excision of damaged bases by four different repair enzymes. All these reactions are concerted SN2reactions characterized by glycosidic bond cleavage coupled with proton transfer. The ABEEM PFF defines a reactivity-valence electronegativity function that accurately simulates the dramatic changes in atomic charges occurring during the reaction—specifically, those involving atoms participating in bond formation and bond breaking, as well as the surrounding atoms near the transition state. By applying the condition of local charge conservation throughout these changes, the simulation successfully reproduced an approximately −1.0e charge transfer process from the active-site residues to the damaged base. The results of the molecular dynamics simulations indicate that a protonated base serves as a superior leaving group. The activation energy for the enzyme-catalyzed excision of a protonated base is generally 10 kcal/mol lower than that for the excision of a neutral base. For damaged guanine［30］and damaged adenine［147］,the repair enzymes employ a bifunctional mechanism to attack the damaged base and cleave the glycosidic bond; whereas for damaged cytosine［148］and damaged thymine［149］,the repair enzymes utilize a monofunctional mechanism to achieve the same purpose. Figure 9illustrates this process. The ABEEM polarizable force field has provided a novel theoretical approach for investigating the reaction mechanisms underlying DNA glycosylase-mediated excision of damaged bases in aqueous environments. We plan to conduct simulations of the entire protein–DNA system in aqueous solution to explore the influence of distal residues on enzyme recognition, thereby revealing the fundamental nature of enzyme-specific recognition.
3.2.6 Mechanism of Deprotonation of Guanine Radical Cation (G·+)
The rapid deprotonation of G·+ impedes long-range hole (positive charge) transport in DNA strands. Investigating its mechanism can facilitate the design of highly conductive DNA structures［150］. In 2024, Wang et al.［150］conducted a detailed analysis using QM and QM/MM (ABEEM) methods within 3H2O and 9H2O models, focusing on charge distribution, charge transfer, geometric structure, and activation energy (Figure 10).
As shown in Figure 10, the proton transfer distances and activation energies obtained from the QM method for the 3H2O and 9H2O models are 0.469 Å and 36.6 kJ/mol, and 0.391 Å and 26.0 kJ/mol, respectively. The relative energy of G·+ in the 3H2O model is significantly higher than that in the 9H2O model, indicating that proton transfer is more facile in the latter. Moreover, the presence of a second hydration shell in the 9H2O model lowers the potential energy barrier, thereby facilitating proton transfer. The microscale deprotonation mechanism of G·+ in the 9H2O model was further simulated using the QM/MM (ABEEM) method. By employing piecewise functions based on valence-state electronegativity (χ*) and incorporating charge local conservation conditions, this method can accurately simulate the effects of polarization and charge transfer. We fitted a linear relationship between the charge distributions calculated by the QM and QM/MM (ABEEM) methods, with a correlation coefficient exceeding 0.97, demonstrating consistency in charge transfer between the two approaches. MD simulations revealed a deprotonation activation energy of (16.3 ± 0.8) kJ/mol, which is consistent with the experimental result of (15.1 ± 1.5) kJ/mol[151]. These findings highlight the successful application of the QM/MM (ABEEM) method in simulating proton transfer. The deprotonation rate of G·+ in free bases (1.5 × 107 s−1) is faster than that in double-stranded DNA (106 to 107 s−1). Furthermore, long-range proton transfer corresponds to a high-energy potential barrier. When designing DNA-based electronic devices, it is advisable to avoid incorporating free bases and instead prioritize molecules that exhibit long-range proton transfer involving G·+. The aforementioned mechanisms provide valuable guidance for the rational design of electronic device architectures.
3.2.7 Molecular Dynamics Simulation of the Lipid Bilayer
The structure of the cell membrane is based on a lipid bilayer, which is crucial for various biological processes. To gain a deeper understanding of its properties and interaction mechanisms, extensive research has been conducted in the academic community. In 2025, Wang et al.［152］developed an ABEEM PFF model for PC lipids (Figure 11),which successfully reproduces key experimental characteristics of the lipid bilayer, including volume, thickness, area, NMR order parameters, and electron density.
For volume and area, the DMPC simulation values are 60.5 Å2and 1230 Å3, respectively, which are consistent with the experimental values［153］of 60.6 Å2and 1229 Å3; compared with the experimental results, the carbon atom deuterium order parameters optimized using ABEEM PFF are significantly superior to those obtained with Drude PFF［154-155］; the electron density distributions (EDPs) on the bilayer exhibit excellent correlation with the experimental X-ray scattering data, with an RMSD of 0.05 e/Å3. The predicted dipole potential in the bulk water region is 0.403 V, which is an important characteristic reflecting lipid polarizability. ABEEM PFF effectively captures charge fluctuations and provides a more accurate representation of charge distribution. Based on this feature, we employed unsupervised machine learning methods to investigate the correlation between the charge distribution of P atoms in lipids and their conformational states, revealing that the orientation of the P-N vector relative to the membrane normal exhibits charge dependence. Furthermore, our analysis indicates that the charge of P atoms in the lipid bilayer influences the dynamic behavior of water molecules: the higher the P charge, the greater the degree of water molecule aggregation. Advancing this model requires concerted efforts in algorithm development, software engineering, and practical applications. Currently, the computational cost of non-parallelized ABEEM PFF simulations is indeed higher than that of classical force fields (FF). However, by employing the ABEEM-DBSS［146］method, for all considered lipid bilayers (ranging from 14,112 to 16,704 atoms, including 72 lipids and 2,160 water molecules), the ABEEM PFF code can generate 1 ns of data on an AMD Ryzen 9 9950X approximately four times faster than AMOBEA 09. This underscores the critical importance of program parallelization for applying ABEEM PFF to the simulation of complex lipid membranes and membrane proteins.
4 Conclusion and Outlook
The study of biomolecular aqueous solutions holds significant importance in life sciences, medicine, and industry. Molecular dynamics simulations serve as a crucial tool for investigating the structure and dynamic behavior of biomolecules in liquid environments. This paper takes ABEEM PFF as the primary subject of study, providing a detailed introduction to the ABEEM-7P water model and the application of ABEEM PFF in biomolecular systems. Through several specific examples, it illustrates the use of ABEEM PFF in molecular dynamics simulations of biomolecular aqueous solutions.
(1) The ABEEM model contains atomic, bond, and lone-pair electron sites. The multi-site model reflects the anisotropy of charge distribution around atoms and particularly considers hydrogen-bonding interaction regions. By solving linear equations related to the molecular geometry, the charge distribution of the target molecule can be determined without a self-consistent iterative process, thereby enabling rapid acquisition of the charge distribution. The ABEEM-7P water model is crucial for understanding hydrogen-bonding interactions among water molecules and thus can accurately describe the electrostatic interactions in aqueous solutions.
(2) ABEEM PFF can not only effectively describe the structure and conformational changes of biomolecules in aqueous solutions, but also calculate their solvation energies and free energies. During simulations, atomic charges in ABEEM PFF fluctuate with the surrounding environment, thereby accurately capturing the polarization effects and charge transfer within molecular systems. It has been successfully applied to aqueous solution systems of biomolecules such as water clusters, proteins, nucleic acids, and lipids.
In summary, this paper has demonstrated through a comprehensive review of extensive related research that the ABEEM-7P water model and ABEEM PFF can be effectively employed in molecular dynamics simulations of biomolecular aqueous solutions, yielding satisfactory results. Nevertheless, several pressing issues remain to be addressed. The development of algorithms and software, as well as the creation of new models and further application studies, are of paramount importance for the advancement of our modeling framework. In future work, integrating machine learning techniques with the ABEEM model will also constitute a significant direction of development.
PDF
