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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

Molecular dynamics simulations on hydrogen adsorption in finite single walled carbon nanotube bundles.

Molecular dynamics simulations of the adsorption of hydrogen molecules in finite single-walled carbon nanotube bundles are presented using a curvature dependent force field. The heat of formation and the effective adsorption capacity are expressed as a function of H(2) distance from adsorbent. The heat of adsorption decreases rapidly with the distance and increasing H(2) loading results in weakening adsorption strength. The effects of nanotube packing and bundle thickness on hydrogen adsorption strength were investigated and the results show that the heat of adsorption can be improved slightly if hydrogen molecules are placed in thicker and inhomogeneously packed nanotube bundles. Only very small diameter nanotube bundles were found to hold promise for significant hydrogen storage for onboard applications.

Calculation of adsorption free energy for solute-surface interactions using biased replica-exchange molecular dynamics.

The adsorption behavior of a biomolecule, such as a peptide or protein, to a functionalized surface is of fundamental importance for a broad range of applications in biotechnology. The adsorption free energy for these types of interactions can be determined from a molecular dynamics simulation using the partitioning between adsorbed and nonadsorbed states, provided that sufficient sampling of both states is obtained. However, if interactions between the solute and the surface are strong, the solute will tend to be trapped near the surface during the simulation, thus preventing the adsorption free energy from being calculated by this method. This situation occurs even when using an advanced sampling algorithm such as replica-exchange molecular dynamics (REMD). In this paper, the authors demonstrate the fundamental basis of this problem using a model system consisting of one sodium ion (Na(+)) as the solute positioned over a surface functionalized with one negatively charged group (COO(-)) in explicit water. With this simple system, the authors show that sufficient sampling in the coordinate normal to the surface cannot be obtained by conventional REMD alone. The authors then present a method to overcome this problem through the use of an adaptive windowed-umbrella sampling technique to develop a biased-energy function that is combined with REMD. This approach provides an effective method for the calculation of adsorption free energy for solute-surface interactions.

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.

An improved replica-exchange sampling method: temperature intervals with global energy reassignment.

In a molecular dynamics (MD) simulation, representative sampling over the entire phase space is desired to obtain an accurate canonical distribution at a given temperature. For large molecules, such as proteins, this is problematic because systems tend to become trapped in local energy minima. The extensively used replica-exchange molecular dynamics (REMD) simulation technique overcomes this kinetic-trapping problem by allowing Boltzmann-weighted configuration exchange processes to occur between numerous thermally adjacent and compositionally identical simulations that are thermostated at sequentially higher temperatures. While the REMD method provides much better sampling than conventional MD, there are two substantial difficulties that are inherent in its application: (1) the large number of replicas that must be used to span a designated temperature range and (2) the subsequent long time required for configurations sampled at high temperatures to exchange down for potential inclusion within the low-temperature ensemble of interest. In this work, a new method based on temperature intervals with global energy reassignment (TIGER) is presented that overcomes both of these problems. A TIGER simulation is conducted as a series of short heating-sampling-quenching cycles. At the end of each cycle, the potential energies of all replicas are simultaneously compared at the same temperature using a Metropolis sampling method and then globally reassigned to the designated temperature levels. TIGER is compared with regular MD and REMD methods for the alanine dipeptide in water. The results indicate that TIGER increases sampling efficiency while substantially reducing the number of central processing units required for a comparable conventional REMD simulation.