Charge Equilibration Research Articles

Development and Application of Non-Additive Force Fields for Molecular Simulations of Lipid Bilayers and Integral Membrane Proteins

Molecular simulations today are applied across many scientific disciplines. Complementing experiment, these tools afford a molecular-level understanding and interpretation of physico-chemical processes at resolutions and timescales difficult or practically inaccessible to experiment. At the heart of such methods is the description of interactions between atoms and molecules, the force field. Traditionally, non-reactive force fields have treated electrostatic interactions using an additive, Coulomb model between fixed partial charges on atomic sites. Though quite successful, there has been conjecture as to the effects of incorporating non-additivity in classical force fields, particularly in biological systems. Over the last several decades, attempts to incorporate electrostatic non-additivity in the form of inducible dipole interactions or dynamically varying partial charges have provided a vast body of knowledge that has aided in the development of a new class of force fields attempting to explicitly account for non-additive effects. We will present our recent work in developing one such class of models, charge equilibration force fields, and applications of such models to aqueous solution interfaces, membrane bilayers and simple integral membrane peptides such as the gramicidin A bacterial channel, and recent work on modeling of protein-ligand interaction free energetics.

The Extended Universal Force Field (XUFF):Theory and Applications

We present a new force field derived from a modification of the universal force field (UFF). Although the UFF is widely used in the field of molecular mechanics and dynamics, it is unsuitable for such applications as correctly estimating the hydrogen bond. We have derived a new variable-charge force field from the UFF by incorporating a modified charge equilibration (MQEq) method for the electrostatic term and a scale factor for the van der Waals term, which we refer to as the extended universal force field (XUFF). The MQEq parameters and the scale factor are optimized. We have implemented the XUFF in the ABINIT-MP program and show that the XUFF is very useful for analyzing the conformation of a molecular system such as hydrogen-bonding molecular pairs, ligands, and proteins.

Merging bond-order potentials with charge equilibration

A method is presented for extending any bond-order potential (BOP) to include charge transfer between atoms through a modification of the split-charge equilibration (SQE) formalism. Variable limits on the maximum allowed charge transfer between atomic pairs are defined by mapping bond order to an amount of shared charge in each bond. Charge transfer is interpreted as an asymmetry in how the shared charge is distributed between the atoms of the bond. Charge equilibration (QE) assesses the asymmetry of the shared charge, while the BOP converts this asymmetry to the actual amount of charge transferred. When applied to large molecules, this BOP/SQE method does not exhibit the unrealistic growth of charges that is often associated with QE models.

Surface and interstitial transition barriers in rutile (110) surface growth

We present calculated surface and interstitial transition barriers for Ti, O, O-2, TiO, and TiO2 atoms and clusters at the rutile (110) surface. Defect structures involving these small clusters, including adcluster and interstitial binding sites, were calculated by energy minimization using density-functional theory (DFT). Transition energies between these defect sites were calculated using the NEB method. Additionally, a modified SMB-Q charge equilibration empirical potential and a fixed-charge empirical potential were used for a comparison of the transition energy barriers. Barriers of 1.2-3.5 eV were found for all studied small cluster transitions upon the surface except for transitions involving O-2. By contrast, the O-2 diffusion barriers along the [001] direction upon the surface are only 0.13 eV. The QEq charge equilibration model gave mixed agreement with the DFT calculations, with the barriers ranging between 0.8 and 5.8 eV.

Revised charge equilibration parameters for more accurate hydration free energies of alkanes

We present a refined alkane charge equilibration (CHEQ) force field, improving our previously reported CHEQ alkane force field [1] to better reproduce experimental hydration free energies. Experimental hydration free energies of ethane, propane, butane, pentane, hexane, and heptane are reproduced to within 3.6% on average. We demonstrate that explicit polarization results in a shift in molecular dipole moment for water molecules associated with the alkane molecule. We also show that our new parameters do not have a significant effect on the alkane–water interactions as measured by the radial distribution function (RDF).

Classical simulation of deposition of thiophene oligomers on TiO 2-anatase: Relevance of long-range electrostatic interactions

We performed classical molecular dynamics simulations of the vapor-deposition of α -T4 oligomers on the TiO 2-anatase (1 0 1) surface, comparing different sets of charges associated with the atoms of the model. The potential energy surfaces for α -T4 and TiO 2 were described by re-parametrizations of the Universal force field with charges given by the charge equilibration (QEq) scheme, or with fixed charges obtained by an ab initio method using the Hirshfeld partition. The two sets of charges lead to completely different results for the interface formation, and for the characteristics of the organic film, with a clearly defined α -T4 contact layer in the QEq case, and a more homogeneous molecular distribution when using Hirshfeld charges. The main reason for the discrepancy was found to be the incorrect charge assignment given by QEq to the sulfur and α -carbon atoms in thiophenes, and highlight the relevance of long-range interactions in the organization of molecular films.

Erratum to: Synthesis of exotic nuclei in heavy-ion collisions at higher energies

The charge equilibration upper energy-limit formula is proposed. It shows not only the basic mechanism of charge equilibration but also does imply the possibility of controlling the exotic nuclear synthesis. To understand the mechanism of charge equilibration, charge and mass distributions of final fragments of the heavy-ion collisions are studied based on three-dimensional time-dependent Hartree-Fock calculations with a Skyrme interaction.

The electronegativity equalization method and the split charge equilibration applied to organic systems: Parametrization, validation, and comparison

An extensive benchmark of the electronegativity equalization method (EEM) and the split charge equilibration (SQE) model on a very diverse set of organic molecules is presented. These models efficiently compute atomic partial charges and are used in the development of polarizable force fields. The predicted partial charges that depend on empirical parameters are calibrated to reproduce results from quantum mechanical calculations. Recently, SQE is presented as an extension of the EEM to obtain the correct size dependence of the molecular polarizability. In this work, 12 parametrization protocols are applied to each model and the optimal parameters are benchmarked systematically. The training data for the empirical parameters comprise of MP2/Aug-CC-pVDZ calculations on 500 organic molecules containing the elements H, C, N, O, F, S, Cl, and Br. These molecules have been selected by an ingenious and autonomous protocol from an initial set of almost 500,000 small organic molecules. It is clear that the SQE model outperforms the EEM in all benchmark assessments. When using Hirshfeld-I charges for the calibration, the SQE model optimally reproduces the molecular electrostatic potential from the ab initio calculations. Applications on chain molecules, i.e., alkanes, alkenes, and alpha alanine helices, confirm that the EEM gives rise to a divergent behavior for the polarizability, while the SQE model shows the correct trends. We conclude that the SQE model is an essential component of a polarizable force field, showing several advantages over the original EEM.

Interfacial structure, thermodynamics, and electrostatics of aqueous methanol solutions via molecular dynamics simulations using charge equilibration models.

We present results from molecular dynamics simulations of methanol-water solutions using charge equilibration force fields to explicitly account for nonadditive electronic interaction contributions to the potential energy. We study solutions across the concentration range from 0.1 to 0.9 methanol mole fraction. At dilute concentrations, methanol density is enhanced at the liquid-vapor interface, consistent with previous molecular dynamics and experimental studies. Interfacial thickness exhibits a monotonic increase with increasing methanol mole fraction, while surface tensions display monotonic decrease with methanol concentration, in qualitative agreement with experimental data and previous molecular dynamics predictions using polarizable force fields. In terms of interfacial structure, in keeping with predictions of traditional force fields, there is a unique preferential orientation of methanol molecules at the interface. Moreover, there is a free energetic preference for methanol molecules at the interface as evidenced by potential of mean force calculations. The pmf calculations suggest an interfacial state with 0.8 kcal/mol stability relative to the bulk, again in qualitative agreement with previous simulation and experimental studies. Interfacial potentials based on double integration of total charge density range from -610 to -330 mV over the dilute to concentrated regimes, respectively. The preponderance of methanol at the interface at all mole fractions gives rise to a dominant methanol contribution to the total interfacial potential. Interestingly, there is a transition of the water surface potential contribution from negative to positive upon the transition from methanol mole fraction of 0.1 to 0.2. The dipole and quadrupole contributions to the water component of the total interfacial potential are effectively of equal magnitude and opposite sign, thus cancelling one another. We compute the in-plane component of the dielectric permittivity along the interface normal. We observe a nonmonotonic behavior of the methanol in-plane dielectric permittivity that tracks the methanol density profiles at low methanol mole fractions. At higher methanol mole fractions, the total in-plane permittivity is dominated by methanol and displays a monotonic decrease from bulk to vapor. We finally probe the nature of hydration of water in the bulk versus interfacial regions for methanol mole fractions of 0.1 and 0.2. In the bulk, methanol perturbs water structure so as to give rise to water hydrogen bond excesses. Moreover, we observe negative hydrogen bond excess in the vicinity of the alkyl group, as reported by Zhong et al. for bulk ethanol-water solutions using charge equilibration force fields, and positive excess in regions hydrogen bonding to nearest-neighbor methanol molecules. Within the interfacial region, water and methanol density reduction lead to concomitant water hydrogen bond deficiencies (negative hydrogen-bond excess).

Charge equilibration force fields for lipid environments: applications to fully hydrated DPPC bilayers and DMPC-embedded gramicidin A.

Polarizable force fields for lipid and solvent environments are used for molecular dynamics simulations of a fully hydrated dipalmitoylphosphatidylcholine (DPPC) bilayer and gramicidin A (gA) dimer embedded in a dimyristoylphosphatidylcholine (DMPC) bilayer. The lipid bilayer is modeled using the CHARMM charge equilibration (CHEQ) polarizable force field for lipids and the TIP4P-FQ force field to represent solvent. For the DPPC bilayer system, results are compared to the same system simulated using the nonpolarizable CHARMM27r (C27r) force field and TIP3P water. Calculated atomic and electron density profiles, head group orientations as measured by the phosphorus-nitrogen vector orientation, and deuterium order parameters are found to be consistent with previous simulations and with experiment. The CHEQ model exhibits greater water penetration into the bilayer interior, as demonstrated by the potential of mean force calculated from the water density profile. This is a result of the variation of the water molecular dipole from 2.55 D in the bulk to 1.88 D in the interior. We discuss this finding in the context of previous studies (both simulation and experiment) that have investigated the extent of penetration of water into DPPC bilayers. We also discuss the effects of including explicit polarization on the water dipole moment variation as a function of distance from the bilayer. We show distributions of atomic charges over the course of the simulation since the CHEQ model allows the charges to fluctuate. We have calculated the interfacial dipole potential, which the CHEQ model predicts to be 0.95 V compared to 0.86 V as predicted by the C27r model. We also discuss dielectric permittivity profiles and the differences arising between the two models. We obtain bulk values of 72.77 for the CHEQ model (TIP4P-FQ water) and 91.22 for C27r (TIP3P), and values approaching unity in the membrane interior. Finally, we present results of simulations of gA embedded in a DMPC bilayer using the CHEQ model and discuss structural properties.

Coulomb explosion as a probe to understand the mechanism of electron stripping from ions interacting with crystalline solids

When an ion impinges on a solid, it rapidly undergoes a process in which its electrons are stripped away provided the velocity of the orbiting electrons is smaller than the projectile speed. Electron stripping determines any posterior behavior of the ions in the solid, and it is assumed that it takes place on the surface of the solid, but no information is available on the details of the process. Here we show, using the Coulomb explosion of ${\text{C}}_{2}^{+}$ ions moving in Si as a tool, that electron stripping takes place in an orderly manner and that the number of electrons stripped, before charge equilibration, depends on a characteristic length. We also propose a relation capable of quantifying this dependence. We foresee these results as a starting point to a more general understanding of ion-solid interaction, with important consequences on ion beam analysis and modification techniques, and special significance in silicon technology.

Valence-Bond/Coherent-States Approach to the Charge Equilibration Model I. Valence-Bond Models for Diatomic Molecules

A valence-bond (VB)/coherent-states (CS) approach to the charge equilibration (CE) model for diatomic molecules is presented in this work (VB part) and its sequel (CS part). By emphasizing theoretical aspects, this approach obtains the classical-electrostatics CE model from a quantum VB model in conjunction with the CS theory. For the VB part, a VB generalized CE (VB/GCE) model, which contains the CE model as a subcase, is derived from a two-electron, three-state VB model via the sequential application of seven approximations. Unlike its CE subcase, the VB/GCE model provides a satisfactory charge-transfer description at dissociation as illustrated with HF((g)) and other molecules. Through the previous derivation, CE charges and CE Coulomb interactions are elucidated in terms of VB Mulliken charges and VB atomic interactions, respectively. Modifications in the above derivation can generate a family of related VB/GCE models that includes the aforesaid VB/GCE model. Despite their classical appearance, all of the VB/GCE and CE models involve VB wave functions corresponding to ground and first-excited states. Moreover, all of the VB/GCE and CE energy and charge optimizations are proven to be equivalent to the variational eigenvector equation procedures of the underlying VB models. The quantum-mechanics/classical-electrostatics connection implicit in this work is further elaborated by means of VB CS sets in the sequel. The VB/CS treatment of polyatomic molecules and additional tests of the present approach will be reported in later papers in this series.

Dielectric properties of solids in the regular and split-charge equilibration formalisms

We investigate the generic dielectric properties of solids in the split-charge equilibration (SQE) formalism, which contains the regular charge equilibration (QE) method as a limiting case, but augments it with a bond hardness term. It is shown that QE always mimics ideal conductors, while any positive bond hardness used in SQE turns the solid into a dielectric. Crystals with simple cubic and rocksalt structure are considered explicitly. For these symmetries, we solve the continuum limit of the SQE formalism analytically. As a result, we provide simple analytical expressions for how dielectric constant and penetration depth of the electrostatic field depend on atomic hardness, bond hardness, and lattice constant. This mapping may prove useful not only for force field parametrization but also for solving dielectric responses on coarse-grained scales. Successful comparison of numerical data to analytical solutions is made, including those containing discretization corrections.

Self-Consistent Charge Equilibration Method and Its Application to Au13Nan(n= 1,10) Clusters

We modified the charge equilibration method (QEq) of Rappe et al. by including the third and fourth order terms in the series expansion for the energy of charged atoms. This leads to a self-consistent scheme for obtaining charges and total energy in a molecule. We combined the modified QEq with a scaled morse potential (SMP) to get a new potential, SMP/QEq. The SMP/QEq method allows to do realistic geometry optimizations and calculate reliable charge distributions and ionization energies. The performance of SMP/QEq was fully tested by comparing to high level theories and experiments. We established a connection between first-principles DFT and SMP/QEq which gives insight into SMP/QEq and suggests ways to improve it.

Molecular geometry-dependent atomic charge calculation with modified charge equilibration method

We have improved a modified charge equilibration (MQEq) method for calculating the geometry-dependent distribution of atomic charges. In this paper, Ohno-Klopman, Ohno and DasGupta-Huzinaga equations are adopted to express the shielding effect, and the calculated atomic charges with these MQEq methods are in good agreement with those by the HF/6-31G(d,p) calculations for several organic molecules. These MQEq methods would be useful to estimate the charge distribution for large molecules.

THE INFLUENCE OF SUBSTITUENT ELECTRONIC EFFECT ON ETHYLENE OLIGOMERIZATION ACTIVITIES OF BIS(IMINO)PYRIDYL Fe(II) CATALYSTS: A COMBINED MOLECULAR MECHANICS AND CHARGE EQUILIBRATION METHOD

Bis(imino)pyridyl Fe(II) complexes are important catalysts in ethylene oligomerization for preparing α-olefins. The metal net charge-activity relationship of bis(imino)pyridyl Fe(II) complexes was investigated by molecular mechanics (MM) and net charge equilibration (QEq) method with modified Dreiding force field. It was found that metal net charge was in reverse ratio to ethylene oligomerization activity. Electron-donor substituents with less steric hindrance to the central metal were favorable to Fe complex activity. Metal net charge-activity relationship could be used to assist the design of new Fe oligomerization catalysts with higher activity.

Electrostatic Polarization Effects and Hydrophobic Hydration in Ethanol−Water Solutions from Molecular Dynamics Simulations

Nonadditive electrostatic force fields based on the charge equilibration formalism coupled with long time-scale molecular dynamics simulations are used to investigate the microscopic structural aspects of hydrophobic hydration in ethanol-water solutions. Employing a combination of polarizable ethanol and water force fields (developed independently), we find that solution properties are satisfactorily reproduced across the ethanol mole fraction range between 0.1 and 0.9. Solution densities are predicted within 3.6% of experimental measurements, while excess mixing enthalpies are overestimated as in earlier studies. The solvation free energy of ethanol in infinite dilution is determined via thermodynamic integration to be 5.70 +/- 0.23 kcal/mol, overestimating the free energetics of solvation relative to experiment (5.01 kcal/mol). Bulk solution dielectric constants and diffusion constants reproduce experimental trends and are in reasonable agreement across the ethanol concentration range studied. Because of explicit accounting of induction effects, ethanol and water exhibit varying molecular dipole moment distributions with concentration. The polarizable ethanol model, possessing higher condensed-phase polarizability relative to the TIP4P-FQ water model (4.54 A(3) versus 1.1 A(3), respectively), displays greater variation upon perturbation by the electric field of water. With regard to hydrophobic hydration, the current force fields indicate positive hydrogen bonding excess for water in the dilute ethanol concentration range, consistent with previous theoretical and experimental studies. Strikingly, we find that there are both positive and negative hydrogen bond excess contributions within the first hydration shell of both the ethanol hydroxyl oxygen and ethylene carbon atoms. The larger positive contributions dominate the overall hydrogen bonding patterns to yield overall net positive excesses. Moreover, we do not find evidence of excess hydrogen bonding vicinal to the nonpolar moieties as has been suggested based on the "iceberg"-like models proposed by Frank and Evans. The present results suggest negative excess in the regions surrounding the alkyl groups that vis-a-vis corresponds to a reduction in the average molecular dipole moment of resident water molecules due to smaller dipole induction in the weaker electrostatic fields of the nonpolar groups.

The Nature of Resonance in Allyl Ions and Radical

A recent valence bond scheme based on Lewis structures, the valence bond BOND (VBB) method (BOND: Breathing Orbitals Naturally Delocalized) method (Linares, M.; Braida, B.; Humbel, S. J. Phys. Chem. A 2006, 110, 2505-2509), is applied to explore the nature of resonance in allyl systems. Whereas allyl radical is correctly described by the resonance between the two traditional Lewis structures, a third "long-bonded" structure, which apparently creates a pi bond between the two distant carbon atoms, appears to plays an important role in allyl ions description. The similar vertical resonance energy (VRE) for both allyl ions is rather moderate (approximately 37 kcal/mol) in the two-structure description but is significantly enhanced when the long-bonded structure is included into the VBB wave function (by up to 20 kcal/mol). The allyl radical is much less resonant and is correctly described by the traditional two-structure picture. The development of VBB Lewis structures into "pure" valence bond determinants enlightens the role of the third structure in the description of allyl ions. The existence of a long bond between the two distant carbon atoms is clearly ruled out. Charge equilibration effect is shown to be a minor factor. The third structure is finally attributed to one- and three-electron bonding character revealed in the pi systems of the cation and anion, respectively. This makes these systems two surprising examples of odd electron bonding within a singlet state. Last, the two-structure description of allyl radical is improved by addition of missing ionic structures.

Catalytic activity correlation of Ni(II), Co(II) and Pd(II) complexes to metal atom net charge

The metal atom net charge correlation (MANCC) method was developed in prediction of catalyst activity of asymmetric late-transition metal complexes, 2-quinoxalinyl-6-iminopyridine Ni (II), 2-imino-1, 10-phenanthroline Co(II) and 2-methoxycarbonyl-6-iminopyridine Pd(II) complexes, from the net charge of the metal atom for ethylene polymerization. Dreiding force field was modified according to the X-ray diffraction data. We found that the asymmetric structure of the complexes resulted in a charge difference between two halogen atoms coordinated to the metal atom. In order to remove such contribution we introduced the effective charge Qeff, which was obtained by the charge equilibration (QEq) approach. The results verified the successful introduction of Qeff and showed that the catalytic activities of different complexes are related to central metal atom effective charge.

Real-Time QCM-D Monitoring of Electrostatically Driven Lipid Transfer between Two Lipid Bilayer Membranes

The lipid exchange/transfer between lipid membranes is important for many biological functions. To learn more about how the dynamics of such processes can be studied, we have investigated the interaction of positively and negatively charged lipid vesicles with supported lipid bilayers (SLBs) of opposite charge. The vesicle-SLB interaction leads initially to adsorption of lipid vesicles on the SLB, as deduced from the mass uptake kinetics and the concerted increase in dissipation, monitored by the quartz crystal microbalance with dissipation (QCM-D) technique. Eventually, however, vesicles (and possibly other lipid structures) desorb from the SLB surface, as judged from the mass loss and the dissipation decrease. The mass loss is approximately as large as the initial mass increase; i.e., at the end of the process the mass load is that of a SLB. We interpret this interesting kinetics in terms of initial strong electrostatic attraction between the added vesicles and the SLB, forming a structure where lipid transfer between the two bilayers occurs on a time scale of 10-40 min. We suggest that this lipid transfer causes a charge equilibration with an accompanying weakening of the attraction, and eventually repulsion, between the SLB and vesicles, leading to desorption of vesicles from the SLB. The composition of the latter has thus been modified compared to the initial one, although no net mass increase or decrease has occurred. Direct evidence for the lipid exchange was obtained by sequential experiments with alternating positive and negative vesicles, as well as by using fluorescently labeled lipids and FRAP. The above interpretation was further strengthened by combined QCM-D and optical reflectometry measurements.