Simulated annealing with adaptive cooling rates.

As one of the most robust global optimization methods, simulated annealing has received considerable attention with many variations that attempt to improve the cooling schedule. This paper introduces a variant of molecular dynamics-based simulated annealing that is useful for optimizing atomistic structures, and makes use of the heat capacity of the system, determined on the fly during optimization, to adaptively control the cooling rate. This adaptive cooling approach is demonstrated to be more computationally efficient than classical simulated annealing when applied to Lennard-Jones clusters. The increase in efficiency is approximately a factor of two for clusters with 25-40 atoms, and improves as the size of the system increases.

Cluster analysis of molecular simulation trajectories for systems where both conformation and orientation of the sampled states are important.

Clustering methods have been widely used to group together similar conformational states from molecular simulations of biomolecules in solution. For applications such as the interaction of a protein with a surface, the orientation of the protein relative to the surface is also an important clustering parameter because of its potential effect on adsorbed-state bioactivity. This study presents cluster analysis methods that are specifically designed for systems where both molecular orientation and conformation are important, and the methods are demonstrated using test cases of adsorbed proteins for validation. Additionally, because cluster analysis can be a very subjective process, an objective procedure for identifying both the optimal number of clusters and the best clustering algorithm to be applied to analyze a given dataset is presented. The method is demonstrated for several agglomerative hierarchical clustering algorithms used in conjunction with three cluster validation techniques. © 2016 Wiley Periodicals, Inc.

TIGER2 with solvent energy averaging (TIGER2A): An accelerated sampling method for large molecular systems with explicit representation of solvent.

The recently developed "temperature intervals with global exchange of replicas" (TIGER2) accelerated sampling method is found to have inaccuracies when applied to systems with explicit solvation. This inaccuracy is due to the energy fluctuations of the solvent, which cause the sampling method to be less sensitive to the energy fluctuations of the solute. In the present work, the problem of the TIGER2 method is addressed in detail and a modification to the sampling method is introduced to correct this problem. The modified method is called "TIGER2 with solvent energy averaging," or TIGER2A. This new method overcomes the sampling problem with the TIGER2 algorithm and is able to closely approximate Boltzmann-weighted sampling of molecular systems with explicit solvation. The difference in performance between the TIGER2 and TIGER2A methods is demonstrated by comparing them against analytical results for simple one-dimensional models, against replica exchange molecular dynamics (REMD) simulations for sampling the conformation of alanine dipeptide and the folding behavior of (AAQAA)3 peptide in aqueous solution, and by comparing their performance in sampling the behavior of hen egg-white lysozyme in aqueous solution. The new TIGER2A method solves the problem caused by solvent energy fluctuations in TIGER2 while maintaining the two important characteristics of TIGER2, i.e., (1) using multiple replicas sampled at different temperature levels to help systems efficiently escape from local potential energy minima and (2) enabling the number of replicas used for a simulation to be independent of the size of the molecular system, thus providing an accelerated sampling method that can be used to efficiently sample systems considered too large for the application of conventional temperature REMD.

Simulation of multiphase systems utilizing independent force fields to control intraphase and interphase behavior.

Fixed-charge empirical force fields have been developed and widely used over the past three decades for all-atom molecular simulations. Most simulation programs providing these methods enable only one set of force field parameters to be used for the entire system. Whereas this is generally suitable for single-phase systems, the molecular environment at the interface between two phases may be sufficiently different from the individual phases to require a different set of parameters to be used to accurately represent the system. Recently published simulations of peptide adsorption to material surfaces using the CHARMM force field have clearly demonstrated this issue by revealing that calculated values of adsorption free energy substantially differ from experimental results. Whereas nonbonded parameters could be adjusted to correct this problem, this cannot be done without also altering the conformational behavior of the peptide in solution, for which CHARMM has been carefully tuned. We have developed a dual-force-field approach (Dual-FF) to address this problem and implemented it in the CHARMM simulation package. This Dual-FF method provides the capability to use two separate sets of nonbonded force field parameters within the same simulation: one set to represent intraphase interactions and a separate set to represent interphase interactions. Using this approach, we show that interfacial parameters can be adjusted to correct errors in peptide adsorption free energy without altering peptide conformational behavior in solution. This program thus provides the capability to enable both intraphase and interphase molecular behavior to be accurately and efficiently modeled in the same simulation.

Modified reactive empirical bond-order potential for heterogeneous bonding environments.

An improvement to the AIREBO potential for hydrocarbons is presented in which contributions to the bond order are determined by the local bonding environment around the bond, rather than the average of the environments around the two constituent atoms. This bond-centric approach decreases the errors by ~80% in the fullerene-type systems for which the original approach leads to the most severe errors. With the newly developed and parameterized method, energy errors are less than 0.7 eV for a collection of hydrocarbon molecules not used in the fitting. This modified AIREBO potential is expected to be more useful not only for the molecular hydrocarbons and fullerene isomers studied here, but also for the full range of carbon and hydrocarbon systems to which the AIREBO potential has been applied.

Bond-order potentials with split-charge equilibration: application to C-, H-, and O-containing systems.

A method for extending charge transfer to bond-order potentials, known as the bond-order potential/split-charge equilibration (BOP/SQE) method [P. T. Mikulski, M. T. Knippenberg, and J. A. Harrison, J. Chem. Phys. 131, 241105 (2009)], is integrated into a new bond-order potential for interactions between oxygen, carbon, and hydrogen. This reactive potential utilizes the formalism of the adaptive intermolecular reactive empirical bond-order potential [S. J. Stuart, A. B. Tutein, and J. A. Harrison, J. Chem. Phys. 112, 6472 (2000)] with additional terms for oxygen and charge interactions. This implementation of the reactive potential is able to model chemical reactions where partial charges change in gas- and condensed-phase systems containing oxygen, carbon, and hydrogen. The BOP/SQE method prevents the unrestricted growth of charges, often observed in charge equilibration methods, without adding significant computational time, because it makes use of a quantity which is calculated as part of the underlying covalent portion of the potential, namely, the bond order. The implementation of this method with the qAIREBO potential is designed to provide a tool that can be used to model dynamics in a wide range of systems without significant computational cost. To demonstrate the usefulness and flexibility of this potential, heats of formation for isolated molecules, radial distribution functions of liquids, and energies of oxygenated diamond surfaces are calculated.

Assessment of the transferability of a protein force field for the simulation of peptide-surface interactions.

In order to evaluate the transferability of existing empirical force fields for all-atom molecular simulations of protein adsorption behavior, we have developed and applied a method to calculate the adsorption free energy (DeltaG(ads)) of model peptides on functionalized surfaces for comparison with available experimental data. Simulations were conducted using the CHARMM program and force field using a host-guest peptide with the sequence TGTG-X-GTGT (where G and T are glycine and threonine amino acid residues, respectively, with X representing valine, threonine, aspartic acid, phenylalanine or lysine) over nine different functionalized alkanethiol self-assembled monolayer (SAM) surfaces with explicitly represented solvent. DeltaG(ads) was calculated using biased-energy replica exchange molecular dynamics to adequately sample the conformational states of the system. The simulation results showed that the CHARMM force-field was able to represent DeltaG(ads) within 1 kcal/mol of the experimental values for most systems, while deviations as large as 4 kcal/mol were found for others. In particular, the simulations reveal that CHARMM underestimates the strength of adsorption on the hydrophobic and positively charged amine surfaces. These results clearly show that improvements in force field parameterization are needed in order to accurately represent interactions between amino acid residues and functional groups of a surface and they provide a means for force field evaluation and modification for the eventual development and validation of an interfacial force field for the accurate simulation of protein adsorption behavior.

TIGER2: an improved algorithm for temperature intervals with global exchange of replicas.

An empirical sampling method for molecular simulation based on "temperature intervals with global exchange of replicas" (TIGER2) has been developed to reduce the high demand for computational resources and the low computational efficiency of the conventional replica-exchange molecular dynamics (REMD) method. This new method overcomes the limitation of its previous version, called TIGER, which requires the assumption of constant heat capacity during quenching of replicas from elevated temperatures to the baseline temperature. The robustness of the TIGER2 method is examined by comparing it against a Metropolis Monte Carlo simulation for sampling the conformational distribution of a single butane molecule in vacuum, a REMD simulation for sampling the behavior of alanine dipeptide in explicit solvent, and REMD simulations for sampling the folding behavior of two peptides, (AAQAA)(3) and chignolin, in implicit solvent. The agreement between the results from these conventional sampling methods and the TIGER2 simulations indicates that the TIGER2 algorithm is able to closely approximate a Boltzmann-weighted ensemble of states for these systems but without the limiting assumptions that were required for the original TIGER algorithm. TIGER2 is an efficient replica-exchange sampling method that enables the number of replicas that are used for a replica-exchange simulation to be substantially reduced compared to the conventional REMD method.

Empirical bond-order potential for hydrocarbons: adaptive treatment of van der Waals interactions.

Bond-order potentials provide a powerful class of models for simulating chemically reactive systems with classical potentials. In these models, the covalent bonding interactions adapt to the environment, allowing bond strength to change in response to local chemical changes. However, the non-bonded interactions should also adapt in response to chemical changes, an effect which is neglected in current bond-order potentials. Here the AIREBO potential is extended to include adaptive Lennard-Jones terms, allowing the van der Waals interactions to vary adaptively with the chemical environment. The resulting potential energy surface and its gradient remain continuous, allowing it to be used for dynamics simulations. This new potential is parameterized for hydrocarbons, and is fit to the energetics and densities of a variety of condensed phase molecular hydrocarbons. The resulting model is more accurate for modeling aromatic and other unsaturated hydrocarbon species, for which the original AIREBO potential had some deficiencies. Testing on compounds not used in the fitting procedure shows that the new model performs substantially better in predicting heats of vaporization and pressures (or densities) of condensed-phase molecular hydrocarbons.

Simulation of carbon nanotube welding through Ar bombardment.

Single-walled carbon nanotubes show promise as nanoscale transistors for nanocomputing applications. This use will require appropriate methods for creating electrical connections between distinct nanotubes, analogous to welding of metallic wires at larger length scales, but methods for performing nanoscale chemical welding are not yet sufficiently understood. This study examines the effect of Ar bombardment on the junction of two crossed single-walled carbon nanotubes, to understand the value and limitations of this method for generating connections between nanotubes. A geometric criterion was used to assess the quality of the junctions formed, with the goal of identifying the most productive conditions for experimental ion bombardment. In particular, the effects of nanotube chirality, Ar impact kinetic energy, impact particle flux and fluence, and annealing temperature were considered. The most productive bombardment conditions, leading to the most crosslinking of the tubes with the smallest loss of graphitic (i.e., conductive) character, were found to be at relatively mild impact energies (100 eV), with low flux and high-temperature (3000 K) annealing. Particularly noteworthy for experimental application, a high junction quality is maintained for a relatively broad range of fluences, from 3 × 1019 m-2 to at least 1 × 1020 m-2.

Application of advanced sampling and analysis methods to predict the structure of adsorbed protein on a material surface.

The use of standard molecular dynamics simulation methods to predict the interactions of a protein with a material surface have the inherent limitations of lacking the ability to determine the most likely conformations and orientations of the adsorbed protein on the surface and to determine the level of convergence attained by the simulation. In addition, standard mixing rules are typically applied to combine the nonbonded force field parameters of the solution and solid phases the system to represent interfacial behavior without validation. As a means to circumvent these problems, the authors demonstrate the application of an efficient advanced sampling method (TIGER2A) for the simulation of the adsorption of hen egg-white lysozyme on a crystalline (110) high-density polyethylene surface plane. Simulations are conducted to generate a Boltzmann-weighted ensemble of sampled states using force field parameters that were validated to represent interfacial behavior for this system. The resulting ensembles of sampled states were then analyzed using an in-house-developed cluster analysis method to predict the most probable orientations and conformations of the protein on the surface based on the amount of sampling performed, from which free energy differences between the adsorbed states were able to be calculated. In addition, by conducting two independent sets of TIGER2A simulations combined with cluster analyses, the authors demonstrate a method to estimate the degree of convergence achieved for a given amount of sampling. The results from these simulations demonstrate that these methods enable the most probable orientations and conformations of an adsorbed protein to be predicted and that the use of our validated interfacial force field parameter set provides closer agreement to available experimental results compared to using standard CHARMM force field parameterization to represent molecular behavior at the interface.

Physical adsorption strength in open systems.

For a physical adsorption system, the distances of adsorbates from the surface of a substrate can vary significantly, depending on particle loading and interatomic interactions. Although the total adsorption energy is quantified easily, the normalized, per-particle adsorption energies are more ambiguous if some of these particles are far away from the surface and are interacting only weakly with the substrate. A simple analytical procedure is proposed to characterize the distance dependence of the physisorption strength and effective adsorption capacity. As an example, the method is utilized to describe H2 physisorption in a finite bundle of single-walled carbon nanotubes.

Molecular dynamics simulations on the effects of diameter and chirality on hydrogen adsorption in single walled carbon nanotubes.

We present systematic molecular dynamics simulation studies of hydrogen storage in single walled carbon nanotubes of various diameters and chiralities using a recently developed curvature-dependent force field. Our main objective is to address the following fundamental issues: 1. For a given H2 loading and nanotube type, what is the H2 distribution in the nanotube bundle? 2. For a given nanotube type, what is the maximal loading (H2 coverage)? 3. What is the diameter range and chirality for which H2 adsorption is most energetically favorable? Our simulation results suggest strong dependence of H2 adsorption energies on the nanotube diameter but less dependence on the chirality. Substantial lattice expansion upon H2 adsorption was found. The average adsorption energy increases with the lowering of nanotube diameter (higher curvature) and decreases with higher H2 loading. The calculated H2 vibrational power spectra and radial distribution functions indicate a strong attractive interaction between H2 and nanotube walls. The calculated diffusion coefficients are much higher than what has been reported for H2 in microporous materials such as zeolites, indicating that diffusivity does not present a problem for hydrogen storage in carbon nanotubes.

Parameterization of an interfacial force field for accurate representation of peptide adsorption free energy on high-density polyethylene.

Interfacial force field (IFF) parameters for use with the CHARMM force field have been developed for interactions between peptides and high-density polyethylene (HDPE). Parameterization of the IFF was performed to achieve agreement between experimental and calculated adsorption free energies of small TGTG-X-GTGT host-guest peptides (T = threonine, G = glycine, and X = variable amino-acid residue) on HDPE, with ±0.5 kcal/mol agreement. This IFF parameter set consists of tuned nonbonded parameters (i.e., partial charges and Lennard-Jones parameters) for use with an in-house-modified CHARMM molecular dynamic program that enables the use of an independent set of force field parameters to control molecular behavior at a solid-liquid interface. The R correlation coefficient between the simulated and experimental peptide adsorption free energies increased from 0.00 for the standard CHARMM force field parameters to 0.88 for the tuned IFF parameters. Subsequent studies are planned to apply the tuned IFF parameter set for the simulation of protein adsorption behavior on an HDPE surface for comparison with experimental values of adsorbed protein orientation and conformation.

A molecular modeling study of the effect of surface chemistry on the adsorption of a fibronectin fragment spanning the 7-10th type III repeats.

Although it is well documented that proteins adsorb onto biomaterial surfaces, relatively little is quantitatively understood about the effects of adsorption on protein orientation and conformation. Because this is the primary determining factor of protein bioactivity, the ability to accurately predict a protein's orientation and conformation following adsorption will be essential for the rational design of biomaterial surfaces to control biological responses. Force field-based computational chemistry methods provide an excellent means to theoretically address this issue, with the nontrivial requirement that the force field must be tailored to appropriately represent protein adsorption behavior. Accordingly, we have modified an existing force field (CHARMm) based on semiempirical quantum-mechanical peptide adsorption data to enable it to simulate protein adsorption behavior in an implicit aqueous environment. This modified force field was then applied to predict the adsorption behavior of the 7-10 type III repeats of fibronectin on functionalized surfaces. Predicted changes in adsorption energy and adsorption-induced conformation as a function of surface chemistry were found to correlate well with experimentally observed trends for these same systems. This work represents a first attempt towards the development of a molecular mechanics force field that is specifically parameterized to accurately simulate protein adsorption to biomaterial surfaces.

Polarizable molecular dynamics simulations of aqueous dipeptides.

Molecular dynamics simulations were carried out for concentrated aqueous solutions of three dipeptides: Gly-Ala, Gly-Pro, and Ala-Pro. The simulations were performed using both polarizable and nonpolarizable force fields, as a method of assessing the effects of polarization in a well-characterized biomolecular system, and to determine whether the models are adequate to reproduce observed aggregation behavior. The structure and dynamics of both solute and solvent were analyzed and the results compared to experiment, including neutron diffraction measurements of peptide aggregation. The polarizable water is depolarized in concentrated peptide solutions, reflecting its ability to adapt to heterogeneous electrostatic environments. Significant differences between the polarizable and nonpolarizable models are found in terms of both the structure and the dynamics of water as a solvent. Although the water shows more realistic structure and dynamics in the polarizable simulations, consistent with enhanced peptide-water interaction, the peptide aggregation behavior agrees less well with the experiment. Neither model successfully reproduces the experimentally observed dipeptide aggregation behavior.

Comparison between empirical protein force fields for the simulation of the adsorption behavior of structured LK peptides on functionalized surfaces.

All-atom empirical molecular mechanics protein force fields, which have been developed to represent the energetics of peptide folding behavior in aqueous solution, have not been parameterized for protein interactions with solid material surfaces. As a result, their applicability for representing the adsorption behavior of proteins with functionalized material surfaces should not be assumed. To address this issue, we conducted replica-exchange molecular dynamics simulations of the adsorption behavior of structured peptides to functionalized surfaces using three protein force fields that are widely used for the simulation of peptide adsorption behavior: CHARMM22, AMBER94, and OPLS-AA. Simulation results for peptide structure both in solution and when adsorbed to the surfaces were compared to experimental results for similar peptide-surface systems to provide a means of evaluating and comparing the performance of these three force fields for this type of application. Substantial differences in both solution and adsorbed peptide conformations were found amongst these three force fields, with the CHARMM22 force field found to most closely match experimental results.

Development of a tuned interfacial force field parameter set for the simulation of protein adsorption to silica glass.

Adsorption free energies for eight host-guest peptides (TGTG-X-GTGT, with X = N, D, G, K, F, T, W, and V) on two different silica surfaces [quartz (100) and silica glass] were calculated using umbrella sampling and replica exchange molecular dynamics and compared with experimental values determined by atomic force microscopy. Using the CHARMM force field, adsorption free energies were found to be overestimated (i.e., too strongly adsorbing) by about 5-9 kcal/mol compared to the experimental data for both types of silica surfaces. Peptide adsorption behavior for the silica glass surface was then adjusted using a modified version of the CHARMM program, which we call dual force-field CHARMM, which allows separate sets of nonbonded parameters (i.e., partial charge and Lennard-Jones parameters) to be used to represent intra-phase and inter-phase interactions within a given molecular system. Using this program, interfacial force field (IFF) parameters for the peptide-silica glass systems were corrected to obtain adsorption free energies within about 0.5 kcal/mol of their respective experimental values, while IFF tuning for the quartz (100) surface remains for future work. The tuned IFF parameter set for silica glass will subsequently be used for simulations of protein adsorption behavior on silica glass with greater confidence in the balance between relative adsorption affinities of amino acid residues and the aqueous solution for the silica glass surface.

Development of molecular simulation methods to accurately represent protein-surface interactions: The effect of pressure and its determination for a system with constrained atoms.

When performing molecular dynamics simulations for a system with constrained (fixed) atoms, traditional isobaric algorithms (e.g., NPT simulation) often cannot be used. In addition, the calculation of the internal pressure of a system with fixed atoms may be highly inaccurate due to the nonphysical nature of the atomic constraints and difficulties in accurately defining the volume occupied by the unconstrained atoms in the system. The inability to properly set and control pressure can result in substantial problems for the accurate simulation of condensed-phase systems if the behavior of the system (e.g., peptide/protein adsorption) is sensitive to pressure. To address this issue, the authors have developed an approach to accurately determine the internal pressure for a system with constrained atoms. As the first step in this method, a periodically extendable portion of the mobile phase of the constrained system (e.g., the solvent atoms) is used to create a separate unconstrained system for which the pressure can be accurately calculated. This model system is then used to create a pressure calibration plot for an intensive local effective virial parameter for a small volume cross section or "slab" of the system. Using this calibration plot, the pressure of the constrained system can then be determined by calculating the virial parameter for a similarly sized slab of mobile atoms. In this article, the authors present the development of this method and demonstrate its application using the CHARMM molecular simulation program to characterize the adsorption behavior of a peptide in explicit water on a hydrophobic surface whose lattice spacing is maintained with atomic constraints. The free energy of adsorption for this system is shown to be dramatically influenced by pressure, thus emphasizing the importance of properly maintaining the pressure of the system for the accurate simulation of protein-surface interactions.

Development of molecular simulation methods to accurately represent protein-surface interactions: Method assessment for the calculation of electrostatic effects.

The simulation of the interactions of proteins with charged surfaces in a condensed-phase aqueous solution containing electrolytes using empirical force field based methods is predominantly governed by nonbonded interactions between the atoms of the protein, surface, and the solvent. Electrostatic effects represent the strongest type of these interactions and the type that is most difficult to accurately represent because of their long-range influence. While many different methods have been developed to represent electrostatic interactions, the particle mesh Ewald summation (PME) method is generally considered to be the most accurate one for calculating these effects. However, the PME method was designed for systems with three-dimensional (3D) periodicity, and not for interfacial systems such as the case of protein adsorption to a charged surface. Interfacial systems such as these have only two-dimensional periodicity, which may not be appropriate for treatment with PME due to the possibility that the presence of multiple charged image surfaces parallel to the primary simulation cell's surface, may introduce nonphysical effects on the behavior of the charged molecules in the system. In an effort to address this issue, the authors have conducted a set of nanosecond-scale molecular dynamics simulations to calculate the equilibrium distribution of Na(+) and Cl(-) ions near a charged surface using PME and a range of radial cutoff methods for treating electrostatic interactions, where the cutoffs prevent interaction with the periodic images of the system. The resulting ion concentration profiles were compared to one another and to a continuum analytical solution of the theoretical ion distribution obtained from the Poisson-Boltzmann equation. Their results show that the PME method does not introduce the suspected nonphysical effects in the ion distributions due to the 3D periodic images of the system, thus indicating that it is appropriate for use for this type of molecular simulation. Although their interest is motivated by protein-surface interactions, the conclusions are applicable for the treatment of electrostatics in other aqueous systems with two-dimensional periodicity.