Effective force fields for condensed phase systems from ab initio molecular dynamics simulation: a new method for force-matching

A novel least-squares fitting approach is presented to obtain classical force fields from trajectory and force databases produced by ab initio (e.g., Car-Parrinello) molecular dynamics (MD) simulations. The method was applied to derive effective nonpolarizable three-site force fields for liquid water at ambient conditions from Car-Parrinello MD simulations in the Becke-Lee-Yang-Parr approximation to the electronic density functional theory. The force-matching procedure includes a fit of short-ranged nonbonded forces, bonded forces, and atomic partial charges. The various parameterizations of the water force field differ by an enforced smooth cut-off applied to the short-ranged interaction term. These were obtained by fitting to the trajectory and force data produced by Car-Parrinello MD simulations of systems of 32 and 64 H(2)O molecules. The new water force fields were developed assuming both flexible or rigid molecular geometry. The simulated structural and self-diffusion properties of liquid water using the fitted force fields are in close agreement with those observed in the underlying Car-Parrinello MD simulations. The resulting empirical models compare to experiment much better than many conventional simple point charge (SPC) models. The fitted potential is also shown to combine well with more sophisticated intramolecular potentials. Importantly, the computational cost of the new models is comparable to that for SPC-like potentials.     

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.     

A study of pore blockage in silicalite zeolite using free energy perturbation calculations

Binary systems consisting of large coadsorbed molecules (n-hexane, cyclohexane, and benzene) with smaller penetrant molecules (methane) were simulated to investigate the mechanisms of pore blockage in the zeolite silicalite. Benzene and cyclohexane trap the methane molecules in the zeolite channels on the time scales of molecular dynamics simulations. Minimum energy paths for methane diffusion past the blocking molecules were determined, and free energy perturbation calculations were carried out along the paths to get the rate constants of methane hopping past coadsorbed benzene and cyclohexane molecules, which adsorb in the channel intersections. Three principal diffusion pathways were found in both the methane/benzene and methane/cyclohexane systems. Minima which were connected by low-energy pathways were grouped together into macrostates. Using the calculated hopping rates between macrostates, kinetic Monte Carlo was then used to obtain the diffusivity of methane with a coadsorbate benzene loading such that all channel intersections are filled by benzene - conditions where molecular dynamics simulations fail. Passage of methane across cyclohexane molecules involved pushing the cyclohexane molecules into the channels from their preferred channel intersection positions.     

Molecular dynamics and X-ray scattering simulations of cyclic siloxane-based liquid crystal mesogens

Molecular dynamics simulations of cyclic siloxane-based liquid crystals offer new insights into the conformational flexibility of these materials. Interdigitation between the cholesteryl-4'-allyloxybenzoate and biphenyl-4'-allyloxybenzoate mesogens pendant on the cyclic siloxane ring is observed in the simulated structures. All molecular models considered viz. disc, cone, and cylinder, display a large conformational flexibility, which is important regarding the liquid crystalline phase behavior. The disc molecular model exhibits the largest flexibility as indicated by mean dihedral angles and their range for certain principal torsions, evaluated from the molecular dynamics simulations. Results from the dynamics simulations of cylinder molecular pairs indicate a significant amount of conformational flexibility in the siloxane rings. The degree of interdigitation between mesogens is dependent on the flexibility of the siloxane rings, as shown by calculations for a fixed ring system resulting in less interdigitation, also reflected in calculated X-ray scattering sections along the starting molecular direction. Weaker molecular transforms were observed for the non-fixed system due to a lack of boundary conditions. In general, the qualitative agreement of the starting structure's reflections and those shown by the experimental data is encouraging.     

Quasi-Hamiltonian equations of motion for internal coordinate molecular dynamics of polymers

Conventional molecular dynamics simulations of macromolecules require long computational times because the most interesting motions are very slow compared to the fast oscillations of bond lengths and bond angles that limit the integration time step. Simulation of dynamics in the space of internal coordinates, that is, with bond lengths, bond angles, and torsions as independent variables, gives a theoretical possibility of eliminating all uninteresting fast degrees of freedom from the system. This article presents a new method for internal coordinate molecular dynamics simulations of macromolecules. Equations of motion are derived that are applicable to branched chain molecules with any number of internal degrees of freedom. Equations use the canonical variables and they are much simpler than existing analogs. In the numerical tests the internal coordinate dynamics are compared with the traditional Cartesian coordinate molecular dynamics in simulations of a 56 residue globular protein. For the first time it was possible to compare the two alternative methods on identical molecular models in conventional quality tests. It is shown that the traditional and internal coordinate dynamics require the same time step size for the same accuracy and that in the standard geometry approximation of amino acids, that is, with fixed bond lengths, bond angles, and rigid aromatic groups, the characteristic step size is 4 fs, which is 2 times higher than with fixed bond lengths only. The step size can be increased up to 11 fs when rotation of hydrogen atoms is suppressed. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18 : 1354–1364, 1997     

A configurational temperature for molecules with hard-core or discontinuous interactions

We extend the usual formula for a configurational temperature so that it applies to condensates in which the molecules interact through hard-core or discontinuous potentials. The new formula involves extra terms which may be calculated during the course of a simulation. The formula is tested by its application to a number of systems with discontinuous or hard-core potentials in thermodynamic equilibrium. Metropolis Monte Carlo simulations were performed on these systems in a canonical ensemble and the configurational temperature is compared with the input temperature. The two are in agreement to within less than 0.1%.     

On using the NMR chemical shift to assess polar–nonpolar cross-intermolecular interactions

The possibility of using the NMR chemical shift to evaluate and develop intermolecular potentials for cross-interactions between polar and nonpolar molecules has been examined using the method of molecular dynamics. Such interaction potential models are known to be notoriously difficult to develop. Our work has shown that chemical shift can be obtained quite efficiently in simulations and converges much faster than other properties traditionally used for such evaluations (for example, the infinite dilution activity coefficients, Henry’s constants or the solubility of solutes in solvents). We have also found chemical shift to be quite sensitive to the intermolecular potentials which makes it a rather promising property to investigate polar–nonpolar interactions in fluids.     

First principles and classical molecular dynamics simulations of solvated benzene

We have performed extensive ab initio and classical molecular dynamics (MD) simulations of benzene in water in order to examine the unique solvation structures that are formed. Qualitative differences between classical and ab initio MD simulations are found and the importance of various technical simulation parameters is examined. Our comparison indicates that nonpolarizable classical models are not capable of describing the solute-water interface correctly if local interactions become energetically comparable to water hydrogen bonds. In addition, a comparison is made between a rigid water model and fully flexible water within ab initio MD simulations which shows that both models agree qualitatively for this challenging system.     

Fluctuating charge normal modes: An algorithm for implementing molecular dynamics simulations with polarizable potentials

A new method for performing molecular dynamics simulations with fluctuating charge polarizable potentials is introduced. In fluctuating charge models, polarizability is treated by allowing the partial charges to be variables, with values that are coupled to charges on the same molecule as well as those on other molecules. The charges can be efficiently propagated in a molecular dynamics simulation using extended Lagrangian dynamics. By making a coordinate change from the charge variables to a set of normal mode charge coordinates for each molecule, a new method is constructed in which the normal mode charge variables uncouple from those on the same molecule. The method is applied to the TIP4P-FQ model of water and compared to other methods for implementing the dynamics. The methods are compared using different molecular dynamics time steps.     

Consequences of Strain for the Structure of Aliphatic Molecules

The chemist is accustomed to deriving structures and preferred conformations of organic compounds from rigid molecular models and standard values for bond lengths, bond angles, and torsional profiles. In the case of strained compounds, this rigid structural model has to be abandoned and replaced by a flexible one which takes individual conditions of strain into consideration. It is shown, on the basis of new experimental structure data, that the force field method is suitable and highly reliable for the calculation of structural parameters and preferred conformations of strained compounds. It is, therefore, capable of replacing the rigid molecular model. Furthermore, the systematic analysis of strain induced angle and bond deformation gives a new pivot for the development of a qualitative discussion of deformation in strained molecules and hence for improved conformational analysis. — In the course of this work we were able to isolate two rotamers of D,L -3,4-di(1-adamantyl)-2,2,5,5-tetramethylhexane; this is the first isolation of a rotamer pair of an aliphatic hydrocarbon.     

Simulation of liquid water using semiempirical Hamiltonians and the divide and conquer approach

This work examines the ability of semiempirical methods to describe the structure of liquid water. Particularly, the standard AM1 and PM3 methods together with recently developed PM3-PIF and PM3-MAIS parametrizations have been considered. We perform molecular dynamics simulations for a system consisting of 64 or 216 water molecules in a periodic cubic box. The whole system is described quantum mechanically. Calculations with 64 molecules have been carried out using standard SCF techniques whereas calculations with 216 molecules have been done using the divide and conquer approach. This method has also been used in one simulation case with 64 molecules for test purposes. Within this scope, the molecular dynamics program ROAR have been coupled with a linear scaling semiempirical code (DivCon) implemented in a parallel MPI version. The predicted liquid water structure using either AM1 or PM3 is shown to be very poor due to well-known limitations of these methods describing hydrogen bonds. In contrast, PM3-PIF and PM3-MAIS calculations lead to results in reasonably good agreement with experimental data. The best results for the heat of vaporization are obtained with the PM3-PIF method. The average induced dipole moment of the water molecule in the liquid is underestimated by all semiempirical techniques, which seems to be related to the NDDO approximation and to the use of minimal basis sets. A brief discussion on charge-transfer effects in liquid water is also presented.     

Coarse-grained molecular dynamics simulations of the sphere to rod transition in surfactant micelles

Surfactant molecules self-assemble in aqueous solutions to form various micellar structures such as spheres, rods, or lamellae. Although phase transitions in surfactant solutions have been studied experimentally, their molecular mechanisms are still not well understood. In this work, we show that molecular dynamics (MD) simulations using the coarse-grained (CG) MARTINI force field and explicit CG solvent, validated against atomistic MD studies, can accurately represent micellar assemblies of cetyltrimethylammonium chloride (CTAC). The effect of salt on micellar structures is studied for aromatic anionic salts, e.g., sodium salicylate, and simple inorganic salts, e.g., sodium chloride. Above a threshold concentration, sodium salicylate induces a sphere to rod transition in the micelle. CG MD simulations are shown to capture the dynamics of this shape transition and support a mechanism based on the reduction in the micelle-water interfacial tension induced by the adsorption of the amphiphilic salicylate ions. At the threshold salt concentration, the interface is nearly saturated with adsorbed salicylate ions. Predictions of the effect of salt on the micelle structure in different CG solvent models, namely, single-site standard water and three-site polarizable water, show qualitative agreement. This suggests that phase transitions in aqueous micelle solutions could be investigated by using standard CG water models which allow for 3 orders of magnitude reduction in the computational time as compared to that required for atomistic MD simulations.     

Quantum effects in liquid water and ice: model dependence

This paper explores the influence of choice of potential model on the quantum effects observed in liquid water and ice. This study utilizes standard rigid models and a more formal context for the rigid-body centroid molecular dynamics methodology used to perform the quantum simulations is provided. Quantum and classical molecular dynamics simulations are carried out for liquid water and ice Ih at 298 and 220 K, respectively, with the simple point charge/extended and TIP4P-Ew water models. The results obtained for equilibrium and dynamical properties are compared with those recently reported on TIP4P [L. Hernandez de la Pena and P. G. Kusalik, J. Chem. Phys. 121, 5992 (2004); L. Hernandez de la Pena et al., J. Chem. Phys 123, 144506 (2005)]. For the liquid, an energy shift of about 8% and an average molecular uncertainty of about 11 degrees were found independently of the water model. The self-diffusion coefficient consistently increases by more than 50% when going from the classical to the quantum system and quantum dynamics are found to reproduce the experimental isotopic shifts with the models examined. The ice results compare remarkably well with those previously reported for the TIP4P water model; they confirm that quantum effects are considerable and that the quantum mechanical uncertainty and the energy shifts due to quantization are smaller in ice than in liquid water. The relevance of these findings in the context of the construction of water models is briefly discussed.     

Computer simulations of molecular motion in liquid crystals by the method of Brownian dynamics

Computer simulations of the molecular motion of polymer chains in the presence of a strong nematic field were carried out by the method of Brownian dynamics. Two models were studied: the first model (linear liquid crystal) is a freely jointed chain with rigid bonds, the second model (comb-like liquid crystal) is a chain with fixed bond angles and rigid side groups. The influence of ordering on chain conformations, orientational and translational mobility and spectra of relaxation times was investigated.     

Carbon nanotubes in benzene: internal and external solvation

The structure and dynamics of benzene inside and outside of single-walled carbon nanotubes (SWNTs) in the (n,n) armchair configuration are studied via molecular dynamics computer simulations. Irrespective of the nanotube diameter, benzene molecules form cylindrical solvation shell structures on the outside of the nanotubes. Their molecular planes near the SWNTs in the first external solvation shell are oriented parallel to the nanotube surface, forming a π-stacked structure between the two. By contrast, the benzene distributions in the interior of the SWNTs are found to vary markedly with the nanotube diameter. In the case of the (7,7) and (8,8) nanotubes, internal benzene forms a single-file distribution, either in a vertex-to-vertex (n = 7) or face-to-face (n = 8) orientation between two neighboring molecules. Inside a slightly wider (9,9) nanotube channel, however, a cylindrical single-shell distribution of benzene arises. A secondary solvation structure, which begins to appear inside (10,10), develops into a full structure separate from the first internal solvation shell in (12,12). The ring orientation of internal benzene is generally parallel to the nanotube wall for n = 9-12, while it becomes either slanted with respect to (n = 7), or perpendicular to (n = 8), the nanotube axis. The confinement inside the small nanotube pores exerts a strong influence on the dynamics of benzene. Both translational and rotational dynamics inside SWNTs are slower and more anisotropic than in liquid benzene. It is also found that reorientational dynamics of internal benzene deviate dramatically from the rotational diffusion regime and change substantially with the nanotube diameter.     

Interior and interfacial aqueous solvation of benzene dicarboxylate dianions and their methylated analogues: A combined molecular dynamics and photoelectron spectroscopy study

Aqueous solvation of benzene dicarboxylate dianions (BCD(2-)) was studied by means of photoelectron spectroscopy and molecular dynamics simulations. Photoelectron spectra of hydrated o- and p-BCD(2-) with up to 25 water molecules were obtained. An even-odd effect was observed for the p-BCD(2-) system as a result of the alternate solvation of the two negative charges. However, the high polarizability of the benzene ring makes the two carboxylate groups interact with each other in p-BCD(2-), suppressing the strength of this even-odd effect compared with the linear dicarboxylate dianions linked by an aliphatic chain. No even-odd effect was observed for the o-BCD(2-) system, because each solvent molecule can interact with the two carboxylate groups at the same time due to their proximity. For large solvated clusters, the spectral features of the solute decreased while the solvent features became dominant, suggesting that both o- and p-BCD(2-) are situated in the center of the solvated clusters. Molecular dynamics simulations with both nonpolarizable and polarizable force fields confirmed that all three isomers (o-, m-, and p-BCD(2-)) solvate in the aqueous bulk. However, upon methylation the hydrophobic forces overwhelm electrostatic interactions and, as a result, the calculations predict that the tetramethyl-o-BCD(2-) is located at the water surface with the carboxylate groups anchored in the liquid and the methylated benzene ring tilted away from the aqueous phase.     

Dynamical simulation of dipolar Janus colloids: Equilibrium structure and thermodynamics

The static microstructures and thermodynamics of a colloidal dispersion of dipolar Janus (DJ) particles-that is, dipolar spheres in which each hemisphere is specified by a different charge interaction-have been investigated through simulation. DJ particles are modeled at a high level of detail with pairwise potentials represented as a sum of a spherically symmetric soft repulsion and an orientation-dependent electrostatic component using continuous potentials. The latter is important because it allows for the use of conventional molecular dynamics simulations, and is in contrast to the patch model and dipolar hard sphere model, which are discontinuous and therefore do not. The electrostatics are represented through a rigorous pointwise (PW) covering of two different hemispheres filled by points of corresponding charge. An isotropic coarse-graining (CG) of the PW models serves as a limit of the structure wherein the orientations of the DJ particles can be pairwise averaged. Over the range of volume fractions and DJ charge densities studied-consistent with reversible structures absent of long-range correlations-the CG model agrees well with the PW model with respect to equilibrium structure (isotropic pair correlation) and ensemble free energy. Time-dependent relaxation simulations of the PW model suggest that chain structures are not expected in liquid phases in contrast to that which has been observed for point dipole models of simple polar fluids.     

Decoupling of exchange and persistence times in atomistic models of glass formers

With molecular dynamics simulations of a fluid mixture of classical particles interacting with pairwise additive Weeks-Chandler-Andersen potentials, we consider the time series of particle displacements and thereby determine the distributions for local persistence times and local exchange times. These basic characterizations of glassy dynamics are studied over a range of supercooled conditions and were shown to have behaviors, most notably decoupling, similar to those found in kinetically constrained lattice models of structural glasses. Implications are noted.