Ammm:Md tp

Ensemble

Typical microscopic ensembles

Microcanonical: $\;\Omega (E,V,N)=e^{S/k}$ $\;\Omega (E,V,N)=e^{S/k}$
- Constant energy; the probability density is proportional to $\delta =\aleph (q,p)-E\,$

Canonical: $\;Z(T,V,N)=e^{-\aleph /kT}$ $\;Z(T,V,N)=e^{-\aleph /kT}$
- the probability density is proportional to $exp\left(-{\frac {\aleph }{kT}}\right)$

Grand canonical: $\;\Xi (T,V,\mu )=\sum _{N}^{}e^{(\mu N/kT)}{Z}_{NVT}=e^{\beta PV}$ $\;\Xi (T,V,\mu )=\sum _{N}^{}e^{(\mu N/kT)}{Z}_{NVT}=e^{\beta PV}$
- th density is proportional to $exp\left(-{\frac {\aleph -\mu N}{kT}}\right)$

isothermal-isobaric: $\Delta (N,P,T)=\int Z(N,V,T)\exp(-\beta PV)CdV.\,\;$ $\Delta (N,P,T)=\int Z(N,V,T)\exp(-\beta PV)CdV.\,\;$
- the density is proportional to $exp\left(-{\frac {\aleph +PV}{kT}}\right)$
- Note that $H=<\aleph >+P<V>\,$

In different situations, different ensembles shall be applied in MD simulations:

phase equilibrium; absorption where particle number N is a variable: Grand canonical

free energy: NPT or NVT. Note that free energy calculated from NPT or NVT is the same but the enthalpy and entropy components aren't necessary the same.

It should be noted that the choice of the ensemble (canonical or isobaric) in eq 1 is somewhat arbitrary. In a constant volume ensemble, eq 1 would correspond to the change in Helmholtz free energy at a specific volume. Not considering finite size effects, this Helmholtz free energy change would be equal to the corresponding Gibbs free energy change at the pressure determined by the equation of state of the system. Both quantities correspond to the change in chemical potential of the solute.19 The choice of the ensemble is, however, critical when considering the entropy and enthalpy changes.20

NPT resembles the experimental condition (open container) most closely; but the converge is slow.

How does it affect the error of sampling?

$\sigma ={\sqrt {{\frac {1}{N}}\sum _{i=1}^{N}(x_{i}-\mu )^{2}}},$

where N is the number of independent samples.

${{N}_{true}}={{N}_{record}}{\frac {{t}_{interval}}{t{}_{correlation}}}$

Comparing Ensemble

The averages of any quantity calculated in, say, the constant-NVE ensemble and the constant-NVT ensemble, will be equal in the thermodynamic limit, as long as we choose E and Г consistently so that $E=<E>_{NVT}$ . In fact, there are some restrictions on the kinds of functions for which this principle holds. The average should be a sum of single-particle functions

$A=\sum _{i=1}^{N}A_{i}$ for example: K, U, P, V, T....

The fluctuation is however, ensemble dependent!

$\delta ^{2}(A)=<\delta A^{2}>=<A^{2}>-<A>^{2}\,$

For example, the heat capacity:

In NVT, $kT^{2}C_{V}=<\delta \aleph ^{2}>_{NVT}=<\delta K^{2}>_{NVT}+<\delta U^{2}>_{NVT}$

and for a system of N partciles,

$<\delta K^{2}>_{NVT}={\frac {3N}{2}}(kT)^{2}$ so that

$kT^{2}(C_{V}-3/2Nk)=<\delta U^{2}>_{NVT}\,$

In NPT, $kT^{2}C_{P}=<\delta (\aleph +PV)^{2}>_{NPT}$

In NVE, the potential energy component of Cv ^[1]

In weak-coupling (when Berendsen thermostat is used)

Mor examples of fluctuation, structure and correlation properties in Allen and Tildesley book.

Temperature control

In Tinker, recommended thermostat is Bussi, used with RESPA integrator, 2fs time step. Recommended pressure control for homogeneous systems (protein in water is OK) is Monte Carlo, used with RESPA integrator, 2fs time step. For inhomogeneous system like lipids, use Langevin piston for pressure/temperature control and its own integrator. See Using Tinker&AMOEBA for details.

In the NVT, NPT, or $\mu$ VT ensemble simulations, the temperature is controlled by scaling atomic velocity in certain frequency. There are several ways to scale the velocity and maintain the T at desired value, however, not all produce a true ensemble (fluctuation).

Reference reading: ^[2]

Direct velocity scaling

or

Berendsen method of temperature-bath coupling

The following methods ensure a correct ensemble fluctuation:

Bussi Parrilleno

Nosé and Nosé–Hoover dynamics

The extended system (model plus a fictitious for temperature) Hamiltonian is

The choice of the fictitious mass Q of that additional degree of freedom is arbitrary but is critical to the success of a run. If Q is too small, the frequency of the harmonic motion of the extended degree of frequency is too high. This forces a smaller timestep to be used in integration. However, if Q is too large, the thermalization process is not efficient—as Q approaches infinity, there is no energy exchange between the heat bath and the model.
Q should be different for different models—Nosé (1991) suggests that Q should be proportional to gkBT, where g is the number of degrees of freedom in the model, kB is the Boltzmann constant, and T is the temperature.
To determine the proportionality constant, studies were done with a model consisting of a box of liquid argon containing 343 atoms at 87 K. Setting Q to 2.5 10-5 kcal mol-1 fs-2 for this model yielded good results. This proportionality constant, together with the gkBT term, was then used to generate Q for a box of water and an amorphous cell of polypropylene, which also yielded satisfactory results.
For the Nosé–Hoover method, the smaller the timestep, the closer it approaches the target temperature.
For energy to be conserved, the Nosé–Hoover method requires an accurate integrator. such as 0.5 fs as the timestep should be used.

Andersen method

One version of the Andersen method of temperature control involves randomizing the velocities of all atoms at a predefined “collision frequency”. (The other version involves choosing one atom at each timestep and changing its velocity according to the Boltzmann distribution.)

$f(v)=4\pi \left({\frac {m}{2\pi kT}}\right)^{3/2}\!\!v^{2}\exp \left[{\frac {-mv^{2}}{2kT}}\right]\qquad$

pressure control

Pressure is another state variable of the system. It is defined as the force per unit area. Standard atmospheric pressure is 1.013 bar, where 1 bar = 105 Pa (roughly 1 atm). Pressure is actually a tensor:

$P={\begin{pmatrix}Pxx&Pxy&Pxz\\Pyx&Pyy&Pyz\\Pzx&Pzy&Pzz\end{pmatrix}}$

In an isotropic situation, where forces are the same in all directions and there is no viscous force, the pressure tensor becomes

$P=p{\begin{pmatrix}1&0&0\\0&1&0\\0&0&1\end{pmatrix}}$

where p is the scalar pressure we are familiar with.

The thermodynamic pressure is defined by: $PV=Nk_{b}T+{\frac {2}{3}}<W>$

where W is the virial

$W={\frac {1}{2}}\sum _{i=1}^{N}{r}_{i}.{f}_{i}$

The instantaneous pressure of a system is then

$P={\frac {2}{3V}}(K+W)={\frac {1}{V}}{\begin{pmatrix}\sum _{i}^{N}miVi.{Vi}^{T}+\sum _{i}^{N}ri.{fi}^{T}\end{pmatrix}}$

$P={\frac {1}{V}}\left[{\begin{matrix}\sum _{i}{{{m}_{i}}{{v}_{ix}}{{v}_{ix}}}&\sum _{i}{{{m}_{i}}{{v}_{ix}}{{v}_{iy}}}&\sum _{i}{{{m}_{i}}{{v}_{ix}}{{v}_{iz}}}\\\sum _{i}{{{m}_{i}}{{v}_{iy}}{{v}_{ix}}}&\sum _{i}{{{m}_{i}}{{v}_{iy}}{{v}_{iy}}}&\sum _{i}{{{m}_{i}}{{v}_{iy}}{{v}_{iz}}}\\\sum _{i}{{{m}_{i}}{{v}_{iz}}{{v}_{ix}}}&\sum _{i}{{{m}_{i}}{{v}_{iz}}{{v}_{iy}}}&\sum _{i}{{{m}_{i}}{{v}_{iz}}{{v}_{iz}}}\\\end{matrix}}\right]+{\frac {1}{V}}\left[{\begin{matrix}\sum _{i}{{{r}_{ix}}{{f}_{ix}}}&\sum _{i}{{{r}_{ix}}{{f}_{iy}}}&\sum _{i}{{{r}_{ix}}{{f}_{iz}}}\\\sum _{i}{{{r}_{iy}}{{f}_{ix}}}&\sum _{i}{{{r}_{iy}}{{f}_{iy}}}&\sum _{i}{{{r}_{iy}}{{f}_{iz}}}\\\sum _{i}{{{r}_{iz}}{{f}_{ix}}}&\sum _{i}{{{r}_{iz}}{{f}_{iy}}}&\sum _{i}{{{r}_{iz}}{{f}_{iz}}}\\\end{matrix}}\right]$

The scalar pressure is the average of the trace:

 $p={\frac {1}{3}}\left(Pxx+Pyy+Pzz\right)$

Berendsen method

At each step, the x, y, and z coordinates of each atom are scaled by the factor:

Note that this method (as implemented) changes the cell uniformly,so that the size of the cell is changed, but not its shape. Therefore, for simulations such as crystal phase transitions, where both the cell size and shape are expected to change, this method is not appropriate.
The pressure fluctuation has been observed to be large during test runs with the Discover program and in studies in the literature (Brown and Clark 1984). Negative pressures sometimes occur because the virial can be negative, even though this defies the usual sense that pressure is a positive number. The calculated pressure depends on the cutoff distances used in the simulation.

Parrinello–Rahman method of pressure and stress control

Both size and shape of cell can change under the stress tensor.
The Lagrangians of the system are modified such that a term representing the kinetic energy of the cell depends on a user-defined mass-like parameter W. An elastic energy term pΩ is related to the pressure and volume or the stress and strain of the system. The equations of motion for the atoms and cell vectors can be derived from this Lagrangian. The motion of the cell vectors, which determines the cell shape and size, is driven by the difference between the target and internal stress.

Note the user-defined variable W, which determines the rate of change of the volume/shape matrix. A large W means a heavy, slow cell. In the limiting case, infinite W reverts to constant-volume dynamics. A small W means fast motion of the cell vectors. Although that may mean that the target stress can be reached faster, there may not be enough time for equilibration.

Types of MD

Langevin

The Langevin dynamics method (McCammon et al. 1976, Levy et al. 1979) approximates a full molecular dynamics simulation of a system by eliminating unimportant degrees of freedom (e.g. solvent). The effects of the eliminated degrees of freedom are simulated by mean and stochastic forces. Friction coefficients must be assigned to selected atoms before starting the dynamics run. Langevin dynamics includes a constant-temperature bath.

For a system of $N$ particles with masses $M$ , with coordinates $X = X (t)$ the resulting Langevin equation is ^[3]

M{\ddot {X}}=-\nabla U(X)-\gamma M{\dot {X}}+{\sqrt {2\gamma k_{B}TM}}R(t)\,,

where $U (X)$ is the potential energy; is the friction coefficent or damping factor; the dot is a time derivative such that ${\dot {X}}$ is the velocity and ${\ddot {X}}$ is the acceleration; and $R (t)$ is a random (Gaussian) force with zero-mean, satisfying

\left\langle R(t)\right\rangle =0

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \left\langle R(t)R(t') \right\rangle = \delta(t-t')}

As $γ$ grows, it spans the inertial all the way to the diffusive (Brownian) regime. The limit of non-inertia is commonly described as Brownian dynamics, or overdamped, creeping motion.

The Hamiltonian of a generalized LE is

Brownian dynamics

In Brownian dynamics, no acceleration is assumed to take place. Thus, the Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle M\ddot{X}(t)} term is neglected, and the sum of these terms is zero.

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0 = - \nabla U(X) - \gamma M\dot{X}+ \sqrt{2 \gamma k_B T M} R(t)}

Defining $ζ = γ M$ , and using the Einstein relation, $D = k B T / ζ$ , it is often convenient to write the equation as,

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \dot{X}(t) = - \nabla U(X)/\zeta + \sqrt{2 D} R(t)}

A free program for brownian dynamics: http://adrik.bchs.uh.edu/uhbd.html

References

↑ Ren and Ponder, J. Phys. Chem. B 2003, 107, 5933-5947
↑ Forcefield-Based Simulations/September 1998, MSI
↑ Template:Cite bookfckLR

[1] Ren and Ponder, J. Phys. Chem. B 2003, 107, 5933-5947

[2] Forcefield-Based Simulations/September 1998, MSI

[3] Template:Cite bookfckLR

[1]

[2]

[3]

Ammm:Md tp

Contents

Ensemble

Comparing Ensemble