Systematic Molecular Dynamics Simulations Research Articles

Linear-Scaling ab-initio Calculations for Large and Complex Systems

A brief review of the SIESTA project is presented in the context of linear-scaling density-functional methods for electronic-structure calculations and molecular-dynamics simulations of systems with a large number of atoms. Applications of the method to different systems are reviewed, including carbon nanotubes, gold nanostructures, adsorbates on silicon surfaces, and nucleic acids. Also, progress in atomic-orbital bases adapted to linear-scaling methodology is presented.

Symplectic integration of closed chain rigid body dynamics with internal coordinate equations of motion

Internal coordinate molecular dynamics (ICMD) is a recent efficient method for modeling polymer molecules which treats them as chains of rigid bodies rather than ensembles of point particles as in Cartesian MD. Unfortunately, it is readily applicable only to linear or tree topologies without closed flexible loops. Important examples violating this condition are sugar rings of nucleic acids, proline residues in proteins, and also disulfide bridges. This paper presents the first complete numerical solution of the chain closure problem within the context of ICMD. The method combines natural implicit fixation of bond lengths and bond angles by the choice of internal coordinates with explicit constraints similar to Cartesian dynamics used to maintain the chain closure. It is affordable for large molecules and makes possible 3–5 times faster dynamics simulations of molecular systems with flexible rings, including important biological objects like nucleic acids and disulfide-bonded proteins.

Dynamic properties of individual water molecules in a hydrophobic pore lined with acyl chains: a molecular dynamics study

Recently, a certain class of synthetic molecules has been shown to form ion channels, the pore of which is lined with hydrophobic acyl chains [M. Sokabe, in: F. Oosawa, H. Hayashi, T. Yoshioka (Eds.), Transmembrane Signaling and Sensation, JSSP/VNU Science Press BV, Tokyo, 1984, p. 119; F. Hayashi, M. Sokabe, M. Takagi, K. Hayashi, U. Kishimoto, Biochim. Biophys. Acta, 510 (1978) 305; M.J. Pregel, L. Jullien, J. Canceill, L. Lacombe, J.M. Lehn, J. Chem. Soc. Perkin Trans., 2 (1995) 417; Y. Tanaka, Y. Kobuke, M. Sokabe, Angew. Chem. Int. Ed. Engl., 34 (1995) 693; M. Sokabe, Z. Qi, K. Donowaki, H. Ishida, K. Okubo, Biophys. J., 70 (1996) A201; H. Ishida, K. Donowaki, Y. Inoue, Z. Qi, M. Sokabe, Chem. Lett. (1997) p. 953]. As an initial step towards understanding the physical mechanisms of ion permeation across such a hydrophobic pore, systematic molecular dynamics simulations were performed to investigate dynamic and energetic properties of water molecules inside the pore using a dimer of alanine- N′-acylated cyclic peptide as a channel model. Dynamic energy profiles for water molecules indicated that the energy barrier at the middle region of the pore is approximately 2–3 kcal/mol higher than that in the cap water region which was defined as a vicinity region of the channel entrance. Energetics analyses demonstrated that the mutual interactions among intrapore water molecules are the major factor to give favorable interaction (negative energy contribution) for themselves. The pore, despite being lined with acyl chains, has a favorable van der Waals interaction with intrapore water molecules. These results may help to explain why water-filled channels can be formed by the hydrophobic helices in natural channels.

A domain decomposition parallel processing algorithm for molecular dynamics simulations of systems of arbitrary connectivity

We describe in this paper methods for applying domain decomposition to a general purpose molecular dynamics program. The algorithm is suitable for either distributed memory parallel computers or shared memory machines with message passing libraries. A method is discussed in detail which allows molecules of arbitrary connectivity to be simulated within the domain decomposition approach. The algorithm also contains techniques to handle both rigid bond constraints and special CH 2 constraints which involve five atoms. Examples are given for molecular systems containing both types of constraints plus two, three- and four-body intramolecular potentials as well as short-range and long-range non-bonded potentials. The algorithm has been implemented on a Silicon Graphics Power Challenge machine using the MPI message passing library.

Atomic stress isobaric scaling for systems subjected to holonomic constraints

We show that in the molecular dynamics simulation of systems where the molecules are subjected to holonomic constraints the isobaric–isothermal (NPT) techniques can be applied by scaling in space the atomic coordinates rather than the molecular center of mass, as has been customary until now. The correct constrained distribution function is obtained by applying the original Andersen transformation or equivalently the Hoover-type NPT equations of motion. The latter set of equations are shown to be derived from an underlying Lagrangian. The atomic scaling is preferred to the molecular one when simulating large molecules since molecular metastabilities are more easily overcome. However, also for systems composed of small molecules, the approach can be interesting since in this way cell rotations are automatically avoided.

Molecular dynamics for nonequilibrium systems in which there are a small number of very hot particles in a cold bath: Reference system propagator methods

In molecular dynamics simulations of systems with very hot particles in contact with a bath of cold particles one must choose a very small time step dictated by the fast (hot) particles. This requires the recalculation of all the forces after each of these small time steps. We show how reference system propagator algorithms (RESPA) can be used to speed up these simulations, and also comment on some interesting physical phenomena involved in these nonequilibrium systems.

Comparison of two methods for solving linear equations occurring in molecular dynamics applications

Consider molecular dynamics simulations of systems containing large molecules such as proteins, where all fast bond vibrations are frozen. An important part of the computational overhead then comes from solving a large set of linear equations for determining the bonded forces (constraint forces). The corresponding matrix is symmetric, positive definite and has a very sparse nature. In this article we compare efficient implementations of the Conjugate Gradient method and an improved version of Cholesky Factorization for solving such types of linear equations. Based on both theoretical predictions and practical simulation tests on three different proteins, we show that the improved Cholesky method is a factor 15 faster than the Conjugate Gradient method.

A new parallel method for molecular dynamics simulation of macromolecular systems

Short-range molecular dynamics simulations of molecular systems are commonly parallelized by replicated-data methods, in which each processor stores a copy of all atom positions. This enables computation of bonded 2-, 3-, and 4-body forces within the molecular topology to be partitioned among processors straightforwardly. A drawback to such methods is that the interprocessor communication scales as N (the number of atoms) independent of P (the number of processors). Thus, their parallel efficiency falls off rapidly when large numbers of processors are used. In this article a new parallel method for simulating macromolecular or small-molecule systems is presented, called force-decomposition. Its memory and communication costs scale as N/√P, allowing larger problems to be run faster on greater numbers of processors. Like replicated-data techniques, and in contrast to spatial-decomposition approaches, the new method can be simply load balanced and performs well even for irregular simulation geometries. The implementation of the algorithm in a prototypical macromolecular simulation code ParBond is also discussed. On a 1024-processor Intel Paragon, ParBond runs a standard benchmark simulation of solvated myoglobin with a parallel efficiency of 61% and at 40 times the speed of a vectorized version of CHARMM running on a single Cray Y-MP processor. © 1996 by John Wiley & Sons, Inc.

A new parallel method for molecular dynamics simulation of macromolecular systems

Short-range molecular dynamics simulations of molecular systems are commonly parallelized by replicated-data methods, in which each processor stores a copy of all atom positions. This enables computation of bonded 2-, 3-, and 4-body forces within the molecular topology to be partitioned among processors straightforwardly. A drawback to such methods is that the interprocessor communication scales as N (the number of atoms) independent of P (the number of processors). Thus, their parallel efficiency falls off rapidly when large numbers of processors are used. In this article a new parallel method for simulating macromolecular or small-molecule systems is presented, called force-decomposition. Its memory and communication costs scale as N/√P, allowing larger problems to be run faster on greater numbers of processors. Like replicated-data techniques, and in contrast to spatial-decomposition approaches, the new method can be simply load balanced and performs well even for irregular simulation geometries. The implementation of the algorithm in a prototypical macromolecular simulation code ParBond is also discussed. On a 1024-processor Intel Paragon, ParBond runs a standard benchmark simulation of solvated myoglobin with a parallel efficiency of 61% and at 40 times the speed of a vectorized version of CHARMM running on a single Cray Y-MP processor. © 1996 by John Wiley & Sons, Inc.

Stability of Cellulose Structures Studied by MD Simulations. Could Mercerized Cellulose II Be Parallel?

By extensive and systematic molecular dynamics simulations eight stable cellulose crystal structures have been identified. We organized them in a coherent scheme, which gives a clue to possible interconversions between the various structures. From this scheme it can be inferred that mercerized cellulose II could in fact be a parallel structure or contain parallel regions. This parallel cellulose II structure fits into the same unit cell dimensions as antiparallel cellulose II and therefore can remain unnoticed. Regenerated cellulose II is likely to have all hydroxymethyl groups in the gt conformation similar to what was found in two recent determinations of the cellotetraose structures. Due to its characteristic layer structure cellulose I could crystallize in several different packing modes. Biosynthesis, however, enforces a parallel packing. It is shown that Iα and Iβ have a similar and specific parallel packing mode.

1-N-PHENYLAMINO-1-PHENYL-METHANEPHOSPHONIC ACID DIETHYL ESTER: QUANTUM MECHANICAL AND FORCE FIELD STUDIES

Abstract Using Molecular Modelling methods like Molecular Mechanics (MM), systematic conformational searches and Molecular Dynamics (MD) simulations (QUANTA 3.3.1/CHARMm 22) in combination with semi-empirical methods (MOPAC 6.0 PM3, VAMP 4.4 PM3) the conformational aspects of 1-N-phenylamino-1-phenylmethanepnosphonic acid diethyl ester were investigated in order to find the most stable con-formers. The results of our quantum mechanical structural analysis are in good agreement with the conformations obtained from X-ray diffraction studies. It is shown that PM3 reproduces quite well experimental geometries including phosphorus atoms in a-amino-phosphonic ester derivatives.

Multiple time scale methods for constant pressure molecular dynamics simulations of molecular systems

An efficient integration scheme is proposed based on the reversible reference system propagator algorithm (r-RESPA) for the molecular dynamics simulation of molecular solids or liquids at constant pressure. Starting from the Parrinello-Rhaman Lagrangian, the algorithm has been obtained using a Trotter factorization for the classical propagator. This r-RESPA scheme allows us to simulate efficiently the dynamics of completely flexible molecules in the condensed phases, under the constraint of constant pressure. The method may find potential application in the study of mesophases such as nematic or plastic crystals, or in studying the effect of pressure on the spectral properties of solids.

Molecular dynamics free energy calculation in four dimensions

A new method for the calculation of free energy differences by molecular dynamics simulations of molecular systems is presented. The main characteristic of the method is that the three-dimensional physical space is augmented by a further dimension. This allows the system to circumvent energy barriers in three-dimensional space. When simulating physical (three-dimensional) states, the dynamics along the 4th coordinate is uncoupled from the dynamics described with the original three coordinates. In the method presented here, free energy differences between physical states are determined using pathways leading over states where the dynamics along all four dimensions are coupled. Due to the barrier lowering in four dimensions, an increase of efficiency can be expected when low free energy physical states are separated by high barriers. The method is tested by application to an atomic liquid.

Efficient linear scaling algorithm for tight-binding molecular dynamics.

A novel formulation for tight-binding total energy calculations and tight-binding molecular dynamics, which scales linearly with the size of the system, is presented. The linear complexity allows us to treat systems of very large size and the algorithm is already faster than the best implementation of classical diagonalization for systems of 64 atoms. In addition, it is naturally parallelizable and it permits us therefore to perform molecular dynamics simulations of systems of unprecedented size. Finite electronic temperatures can also be taken into account. We illustrate this method by investigating structural and dynamical properties of solid and liquid carbon at different densities.

Molecular dynamics study of pressure in molecular systems

In molecular dynamics simulation of molecular systems, an atomistic model is needed to describe the intramolecular effects on system properties such as pressure. An expression for computing the pressure is derived based on the virial theorem, with explicit kinetic, intra-, and intermolecular contributions. It is shown that the virial terms arising from three- or four-body forces that depend only on internal angles are zero by using the Wilson S vector technique, and that only the two-body forces appear in the pressure expression. Molecular dynamics simulations are carried out for an atomistic model of benzene with intramolecular interactions based on ab initio harmonic potentials and intermolecular interactions given by semiempirical atom–atom potentials. Calculated total pressure depends on parameters of intermolecular potentials. An intramolecular contribution to the pressure is studied for both solid and liquid phases. In solids where the molecules are compressed, the intramolecular pressure contributions are appreciable and positive, while in liquids where molecular deformations are relatively small, the contributions are small and negative. A possible further improvement of the atomistic model of benzene is discussed.

An algorithm for calculating intramolecular angle-dependent forces on vector computers

We describe an approach based on projection methods for the calculation of angle-bending and torsional forces in molecular dynamics simulations. These forces are important in molecular dynamics simulations of systems containing polyatomic molecules. A significant speedup can be achieved using projection methods, because they require fewer high-cost operations than traditional crossproduct methods. Initial tests on a Cray X-M P show factors of 7 and 2.5 increase in speed for the calculation of angle-bending and torsional forces, respectively, relative to a comparable cross-product formulation. Our analysis of projection methods for calculating intramolecular angle-dependent forces provides a framework for the development of efficient programming structures.

Computer simulation and perturbation theory of fluids modelled using three- and six-site Lennard-Jones potentials

Systematic NVT molecular dynamics simulations have been performed for pure fluids of nonspherical molecules modelled using a triangular three-site and a hexagonal six-site Lennard-Jones (12:6) potential in order to obtain thermodynamic properties of these fluids, particularly the residual Helmholtz free energy. Molecular geometries, shown earlier to work well for propane and benzene, have been utilized in this study. Free energy values obtained from these simulations have been compared with predictions from the reference fluid simulations in a Weeks-Chandler-Andersen type site-frame perturbation expansion to the first order. Comparisons clearly show this site-frame expansion to work well at high densities and no significant effect of molecular anisotropy has been found on the deviations of these predictions, from the exact simulation values, vis-a-vis the WCA theory for the Lennard-Jones (12 : 6) fluid.

Treatment of rotational isomers in free energy evaluations. Analysis of the evaluation of free energy differences by molecular dynamics simulations of systems with rotational isomeric states

Using 1,2-dimethoxy ethane in aqueous solution as an example, the applicability of the perturbation method and the thermodynamic integration technique to evaluate free energy differences is considered for systems with multiple rotational isomeric states. It is shown that the naive application of these methods to evaluate free energy differences for such systems, even in a simple case such as the free energy of hydration of 1,2-dimethoxy ethane, may lead to unreasonable results. This problem is due to the fact that, in a conventional simulation, it is unlikely that all isomeric states will be sampled with the appropriate equilibrium probabilities. A procedure is proposed to estimate the contributions of isomeric states in free energy difference calculations.

Molecular Graphics and the Computer Simulation of Liquid Crystals

Abstract Molecular dynamics simulations of systems of hard ellipsoids of revolution have been carried out. Using configurations generated in these simulations, it is shown how molecular graphics can be used to examine the structure within the nematic liquid crystal phase, and close to the isotropic-nematic phase transition. A way of visualizing the degree of orientational order in these systems is presented.

Molecular dynamics for discontinuous potentials

A method to perform molecular dynamics simulations for systems of particles interacting with discontinuous potentials is presented. It is a generalization of Alder and Wainwright's algorithm for potentials represented by an arbitrary number of discontinuous horizontal line segments. The method is applied to several hard molecular fluids of various shapes, discrete Lennard-Jonesium, surfaces and mixtures in order to test its generality and flexibility. Results are compared, when possible, with previously published material, which includes numerical, theoretical and experimental results. These comparisons show that the method is a powerful tool in the study of molecular fluids and looks like a promising alternative approach to perform molecular dynamics for continuous potentials represented by a series of discontinuous line segments. The flexibility of the method will allow the testing of theories based on integral equations, like RISM and RAM, and to develop equations of state for complex hard molecules to...