Molecular Mechanics Force Field Research Articles

Reweighting from Molecular Mechanics Force Fields to the ANI-2x Neural Network Potential.

To achieve chemical accuracy in free energy calculations, it is necessary to accurately describe the system's potential energy surface and efficiently sample configurations from its Boltzmann distribution. While neural network potentials (NNPs) have shown significantly higher accuracy than classical molecular mechanics (MM) force fields, they have a limited range of applicability and are considerably slower than MM potentials, often by orders of magnitude. To address this challenge, Rufa et al. [Rufa et al. bioRxiv 2020, 10.1101/2020.07.29.227959.] suggested a two-stage approach that uses a fast and established MM alchemical energy protocol, followed by reweighting the results using NNPs, known as endstate correction or indirect free energy calculation. This study systematically investigates the accuracy and robustness of reweighting from an MM reference to a neural network target potential (ANI-2x) for an established data set in vacuum, using single-step free-energy perturbation (FEP) and nonequilibrium (NEQ) switching simulation. We assess the influence of longer switching lengths and the impact of slow degrees of freedom on outliers in the work distribution and compare the results to those of multistate equilibrium free energy simulations. Our results demonstrate that free energy calculations between NNPs and MM potentials should be preferably performed using NEQ switching simulations to obtain accurate free energy estimates. NEQ switching simulations between the MM potentials and NNPs are efficient, robust, and trivial to implement.

PoseBusters: AI-based docking methods fail to generate physically valid poses or generalise to novel sequences.

The last few years have seen the development of numerous deep learning-based protein-ligand docking methods. They offer huge promise in terms of speed and accuracy. However, despite claims of state-of-the-art performance in terms of crystallographic root-mean-square deviation (RMSD), upon closer inspection, it has become apparent that they often produce physically implausible molecular structures. It is therefore not sufficient to evaluate these methods solely by RMSD to a native binding mode. It is vital, particularly for deep learning-based methods, that they are also evaluated on steric and energetic criteria. We present PoseBusters, a Python package that performs a series of standard quality checks using the well-established cheminformatics toolkit RDKit. The PoseBusters test suite validates chemical and geometric consistency of a ligand including its stereochemistry, and the physical plausibility of intra- and intermolecular measurements such as the planarity of aromatic rings, standard bond lengths, and protein-ligand clashes. Only methods that both pass these checks and predict native-like binding modes should be classed as having "state-of-the-art" performance. We use PoseBusters to compare five deep learning-based docking methods (DeepDock, DiffDock, EquiBind, TankBind, and Uni-Mol) and two well-established standard docking methods (AutoDock Vina and CCDC Gold) with and without an additional post-prediction energy minimisation step using a molecular mechanics force field. We show that both in terms of physical plausibility and the ability to generalise to examples that are distinct from the training data, no deep learning-based method yet outperforms classical docking tools. In addition, we find that molecular mechanics force fields contain docking-relevant physics missing from deep-learning methods. PoseBusters allows practitioners to assess docking and molecular generation methods and may inspire new inductive biases still required to improve deep learning-based methods, which will help drive the development of more accurate and more realistic predictions.

Development and Comprehensive Benchmark of a High-Quality AMBER-Consistent Small Molecule Force Field with Broad Chemical Space Coverage for Molecular Modeling and Free Energy Calculation.

Biomolecular simulations have become an essential tool in contemporary drug discovery, and molecular mechanics force fields (FFs) constitute its cornerstone. Developing a high quality and broad coverage general FF is a significant undertaking that requires substantial expert knowledge and computing resources, which is beyond the scope of general practitioners. Existing FFs originate from only a limited number of groups and organizations, and they either suffer from limited numbers of training sets, lower than desired quality because of oversimplified representations, or are costly for the molecular modeling community to access. To address these issues, in this work, we developed an AMBER-consistent small molecule FF with extensive chemical space coverage, and we provide Open Access parameters for the entire modeling community. To validate our FF, we carried out benchmarks of quantum mechanics (QM)/molecular mechanics conformer comparison and free energy perturbation calculations on several benchmark data sets. Our FF achieves a higher level of performance at reproducing QM energies and geometries than two popular open-source FFs, OpenFF2 and GAFF2. In relative binding free energy calculations for 31 protein-ligand data sets, comprising 1079 pairs of ligands, the new FF achieves an overall root-mean-square error of 1.19 kcal/mol for ΔΔG and 0.92 kcal/mol for ΔG on a subset of 463 ligands without bespoke fitting to the data sets. The results are on par with those of the leading commercial series of OPLS FFs.

Electron Spin Resonance Spectrum Simulations of Possible Radicals of the Ketoprofen Molecule

In this study, the molecular structure of Ketoprofen molecule, which is the drug active ingredient , was revealed by using the Density Functional Theory method and the Molecular Mechanical Force Field method. First of all, conformational space scanning in the Ketoprofen molecule was performed by the Molecular Mechanical Force Field method. The most stable structure of the ketoprofen molecule was found with the help of the Density Functional Theory method. The geometry parameters of the ketoprofen molecule and the Nuclear Magnetic Resonance parameters were calculated with the help of the Density Functional Theory method. Possible radicals were modeled using the most stable structure of the ketoprofen molecule. Electron Spin Resonance parameters of these possible radicals were calculated with the Density Functions Theory method. The calculated Electron Spin Resonance parameters were used in the JEOL IsoSimu/Fa Version 2.2.0 simulation program and theoretical Electron Spin Resonance spectra of possible radicals were obtained.

Simulating the Fluorescence of the Locally Excited State of DMABN in Solvents of Different Polarities.

4-(N,N-Dimethylamino)benzonitrile (DMABN) is a luminescent probe that can be used for tracking changes in the surrounding solvent due to the large change in polarity between its ground and excited states. An important characteristic of DMABN is that it exhibits dual fluorescence with two different emission energies that can be monitored, allowing for better characterization of the surrounding system. The first excited state is called the locally excited (LE) state and is characterized by the movement of charge over the conjugated ring structure. In nonpolar solvents and in the gas phase, the fluorescence of DMABN is entirely attributed to the transition from the near-planar LE state. In more polar environments, emission occurs from both the LE and a second excited state, corresponding to a twisted intramolecular charge-transfer (ICT) structure. For the sake of simplicity, this work considers transitions between only the ground and LE state. Molecular mechanical force field models of DMABN in its ground and LE states have been developed to investigate the sensitivity of the LE state to the polarity of the solvent. Both nonpolarizable and polarizable force fields were developed to simulate the molecule in a series of 10 different solvents of different polarities. The calculated Stokes shift of DMABN increases with increasing orientation polarizability of the surrounding solvent, which is the expected trend, as seen in experimental studies.

Simulating the Fluorescence of Di-8-ANEPPS in Solvents of Different Polarity.

Voltage-sensitive fluorescent probes, such as di-8-amino-naphthyl-ethylene-pyridinium-propyl-sulfonate (di-8-ANEPPS), are extremely useful for monitoring the membrane potential in biological and biophysical studies. However, because di-8-ANEPPS is very sensitive to its environment, it can be difficult to distinguish the degree to which a given external factor affects the observed fluorescence. Molecular dynamics simulations based on detailed atomic models make it possible to examine the particular characteristics of the system and predict the effects of the surroundings. Here, the sensitivity of the spectra of di-8-ANEPPS to solvent polarity is investigated by modeling the electronic transition between the ground and excited states using classical molecular mechanical force fields. The absorption and emission of di-8-ANEPPS were simulated in 12 solvents of increasing polarity using nonpolarizable and polarizable force fields. While the computational results and experimental data do not match perfectly, classical Lippert plots of both models show the expected increase of the Stokes shift of di-8-ANEPPS with the orientation polarizability of the surrounding solvent.

In silico pharmacological study of lacourtianal, a new terpenoid isolated from the stem bark of Chrysophyllum lacourtianum De Wild (Sapotaceae) against Alzheimer's disease

Urgent need to treat, prevent, delay the onset, slow the progression and reduce the symptoms of Alzheimer's disease (AD) is dictated by the growing number of people affected by the disease. This research requires an in silico approach. The aim of this work was to carry out an in silico pharmacological study of lacourtianal, a new terpenoid against Alzheimer's disease. The structure of lacourtianal and reference drugs were drawn on chemdraw 2D. Energy minimization was performed using the molecular mechanics force field (MM2) and saved in the PDB using chemdraw 3D. The compound's Smiles format was entered into the pk-CSM web servers for ADME/Tox parameter prediction. Conformational site analysis and docking parameters such as binding energy, interaction profiles with AD target residues (AchE, BuchE, β-secretase and GSK-3β) were determined using AutoDock 4.2 and Discovery Studio visualizer. Lacourtianal is a valuable oral drug candidate. It crosses the Blood-brain barrier (BBB) and shows considerable docking scores for AD targets. These scores were higher for BuchE and β-secretase compared with reference compounds. It binds to residues Phe297, Phe 338 and Tyr31, Trp286 in the acyl pocket and peripheral site of AchE, respectively. Importantly, it establishes low-energy interactions (hydrogen and pi) with the His438 residue of the Butyrylcholinesterase catalytic triad. It also establishes favorable interactions with residues Tyr71 and Tyr216, which are residues controlling substrate accessibility to the active site of the BACE1 and GSK 3β enzymes respectively. These results show that lacourtianal has promising therapeutic potential for the treatment of Alzheimer's disease.
 Keywords: Lacourtianal, in silico, AchE, BuchE, β-secretase, GSK-3β and AD

Boosting Graph Neural Networks with Molecular Mechanics: A Case Study of Sigma Profile Prediction.

Sigma profiles are quantum-chemistry-derived molecular descriptors that encode the polarity of molecules. They have shown great performance when used as a feature in machine learning applications. To accelerate the development of these models and the construction of large sigma profile databases, this work proposes a graph convolutional network (GCN) architecture to predict sigma profiles from molecule structures. To do so, the usage of molecular mechanics (force field atom types) is explored as a computationally inexpensive node-level featurization technique to encode the local and global chemical environments of atoms in molecules. The GCN models developed in this work accurately predict the sigma profiles of assorted organic and inorganic compounds. The best GCN model here reported, obtained using Merck molecular force field (MMFF) atom types, displayed training and testing set coefficients of determination of 0.98 and 0.96, respectively, which are superior to previous methodologies reported in the literature. This performance boost is shown to be due to both the usage of a convolutional architecture and node-level features based on force field atom types. Finally, to demonstrate their practical applicability, we used GCN-predicted sigma profiles as the input to machine learning models previously developed in the literature that predict boiling temperatures and aqueous solubilities. Using the predicted sigma profiles as input, these models were able to compute both physicochemical properties using significantly less computational resources and displayed only a slight decrease in performance when compared with sigma profiles obtained from quantum chemistry methods.

Adaptive-Partitioning Multilayer Dynamics Simulations: 2. Implementations of the Permuted and Interpolated Adaptive-Partitioning Gradients.

Recently, an adaptive-partitioning multilayer Q1/Q2/MM method was proposed, where Q1 and Q2 denote, respectively, two distinct quantum-mechanical levels of theory and MM, the molecular-mechanical force fields. Such a multilayer model resembles the ONIOM (our own N-layered integrated molecular orbital and molecular mechanics) model by Morokuma and co-workers, but it is distinguished by on-the-fly reclassifying atoms to be Q1, Q2, or MM in dynamics simulations. To smoothly blend the levels of descriptions of the atoms, buffer zones are introduced between adjacent layers, and the energy is smoothly interpolated. In particular, the Q1/Q2 interaction energy was expressed in two different formalisms: permuted and interpolated adaptive-partitioning (PAP and IAP), respectively. While the PAP energy is based on a weighted many-body expansion, the IAP energy is derived via alchemical quantum calculations with interpolated Fock and overlap matrices. In this article, we examine in-depth the irregularities in the IAP energy near the boundary between the buffer and Q2 zones, which were found prominent in some calculations. These irregularities are due to basis-set linear dependencies, which can be effectively suppressed using a cutoff for the weighted atomic orbital coefficients. Furthermore, we derived and implemented the gradients for both PAP and IAP. Test calculations on a series of water cluster models show perfectly smooth gradients in PAP, while a minor discontinuity occurs in IAP gradients at the buffer/Q2 boundary. The energy and gradient discontinuities in IAP become smaller when moving the buffer/Q2 boundary further away from the Q1 center and when increasing the size of the basis sets used. Overall, those discontinuities are controllable, and possible ways to further diminish them are discussed.

Simple, near-universal relationships between bond lengths, strengths, and anharmonicities

Harmonic bond force constants and bond lengths are shown to generally obey the simple relationships, ke=ζ2Re−3 (hydrides) and ke=10ζ1/2Re−4 (all other bond types), where ζ is the reduced nuclear charge and Re is the equilibrium bond length. Equally simple power-law relationships are found for higher-order bond force constants. Although not spectroscopically accurate, these models are nonetheless of significant heuristic value for identifying strongly multireference states of diatomic molecules (including electronically coupled excited states ill-suited for inclusion in laser-cooling schemes), rationalizing the observed trends in vibrational frequencies for diatomics and/or local mode oscillators within molecules or complexes and estimating and/or validating covalent bonding parameters within molecular mechanics force fields. Particular advantages of our approach over other bond length-strength scaling relationships proposed in the literature include its simplicity and generality and its appropriate asymptotic behavior. Notably, the relationships derived in this work can be used to predict harmonic and higher-order force constant bonds between any pair of atoms in the Periodic Table (including transition metals and lanthanides) without requiring row- or column-dependent parameterization, to accuracies commensurate with conventional force field transferability errors. We therefore anticipate that they will expedite force field development for metal-containing complexes and materials, which are structurally well-characterized but challenging to parameterize ab initio.

Structural characterization of an intrinsically disordered protein complex using integrated small-angle neutron scattering and computing.

Characterizing structural ensembles of intrinsically disordered proteins (IDPs) and intrinsically disordered regions (IDRs) of proteins is essential for studying structure-function relationships. Due to different neutron scattering lengths of hydrogen and deuterium, selective labeling and contrast matching in small-angle neutron scattering (SANS) becomes an effective tool to study dynamic structures of disordered systems. However, experimental timescales typically capture measurements averaged over multiple conformations, leaving complex SANS data for disentanglement. We hereby demonstrate an integrated method to elucidate the structural ensemble of a complex formed by two IDRs. We use data from both full contrast and contrast matching with residue-specific deuterium labeling SANS experiments, microsecond all-atom molecular dynamics (MD) simulations with four molecular mechanics force fields, and an autoencoder-based deep learning (DL) algorithm. From our combined approach, we show that selective deuteration provides additional information that helps characterize structural ensembles. We find that among the four force fields, a99SB-disp and CHARMM36m show the strongest agreement with SANS and NMR experiments. In addition, our DL algorithm not only complements conventional structural analysis methods but also successfully differentiates NMR and MD structures which are indistinguishable on the free energy surface. Lastly, we present an ensemble that describes experimental SANS and NMR data better than MD ensembles generated by one single force field and reveal three clusters of distinct conformations. Our results demonstrate a new integrated approach for characterizing structural ensembles of IDPs. This article is protected by copyright. All rights reserved.

Neural potentials of proteins extrapolate beyond training data.

We evaluate neural network (NN) coarse-grained (CG) force fields compared to traditional CG molecular mechanics force fields. We conclude that NN force fields are able to extrapolate and sample from unseen regions of the free energy surface when trained with limited data. Our results come from 88 NN force fields trained on different combinations of clustered free energy surfaces from four protein mapped trajectories. We used a statistical measure named total variation similarity to assess the agreement between reference free energy surfaces from mapped atomistic simulations and CG simulations from trained NN force fields. Our conclusions support the hypothesis that NN CG force fields trained with samples from one region of the proteins' free energy surface can, indeed, extrapolate to unseen regions. Additionally, the force matching error was found to only be weakly correlated with a force field's ability to reconstruct the correct free energy surface.

Development and application of fragment-based de novo inhibitor design approaches against Plasmodium falciparum GST.

Modulation of disease progression is frequently started by identifying biochemical pathway catalyzed by biomolecule that is prone to inhibition by small molecular weight ligands. Such ligands (leads) can be obtained from natural resources or synthetic libraries. However, de novo design based on fragments assembly and optimization is showing increasing success. Plasmodium falciparum parasite depends on glutathione-S-transferase (PfGST) in buffering oxidative heme as an approach to resist some antimalarials. Therefore, PfGST is considered an attractive target for drug development. In this research, fragment-based approaches were used to design molecules that can fit to glutathione (GSH) binding site (G-site) of PfGST. The involved approaches build molecules from fragments that are either isosteric to GSH sub-moieties (ligand-based) or successfully docked to GSH binding sub-pockets (structure-based). Compared to reference GST inhibitor of S-hexyl GSH, ligands with improved rigidity, synthetic accessibility, and affinity to receptor were successfully designed. The method involves joining fragments to create ligands. The ligands were thenexplored using molecular docking, Cartesian coordinate's optimization, and simplified free energy determination as well as MD simulation and MMPBSA calculations. Several tools were used which includeOPENEYE toolkit, Open Babel, Autodock Vina, Gromacs, and SwissParam server, and molecular mechanics force field of MMFF94 for optimization and CHARMM27 for MD simulation. In addition,in-house scripts written in Matlab were used to control fragments connectionand automation of the tools.

Protein-Ligand CH-π Interactions: Structural Informatics, Energy Function Development, and Docking Implementation.

Here, we develop an empirical energy function based on quantum mechanical data for the interaction between methane and benzene that captures the contribution from CH-π interactions. Such interactions are frequently observed in protein-ligand crystal structures, particularly for carbohydrate ligands, but have been hard to quantify due to the absence of a model for CH-π interactions in typical molecular mechanical force fields or docking scoring functions. The CH-π term was added to the AutoDock Vina (AD VINA) scoring function enabling its performance to be evaluated against a cohort of more than 1600 occurrences in 496 experimental structures of protein-ligand complexes. By employing a conformational grid search algorithm, inclusion of the CH-π term was shown to improve the prediction of the preferred orientation of flexible ligands in protein-binding sites and to enhance the detection of carbohydrate-binding sites that display CH-π interactions. Last but not least, this term was also shown to improve docking performance for the CASF-2016 benchmark set and a carbohydrate set.

The application of molecular dynamic simulations

Molecular dynamic(MD) simulations have been developed into a mature tool for researchers and have had a huge impact on the aspects of drug discovery, molecular biology, and prediction of specific structures and related functions of macromolecules in recent years. The aim of this article is to provide some essential knowledge about the theories and applications of this useful technology. This introduction to molecular dynamics starts from the definition of some basic concepts, including the molecular mechanics force fields and how they are applied to the motion systems of molecular dynamics simulations and their developments recently. Also, the outstanding strengths and limitations of the MD simulations are discussed, as are some possible solutions to them. This article shows the wide applications of MD simulations, including drug discovery, prediction of structure, and allostery. Due to the growing power of computational calculations and its important role in simulating the biological process between receptors and ligands like allosteric binding sites, MD simulations are thought to be more and more important in the near future.

Toward a general neural network force field for protein simulations: Refining the intramolecular interaction in protein.

Molecular dynamics (MD) is an extremely powerful, highly effective, and widely used approach to understanding the nature of chemical processes in atomic details for proteins. The accuracy of results from MD simulations is highly dependent on force fields. Currently, molecular mechanical (MM) force fields are mainly utilized in MD simulations because of their low computational cost. Quantum mechanical (QM) calculation has high accuracy, but it is exceedingly time consuming for protein simulations. Machine learning (ML) provides the capability for generating accurate potential at the QM level without increasing much computational effort for specific systems that can be studied at the QM level. However, the construction of general machine learned force fields, needed for broad applications and large and complex systems, is still challenging. Here, general and transferable neural network (NN) force fields based on CHARMM force fields, named CHARMM-NN, are constructed for proteins by training NN models on 27 fragments partitioned from the residue-based systematic molecular fragmentation (rSMF) method. The NN for each fragment is based on atom types and uses new input features that are similar to MM inputs, including bonds, angles, dihedrals, and non-bonded terms, which enhance the compatibility of CHARMM-NN to MM MD and enable the implementation of CHARMM-NN force fields in different MD programs. While the main part of the energy of the protein is based on rSMF and NN, the nonbonded interactions between the fragments and with water are taken from the CHARMM force field through mechanical embedding. The validations of the method for dipeptides on geometric data, relative potential energies, and structural reorganization energies demonstrate that the CHARMM-NN local minima on the potential energy surface are very accurate approximations to QM, showing the success of CHARMM-NN for bonded interactions. However, the MD simulations on peptides and proteins indicate that more accurate methods to represent protein-water interactions in fragments and non-bonded interactions between fragments should be considered in the future improvement of CHARMM-NN, which can increase the accuracy of approximation beyond the current mechanical embedding QM/MM level.

Comparative Study of Molecular Mechanics Force Fields for β-Peptidic Foldamers: Folding and Self-Association.

Computer-assisted study and design of non-natural peptidomimetics is increasingly important in the development of novel constructs with widespread applicability. Among these methods, molecular dynamics can accurately describe monomeric as well as oligomeric states of these compounds. We studied seven different sequences composed of cyclic and acyclic β-amino acids, the closest homologues of natural peptides, and compared the performance on them of three force field families in which specific modifications were made to improve reproduction of β-peptide structures. Altogether 17 systems were simulated, each for 500 ns, testing multiple starting conformations and in three cases also oligomer formation and stability from eight β-peptide monomers. The results indicated that our recently developed CHARMM force field extension, based on torsional energy path matching of the β-peptide backbone against quantum-chemical calculations, performs best overall, reproducing the experimental structures accurately in all monomeric simulations and correctly describing all the oligomeric examples. The Amber and GROMOS force fields could only treat some of the seven peptides (four in each case) without further parametrization. Amber was able to reproduce the experimental secondary structure of those β-peptides which contained cyclic β-amino acids, while the GROMOS force field had the lowest performance in this sense. From the latter two, Amber was able to hold together already formed associates in the prepared state but was not able to yield spontaneous oligomer formation in the simulations.

Multiscale modeling and simulation of surface‐enhanced spectroscopy and plasmonic photocatalysis

AbstractPlasmonic metal nanoparticles (PMNPs) are capable of localized surface plasmon resonance (LSPR) and have become an important component in many experimental settings, such as the surface‐enhanced spectroscopy and plasmonic photocatalysts, in which PMNPs are used to regulate the nearby molecular photophysical and photochemical behaviors by means of the complex interplay between the plasmon and molecular quantum transitions. Building computational models of these coupled plasmon‐molecule systems can help us better understand the bound molecular properties and reactivity, and make better decisions to design and control such systems. Ab initio modeling the nanosystem remains highly challenging. Many hybrid quantum‐classical (or ‐quantum) computing models have thus been developed to model the coupled systems, in which the molecular system of interest is designated as the quantum mechanical (QM) sub‐region and treated by the excited‐state electronic structure approaches such as the time‐dependent density functional theory (TDDFT), while the electromagnetic response of PMNPs is usually described using either a computational/classical electrodynamic (CED) model, polarizable continuum model(PCM), a polarizable molecular mechanics (MM) force field, or a collective of optical oscillators in QED model, leading to many hybrid approaches, such as QM/CED, QM/PCM, QM/MM or ab initio QED. In this review, we summarize recent advances in the development of these hybrid models as well as their advantages and limitations, with a specific emphasis on the TDDFT‐based approaches. Some numerical simulations on the plasmon‐enhanced absorption and Raman spectroscopy, plasmon‐driven water splitting reaction and interfacial electronic injection dynamics in dye‐sensitized solar cell are demonstrated.This article is categorized under: Electronic Structure Theory > Ab Initio Electronic Structure Methods Theoretical and Physical Chemistry > Spectroscopy Software > Quantum Chemistry Electronic Structure Theory > Combined QM/MM Methods

GAFF-AIC: reoptimisation of the GAFF force field for realistic densities and viscosities in aromatic isocyanates

ABSTRACT Molecular mechanical force field methods, which help to investigate material properties beyond molecular scales, describe the intra- and intermolecular interactions of molecules assuming a classical force field comprised of several terms. Force fields should accurately describe the energy interactions of all moleculer elements to, eventually, accurately reproduce macromolecular properties of that materials. In the case of aromatic isocyanates, general purpose force fields failed to properly account for molecular interactions and, consequently, at the prediction of macromolecular properties of isocyanate-based materials. The present study presents a new force field based on the General AMBER Force Field (GAFF), namely GAFF-AIC, which shows improved results for the prediction of macromolecular properties, in particular density and viscosities, of materials based on aromatic isocyanates. For that, a new set of parameters for intramolecular and intermolecular interactions has been optimised by means of a least square fitting. It is shown that GAFF-AIC is able to accurately reproduce the experimental density and viscosity of liquid 1-naphtyl isocyanate. Moreover, these newly proposed force field parameters were evaluated to see how efficient they were at predicting the liquid densities and viscosities for several considered aromatic isocyanates, resulting in an accurate estimation that is in excellent agreement with experimental results.

On hidden anisotropy of formally charged fragments.

The proper and precise reproduction of the molecular electrostatic potential (MEP) is crucial to describe correctly electrostatic interactions in molecular modeling. Most of the classical molecular mechanics force fields for biomolecules and drug-like molecules use the atom-centered point charges to describe MEP. However, it has been systematically pointed out in literature that such an approximation is not always enough, and some groups, like amino group or heavy halogens, require the use of anisotropic model for better description of their MEP. At the same time, the formally charged groups have not been as extensively and systematically studied as their neutral counterparts. In this report, we demonstrate that the anisotropic models for formally charged groups do bring improvements in the reference MEP reproduction, that are comparable in magnitude to those for neutral groups.