Molecular dynamics implementation cheat sheet

A cheat-sheet for implementing a molecular dynamics simulation
Source: "Nanosystems: Molecular Machinery, Manufacturing, and Computation"

The MM2 model (3.3.2.)

<math display="inline">\mathcal{V}_{s} = \frac{1}{2}k_{s}(r - r_{0})^{2}\lbrack 1 - k_{cubic}(r - r_{0})\rbrack</math> … stretching (3.4)
<math display="inline">\mathcal{V}_{\theta} = \frac{1}{2}k_{\theta}(\theta - \theta_{0})^{2}\lbrack 1 + k_{sextic}(\theta - \theta_{0})^{4}\rbrack</math> … bending (3.5)
<math display="inline">k_{s\bot} = k_{\theta}/r_{0}^{2}</math> … per length angular stiffness <math display="inline">k_{s\bot} \approx k_{s}/20</math> at <math display="inline">(r_{0},\theta_{0})</math> (3.6)
<math display="inline">\mathcal{V}_{\omega} = \frac{1}{2}\lbrack V_{1}(1 - \cos(1\omega)) + V_{2}(1 - \cos(2\omega)) + V_{3}(1 + \cos(3\omega))\rbrack</math> … torsion (3.7)
<math display="inline">\mathcal{V}_{vdw} = \epsilon_{vdw}\lbrack 2.48 \times 10^{5}\exp( - 12.5\frac{r}{r_{vdw0}}) - 1.924(\frac{r}{r_{vdw0}})^{- 6}\rbrack</math> … vdW attraction & non-bonded Pauli repulsion (3.8)
<math display="inline">\mathcal{V}_{s\theta} = k_{s\theta}(\theta - \theta_{0})\lbrack(r_{A} - r_{A0}) + (r_{B} - r_{B0})\rbrack</math> … stretch-bend interaction (3.9)

---

<math display="inline">k_{s,C - C} = 440N/m</math> — <math display="inline">r_{0,C - C} = 111.3pm</math> … (Table 3.2.)
<math display="inline">k_{\theta,C - C - C} = 450zJ/rad^{2}</math> — <math display="inline">\theta_{0,C - C - C} = 1.911rad = 109.47{^\circ}</math> … (Table 3.3.)
<math display="inline">\epsilon_{vdw,Csp^{3}} = 357yJ</math> — <math display="inline">r_{vdw,Csp^{3}} = 190pm</math> … well depth <math display="inline">\epsilon_{vdw}</math> varies widely for other atoms (~10x) (Table 3.1.)
<math display="inline">V_{1,C - C - C - C} = 1.39zJ,\ V_{2,C - C - C - C} = 1.88zJ,\ V_{3,C - C - C - C} = 0.65zJ</math> … (Table 3.5.)
<math display="inline">k_{s\theta,C - C - C} = 1.2nN/rad</math> … (Table 3.6.)
<math display="inline">k_{cubic} = ?m^{- 1}</math> — <math display="inline">k_{sextic} = ?rad^{- 4}</math>

<math>n … nano = 10^{-9}</math> — <math>p … pico = 10^{-12}</math>
<math>z … zepto = 10^{-21}</math> — <math>y … yocto = 10^{-24}</math>

Bonds under large loads (beyond MM2) (3.3.3.)

<math display="inline">\mathcal{V}_{morse} = D_{e}({1 - \exp\lbrack - \beta(r - r_{0})\rbrack}^{2} - 1)</math> … tensile (3.10)
<math display="inline">\mathcal{V}_{lippincott} = D_{e}\lbrack 1 - \exp( - \frac{k_{s}r_{0}(r - r_{0})^{2}}{2D_{e}r})\rbrack,\ \ r \geq r_{0}</math> … tensile – more accurate for larger distances (3.15)
For compressive loads <math display="inline">\mathcal{V}_{vdw}</math> see (3.8) above <math display="inline">F_{vdw} = - \frac{\partial}{\partial r}\mathcal{V}_{vdw}</math> … (3.16)
Note: Popular Lennard-Jones 6-12 potential is too steep in repulsive regime!
MM2 potential within 10% of experiment for <math display="inline">0.5r_{vdw}</math> (<math display="inline">> 100yJ</math> repulsion) 🙂

<math display="inline">\beta = \sqrt{k_{s}/(2D_{e})}</math> … (3.13)
<math display="inline">D_{e} \approx D_{0} + \frac{\hslash}{2}\sqrt{k_{s}/\mu} = D_{0} + \frac{\hslash\omega}{2};\ \ \mu = \frac{m_{1}m_{2}}{m_{1} + m_{2}}</math> … (3.14)
<math display="inline">D_{0}</math> … potential well depth
<math display="inline">D_{e,C - C} = 556yJ</math> — <math display="inline">k_{s,C - C} = 440N/m</math> — <math display="inline">r_{0,C - C} = 152.3pm</math> … (Table 3.8)

Verlet Integration (not in Nanosystems)

<math display="inline">\overset{\rightarrow}{x}(t + \Delta t) = 2\overset{\rightarrow}{x}(t) - \overset{\rightarrow}{x}(t - \Delta t) + \overset{\rightarrow}{a}(t)\Delta t^{2} + O(\Delta t^{4})</math>
<math display="inline">\overset{\rightarrow}{a}(t) = m\nabla_{\overset{\rightarrow}{x}}\mathcal{V}(t)</math>

MM2 vs MM3 (Chapter 3.3.2.g.)

★ MM2 prioritizes accuracy in energy and geometry
★ MM3 prioritizes accuracy of vibration frequencies
★ MM3 predicts greater angle-bending stiffness by ~1.5x or more
★ MM3 has a more complex functional form – e.g. it adds:
– stretch-torsion interaction, cubic bending, quartic stretching
★ MM3 has lower energies forces and stiffness in deep repulsive regime by ~10%
★ MM2 overall seems like a more conservative (safe side wrong) choice for diamondoid nanomachinery

Bond cleavage & radical coupling (3.4.2)

The Morse potential <math display="inline">\mathcal{V}_{morse}</math> can approximate homolytic reactions
<math display="inline">\mathcal{V}_{anti - morse} = \frac{1}{2}D_{e}({1 + \exp\lbrack - \beta(r - r_{0})\rbrack}^{2} - 1)</math> … unpaired spins (3.21)
But unpaired spins being the limiting factor in reaction rate should
typically be avoidable in piezomechanosynthesis.
See Nanosystems Chapter 8.4.3.b. Radical coupling and intersystem crossing.

Nonbonded & large compression limits (3.3.3.b.)

<math display="inline">k_{s,vdw} \approx \frac{12.5}{r_{vdw0}}F_{vdw} \approx 3.5 \times 10^{10}m^{- 1} \cdot F_{vdw}</math> … (3.18)
<math display="inline">\mathcal{V}_{vdw} \approx 0.08r_{vdw0}F_{vdw} \approx 2.9 \times 10^{- 11}m \cdot F_{vdw}</math> …(3.19)

Abstraction reactions (Chapter 3.4.3.)

<math>\mathcal{V}_{LEPS}</math> extended London-Eyring-Polyano-Sato potential – details omitted here

Continuum models of van der Waals attraction (Chapter 3.5.1.)

Hamaker constant – details omitted here

Related

• Snapback – sideward bond bending stiffness is only 1/20th of radial bond stretching/compressing stiffness