Set Of Atomic Charges Research Articles

SAMPL9 blind predictions for toluene/water partition coefficients using nonequilibrium alchemical approaches.

We present our blind prediction of the toluene-water partition coefficients in the context of the SAMPL9 challenge. For the calculation of the solvation free energies in water, toluene, and 1-octanol, we used an efficient MD-based nonequilibrium alchemical technique relying on the GAFF2 non-polarizable force field. The method is based on the fast-growth of an initially decoupled solute. Canonical sampling of the associated end-state is efficiently obtained by performing a Hamiltonian replica exchange simulation of the gas-phase solute molecule alone, combined with equilibrium configurations of the solvent. Before submitting the prediction, a pre-assessment of the method and of the force field was made by comparing with the known experimental counterpart the calculated octanol-water partition coefficients using different set of atomic charges. The analysis allowed to optimize our blind prediction for the toluene-water partition coefficients, providing at the same time valid clues for improving the performance and reliability of the non-polarizable force field in free energy calculations of drug-receptor systems.

Revisiting OPLS-AA Force Field for the Simulation of Anionic Surfactants in Concentrated Electrolyte Solutions.

Hereby, we developed a set of nonbonded parameters within all-atom optimized potentials for liquid simulations (OPLS-AA) force field for the simulation of concentrated electrolyte solutions of anionic surfactants. More specifically, the aim of this paper is to assess the performance of five sets of atomic charges calculated using different population analyses (DDEC6, CHelpG, CHelpG-SMD, RESP, and CM5), as well as the original set of charges used in the literature for sodium dodecyl sulfate (SDS) simulation. Recently, Farafonov et al. have revised the SDS OPLS-AA force field; however, we were unable to obtain the experimental rodlike micelles using this parameter set on long time scale. In fact, the initial SDS bilayer micelle adopted a rodlike shape transiently and then broke down into spherical micelles. Updating OPLS-AA force field with DDEC6, CHelpG, and CHelpG-SMD charges resulted in stable rod micelles for a long simulation time (1 μs). The atomic charges of Farafonov (taken from Shelley et al.), RESP, and CM5 could not correctly describe SDS in concentrated electrolyte solutions. Analysis of the interaction of SDS with the counterions and solvent highlights the role of a balance of the intermolecular forces that must be met to describe adequately the anionic surfactant electrolyte solutions. Further, the optimization of the SDS Lennard-Jones parameters enabled the Farafonov set to properly reproduce the experimental rod micelle. In addition, we have examined the performance of different parameters of sodium ions: the first developed based on the Kirkwood-Buff integrals (KBI) and the second developed by Joung et al. The excessive ion pairing caused by KBI parameters screens significantly SDS-water interactions, which stabilize the rod micelle. Further, a tight interaction of the Na+-SDS head group resulted in stabilization of the bilayer micelle as observed in the case of Na+ parameters developed by Joung et al.

Effect of Atomic Charges on Octanol-Water Partition Coefficient Using Alchemical Free Energy Calculation.

The octanol–water partition coefficient (logPow) is an important index for measuring solubility, membrane permeability, and bioavailability in the drug discovery field. In this paper, the logPow values of 58 compounds were predicted by alchemical free energy calculation using molecular dynamics simulation. In free energy calculations, the atomic charges of the compounds are always fixed. However, they must be recalculated for each solvent. Therefore, three different sets of atomic charges were tested using quantum chemical calculations, taking into account vacuum, octanol, and water environments. The calculated atomic charges in the different environments do not necessarily influence the correlation between calculated and experimentally measured ∆Gwater values. The largest correlation coefficient values of the solvation free energy in water and octanol were 0.93 and 0.90, respectively. On the other hand, the correlation coefficient of logPow values calculated from free energies, the largest of which was 0.92, was sensitive to the combination of the solvation free energies calculated from the calculated atomic charges. These results reveal that the solvent assumed in the atomic charge calculation is an important factor determining the accuracy of predicted logPow values.

Atomic charges for conformationally rich molecules obtained through a modified principal component regression.

A modification of the principal component regression model is proposed for obtaining a fixed set of atomic charges (referred to as dipole-derived charges) optimized for reproducing the dipole moment of a conformationally rich molecule, i.e., a molecule with multiple local minima on the potential energy surface. The method does not require any adjustable parameters and requires the geometries of conformers, their dipole moments and atomic polar tensor (APT) charges as the only input data. The fixed atomic charges generated by the method not only reproduce the molecular dipole moment in all the conformers accurately, but are also numerically close to the APT charges, thereby ensuring accurate reproduction of the dipole moment variations caused by small geometrical distortions (e.g., by vibrations) of the conformers. The proposed method has been applied to canonical 2'-deoxyribonucleotides, the model DNA monomers, and the dipole-derived charges have been shown to outperform both the averaged APT and RESP charges in reproducing the dipole moments of large sets of conformers, thus demonstrating a potential usefulness of the dipole-derived charges as a 'reference point' for modeling polarization effects in conformationally rich molecules.

LigParGen web server: an automatic OPLS-AA parameter generator for organic ligands.

The accurate calculation of protein/nucleic acid–ligand interactions or condensed phase properties by force field-based methods require a precise description of the energetics of intermolecular interactions. Despite the progress made in force fields, small molecule parameterization remains an open problem due to the magnitude of the chemical space; the most critical issue is the estimation of a balanced set of atomic charges with the ability to reproduce experimental properties. The LigParGen web server provides an intuitive interface for generating OPLS-AA/1.14*CM1A(-LBCC) force field parameters for organic ligands, in the formats of commonly used molecular dynamics and Monte Carlo simulation packages. This server has high value for researchers interested in studying any phenomena based on intermolecular interactions with ligands via molecular mechanics simulations. It is free and open to all at jorgensenresearch.com/ligpargen, and has no login requirements.

Implicit Solvent Model for Million-Atom Atomistic Simulations: Insights into the Organization of 30-nm Chromatin Fiber

Molecular dynamics (MD) simulations based on the implicit solvent generalized Born (GB) models can provide significant computational advantages over the traditional explicit solvent simulations. However, the standard GB becomes prohibitively expensive for all-atom simulations of large structures; the model scales poorly, ∼n2, with the number of solute atoms. Here we combine our recently developed optimal point charge approximation (OPCA) with the hierarchical charge partitioning (HCP) approximation to present an ∼n log n multiscale, yet fully atomistic, GB model (GB-HCPO). The HCP approximation exploits the natural organization of biomolecules (atoms, groups, chains, and complexes) to partition the structure into multiple hierarchical levels of components. OPCA approximates the charge distribution for each of these components by a small number of point charges so that the low order multipole moments of these components are optimally reproduced. The approximate charges are then used for computing electrostatic interactions with distant components, while the full set of atomic charges are used for nearby components. We show that GB-HCPO can deliver up to 2 orders of magnitude speedup compared to the standard GB, with minimal impact on its accuracy. For large structures, GB-HCPO can approach the same nominal speed, as in nanoseconds per day, as the highly optimized explicit-solvent simulation based on particle mesh Ewald (PME). The increase in the nominal simulation speed, relative to the standard GB, coupled with substantially faster sampling of conformational space, relative to the explicit solvent, makes GB-HCPO a suitable candidate for MD simulation of large atomistic systems in implicit solvent. As a practical demonstration, we use GB-HCPO simulation to refine a ∼1.16 million atom structure of 30 nm chromatin fiber (40 nucleosomes). The refined structure suggests important details about spatial organization of the linker DNA and the histone tails in the fiber: (1) the linker DNA fills the core region, allowing the H3 histone tails to interact with the linker DNA, which is consistent with experiment; (2) H3 and H4 tails are found mostly in the core of the structure, closer to the helical axis of the fiber, while H2A and H2B are mostly solvent exposed. Potential functional consequences of these findings are discussed. GB-HCPO is implemented in the open source MD software NAB in Amber 2016.

Searching the Force Field Electrostatic Multipole Parameter Space.

We show by tensor decomposition analyses that the molecular electrostatic potential for amino acid peptide models has an effective rank less than twice the number of atoms. This rank indicates the number of parameters that can be derived from the electrostatic potential in a statistically significant way. Using this as a guideline, we investigate different strategies for deriving a reduced set of atomic charges, dipoles, and quadrupoles capable of reproducing the reference electrostatic potential with a low error. A full combinatorial search of selected parameter subspaces for N-methylacetamide and a cysteine peptide model indicates that there are many different parameter sets capable of providing errors close to that of the global minimum. Among the different reduced multipole parameter sets that have low errors, there is consensus that atoms involved in π-bonding require higher order multipole moments. The possible correlation between multipole parameters is investigated by exhaustive searches of combinations of up to four parameters distributed in all possible ways on all possible atomic sites. These analyses show that there is no advantage in considering combinations of multipoles compared to a simple approach where the importance of each multipole moment is evaluated sequentially. When combined with possible weighting factors related to the computational efficiency of each type of multipole moment, this may provide a systematic strategy for determining a computational efficient representation of the electrostatic component in force field calculations.

Activity prediction of substrates in NADH-dependent carbonyl reductase by docking requires catalytic constraints and charge parameterization of catalytic zinc environment.

Molecular docking of substrates is more challenging compared to inhibitors as the reaction mechanism has to be considered. This becomes more pronounced for zinc-dependent enzymes since the coordination state of the catalytic zinc ion is of greater importance. In order to develop a predictive substrate docking protocol, we have performed molecular docking studies of diketone substrates using the catalytic state of carbonyl reductase 2 from Candida parapsilosis (CPCR2). Different docking protocols using two docking methods (AutoDock Vina and AutoDock4.2) with two different sets of atomic charges (AM1-BCC and HF-RESP) for catalytic zinc environment and substrates as well as two sets of vdW parameters for zinc ion were examined. We have selected the catalytic binding pose of each substrate by applying mechanism based distance criteria. To compare the performance of the docking protocols, the correlation plots for the binding energies of these catalytic poses were obtained against experimental Vmax values of the 11 diketone substrates for CPCR2. The best correlation of 0.73 was achieved with AutoDock4.2 while treating catalytic zinc ion in optimized non-bonded (NBopt) state with +1.01 charge on the zinc ion, compared to 0.36 in non-bonded (+2.00 charge on the zinc ion) state. These results indicate the importance of catalytic constraints and charge parameterization of catalytic zinc environment for the prediction of substrate activity in zinc-dependent enzymes by molecular docking. The developed predictive docking protocol described here is in principle generally applicable for the efficient in silico substrate spectra characterization of zinc-dependent ADH.

Charge density distributions derived from smoothed electrostatic potential functions: design of protein reduced point charge models

To generate reduced point charge models of proteins, we developed an original approach to hierarchically locate extrema in charge density distribution functions built from the Poisson equation applied to smoothed molecular electrostatic potential (MEP) functions. A charge fitting program was used to assign charge values to the so-obtained reduced representations. In continuation to a previous work, the Amber99 force field was selected. To easily generate reduced point charge models for protein structures, a library of amino acid templates was designed. Applications to four small peptides, a set of 53 protein structures, and four KcsA ion channel models, are presented. Electrostatic potential and solvation free energy values generated by the reduced models are compared with the corresponding values obtained using the original set of atomic charges. Results are in closer agreement with the original all-atom electrostatic properties than those obtained with a previous reduced model that was directly built from the smoothed MEP functions [Leherte and Vercauteren in J Chem Theory Comput 5:3279-3298, 2009].

An n log n Generalized Born Approximation.

Molecular dynamics (MD) simulations based on the generalized Born (GB) model of implicit solvation offer a number of important advantages over the traditional explicit solvent based simulations. Yet, in MD simulations, the GB model has not been able to reach its full potential partly due to its computational cost, which scales as ∼n(2), where n is the number of solute atoms. We present here an ∼n log n approximation for the generalized Born (GB) implicit solvent model. The approximation is based on the hierarchical charge partitioning (HCP) method (Anandakrishnan and Onufriev J. Comput. Chem. 2010 , 31 , 691 - 706 ) previously developed and tested for electrostatic computations in gas-phase and distant dependent dielectric models. The HCP uses the natural organization of biomolecular structures to partition the structures into multiple hierarchical levels of components. The charge distribution for each of these components is approximated by a much smaller number of charges. The approximate charges are then used for computing electrostatic interactions with distant components, while the full set of atomic charges are used for nearby components. To apply the HCP concept to the GB model, we define the equivalent of the effective Born radius for components. The component effective Born radius is then used in GB computations for points that are distant from the component. This HCP approximation for GB (HCP-GB) is implemented in the open source MD software, NAB in AmberTools, and tested on a set of representative biomolecular structures ranging in size from 632 atoms to ∼3 million atoms. For this set of test structures, the HCP-GB method is 1.1-390 times faster than the GB computation without additional approximations (the reference GB computation), depending on the size of the structure. Similar to the spherical cutoff method with GB (cutoff-GB), which also scales as ∼n log n, the HCP-GB is relatively simple. However, for the structures considered here, we show that the HCP-GB method is more accurate than the cutoff-GB method as measured by relative RMS error in electrostatic force compared to the reference (no cutoff) GB computation. MD simulations of four biomolecular structures on 50 ns time scales show that the backbone RMS deviation for the HCP-GB method is in reasonable agreement with the reference GB simulation. A critical difference between the cutoff-GB and HCP-GB methods is that the cutoff-GB method completely ignores interactions due to atoms beyond the cutoff distance, whereas the HCP-GB method uses an approximation for interactions due to distant atoms. Our testing suggests that completely ignoring distant interactions, as the cutoff-GB does, can lead to qualitatively incorrect results. In general, we found that the HCP-GB method reproduces key characteristics of dynamics, such as residue fluctuation, χ1/χ2 flips, and DNA flexibility, more accurately than the cutoff-GB method. As a practical demonstration, the HCP-GB simulation of a 348 000 atom chromatin fiber was used to refine the starting structure. Our findings suggest that the HCP-GB method is preferable to the cutoff-GB method for molecular dynamics based on pairwise implicit solvent GB models.

Molecular Modeling of Crystalline Alkylthiophene Oligomers and Polymers

We present the results of a thorough molecular modeling study of several alkylthiophene-based oligomers and polymers. In particular, we consider two polymers whose limit-ordered crystal structures have been recently reported by our group, on the basis of powder X-ray data analysis: poly(3-(S)-2-methylbutylthiophene) (P3MBT) and form I' of poly(3-butylthiophene) (P3BT). We first describe the development of a series general purpose force fields for the simulation of these and related systems. The force fields incorporate the results of ab initio calculations of the bond torsion energies of selected oligomers and differ in the set of atomic charges used to represent the electrostatic interactions. We then present the results of an extensive validation of these force fields, by means of molecular mechanics (MM) energy minimizations and molecular dynamics (MD) simulations of the crystal structures of these oligomers and polymers. While our "best" force field does not outperform the others on each of the investigated systems, it provides a balanced description of their overall structure and energetics. Finally, our MM minimizations and MD simulations confirm that the reported crystal structures of P3MBT and P3BT are stable and correspond to well-defined energetic minima. The room-temperature MD simulations reveal a certain degree of side-chain disorder, even in our virtually defect-free polymer crystal models.

Coarse Point Charge Models For Proteins From Smoothed Molecular Electrostatic Potentials

To generate coarse electrostatic models of proteins, we developed an original approach to hierarchically locate maxima and minima in smoothed molecular electrostatic potentials. A charge-fitting program was used to assign charges to the so-obtained reduced representations. Templates are defined to easily generate coarse point charge models for protein structures, in the particular cases of the Amber99 and Gromos43A1 force fields. Applications to four small peptides and to the ion channel KcsA are presented. Electrostatic potential values generated by the reduced models are compared with the corresponding values obtained using the original sets of atomic charges.

Intra-protein hydrogen bonding is dynamically stabilized by electronic polarization

Molecular dynamics (MD) simulation has been carried out to study dynamical stability of intra-protein hydrogen bonds based on two set of atomic charges, the standard AMBER charge and the polarized protein-specific charge (PPC). The latter is derived from quantum mechanical calculation for protein in solution using a recently developed molecular fractionation with conjugate caps-Poisson-Boltzmann (MFCC-PB) approach and therefore includes electronic polarization effect of the protein at native structure. MD simulations are performed for a number of benchmark proteins containing helix and/or beta sheet secondary structures. The computational result shows that occupancy percentage of hydrogen bonds averaged over simulation time, as well as the number of hydrogen bonds as a function of simulation time, is consistently higher under PPC than AMBER charge. In particular, some intra-protein hydrogen bonds are found broken during MD simulation using AMBER charge but they are stable using PPC. The breaking of some intra-protein hydrogen bonds in AMBER simulation is responsible for deformation or denaturing of some local structures of proteins during MD simulation. The current study provides strong evidence that hydrogen bonding is dynamically more stable using PPC than AMBER charge, highlighting the stabilizing effect of electronic polarization on protein structure.

Continuum treatment of electronic polarization effect

A continuum treatment of electronic polarization has been explored for in molecular mechanics simulations in implicit solvents. The dielectric constant for molecule interior is the only parameter in the continuum polarizable model. A value of 4 is found to yield optimal agreement with high-level ab initio quantum mechanical calculations for the tested molecular systems. Interestingly, its performance is not sensitive to the definition of molecular volume, in which the continuum electronic polarization is defined. In this model, quantum mechanical electrostatic field in different dielectric environments from vacuum, low-dielectric organic solvent, and water can be used simultaneously in atomic charge fitting to achieve consistent treatment of electrostatic interactions. The tests show that a single set of atomic charges can be used consistently in different dielectric environments and different molecular conformations, and the atomic charges transfer well from training monomers to tested dimers. The preliminary study gives us the hope of developing a continuum polarizable force field for more consistent simulations of proteins and nucleic acids in implicit solvents.

Complementary Use of Simulations and Molecular-Thermodynamic Theory to Model Micellization

Molecular-thermodynamic descriptions of micellization in aqueous media can be utilized to model the self-assembly of surfactants possessing relatively simple chemical structures, where it is possible to identify a priori what equilibrium position they will adopt in the resulting micellar aggregate. For such chemical structures, the portion of the surfactant molecule that is expected to be exposed to water upon aggregate self-assembly can be identified and used as an input to the molecular-thermodynamic description. Unfortunately, for many surfactants possessing more complex chemical structures, it is not clear a priori how they will orient themselves within a micellar aggregate. In this paper, we present a computational approach to identify what portions of a surfactant molecule are hydrated in a micellar environment through the use of molecular dynamics simulations of such molecules at an oil/water interface (modeling the micelle core/water interface). The local environment of each surfactant segment is determined by counting the number of contacts of each segment with the water and oil molecules. After identifying the hydrated and the unhydrated segments of the surfactant molecule, molecular-thermodynamic modeling can be performed to predict: (i) the free-energy change associated with forming a micellar aggregate, (ii) the critical micelle concentration (CMC), and (iii) the optimal shape and size of the micellar aggregate. The computer simulation results were found to be sensitive to the atomic charge parameters utilized during the simulation runs. Two different methods of assigning atomic charges were tested, and the computer simulation and molecular-thermodynamic modeling results obtained using both sets of atomic charges are presented and compared. The combined computer simulation/molecular-thermodynamic modeling approach presented here is validated first by implementing it in the case of anionic (sodium dodecyl sulfate, SDS), cationic (cetyltrimethylammonium bromide, CTAB), zwitterionic (dodecylphosphocholine, DPC), and nonionic (dodecyl poly(ethylene oxide), C12E8) surfactants possessing relatively simple chemical structures and verifying that good predictions of CMCs and micelle aggregation numbers are obtained. In the case of C12E8, the challenges and limitations associated with simulating a single, polymeric E8 moiety at the oil/water interface to model its behavior at the micelle/water interface are discussed. Subsequently, the combined modeling approach is implemented in the case of the anionic surfactant 3-hydroxy sulfonate (AOS) and of the nonionic surfactant decanoyl-n-methylglucamide (MEGA-10), which possess significantly more complex chemical structures. The good predictions obtained for these two surfactants indicate that the combined computer simulation/molecular-thermodynamic modeling approach presented here extends the range of applicability of molecular-thermodynamic theory to allow modeling of the micellization behavior of surfactants possessing more complex chemical structures.

Semiempirical Calculations on the Dipole Moment Enhancement in the Solid State

The enhancement of the molecular dipole moment under the influence of the crystal lattice relative to the value for the corresponding isolated molecule has been calculated for the three low-pressure polymorphs of glycine using the semiempirical MINDO/3 Hamiltonian. In these calculations the influence of the lattice on an arbitrarily chosen reference molecule has been included by modifying the diagonal elements of the Hcore matrix. Starting from a set of atomic charges calculated for the isolated molecule an initial crystal lattice had been set up and included into the calculation of new self-consistent (SCF) charges for the reference molecule. The resulting set of new SCF charges was then used to replace the lattice charges from the preceding run, and a new SCF calculation was performed. This iterative process was continued until the difference between the electronic energies and the sum of the absolute values of the differences between the atomic charges from two subsequent SCF calculations was lower than 10−5 eV and 10−5 e0, respectively. The results agree nicely with the experimental data obtained by other authors for α-glycine.

Collision-induced absorption in liquid carbon tetrachloride

Molecular dynamics (MD) simulation has been used for the computation of the far infrared absorption coefficient of CCl 4, with an excellent accord with experimental results at all frequencies. Nonadditive polarization interactions have been implemented with the chemical potential equalization method, which requires the gas phase values of the octupole moment and mean polarizability tensor. It is shown that the absorption coefficient constitutes an exacting probe of interaction potentials, strongly favoring a physically reasonable set of atomic charges.

The energetic and structural effects of steric crowding in phosphate and dithiophosphinate complexes of lanthanide cations M3+: a computational study.

Metal-ligand binding strength and selectivity result from antagonistic metal-ligand M-L attractions and ligand-ligand L-L repulsions. On the basis of quantum-mechanical (QM) calculations on lanthanide complexes, we show that this interplay determines the binding affinities in the gas phase. In the series of [ML3] complexes (M = La, Eu, and Yb) with negatively charged phosphoryl ligands L- = (MeO)2PO2- and Me2PS2-, the binding energies follow the order Yb3+ > Eu3+ > La3- for a given ligand, and (MeO)2PO2- > Me2PS2- for a given cation. However, adding a neutral LH ligand to [ML3] changes the order to Eu3+ > Yb3+ > La3+ for the oxygen ligand and La3+ > Eu3- > Yb3+ for the sulfur ligand, indicating that steric strain in the first coordination sphere is largest for the smallest cation and for sulfur binding sites. We investigated the question of additional hydration of the [ML3LH] complexes in aqueous solution by molecular dynamics (MD) simulations, using two sets of atomic charges. It was found that pairwise additive potentials overestimate the coordination and hydration numbers of the cations, while adding polarization energy terms for the ligands yields better agreement between QM and MD results and supports the concept of steric strain in the first coordination sphere.

Derivation of an Optimized Potential Model for Phase Equilibria (OPPE) for Sulfides and Thiols

An extension of the anisotropic united atoms model is proposed for thiols and sulfides. A complete derivation of the nonbonded parameters is performed with the aim of obtaining a transferable force field. The electrostatic part of the intermolecular potential is represented by a set of atomic charges. These charges are extracted from quantum chemical calculations thanks to a method developed recently. These charges are shown to depend very weakly on the conformation of the molecule. The repulsion−dispersion interactions are described by an anisotropic united atom Lennard-Jones potential: parameters for methyl and methylene groups are directly taken from a previous study on alkanes, whereas parameters for −S and −SH groups are fitted to experimental data. The resulting potential, which we will term OPPE (optimized potential model for phase equilibria), is tested against liquid properties at atmospheric pressure and vapor−liquid-phase equilibria of various sulfides and thiols to prove its transferability. Excellent agreement is obtained with experimental data.

X-ray Crystal Structures and Molecular Modelling Studies of Calix[4]dibenzocrowns-6 and Their Alkali Metal Cation Complexes

We present X-ray structures of complexes of KClO4·H2O and CsPic·H2O with a new calix[4]arene bearing a dibenzocrown-6 ether bridge. Alkali metal cation complexes with this macrocycle have been studied by molecular dynamics simulations in vacuo, taking account of the counterion, and in an explicit water phase, in order to compare the computed structures to those determined by X-ray crystallography. Two different sets of atomic charges were used in order to evaluate the influence of the electrostatic representation on the computed structures. The results obtained show that MNDO scaled charges give simulation results in better agreement with the available experimental data than MNDO/ESP.