Molecular dynamics scheme for precise estimation of electrostatic interaction via zero-dipole summation principle

We propose a novel idea, zero-dipole summation, for evaluating the electrostatic energy of a classical particle system, and have composed an algorithm for effectively utilizing the idea for molecular dynamics. It conceptually prevents the nonzero-charge and nonzero-dipole states artificially generated by a simple cutoff truncation. The resulting energy formula is nevertheless represented by a simple pairwise function sum, which enables facile application to high-performance computation. By following a heuristic approach to derive the current electrostatic energy formula, we developed an axiomatic approach to construct the method consistently. Explorations of the theoretical details of our method revealed the structure of the generated error, and we analyzed it by comparisons with other methods. A numerical simulation using liquid sodium chloride confirmed that the current method with a small damping factor yielded sufficient accuracy with a practical cutoff distance region. The current energy function also conducts stable numerical integration in a liquid MD simulation. Our method is an extension of the charge neutralized summation developed by Wolf et al. [J. Chem. Phys. 110, 8254 (1999)]. Furthermore, we found that the current method becomes a generalization of the preaveraged potential method proposed by Yakub and Ronchi [J. Chem. Phys. 119, 11556 (2003)], which is based on a viewpoint different from the neutrality. The current study presents these relationships and suggests possibilities for their further applications.     

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.     

On the use of AM1 and MNDO wave functions to compute accurate electrostatic charges

A new strategy to evaluate accurate electrostatic charges from semiempirical wave functions is reported. The rigorous quantum mechanical molecular electrostatic potentials computed from both MNDO and AM1 wave functions are fitted to the point-charge molecular electrostatic potential to obtain the electrostatic charges. The reliability of this strategy is tested by comparing the semiempirical electrostatic charges for 21 molecules with the semiempirical Mulliken charges and with the ab initio STO-3G and 6-31G* electrostatic charges. The ability of the dipoles derived from the semiempirical electrostatic and Mulliken charges as well as from the SCF charge distributions to reproduce the ab initio 6-31G* electrostatic dipoles and the gas phase experimental values is determined. The statistical analysis clearly point out the goodness of the semiempirical electrostatic charges, specially when the MNDO method is used. The excellent relationships found between the MNDO and 6-31G* electrostatic charges permit to define a scaling factor which allows to accurately reproduce the 6-31G* electrostatic charge distribution as well as the experimental dipoles from the semiempirical electrostatic charges.     

Comparison of NDDO and quasi-ab initio approaches to compute semiempirical molecular electrostatic potentials

The suitability of the two most widely used strategies to compute semiempirical MEPs is examined. For this purpose, MEP minima, electrostatic charges, and dipoles for a large number of molecules were computed at the AM1, MNDO, and PM3 levels using both the NDDO strategy developed by Ferenczy, Reynolds, and Richards and our own quasi-ab initio method. Results demonstrate that the quasi-ab initio is preferred over the NDDO method for the computation of MEP minima. It is also found that the best set of semiempirical charges and dipoles are obtained using either the AM1 NDDO or the MNDO quasi-ab initio methods. In these two cases, the quality of the results is fully comparable with 6-31G* values. © 1994 by John Wiley & Sons, Inc.     

A sum-over-state scheme of analysis of hyperpolarizabilities and its application to spiroconjugated molecular system

The static first and second hyperpolarizabilities of a number of spiromolecules with varying degree of polarity have been calculated at the HF and MP2 level using the 6-31+G* basis set and the B3LYP/6-31+G* optimized geometry. The variation of mean second hyperpolarizability in these molecular systems has been explained in terms of the ground state dipole moment, mean linear polarizability and second-order polarizability. A number of relationships among these quantities have been derived in the framework of the sum-over-state scheme and the generalized Thomas–Kuhn sum rule. The spiroconjugation results in the significant increase of the mean polarizability. The appreciable enhancement of first hyperpolarizability due to the spiroconjugation between two dipolar monomer units has been accounted for the rather significant increase of the mean polarizability tensor and the ground state dipole moment. The relatively larger value of the average second hyperpolarizability of the spiroconjugated molecules compared to that of the corresponding monomers arises from the rather significant increase of the nonaxial component γ _xxyy . The replacement of spirocarbon by spirosilicon results in the enhancement of the cubic polarizability manifold. The donor–acceptor substituted spirocompounds are predicted to be the superior third-order nonlinear optical (NLO) phores. The nature of π-conjugation in the monomer units around the spirocenter shows a strong modulation of the NLO properties of spirocompounds. The influence of electron correlation on the NLO properties at the MP2 level has been found to be rather significant.     

An efficient parallelization scheme for molecular dynamics simulations with many-body, flexible, polarizable empirical potentials: application to water

An efficient parallelization scheme for classical molecular dynamics simulations with flexible, polarizable empirical potentials is presented. It is based on the standard Ewald summation technique to handle the long-range electrostatic and induction interactions. The algorithm for this parallelization scheme is designed for systems containing several thousands of polarizable sites in the simulation box. Its performance is evaluated during molecular dynamics simulations under periodic boundary conditions with unit cell sizes ranging from 128 to 512 molecules employing two flexible polarizable water models [DC(F) and TTM2.1-F] containing 1 and 3 polarizable sites, respectively. The time-to-solution for these two polarizable models is compared with the one for a flexible, pairwise-additive water model (TIP4F). The benchmarks were performed on both shared and distributed memory platforms. As a result of the efficient calculation of the induced dipole moments, a superlinear scaling as a function of the number of the processors is observed. To the best of our knowledge, this is the first reported results of parallel scaling and performance for simulations of liquid water with a polarizable potential under periodic boundary conditions.     

Pair interaction potentials with explicit polarization for molecular dynamics simulations of La(3+) in bulk water

Pair interaction potentials (IPs) were defined to describe the La(3+)-OH(2) interaction for simulating the La(3+) hydration in aqueous solution. La(3+)-OH(2) IPs are taken from the literature or parametrized essentially to reproduce ab initio calculations at the second-order Moller-Plesset level of theory on La(H(2)O)(8) (3+). The IPs are compared and used with molecular dynamics (MD) including explicit polarization, periodic boundary conditions of La(H(2)O)(216) (3+) boxes, and TIP3P water model modified to include explicit polarization. As expected, explicit polarization is crucial for obtaining both correct La-O distances (r(La-O)) and La(3+) coordination number (CN). Including polarization also modifies hydration structure up to the second hydration shell and decreases the number of water exchanges between the La(3+) first and second hydration shells. r(La-O) ((1))=2.52 A and CN((1))=9.02 are obtained here for our best potential. These values are in good agreement with experimental data. The tested La-O IPs appear to essentially account for the La-O short distance repulsion. As a consequence, we propose that most of the multibody effects are correctly described by the explicit polarization contributions even in the first La(3+) hydration shell. The MD simulation results are slightly improved by adding a-typically negative 1r(6)-slightly attractive contribution to the-typically exponential-repulsive term of the La-O IP. Mean residence times are obtained from MD simulations for a water molecule in the first (1082 ps) and second (7.6 ps) hydration shells of La(3+). The corresponding water exchange is a concerted mechanism: a water molecule leaving La(H(2)O)(9) (3+) in the opposite direction to the incoming water molecule. La(H(2)O)(9) (3+) has a slightly distorded "6+3" tricapped trigonal prism D(3h) structure, and the weakest bonding is in the medium triangle, where water exchanges take place.     

A divide and conquer real-space approach for all-electron molecular electrostatic potentials and interaction energies

A computational scheme to perform accurate numerical calculations of electrostatic potentials and interaction energies for molecular systems has been developed and implemented. Molecular electron and energy densities are divided into overlapping atom-centered atomic contributions and a three-dimensional molecular remainder. The steep nuclear cusps are included in the atom-centered functions making the three-dimensional remainder smooth enough to be accurately represented with a tractable amount of grid points. The one-dimensional radial functions of the atom-centered contributions as well as the three-dimensional remainder are expanded using finite element functions. The electrostatic potential is calculated by integrating the Coulomb potential for each separate density contribution, using our tensorial finite element method for the three-dimensional remainder. We also provide algorithms to compute accurate electron-electron and electron-nuclear interactions numerically using the proposed partitioning. The methods have been tested on all-electron densities of 18 reasonable large molecules containing elements up to Zn. The accuracy of the calculated Coulomb interaction energies is in the range of 10(-3) to 10(-6) E(h) when using an equidistant grid with a step length of 0.05 a(0).     

Many-body interaction analysis: algorithm development and application to large molecular clusters

A completely automated algorithm for performing many-body interaction energy analysis of clusters (MBAC) [M. J. Elrodt and R. J. Saykally, Chem. Rev. 94, 1975 (1994); S. S. Xantheas, J. Chem. Phys. 104, 8821 (1996)] at restricted Hartree-Fock (RHF)/MA Plesset 2nd order perturbation theory (MP2)/density functional theory (DFT) level of theory is reported. Use of superior guess density matrices (DM's) for smaller fragments generated from DM of the parent system and elimination of energetically insignificant higher-body combinations, leads to a more efficient performance (speed-up up to 2) compared to the conventional procedure. MBAC approach has been tested out on several large-sized weakly bound molecular clusters such as (H(2)O)(n), n=8, 12, 16, 20 and hydrated clusters of amides and aldehydes. The MBAC results indicate that the amides interact more strongly with water than aldehydes in these clusters. It also reconfirms minimization of the basis set superposition error for large cluster on using superior quality basis set. In case of larger weakly bound clusters, the contributions higher than four body are found to be repulsive in nature and smaller in magnitude. The reason for this may be attributed to the increased random orientations of the interacting molecules separated from each other by large distances.     

FAMBE-pH: a fast and accurate method to compute the total solvation free energies of proteins

A fast and accurate method to compute the total solvation free energies of proteins as a function of pH is presented. The method makes use of a combination of approaches, some of which have already appeared in the literature; (i) the Poisson equation is solved with an optimized fast adaptive multigrid boundary element (FAMBE) method; (ii) the electrostatic free energies of the ionizable sites are calculated for their neutral and charged states by using a detailed model of atomic charges; (iii) a set of optimal atomic radii is used to define a precise dielectric surface interface; (iv) a multilevel adaptive tessellation of this dielectric surface interface is achieved by using multisized boundary elements; and (v) 1:1 salt effects are included. The equilibrium proton binding/release is calculated with the Tanford-Schellman integral if the proteins contain more than approximately 20-25 ionizable groups; for a smaller number of ionizable groups, the ionization partition function is calculated directly. The FAMBE method is tested as a function of pH (FAMBE-pH) with three proteins, namely, bovine pancreatic trypsin inhibitor (BPTI), hen egg white lysozyme (HEWL), and bovine pancreatic ribonuclease A (RNaseA). The results are (a) the FAMBE-pH method reproduces the observed pK a's of the ionizable groups of these proteins within an average absolute value of 0.4 p K units and a maximum error of 1.2 p K units and (b) comparison of the calculated total pH-dependent solvation free energy for BPTI, between the exact calculation of the ionization partition function and the Tanford-Schellman integral method, shows agreement within 1.2 kcal/mol. These results indicate that calculation of total solvation free energies with the FAMBE-pH method can provide an accurate prediction of protein conformational stability at a given fixed pH and, if coupled with molecular mechanics or molecular dynamics methods, can also be used for more realistic studies of protein folding, unfolding, and dynamics, as a function of pH.     

A hierarchical construction scheme for accurate potential energy surface generation: an application to the F+H2 reaction

We present a hierarchical construction scheme for accurate ab initio potential energy surface generation. The scheme is based on the observation that when molecular configuration changes, the variation in the potential energy difference between different ab initio methods is much smaller than the variation for potential energy itself. This means that it is easier to numerically represent energy difference to achieve a desired accuracy. Because the computational cost for ab initio calculations increases very rapidly with the accuracy, one can gain substantial saving in computational time by constructing a high accurate potential energy surface as a sum of a low accurate surface based on extensive ab initio data points and an energy difference surface for high and low accuracy ab initio methods based on much fewer data points. The new scheme was applied to construct an accurate ground potential energy surface for the FH(2) system using the coupled-cluster method and a very large basis set. The constructed potential energy surface is found to be more accurate on describing the resonance states in the FH(2) and FHD systems than the existing surfaces.     

Simple and accurate method for determining dissolved inorganic carbon in environmental water by reaction headspace gas chromatography

We investigate a simple and accurate method for quantitatively analyzing dissolved inorganic carbon in environmental water by reaction headspace gas chromatography. The neutralization reaction between the inorganic carbon species (i.e. bicarbonate ions and carbonate ions) in environmental water and hydrochloric acid is carried out in a sealed headspace vial, and the carbon dioxide formed from the neutralization reaction, the self‐decomposition of carbonic acid, and dissolved carbon dioxide in environmental water is then analyzed by headspace gas chromatography. The data show that the headspace gas chromatography method has good precision (relative standard deviation ≤ 1.63%) and accuracy (relative differences ≤ 5.81% compared with the coulometric titration technique). The headspace gas chromatography method is simple, reliable, and can be well applied in the dissolved inorganic carbon detection in environmental water.     

Simple and effective application of the Wang-Landau method for multicanonical molecular dynamics simulation

We propose a novel application of the Wang-Landau method (WLM) for multicanonical molecular dynamics (McMD) simulations. Originally, WLM was developed for Monte Carlo (MC) simulations. Fundamentally, WLM remarkably reduces simulation efforts because it estimates the optimal multicanonical energy function automatically. When WLM is applied to McMD, not only the multicanonical energy but also energy gradient must be estimated adequately. However, because of the rugged multicanonical energy function at the early simulation stage, applications of WLM for MD simulations are difficult and require a smoothing procedure: simulation efforts such as cubic-spline extrapolation and gathering multiple preruns are utilized for smoothing. We propose a simple and effective smoothing method that requires only one additional equation and two time-dependent parameters. As a result, our method produced the correct multicanonical energy function and succeeded in the flat sampling of a small biomolecule with reduced simulation effort.     

Constructing simple yet accurate potentials for describing the solvation of HCl/water clusters in bulk helium and nanodroplets

The infrared spectroscopy of molecules, complexes, and molecular aggregates dissolved in superfluid helium clusters, commonly called HElium NanoDroplet Isolation (HENDI) spectroscopy, is an established, powerful experimental technique for extracting high resolution ro-vibrational spectra at ultra-low temperatures. Realistic quantum simulations of such systems, in particular in cases where the solute is undergoing a chemical reaction, require accurate solute-helium potentials which are also simple enough to be efficiently evaluated over the vast number of steps required in typical Monte Carlo or molecular dynamics sampling. This precludes using global potential energy surfaces as often parameterized for small complexes in the realm of high-resolution spectroscopic investigations that, in view of the computational effort imposed, are focused on the intermolecular interaction of rigid molecules with helium. Simple Lennard-Jones-like pair potentials, on the other hand, fall short in providing the required flexibility and accuracy in order to account for chemical reactions of the solute molecule. Here, a general scheme of constructing sufficiently accurate site-site potentials for use in typical quantum simulations is presented. This scheme employs atom-based grids, accounts for local and global minima, and is applied to the special case of a HCl(H(2)O)(4) cluster solvated by helium. As a first step, accurate interaction energies of a helium atom with a set of representative configurations sampled from a trajectory following the dissociation of the HCl(H(2)O)(4) cluster were computed using an efficient combination of density functional theory and symmetry-adapted perturbation theory, i.e. the DFT-SAPT approach. For each of the sampled cluster configurations, a helium atom was placed at several hundred positions distributed in space, leading to an overall number of about 400,000 such quantum chemical calculations. The resulting total interaction energies, decomposed into several energetic contributions, served to fit a site-site potential, where the sites are located at the atomic positions and, additionally, pseudo-sites are distributed along the lines joining pairs of atom sites within the molecular cluster. This approach ensures that this solute-helium potential is able to describe both undissociated molecular and dissociated (zwitter-) ionic configurations, as well as the interconnecting reaction pathway without re-adjusting partial charges or other parameters depending on the particular configuration. Test calculations of the larger HCl(H(2)O)(5) cluster interacting with helium demonstrate the transferability of the derived site-site potential. This specific potential can be readily used in quantum simulations of such HCl/water clusters in bulk helium or helium nanodroplets, whereas the underlying construction procedure can be generalized to other molecular solutes in other atomic solvents such as those encountered in rare gas matrix isolation spectroscopy.     

Local control of antibody binding to hapten-presenting interfaces: Steric and electrostatic interaction

The binding of labeled antibodies to hapten substituted monolayers at the air/water interface has been studied by means of fluorescence microscopy. Haptens with various spacer lengths between the epitope and a hydrocarbon chain, anchoring the molecule to the interface, have been synthesized. With DMPC,^a unspecific binding has been shown to predominate over specific binding due to electrostatic interactions. At high surface pressures the bound antibody is detached because of steric interference with the lipid head groups. Due to a reduction of electrostatic interactions, no unspecific binding is observed to monolayers of cholesterol, which carries a small dipole moment. Mixed monolayers of cholesterol and DMPC separate into two fluid phases, with preferential antibody binding to the cholesterol-enriched phase.     

Elastic bag model for molecular dynamics simulations of solvated systems: application to liquid water and solvated peptides

The fluctuating elastic boundary (FEB) model for molecular dynamics has recently been developed and validated through simulations of liquid argon. In the FEB model, a flexible boundary which consists of particles connected by springs is used to confine the solvated system, thereby eliminating the need for periodic boundary conditions. In this study, we extend this model to the simulation of bulk water and solvated alanine dipeptide. Both the confining potential and boundary particle interaction functions are modified to preserve the structural integrity of the boundary and prevent the leakage of the solute-solvent system through the boundary. A broad spectrum of structural and dynamic properties of liquid water are computed and compared with those obtained from conventional periodic boundary condition simulations. The applicability of the model to biomolecular simulations is investigated through the analysis of conformational population distribution of solvated alanine dipeptide. In most cases we find remarkable agreement between the two simulation approaches.     

Energy screening for the incremental scheme: application to intermolecular interactions

A systematic screening procedure for small contributions in the incremental expansion of the correlation energy is presented. The performance of the proposed scheme is checked for the calculation of intermolecular interactions in realistic test systems as large as a guanine-cytosine base pair. It is found that the computational cost for the incremental expansion can be reduced considerably without significant loss of accuracy. Typically, the errors of the systems investigated here amount to <3.4, 0.22, and 0.06% for second-, third-, and fourth-order expansions, respectively. For almost all cases, the error in the total correlation energy can be kept below 1 kcal/mol with respect to the canonical CCSD result if the incremental series is truncated in a proper way.     

Noise reduction method for molecular interaction energy: application to in silico drug screening and in silico target protein screening

We developed a new method to improve the accuracy of molecular interaction data using a molecular interaction matrix. This method was applied to enhance the database enrichment of in silico drug screening and in silico target protein screening using a protein-compound affinity matrix calculated by a protein-compound docking software. Our assumption was that the protein-compound binding free energy of a compound could be improved by a linear combination of its docking scores with many different proteins. We proposed two approaches to determine the coefficients of the linear combination. The first approach is based on similarity among the proteins, and the second is a machine-learning approach based on the known active compounds. These methods were applied to in silico screening of the active compounds of several target proteins and in silico target protein screening.     

From cluster to bulk: size dependent energetics of silica and silica-water interaction

We present our computational investigations on the energetics of clusters that consist of H2O and SiO2 using first-principles Born-Oppenheimer molecular dynamics method. Cohesive energy and hydration energy of both pure (or dry) and hydroxylated (or wet) ring-structured clusters have been investigated as functions of system size. We have found clear trends of energy as the cluster size increases. Energetics of a small silica nano-rod that contains 108 atoms is also obtained as a middle reference point for size evolution. Results from cluster and nano-rod calculations are compared with values from bulk quartz calculations using the same level of theoretical treatments.