:orphan:

Examples
========

This gallery walks through every public feature of ``chemistrykit.md``:
the Lennard-Jones fluid in reduced units (energy conservation, pressure,
and the radial distribution function g(r)); Morse/Buckingham/harmonic
bonded potentials and small bonded clusters; periodic boundaries; SHAKE
bond constraints; thermostats; and self-diffusion coefficients.

Each script in this gallery is self-contained and can be run directly with
``python examples/md/<section>/<script>.py``. Every script also carries an
RST module docstring as its title/description and uses ``# %%`` markers to
split narrative text from code, which is exactly what Sphinx-Gallery
renders into the pages below -- the script *is* the source of truth for
what you see, not a copy of it.

Sections
--------

- **lj_fluid** -- Lennard-Jones fluid simulations: Rahman's liquid argon
  g(r), a cross-check of the simulated speed distribution against
  :mod:`chemistrykit.statmech`'s Maxwell-Boltzmann distribution,
  hard-sphere-like packing order, the Weeks-Chandler-Andersen repulsive
  reference fluid, Verlet's integrator and neighbor list, and Einstein
  and Green-Kubo self-diffusion coefficients.
- **pair_potentials** -- the Lennard-Jones potential, Morse vs. harmonic
  bond potentials and the classical vibration of a two-body
  :class:`~chemistrykit.md.systems.pair_potentials.DiatomicOscillator`,
  the Buckingham potential, and a small bonded molecule vibrating under
  coupled harmonic bond-stretch and angle-bend terms.
- **periodic_boundaries** -- periodic boundary conditions and the
  minimum-image convention.
- **constraints** -- rigid bonds enforced with the SHAKE algorithm.
- **thermostats** -- Nose-Hoover, Berendsen weak-coupling, and stochastic
  velocity-rescaling temperature control.


.. raw:: html

  <div id='sg-tag-list' class='sphx-glr-tag-list'></div>


.. raw:: html

    <div class="sphx-glr-thumbnails">

.. thumbnail-parent-div-open

.. thumbnail-parent-div-close

.. raw:: html

    </div>

Constraints
------------

Rigid bond-length constraints enforced with the SHAKE algorithm, removing
fast bond vibrations from a molecular-dynamics trajectory.


.. raw:: html

  <div id='sg-tag-list' class='sphx-glr-tag-list'></div>


.. raw:: html

    <div class="sphx-glr-thumbnails">

.. thumbnail-parent-div-open

.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Ryckaert, Ciccotti, and Berendsen&#x27;s SHAKE algorithm (1977) replaces stiff bond springs by exact length constraints, restored after every Verlet position update by iterative corrections along the old bond vectors (~chemistrykit.md.systems.constraints.shake). A water-like triatomic is simulated twice: with stiff harmonic O-H springs (~chemistrykit.md.systems.pair_potentials.HarmonicMolecule), and with the two bonds held rigid by SHAKE (~chemistrykit.md.systems.constraints.ShakeMolecule), keeping the same flexible angle-bending term. Both runs use the same time step, chosen comfortably for the slow bending motion but too long for the fast bond stretch.">

.. only:: html

  .. image:: /api/gallery/md/constraints/images/thumb/sphx_glr_plot_01_shake_rigid_bonds_thumb.png
    :alt:

  :doc:`/api/gallery/md/constraints/plot_01_shake_rigid_bonds`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">SHAKE: rigid bonds let a bending molecule take larger time steps</div>
    </div>


.. thumbnail-parent-div-close

.. raw:: html

    </div>

Lennard-Jones fluid
--------------------

Simulations of the Lennard-Jones fluid: Rahman's liquid argon, the
Maxwell-Boltzmann speed distribution, hard-sphere-like packing order,
the Weeks-Chandler-Andersen repulsive reference fluid, Verlet's
integrator and neighbor list, and self-diffusion from the Einstein
mean-squared displacement and the Green-Kubo velocity autocorrelation.


.. raw:: html

  <div id='sg-tag-list' class='sphx-glr-tag-list'></div>


.. raw:: html

    <div class="sphx-glr-thumbnails">

.. thumbnail-parent-div-open

.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Rahman&#x27;s 1964 simulation followed 864 argon atoms interacting through a Lennard-Jones potential with \epsilon/k_B = 120\ \text{K} and \sigma = 3.4\ \text{\AA}, at T = 94.4\ \text{K} and a mass density of 1.374\ \text{g/cm}^3 -- liquid argon just above its triple point. Here the same state point is simulated with LJFluid (512 atoms, reduced units), and the radial distribution function (~chemistrykit.md.systems.lj_fluid.LJFluid.radial_distribution_function) is converted back to ångströms, as Rahman compared it with neutron- and X-ray-scattering data on real liquid argon.">

.. only:: html

  .. image:: /api/gallery/md/lj_fluid/images/thumb/sphx_glr_plot_01_rahman_liquid_argon_thumb.png
    :alt:

  :doc:`/api/gallery/md/lj_fluid/plot_01_rahman_liquid_argon`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Rahman's liquid argon: g(r) of a Lennard-Jones liquid in real units</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="A direct link between chemistrykit.md and chemistrykit.statmech: molecular dynamics produces an actual ensemble of particle velocities, and statistical mechanics predicts what their distribution should look like at equilibrium. After a short equilibration run, an LJFluid&#x27;s particle speeds are histogrammed and compared against chemistrykit.statmech.MaxwellBoltzmannSpeedDistribution evaluated at the simulation&#x27;s own (instantaneous, measured) temperature -- no fitting involved.">

.. only:: html

  .. image:: /api/gallery/md/lj_fluid/images/thumb/sphx_glr_plot_02_maxwell_boltzmann_check_thumb.png
    :alt:

  :doc:`/api/gallery/md/lj_fluid/plot_02_maxwell_boltzmann_check`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Cross-checking an MD trajectory against the Maxwell-Boltzmann distribution</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Alder and Wainwright&#x27;s 1957 computer experiments showed that particles with no attraction at all -- hard spheres -- order into a crystal-like arrangement once they are packed densely enough: freezing driven by excluded volume (entropy), not by energy. This package integrates smooth forces rather than hard-sphere collision events, so the hard spheres are approximated by steeply repulsive particles: the Lennard-Jones potential cut off at its own minimum r=2^{1/6}\sigma (~chemistrykit.md.systems.lj_fluid.LennardJones.r_min) and shifted, which leaves only a repulsive wall. Comparing the radial distribution function g(r) of this purely repulsive LJFluid at low density and at a density near hard-sphere freezing shows sharp positional order appearing from packing alone.">

.. only:: html

  .. image:: /api/gallery/md/lj_fluid/images/thumb/sphx_glr_plot_03_hard_sphere_freezing_thumb.png
    :alt:

  :doc:`/api/gallery/md/lj_fluid/plot_03_hard_sphere_freezing`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Alder and Wainwright's hard spheres: packing-driven order with no attraction</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Weeks, Chandler, and Andersen split the Lennard-Jones potential at its minimum r_m=2^{1/6}\sigma into a purely repulsive reference part u_0 and a smooth attractive perturbation w:">

.. only:: html

  .. image:: /api/gallery/md/lj_fluid/images/thumb/sphx_glr_plot_04_wca_repulsive_reference_thumb.png
    :alt:

  :doc:`/api/gallery/md/lj_fluid/plot_04_wca_repulsive_reference`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Weeks-Chandler-Andersen: repulsive forces set the structure of a dense liquid</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Verlet&#x27;s 1967 algorithm advances positions from the two previous positions and the current force,">

.. only:: html

  .. image:: /api/gallery/md/lj_fluid/images/thumb/sphx_glr_plot_05_verlet_integrator_thumb.png
    :alt:

  :doc:`/api/gallery/md/lj_fluid/plot_05_verlet_integrator`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Verlet's integrator: bounded energy error instead of drift</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Checking every pair of N particles costs N(N-1)/2 distance evaluations per step. Verlet&#x27;s 1967 trick: list the pairs within r_c + r_\text{skin} and reuse that list for several steps; a pair outside the list can only come inside the true cutoff r_c after the two particles have closed the skin distance between them. VerletNeighborList implements this. Here a Lennard-Jones fluid is run with a list that is never rebuilt, and at every step the listed pairs are checked against the exact set of pairs inside the cutoff (~chemistrykit.md.utils.neighbor_list.build_neighbor_list), for lists built with no skin and with a skin of 0.3\sigma.">

.. only:: html

  .. image:: /api/gallery/md/lj_fluid/images/thumb/sphx_glr_plot_06_verlet_neighbor_list_thumb.png
    :alt:

  :doc:`/api/gallery/md/lj_fluid/plot_06_verlet_neighbor_list`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Verlet's neighbor list: a skin that lets one pair list serve many steps</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Einstein&#x27;s 1905 analysis of Brownian motion showed that a particle performing a random walk wanders off, on average, not in proportion to time but to its square root: the mean-squared displacement grows linearly,">

.. only:: html

  .. image:: /api/gallery/md/lj_fluid/images/thumb/sphx_glr_plot_07_einstein_mean_squared_displacement_thumb.png
    :alt:

  :doc:`/api/gallery/md/lj_fluid/plot_07_einstein_mean_squared_displacement`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Einstein's diffusion law: mean-squared displacement grows linearly in time</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Green (1954) and Kubo (1957) showed that transport coefficients are time integrals of equilibrium correlation functions. For self-diffusion the relevant function is the velocity autocorrelation C(t)=\langle\vec v(0)\cdot\vec v(t)\rangle, and">

.. only:: html

  .. image:: /api/gallery/md/lj_fluid/images/thumb/sphx_glr_plot_08_green_kubo_velocity_autocorrelation_thumb.png
    :alt:

  :doc:`/api/gallery/md/lj_fluid/plot_08_green_kubo_velocity_autocorrelation`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Green-Kubo: the diffusion coefficient as the integral of the velocity autocorrelation</div>
    </div>


.. thumbnail-parent-div-close

.. raw:: html

    </div>

Pair potentials
----------------

The Lennard-Jones 12-6 potential, Morse vs. harmonic bond potentials and
the classical vibration of a two-body diatomic oscillator, the Buckingham
exp-6 potential, and a small bonded molecule vibrating under coupled
harmonic bond-stretch and angle-bend terms.


.. raw:: html

  <div id='sg-tag-list' class='sphx-glr-tag-list'></div>


.. raw:: html

    <div class="sphx-glr-thumbnails">

.. thumbnail-parent-div-open

.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="The harmonic bond potential (~chemistrykit.md.systems.pair_potentials.HarmonicBond) is the small-oscillation limit of the anharmonic Morse potential -- they agree near the equilibrium separation but diverge away from it (Morse correctly flattens out toward a finite dissociation energy; harmonic does not). DiatomicOscillator integrates the classical two-body bond-vibration problem via the shared velocity-Verlet integrator, for either potential.">

.. only:: html

  .. image:: /api/gallery/md/pair_potentials/images/thumb/sphx_glr_plot_01_morse_vs_harmonic_thumb.png
    :alt:

  :doc:`/api/gallery/md/pair_potentials/plot_01_morse_vs_harmonic`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Morse vs. harmonic bond potentials, and diatomic bond vibration</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="HarmonicMolecule builds a small molecular-mechanics &quot;force field&quot; out of nothing but harmonic bond (~chemistrykit.md.systems.pair_potentials.HarmonicBond) and harmonic angle-bend (~chemistrykit.md.systems.pair_potentials.HarmonicAngle) terms -- the classical, spring-and-hinge picture of molecular structure that founded molecular mechanics. A bent, water-like triatomic (one heavy vertex atom, two light outer atoms, equilibrium angle 104.5 degrees) is displaced away from its equilibrium geometry and left to vibrate under velocity-Verlet integration alone (no thermostat): both bond lengths and the bond angle oscillate indefinitely about their equilibrium values, and the total energy is conserved to within integration error.">

.. only:: html

  .. image:: /api/gallery/md/pair_potentials/images/thumb/sphx_glr_plot_02_harmonic_molecule_thumb.png
    :alt:

  :doc:`/api/gallery/md/pair_potentials/plot_02_harmonic_molecule`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">A bent triatomic molecule: coupled bond-stretch and angle-bend vibrations</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="The Buckingham (exp-6) potential replaces the Lennard-Jones r^{-12} repulsive wall with an exponential -- closer to the true quantum-mechanical electron-overlap repulsion -- while keeping the same r^{-6} dispersion attraction. Matched to a Lennard-Jones potential&#x27;s well depth, position, and long-range tail, the two agree closely everywhere near and beyond the minimum. But the exponential repulsion is eventually overwhelmed by the unbounded -C/r^6 attractive term at short enough range, so the Buckingham potential has a second, spurious turnover well inside the true minimum: past that point it falls back toward -\infty as r\to0 instead of diverging to +\infty the way Lennard-Jones&#x27;s pure r^{-12} repulsion does.">

.. only:: html

  .. image:: /api/gallery/md/pair_potentials/images/thumb/sphx_glr_plot_03_buckingham_vs_lj_thumb.png
    :alt:

  :doc:`/api/gallery/md/pair_potentials/plot_03_buckingham_vs_lj`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Buckingham vs. Lennard-Jones: exponential repulsion and its inner turnover</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Lennard-Jones&#x27;s 1924 pair potential combines a steep r^{-12} repulsion with the r^{-6} London-dispersion attraction:">

.. only:: html

  .. image:: /api/gallery/md/pair_potentials/images/thumb/sphx_glr_plot_04_lennard_jones_potential_thumb.png
    :alt:

  :doc:`/api/gallery/md/pair_potentials/plot_04_lennard_jones_potential`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">The Lennard-Jones 12-6 potential: repulsion, dispersion, and the minimum</div>
    </div>


.. thumbnail-parent-div-close

.. raw:: html

    </div>

Periodic boundaries
--------------------

Periodic boundary conditions and the minimum-image convention: how a
small simulation box stands in for an effectively infinite bulk fluid.


.. raw:: html

  <div id='sg-tag-list' class='sphx-glr-tag-list'></div>


.. raw:: html

    <div class="sphx-glr-thumbnails">

.. thumbnail-parent-div-open

.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Metropolis and coworkers (1953) avoided the walls of a finite box by surrounding it with periodic copies of itself: a particle leaving one face re-enters through the opposite one (~chemistrykit.md.utils.pbc.wrap_positions), and each pair interacts through whichever periodic image of the partner is nearest (~chemistrykit.md.utils.pbc.minimum_image_displacement):">

.. only:: html

  .. image:: /api/gallery/md/periodic_boundaries/images/thumb/sphx_glr_plot_01_minimum_image_convention_thumb.png
    :alt:

  :doc:`/api/gallery/md/periodic_boundaries/plot_01_minimum_image_convention`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Periodic boundary conditions and the minimum-image convention</div>
    </div>


.. thumbnail-parent-div-close

.. raw:: html

    </div>

Thermostats
------------

Temperature control for a Lennard-Jones fluid: the Nose-Hoover
extended-system thermostat, Berendsen weak coupling, and
Bussi-Donadio-Parrinello stochastic velocity rescaling.


.. raw:: html

  <div id='sg-tag-list' class='sphx-glr-tag-list'></div>


.. raw:: html

    <div class="sphx-glr-thumbnails">

.. thumbnail-parent-div-open

.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Nose&#x27;s extended system, in Hoover&#x27;s friction-coefficient form, adds a single dynamical variable \xi to Newton&#x27;s equations:">

.. only:: html

  .. image:: /api/gallery/md/thermostats/images/thumb/sphx_glr_plot_01_nose_hoover_thermostat_thumb.png
    :alt:

  :doc:`/api/gallery/md/thermostats/plot_01_nose_hoover_thermostat`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">The Nose-Hoover thermostat: a friction variable that steers the temperature</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Berendsen and coworkers (1984) couple the system to a heat bath at T_0 by rescaling all velocities each step by">

.. only:: html

  .. image:: /api/gallery/md/thermostats/images/thumb/sphx_glr_plot_02_berendsen_weak_coupling_thumb.png
    :alt:

  :doc:`/api/gallery/md/thermostats/plot_02_berendsen_weak_coupling`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Berendsen's weak coupling: exponential relaxation to the bath temperature</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Bussi, Donadio, and Parrinello (2007) kept Berendsen&#x27;s exponential relaxation of the kinetic energy K but added a stochastic term of exactly the right size,">

.. only:: html

  .. image:: /api/gallery/md/thermostats/images/thumb/sphx_glr_plot_03_stochastic_velocity_rescaling_thumb.png
    :alt:

  :doc:`/api/gallery/md/thermostats/plot_03_stochastic_velocity_rescaling`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Bussi-Donadio-Parrinello: stochastic velocity rescaling samples the canonical ensemble</div>
    </div>


.. thumbnail-parent-div-close

.. raw:: html

    </div>


.. toctree::
   :hidden:
   :includehidden:


   /api/gallery/md/constraints/index.rst
   /api/gallery/md/lj_fluid/index.rst
   /api/gallery/md/pair_potentials/index.rst
   /api/gallery/md/periodic_boundaries/index.rst
   /api/gallery/md/thermostats/index.rst



.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_