Simple and accurate scheme to compute electrostatic interaction: Zero-dipole summation technique for molecular system and application to bulk water

The zero-dipole summation method was extended to general molecular systems, and then applied to molecular dynamics simulations of an isotropic water system. In our previous paper [I. Fukuda, Y. Yonezawa, and H. Nakamura, J. Chem. Phys. 134, 164107 (2011)], for evaluating the electrostatic energy of a classical particle system, we proposed the zero-dipole summation method, which conceptually prevents the nonzero-charge and nonzero-dipole states artificially generated by a simple cutoff truncation. Here, we consider the application of this scheme to molecular systems, as well as some fundamental aspects of general cutoff truncation protocols. Introducing an idea to harmonize the bonding interactions and the electrostatic interactions in the scheme, we develop a specific algorithm. As in the previous study, the resulting energy formula is represented by a simple pairwise function sum, enabling facile applications to high-performance computation. The accuracy of the electrostatic energies calculated by the zero-dipole summation method with the atom-based cutoff was numerically investigated, by comparison with those generated by the Ewald method. We obtained an electrostatic energy error of less than 0.01% at a cutoff length longer than 13 A for a TIP3P isotropic water system, and the errors were quite small, as compared to those obtained by conventional truncation methods. The static property and the stability in an MD simulation were also satisfactory. In addition, the dielectric constants and the distance-dependent Kirkwood factors were measured, and their coincidences with those calculated by the particle mesh Ewald method were confirmed, although such coincidences are not easily attained by truncation methods. We found that the zero damping-factor gave the best results in a practical cutoff distance region. In fact, in contrast to the zero-charge scheme, the damping effect was insensitive in the zero-charge and zero-dipole scheme, in the molecular system we treated. We discussed the origin of this difference between the two schemes and the dependence of this fact on the physical system. The use of the zero damping-factor will enhance the efficiency of practical computations, since the complementary error function is not employed. In addition, utilizing the zero damping-factor provides freedom from the parameter choice, which is not trivial in the zero-charge scheme, and eliminates the error function term, which corresponds to the time-consuming Fourier part under the periodic boundary conditions.     

A nano-confined charged layer defies the principle of electrostatic interaction

The reactivity between two charged molecules and the activity of charged biomolecules are mainly governed by the principle of electrostatic interaction, i.e., like charges repel and opposite charges attract. In the present study it is shown that the principle of electrostatic interaction is violated in the nano-confined biomimetic environment. Thus a positively charged molecule shows more preference to a positively charged surface compared to a negatively charged surface.     

Long-range electrostatic interactions in hybrid quantum and molecular mechanical dynamics using a lattice summation approach

A method based on a lattice summation technique for treating long-range electrostatic interactions in hybrid quantum mechanics/molecular mechanics simulations is presented in this article. The quantum subsystem is studied at the semiempirical level, whereas the solvent is described by a two-body potential of molecular mechanics. Molecular dynamics simulations of a (quantum) chloride ion in (classical) water have been performed to test this technique. It is observed that the application of the lattice summations to solvent-solvent interactions as well as on solute-solvent ones has a significant effect on solvation energy and diffusion coefficient. Moreover, two schemes for the computation of the long-range contribution to the electrostatic interaction energy are investigated. The first one replaces the exact charge distribution of the quantum solute by a Mulliken charge distribution. The long-range electrostatic interactions are then calculated for this charge distribution that interacts with the solvent molecule charges. The second one is more accurate and involves a modified Fock operator containing long-range electron-charge interactions. It is shown here that both schemes lead to similar results, the method using Mulliken charges for the evaluation of long-range interactions being, however, much more computationally efficient.     

Efficient and precise solvation free energies via alchemical adiabatic molecular dynamics

A new molecular dynamics method for calculating free energies associated with transformations of the thermodynamic state or chemical composition of a system (also known as alchemical transformations) is presented. The new method extends the adiabatic dynamics approach recently introduced by Rosso et al. [J. Chem. Phys. 116, 4389 (2002)] and is based on the use of an additional degree of freedom, lambda, that is used as a switching parameter between the potential energy functions that characterize the two states. In the new method, the coupling parameter lambda is introduced as a fictitious dynamical variable in the Hamiltonian, and a system of switching functions is employed that leads to a barrier in the lambda free energy profile between the relevant thermodynamic end points. The presence of such a barrier, therefore, enhances sampling in the end point (lambda = 0 and lambda = 1) regions which are most important for computing relevant free energy differences. In order to ensure efficient barrier crossing, a high temperature T(lambda) is assigned to lambda and a fictitious mass m(lambda) is introduced as a means of creating an adiabatic separation between lambda and the rest of the system. Under these conditions, it is shown that the lambda free energy profile can be directly computed from the adiabatic probability distribution function of lambda without any postprocessing or unbiasing of the output data. The new method is illustrated on two model problems and in the calculation of the solvation free energy of amino acid side-chain analogs in TIP3P water. Comparisons to previous work using thermodynamic integration and free energy perturbation show that the new lambda adiabatic free energy dynamics method results in very precise free energy calculations using significantly shorter trajectories.     

Fast and spectrally accurate Ewald summation for 2-periodic electrostatic systems

A new method for Ewald summation in planar/slablike geometry, i.e., systems where periodicity applies in two dimensions and the last dimension is "free" (2P), is presented. We employ a spectral representation in terms of both Fourier series and integrals. This allows us to concisely derive both the 2P Ewald sum and a fast particle mesh Ewald (PME)-type method suitable for large-scale computations. The primary results are: (i) close and illuminating connections between the 2P problem and the standard Ewald sum and associated fast methods for full periodicity; (ii) a fast, O(N log N), and spectrally accurate PME-type method for the 2P k-space Ewald sum that uses vastly less memory than traditional PME methods; (iii) errors that decouple, such that parameter selection is simplified. We give analytical and numerical results to support this.     

Molecular dynamics simulation with a continuum electrostatic model of the solvent

The accuracy and simplicity of the Poisson-Boltzmann electrostatics model has led to the suggestion that it might offer an efficient solvent model for use in molecular mechanics calculations on biomolecules. We report a successful merger of the Poisson-Boltzmann and molecular dynamics approaches, with illustrative calculations on the small solutes dichloroethane and alanine dipeptide. The algorithm is implemented within the program UHBD. Computational efficiency is achieved by the use of rather coarse finite difference grids to solve the Poisson-Boltzmann equation. Nonetheless, the conformational distributions generated by the new method agree well with reference distributions obtained as Boltzmann distributions from energies computed with fine finite difference grids. The conformational distributions also agree well with the results of experimental measurements and conformational analyses using more detailed solvent models. We project that when multigrid methods are used to solve the finite difference problem and the algorithm is implemented on a vector supercomputer, the computation of solvent electrostatic forces for a protein of modest size will add only about 0.1 s computer time per simulation step relative to a vacuum calculation. © 1995 by John Wiley & Sons, Inc.     

Ewald summation versus direct summation of shifted-force potentials for the calculation of electrostatic interactions in solids: A quantitative study

The calculated Madelung energies and Madelung forces of the electrostatic interaction for nine crystal structures are reported. The method of direct summation with two different shifted-force potentials is compared to the Ewald summation. There is a considerable difference in the convergence of the energy and the force for the two shifted-force potentials regarding the cutoff radius. The convergence depends not only on the potential itself, but also on the crystal structure. One of the shifted-force potentials used is implemented in the CHARMM force field. The energy calculated with this potential shows a good convergence for small cutoff radii. With the other shifted-force potential, the force shows a better convergence for small cutoff radii. The number of pair interactions for obtaining the Madelung limit using the Ewald summation and the direct summation of a shifted-force potential is also reported. For complex structures like zeolites, the number of relevant pair interactions is smaller using the direct summation of a shifted-force potential. For simple structures such as cesium chloride, the number of significant pair interactions is smaller using the Ewald summation. © 1997 by John Wiley & Sons, Inc.     

Computer simulation of the assembly of gold nanoparticles on DNA fragments via electrostatic interaction

Using Monte Carlo simulation, we study the metallization of DNA fragments via the templating of gold nanoparticles. To represent the interaction between metal entities, a nanoparticle-nanoparticle interaction potential was derived on the basis of the many-body Gupta potential. The aggregation of the nanoparticles on the template surface is due to the additive effect of electrostatic attraction between the positive charges on the Au particles and the negative charges of the phosphate groups of DNA molecule and the short-range attraction between the metallic nanoparticles. As a result, the assembly of a continuous nanowire can be templated. Depending on the nanoparticle size and charge, the metallic covering can be both continuous and discontinuous. The question of how size and charge of Au nanoparticles influence the structure of metallic coat is discussed in detail. Both monodisperse and polydisperse nanoparticles are considered. Dispersion in the nanoparticle size was found to have little effect on the calculated characteristics of the aggregate.     

Molecular dynamics simulations of multilayer polyelectrolyte films: effect of electrostatic and short-range interactions

The effect of the strength of electrostatic and short-range interactions on the multilayer assembly of oppositely charged polyelectrolytes at a charged substrate was studied by molecular dynamics simulations. The multilayer buildup was achieved through sequential adsorption of charged polymers in a layer-by-layer fashion from dilute polyelectrolyte solutions. The strong electrostatic attraction between oppositely charged polyelectrolytes at each deposition step is a driving force behind the multilayer growth. Our simulations have shown that a charge reversal after each deposition step is critical for steady multilayer growth and that there is a linear increase in polymer surface coverage after the first few deposition steps. Furthermore, there is substantial intermixing between chains adsorbed during different deposition steps. We show that the polymer surface coverage and multilayer structure are each strongly influenced by the strength of electrostatic and short-range interactions.     

Quantum monte carlo vibrational dynamics in a property-based interaction potential scheme for weakly bound clusters

The typical shallowness of the potential surfaces of weakly bound clusters implies sizable ground-state vibrational excursions in the weak modes, a feature often complicated by considerable anharmonicity. The difficulties of vibrational analysis are exacerbated as the number of weak modes increases with the number of molecules in a cluster. Quantum Monte Carlo (QMC) approaches offer a general suitability to the problem of vibrational dynamics of weakly bound clusters in that they can fully account for anharmonicity and large amplitude motions. We report on the effectiveness and convergence behavior of diffusion quantum Monte Carlo for both energies and the key spectroscopic values of vibrationally averaged rotational constants. QMC involves recurring evaluations of the interaction potential, and we show how property-based, two-and three-body potentials (e.g., those involving intrinsic molecular tensor properties) may be carefully linked to the QMC propagation steps. © 1997 by John Wiley & Sons, Inc.     

Dynamic viscosity estimation of hydrogen sulfide using a predictive scheme based on molecular dynamics

An approach based on molecular dynamics results on Lennard–Jones spheres is proposed to model the viscosity of hydrogen sulfide, H₂S. The molecular parameters, that have a strong physical meaning, are the depth of the potential, and the length at which the potential is null (the “molecular diameter”), which take into account the dipolar moment of the hydrogen sulfide through an isotropic dipolar approximation. The interest of the method is that the adjustment does not involve any viscosity data because only density values have been used in order to estimate the molecular parameters. Consequently, the model is entirely predictive. A comparison between the data generated by our model, REFPROP7 and REFPROP8 database and the few available experimental viscosity data (dilute gas and saturated liquid) is performed and it clearly demonstrates the performance of this predictive model. It is even shown that this model is, without fitting, slightly better than REFPROP7 and REFPROP8 which uses viscosity experimental database to adjust their parameters. In addition, in typical petroleum reservoirs conditions, it is shown that non-negligible deviations appear when comparing results predicted by REFPROP7, REFPROP8 and the model proposed. Due to its predictive nature, we believe that the values evaluated by the proposed model make sense in such reservoir conditions, at least for industrial purposes. Moreover, the scheme proposed is shown to be very easily extended to deal with mixtures involving H₂S with the limit that the Lennard–Jones fluid model is appropriate for the other species of the mixtures.     

Molecular dynamics study of perovskite structures with modified interatomic interaction potentials

The structure of compounds with the perovskite structure ABX₃ (A and B are cations, X are anions O^2—, F^—, Cl^—, Br^—, and I^—), which are widely used in engineering due to unique electrical, optical, and photovoltaic properties, has been considered. Hybrid organic—inorganic halide perovskites important for photovoltaics of a new generation are worth mentioning; they contain cations of organic nitrogen bases as monovalent cations. A molecular dynamics (MD) study of the CaTiO₃ base structure (Ca²⁺, Ti⁴⁺, and O^2—) has been performed in order to develop the methodology of computer simulation and optimization of the shape and parameters of atomic potentials for perovskite systems.     

Molecular dynamics simulation study of interaction between model rough hydrophobic surfaces

We study some aspects of hydrophobic interaction between molecular rough and flexible model surfaces. The model we use in this work is based on a model we used previously (Eun, C.; Berkowitz, M. L. J. Phys. Chem. B 2009, 113, 13222-13228), when we studied the interaction between model patches of lipid membranes. Our original model consisted of two graphene plates with attached polar headgroups; the plates were immersed in a water bath. The interaction between such plates can be considered as an example of a hydrophilic interaction. In the present work, we modify our previous model by removing the charge from the zwitterionic headgroups. As a result of this procedure, the plate character changes: it becomes hydrophobic. By separating the total interaction (or potential of mean force, PMF) between plates into the direct and the water-mediated interactions, we observe that the latter changes from repulsive to attractive, clearly emphasizing the important role of water as a medium. We also investigate the effect of roughness and flexibility of the headgroups on the interaction between plates and observe that roughness enhances the character of the hydrophobic interaction. The presence of a dewetting transition in a confined space between charge-removed plates confirms that the interaction between plates is strongly hydrophobic. In addition, we notice that there is a shallow local minimum in the PMF in the case of the charge-removed plates. We find that this minimum is associated with the configurational changes that flexible headgroups undergo as the two plates are brought together.     

Fast calculation of the electrostatic potential in ionic crystals by direct summation method

An efficient real space method is derived for the evaluation of the Madelung's potential of ionic crystals. The proposed method is an extension of Evjen's method. It takes advantage of a general analysis of the potential convergence in real space. Indeed, we show that the series convergence is exponential as a function of the number of canceled multipolar moments in the unit cell. The method proposed in this work reaches such an exponential convergence rate. Its efficiency is comparable to Ewald's method. However, unlike the latter, it uses only simple algebraic functions.     

ORAC: A Molecular dynamics program to simulate complex molecular systems with realistic electrostatic interactions

In this study, we present a new molecular dynamics program for simulation of complex molecular systems. The program, named ORAC, combines state-of-the-art molecular dynamics (MD) algorithms with flexibility in handling different types and sizes of molecules. ORAC is intended for simulations of molecular systems and is specifically designed to treat biomolecules efficiently and effectively in solution or in a crystalline environment. Among its unique features are: (i) implementation of reversible and symplectic multiple time step algorithms (or r-RESPA, reversible reference system propagation algorithm) specifically designed and tuned for biological systems with periodic boundary conditions; (ii) availability for simulations with multiple or single time steps of standard Ewald or smooth particle mesh Ewald (SPME) for computation of electrostatic interactions; and (iii) possibility of simulating molecular systems in a variety of thermodynamic ensembles. We believe that the combination of these algorithms makes ORAC more advanced than other MD programs using standard simulation algorithms. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1848–1862, 1997     

Extending electrochemical activity of polyaniline to alkaline media via electrostatic interaction and sol–gel route

The electroactivity of polyaniline was extended to pH = 14 alkaline media by preparation of a novel electrostatic interaction conductive hybrid from water-borne conductive polyaniline and silica network containing carboxyl groups via sol–gel process. In addition, the obtained conductive polyaniline hybrid film displayed very low conductivity threshold percolation and demonstrated excellent stability upon cycling.     

Molecular dynamics simulation and binding energy calculation for estimation of oligonucleotide duplex thermostability in RNA-based therapeutics

For oligonucleotide-based therapeutics, a thorough understanding of the thermodynamic properties of duplex formation is critical to developing stable and potent drugs. For unmodified small interfering RNA (siRNA), DNA antisense oligonucleotide (AON) and locked nucleic acid (LNA), DNA/LNA modified oligonucleotides, nearest neighbor (NN) methods can be effectively used to quickly and accurately predict duplex thermodynamic properties such as melting point. Unfortunately, for chemically modified olignonucleotides, there has been no accurate prediction method available. Here we describe the potential of estimating melting temperature (T(m)) for nonstandard oligonucleotides by using the correlation of the experimental T(m) with the calculated duplex binding energy (BE) for oligonucleotides of a given length. This method has been automated into a standardized molecular dynamics (MD) protocol through Pipeline Pilot (PP) using the CHARMm component in Discovery Studio (DS). Results will be presented showing the correlation of the predicted data with experiment for both standard and chemically modified siRNA and AON.     

Errors in the Ewald summation method with electrostatic layer correction for Coulomb interactions in slab geometry

New theoretical expressions for cut-off errors in 2D reciprocal-space summation of the electrostatic layer correction (ELC) term in energy and forces are derived, and a procedure to determine optimal parameters of the method is proposed. The procedure is tested in numerical calculations for charges distributed uniformly in a cubic box and charges located in two layers near the box basis. The summation errors for conventional Ewald method can be used to find out optimal values of the convergence parameter, and real- and reciprocal-space cut-off radii, whereas the ELC errors give possibility to choose an optimal value of an empty space gap in the simulation box.     

Molecular dynamics simulation study of the interaction of trehalose with lipid membranes

The interactions of the cryoprotective agent trehalose with a lipid membrane made of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine at 323 K were studied by means of molecular dynamics simulations. It was observed that trehalose binds to the phospholipid headgroups with its main axis parallel to the membrane normal. Trehalose establishes hydrogen bonds with the carbonyl and phosphate groups and replaces water molecules from the lipid headgroup. Notably, the number of hydrogen bonds (HBs) that the membrane made with its environment was conserved after trehalose binding. The HBs between lipid and trehalose have a longer lifetime than those established between lipid and water. The binding of the sugar does not produce changes either in the lipid area or in the lipid order parameter. The effect of trehalose on the dipole potential is in agreement with experimental results. The contribution of the different components to the membrane dipole potential was analyzed. It was observed that the binding of trehalose produces changes in the different components and the sugar itself contributes to the surface potential due to the polarization of its hydroxyl in the interface.