Capabilities

Simulation program ms2 - Techniques and Assumptions
The model class that is supported by ms2 covers rigid multi-center Lennard-Jones (LJ) 12-6 interaction sites with an arbitrary number of superimposed electrostatic sites. The supported electrostatic models are point charges, point dipoles and point quadrupoles, which can be positioned anywhere within the molecule. Currently, ms2 is designed for electro-neutral species.

Basic techniques
The simulation tool ms2 is capable of sampling phase space, applying the two most fundamental molecular simulation techniques by molecular dynamics (MD) and Monte-Carlo (MC). MD simulations rely on the numerical solution of Newton's equations of motion: for a point in time, the intermolecular interactions, in particular the resulting forces and torques, are evaluated and treated as constant for a specified time step. These are the driving forces of the molecular motions. The displacements for the time step are calculated on that basis, resulting in a new configuration. This process is repeated in a loop. The chronologically ordered configurations are a time discretized approximation of a real-world molecular process. From this, the thermodynamic properties are derived.

MC simulation explores the phase space stochastically for a given molecular system. Molecules in the simulation volume are displaced randomly. The probability of accepting the displacement decays exponentially with increasing energetic cost, analogous to real-world systems. Therefore, the system is driven to yield a representative set of configurations. This Markov chain of configurations allows for a rigorous calculation of static thermodynamic data.

Simulation assumptions
Simulations with ms2 are performed in quasi-infinite Cartesian space using periodic boundary conditions [1] and the minimum image convention [2]. The molecular interactions are considered to be pairwise additive, like in most other simulation packages. Including three- and many-body interactions leads to a significant increase in computational effort, so that they are omitted in most force field approaches.

The computational effort in ms2 is reduced by the introduction of a cut-off radius, up to which the intermolecular interactions are explicitly evaluated. The contribution of the interactions with molecules beyond the cut-off radius are accounted for by correction schemes.

The basic molecular simulation techniques employed in ms2 are well described in the literature [3, 4, 5] and are not repeated here. All non-standard methods that are implemented in ms2 are introduced in our publication.

Overview over accessible thermodynamic quantities
The simulation tool ms2 allows for the determination of static and dynamic thermodynamic properties in equilibrium. The implemented static properties are:

• Vapor-liquid equilibria
  • Density
  • Saturated pressure
  • Heat of vaporization

• Thermal, caloric and entropic properties for homogeneous states
  • Internal energy
  • Enthalpy
  • Chemical potential

• Second derivatives of the partition function
  • Heat capacities c_v, c_p
  • (dU / dV)_p
  • Isothermal compressibility b_T
    
  • Volume expansivity a_p
    
  • (dH / dp)_T
  • Speed of sound

• Henry's law constant
• Second virial coefficient

The transport properties can be calculated on the fly during a MD simulation with a reasonable additional computational effort using the Green-Kubo formalism [6, 7]. The implemented dynamic properties are:

• Diffusion coefficient
  • Self-diffusion coefficient
  • Maxwell-Stefan diffusion coefficient

• Viscosities
  • Shear viscosity
  • Bulk viscosity

• Thermal conductivity

The quality of thermodynamic properties calculated by molecular simulation is basically determined by two factors: firstly, the employed molecular model, which fully defines the thermodynamics and deviates to some degree from the behavior of the real fluid; secondly, the sampling of the phase space during simulation, which is associated with finite statistical uncertainties. The first factor indicates that optimized molecular models are very important, which have been investigated in numerous cases in the past with ms2 [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. These molecular models, combined with ms2 and its analysis methods, allow for time-efficient high quality simulations. However, ms2 is capable of using any rigid Lennard-Jones based molecular model with superimposed electrostatic sites. The second factor can be controlled for a given simulation method by the number of molecules and the number of sampled configurations. The more data, the lower the statistical uncertainties.

Methods

Ensembles
The calculation of different thermodynamic properties requires different simulation conditions. One set of conditions usually corresponds to one ensemble. The following are currently supported by ms2:

• canonical ensemble (NVT)   - MD and MC
• micro-canonical ensemble (NVE)   - MD
• isobaric-isothermal ensemble (NpT)   - MD and MC
• grand equilibrium method (pseudo-\mu VT)   - MC

Most of these are well known and widely in use. For a detailed description and discussion of the different ensembles, we refer to the literature.

Integrators and thermostat
In ms2, two integrators are implemented to solve Newton's equations of motion during MD simulation:

• Leapfrog
• Gear predictor-corrector

These integrators are well known and described in literature.
The computational demand for both integration schemes is similar in ms2. The Gear integration, though being of higher order, is only 0.3\% slower for the calculation of one MD time step than the Leapfrog algorithm for a molecular model composed of three Lennard-Jones sites having three rotational degrees of freedom.
The thermostat incorporated in ms2 is velocity scaling. Here, the velocities are scaled, such that the actual kinetic energy matches the specified temperature. The scaling is applied equally

Cutoff
The program ms2 supports two different cut-off methods for the dispersive and repulsive interactions, point dipole and point quadrupole interactions:

• Site-site cut-off
• Center of mass cut-off

Simulation of point charges are restricted to the center of mass cut-off.
The dispersive and repulsive contributions of molecules beyond the cut-off are accounted for implicitly.
The long range electrostatic contributions to the internal energy and the pressure in the system respectively, are introduced using the reaction field method.

Literature
[1] M. Born, T. Karman; On fluctuations in spatial grids, Phys. Zeitschrift (1912), 13, 297 - 309
[2] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.T. Teller, E. Teller; Equation of State Calculations By Fast Computing Machines, J. Chem. Phys. (1953), 1087 - 1092
[3] M.P Allen, D.J. Tildesley; Computer Simulation of Liquids; Clarendon Press (1987)
[4] D. Frenkel, B. Smith; Understanding Molecular Simulation; Academic Press, Elsevier (1993)
[5] D.C. Rappaport; The Art of Molecular Dynamics Simulation; Cambridge University Press (2004)
[6] M.S. Green; Markoff Random Processes And The Statistical Mechanics of Time-Dependent Phenomena. 2. Irreversible Processes In Fluids; J. Chem. Phys (1954); 22; 398 - 413
[7] R. Kubo; Statistical-Mechanical Theory of Irreversible Processes .1. General Theory And Simple Applications To Magnetic And Conduction Problems; J. Phys. Soc. Japan (1957); 12; 570 - 586
[8] J. Stoll, J. Vrabec, H. Hasse, J. Fischer; Comprehensive study of the vapour-liquid equilibria of the pure two-centre Lennard-Jones plus pointquadrupole fluid; Fluid Phase Equil. (2001); 179; 339 - 362
[9] J. Stoll, J. Vrabec, H. Hasse; Comprehensive study of the vapour-liquid equilibria of the pure two-centre Lennard-Jones plus pointdipole fluid; Fluid Phase Equil. (2003); 209; 29 - 53
[10] T. Schnabel, J. Vrabec, H. Hasse;Henry's law constants of methane, nitrogen, oxygen and carbon dioxide in ethanol from 273 to 498 K: Prediction from molecular simulation; Fluid Phase Equil. (2005); 233; 134 - 143
[11] T. Schnabel, A. Srivastaca, J. Vrabec, H. Hasse; Hydrogen bonding of methanol in supercritical CO2: Comparison between H-1 NMR spectroscopic data and molecular simulation results; J. Phys. Chem. B (2007); 111; 9871 - 9878
[12] T. Schnabel, M. Coratda, J. Vrabec, S. Lago, H. Hasse; Molecular model for formic acid adjusted to vapor-liquid equilibria; Chem. Phys. Letters (2007); 435; 268 - 272
[13] T. Schnabel, J. Vrabec, H. Hasse;Molecular simulation study of hydrogen bonding mixtures and new molecular models for mono- and dimethylamine; Fluid Phase Equil. (2008); 263; 144 - 159
[14] B. Eckl, J. Vrabec, H. Hasse; Molecular modelling and simulation for the process design; Chem. Ing. Tech. (2008); 80; 25 - 33
[15] B. Eckl, J. Vrabec, H. Hasse; An Optimized Model for Ammonia; Mol. Phys. (2008); 106; 1039 - 1046
[16] B. Eckl, J. Vrabec, H. Hasse; On the Application of Force Fields for Predicting a wide Variety of Properties: Ethylene Oxide as an Example ; Fluid Phase Equil. (2008); 274; 16 - 26
[17] Y.L. Huang, J. Vrabec, H. Hasse; Systematic investigation of vapour-liquid equilibria of binary mixtures on the basis of polar two-center Lennard-Jones models; Fluid Phase Equil. (2009); 287; 62 - 69