Paramfit: Automated optimization of force field parameters for molecular dynamics simulations

The generation of bond, angle, and torsion parameters for classical molecular dynamics force fields typically requires fitting parameters such that classical properties such as energies and gradients match precalculated quantum data for structures that scan the value of interest. We present a program, Paramfit, distributed as part of the AmberTools software package that automates and extends this fitting process, allowing for simplified parameter generation for applications ranging from single molecules to entire force fields. Paramfit implements a novel combination of a genetic and simplex algorithm to find the optimal set of parameters that replicate either quantum energy or force data. The program allows for the derivation of multiple parameters simultaneously using significantly fewer quantum calculations than previous methods, and can also fit parameters across multiple molecules with applications to force field development. Paramfit has been applied successfully to systems with a sparse number of structures, and has already proven crucial in the development of the Assisted Model Building with Energy Refinement Lipid14 force field. © 2014 Wiley Periodicals, Inc.     

Investigation of chiral recognition by molecular micelles with molecular dynamics simulations

Molecular dynamics simulations were used to characterize the binding of the chiral drugs chlorthalidone and lorazepam to the molecular micelle poly-(sodium undecyl-(L)-leucine-valine). The project’s goal was to characterize the nature of chiral recognition in capillary electrophoresis separations that use molecular micelles as the chiral selector. The shapes and charge distributions of the chiral molecules investigated, their orientations within the molecular micelle chiral binding pockets, and the formation of stereoselective intermolecular hydrogen bonds with the molecular micelle were all found to play key roles in determining where and how lorazepam and chlorthalidone enantiomers interacted with the molecular micelle.     

Force-field parametrization and molecular dynamics simulations of Congo red

Congo red, a diazo dye widely used in medical diagnosis, is known to form supramolecular systems in solution. Such a supramolecular system may interact with various proteins. In order to examine the nature of such complexes empirical force field parameters for the Congo red molecule were developed. The parametrization of bonding terms closely followed the methodology used in the development of the charmm22 force field, except for the calculation of charges. Point charges were calculated from a fit to a quantum mechanically derived electrostatic potential using the CHELP-BOW method. Obtained parameters were tested in a series of molecular dynamics simulations of both a single molecule and a micelle composed of Congo red molecules. It is shown that newly developed parameters define a stable minimum on the hypersurface of the potential energy and crystal and ab initio geometries and rotational barriers are well reproduced. Furthermore, rotations around C-N bonds are similar to torsional vibrations observed in crystals of diphenyl-diazene, which confirms that the flexibility of the molecule is correct. Comparison of results obtained from micelles molecular dynamics simulations with experimental data shows that the thermal dependence of micelle creation is well reproduced.     

Steered molecular dynamics simulations of protein-ligand interactions

Molecular recognition and specific protein-ligandinteractions are central to many biochemical processes,such as enzyme catalysis, assembly of organelles, en-ergy transduction, signaling, diverse control functions,and replication, expression and storage of the geneticmaterial[1]. Moreover, protein-ligand interactions pro-vide the mechanism of many drug therapies and un-derstanding of such interactions is thus significant forrational drug design[1,2]. For the experimental studiesof protein-ligan…     

Free energy calculations using flexible-constrained, hard-constrained and non-constrained molecular dynamics simulations.

A comparison of different treatments of bond-stretching interactions in molecular dynamics simulation is presented. Relative free energies from simulations using rigid bonds maintained with the SHAKE algorithm, using partially rigid bonds maintained with a recently introduced flexible constraints algorithm, and using fully flexible bonds are compared in a multi-configurational thermodynamic integration calculation of changing liquid water into liquid methanol. The formula for the free energy change due to a changing flexible constraint in a flexible constraint simulation is derived. To allow for a more direct comparison between these three methods, three different pairs of models for water and methanol were used: a flexible model (simulated without constraints and with flexible constraints), a rigid model (simulated with standard hard constraints), and an alternative flexible model (simulated with flexible constraints and standard hard constraints) in which the ideal or constrained bond lengths correspond to the average bond lengths obtained from a short simulation of the unconstrained flexible model. The particular treatment of the bonds induces differences of up to 2 % in the liquid densities, whereas (excess) free energy differences of up to 5.7 (4.3) kJ mol(-1) are observed. These values are smaller than the differences observed between the three different pairs of methanol/water models: up to 5 % in density and up to 8.5 kJ mol(-1) in (excess) free energy.     

Finite element interpolation for combined classical/quantum mechanical molecular dynamics simulations

A method is presented to interpolate the potential energy function for a part of a system consisting of a few degrees of freedom, such as a molecule in solution. The method is based on a modified finite element (FE) interpolation scheme. The aim is to save computer time when expensive methods such as quantum-chemical calculations are used to determine the potential energy function. The expensive calculations are only carried out if the molecule explores new unknown regions of the conformation space. If the molecule resides in regions previously explored, a cheap interpolation is performed instead of an expensive calculation, using known neighboring points. We report the interpolation techniques for the energies and the forces of the molecule, the handling of the FE mesh, and an application to a simple test example in molecular dynamics (MD) simulations. Good performance of the method was obtained (especially for MD simulations with a preceding Monte Carlo mesh generation) without losing accuracy. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1484–1495, 1997     

Ramachandran-type plots for glycosidic linkages: Examples from molecular dynamic simulations using the Glycam06 force field

The goals of this article are to (1) provide further validation of the Glycam06 force field, specifically for its use in implicit solvent molecular dynamic (MD) simulations, and (2) to present the extension of G.N. Ramachandran's idea of plotting amino acid phi and psi angles to the glycosidic phi, psi, and omega angles formed between carbohydrates. As in traditional Ramachandran plots, these carbohydrate Ramachandran-type (carb-Rama) plots reveal the coupling between the glycosidic angles by displaying the allowed and disallowed conformational space. Considering two-bond glycosidic linkages, there are 18 possible conformational regions that can be defined by (alpha, phi, psi) and (beta, phi, psi), whereas for three-bond linkages, there are 54 possible regions that can be defined by (alpha, phi, psi, omega) and (beta, phi, psi, omega). Illustrating these ideas are molecular dynamic simulations on an implicitly hydrated oligosaccharide (700 ns) and its eight constituent disaccharides (50 ns/disaccharide). For each linkage, we compare and contrast the oligosaccharide and respective disaccharide carb-Rama plots, validate the simulations and the Glycam06 force field through comparison to experimental data, and discuss the general trends observed in the plots.     

Force-field development and molecular dynamics simulations of ferrocene-peptide conjugates as a scaffold for hydrogenase mimics

The increasing importance of hydrogenase enzymes in the new energy research field has led us to examine the structure and dynamics of potential hydrogenase mimics, based on a ferrocene-peptide scaffold, using molecular dynamics (MD) simulations. To enable this MD study, a molecular mechanics force field for ferrocene-bearing peptides was developed and implemented in the CHARMM simulation package, thus extending the usefulness of the package into peptide-bioorganometallic chemistry. Using the automated frequency-matching method (AFMM), optimized intramolecular force-field parameters were generated through quantum chemical reference normal modes. The partial charges for ferrocene were derived by fitting point charges to quantum-chemically computed electrostatic potentials. The force field was tested against experimental X-ray crystal structures of dipeptide derivatives of ferrocene-1,1'-dicarboxylic acid. The calculations reproduce accurately the molecular geometries, including the characteristic C2-symmetrical intramolecular hydrogen-bonding pattern, that were stable over 0.1 micros MD simulations. The crystal packing properties of ferrocene-1-(D)alanine-(D)proline-1'-(D)alanine-(D)proline were also accurately reproduced. The lattice parameters of this crystal were conserved during a 0.1 micros MD simulation and match the experimental values almost exactly. Simulations of the peptides in dichloromethane are also in good agreement with experimental NMR and circular dichroism (CD) data in solution. The developed force field was used to perform MD simulations on novel, as yet unsynthesized peptide fragments that surround the active site of [Ni-Fe] hydrogenase. The results of this simulation lead us to propose an improved design for synthetic peptide-based hydrogenase models. The presented MD simulation results of metallocenes thereby provide a convincing validation of our proposal to use ferrocene-peptides as minimal enzyme mimics.     

Molecular dynamics simulations for the growth of CH4-CO2 mixed hydrate

Molecular dynamics simulations are performed to study the growth mechanism of CH4-CO2 mixed hydrate in xco2 = 75%, xco2 = 50%, and zco2 = 25% systems at T = 250 K, 255 K and 260 K, respectively. Our simulation results show that the growth rate of CH4-CO2 mixed hydrate increases as the CO2 concentration in the initial solution phase increases and the temperature decreases. Via hydrate formation, the composition of CO2 in hydrate phase is higher than that in initial solution phase and the encaging capacity of CO2 in hydrates increases with the decrease in temperature. By analysis of the cage occupancy ratio of CH4 molecules and CO2 molecules in large cages to small cages, we find that CO2 molecules are preferably encaged into the large cages of the hydrate crystal as compared with CH4 molecules. Interestingly, CH4 molecules and CO2 molecules frequently replace with each other in some particular cage sites adjacent to hydrate/solution interface during the crystal growth process. These two species of guest molecules eventually act to stabilize the newly formed hydrates, with CO2 molecules occupying large cages and CH4 molecules occupying small cages in hydrate.     

Effect of acetone accumulation on structure and dynamics of lipid membranes studied by molecular dynamics simulations

The modulation of the properties and function of cell membranes by small volatile substances is important for many biomedical applications. Despite available experimental results, molecular mechanisms of action of inhalants and organic solvents, such as acetone, on lipid membranes remain not well understood. To gain a better understanding of how acetone interacts with membranes, we have performed a series of molecular dynamics (MD) simulations of a POPC bilayer in aqueous solution in the presence of acetone, whose concentration was varied from 2.8 to 11.2 mol%. The MD simulations of passive distribution of acetone between a bulk water phase and a lipid bilayer show that acetone favors partitioning into the water-free region of the bilayer, located near the carbonyl groups of the phospholipids and at the beginning of the hydrocarbon core of the lipid membrane. Using MD umbrella sampling, we found that the permeability barrier of ∼0.5 kcal/mol exists for acetone partitioning into the membrane. In addition, a Gibbs free energy profile of the acetone penetration across a bilayer demonstrates a favorable potential energy well of −3.6 kcal/mol, located at 15–16 Å from the bilayer center. The analysis of the structural and dynamics properties of the model membrane revealed that the POPC bilayer can tolerate the presence of acetone in the concentration range of 2.8–5.6 mol%. The accumulation of the higher acetone concentration of 11.2 mol% results, however, in drastic disordering of phospholipid packing and the increase in the membrane fluidity. The acetone molecules push the lipid heads apart and, hence, act as spacers in the headgroup region. This effect leads to the increase in the average headgroup area per molecule. In addition, the acyl tail region of the membrane also becomes less dense. We suggest, therefore, that the molecular mechanism of acetone action on the phospholipid bilayer has many common features with the effects of short chain alcohols, DMSO, and chloroform.     

The Gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations

A new algorithm for density-functional-theory-based ab initio molecular dynamics simulations is presented. The Kohn–Sham orbitals are expanded in Gaussian-type functions and an augmented-plane-wave-type approach is used to represent the electronic density. This extends previous work of ours where the density was expanded only in plane waves. We describe the total density in a smooth extended part which we represent in plane waves as in our previous work and parts localised close to the nuclei which are expanded in Gaussians. Using this representation of the charge we show how the localised and extended part can be treated separately, achieving a computational cost for the calculation of the Kohn–Sham matrix that scales with the system size N as O(NlogN). Furthermore, we are able to reduce drastically the size of the plane-wave basis. In addition, we introduce a multiple-cutoff method that improves considerably the performance of this approach. Finally, we demonstrate with a series of numerical examples the accuracy and efficiency of the new algorithm, both for electronic structure calculations and for ab initio molecular dynamics simulations. Received: 15 December 1998 /Accepted: 18 February 1999 /Published online: 14 July 1999     

Using molecular dynamics simulations to identify the key factors responsible for chiral recognition by an amino acid-based molecular micelle

Molecular dynamics (MD) simulations were used to investigate the binding of six chiral compounds to the amino acid-based molecular micelle (MM) poly-(sodium undecyl-(L)-leucine-leucine) or poly(SULL). The MM investigated is used as a chiral selector in capillary electrophoresis. The project goal was to characterize the chiral recognition mechanism in these separations and to move toward predictive models to identify the best amino acid-based MM for a given separation. Poly(SULL) was found to contain six binding sites into which chiral compounds could insert. Four sites had similar sizes, shapes, and electrostatic properties. Enantiomers of alprenolol, propranolol, 1,1′-bi-2-naphthyl-2,2′-diyl hydrogen phosphate, 1,1′-bi-2-naphthol, chlorthalidone, or lorazepam were separately docked into each binding pocket and MD simulations with the resulting intermolecular complexes were performed. Solvent-accessible surface area calculations showed the compounds preferentially associated with binding sites where they penetrated into the MM core and shielded their non-polar atoms from solvent. Furthermore, with five of the six compounds the enantiomer with the most favorable free energy of MM association also experienced the most favorable intermolecular hydrogen bonding interactions with the MM. This result suggests that stereoselective intermolecular hydrogen bonds play an important role in chiral discrimination in separations using amino acid-based MMs.GRAPHICAL ABSTRACT     

Molecular dynamics method in studies of molecular film growth processes

This review is an attempt to analyze some of the experimental problems arising in the course of growth of orientated molecular films using metal monophthalocyanine (MPc) films as an example and to demonstrate the possibilities of molecular dynamics simulation of these processes for solving experimental problems. Examples of theoretical simulation of adsorption processes are given; formation of a molecular monolayer is considered for copper phthalocyanine films as an example.     

Development of a parallel molecular dynamics code on SIMD computers: Algorithm for use of pair list criterion

In recent years several implementations of molecular dynamics (MD) codes have been reported on multiple instruction multiple data (MIMD) machines. However, very few implementations of MD codes on single instruction multiple data (SIMD) machines have been reported. The difficulty in using pair lists of nonbonded interactions is the major problem with MD codes for SIMD machines, such that, generally, the full connectivity computation has been used. We present an algorithm, the global cut-off algorithm (GCA), which permits the use of pair lists on SIMD machines. GCA is based on a probabilistic approach and requires the cut-off condition to be simultaneously verified on all nodes of the machine. The MD code used was taken from the GROMOS package; only the routines involved in the pair lists and in the computation of nonbonded interactions were rewritten for a parallel architecture. The remaining calculations were performed on the host computer. The algorithm has been tested on Quadrics computers for configurations of 32, 128, and 512 processors and for systems of 4000, 8000, 15,000, and 30,000 particles. Quadrics was developed by Istituto Nazionale di Fisica Nucleare (INFN) and marketed by Alenia Spazio. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 685–694, 1998     

Evaluating substrate specificity of glutathione peroxidase mimic by molecular dynamics simulations and kinetics

The substrate specificities of glutathione peroxidase (GPX) mimic, 6,6′-ditellurobis(6-deoxy-β-cyclodextrin) (6-TeCD), for three hydroperoxides (ROOH), H₂O₂, tert-butyl hydroperoxide (t-BuOOH) and cumene hydroperoxide (CuOOH), are investigated through molecular dynamics (MD) simulations. The most stable conformations and the total interaction energies of complex of 6-TeCD with ROOH are used to evaluate the substrate specificity of 6-TeCD. The steady-state kinetics of 6-TeCD is studied and the Michaelis-Menten constant (K _m) and second-order rate constant k _max/K _ROOH show that 6-TeCD displays different affinity and specificity to ROOH. These results of experiments are well consistent with ones obtained by MD simulations, indicating that MD simulations could be applied to evaluation substrate specificity of small-molecule enzyme mimics.     

LINCS: A linear constraint solver for molecular simulations

In this article, we present a new LINear Constraint Solver (LINCS) for molecular simulations with bond constraints. The algorithm is inherently stable, as the constraints themselves are reset instead of derivatives of the constraints, thereby eliminating drift. Although the derivation of the algorithm is presented in terms of matrices, no matrix matrix multiplications are needed and only the nonzero matrix elements have to be stored, making the method useful for very large molecules. At the same accuracy, the LINCS algorithm is three to four times faster than the SHAKE algorithm. Parallelization of the algorithm is straightforward. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1463–1472, 1997     

Nuclear quantum effects at aqueous metal interfaces captured by molecular dynamics simulations

In this review, we summarize the recent development in modeling nuclear quantum effects at aqueous metal interfaces. First, we review the nuclear quantum effects on the water-metal interface at ultrahigh vacuum. Then, we illustrate the nuclear quantum effects at the potential of zero charge conditions. At last, we give some outlook for the perspective work in modeling the nuclear quantum effects at electrochemical interfaces and some practical simulation strategies.