chemistrykit.md#

chemistrykit.md hub · History · Examples · API reference

chemistrykit.md: molecular dynamics and force fields.

Lennard-Jones fluid simulation in reduced units (energy, pressure, radial distribution function g(r)); pairwise potentials (Morse, Buckingham, harmonic bond/angle terms); SHAKE bond constraints; velocity-rescaling, Berendsen, Nose-Hoover, and stochastic velocity-rescaling thermostats; Einstein and Green-Kubo self-diffusion coefficients; periodic boundary conditions with the minimum-image convention and a Verlet neighbor list. Integrated via chemistrykit.integrators’s velocity-Verlet (chemistrykit.integrators.velocity_verlet_step()) – not a reimplementation, see chemistrykit.md.core.base_system.MolecularDynamicsSystem.step().

MD’s trajectories are a natural cross-check for chemistrykit.statmech’s distribution functions: the speed distribution of an equilibrated LJFluid should follow chemistrykit.statmech.MaxwellBoltzmannSpeedDistribution at the simulation’s instantaneous temperature (see examples/md/lj_fluid/plot_02_maxwell_boltzmann_check.py).

class chemistrykit.md.AnglePotential[source]#

Bases: ABC

Common interface for a three-body bond-angle potential \(U(\theta)\).

Used for angle-bending terms in HarmonicMolecule, where \(\theta\) is the angle at a vertex atom j between the two bonds i-j and k-j.

abstractmethod energy(theta)[source]#

Return \(U(\theta)\) for angle(s) theta, in radians.

abstractmethod torque(theta)[source]#

Return \(-dU/d\theta\) for angle(s) theta, in radians.

class chemistrykit.md.BerendsenThermostat(target_temperature, tau, k_b=1.0)[source]#

Bases: object

The Berendsen weak-coupling thermostat.

Rescales every velocity by a common factor after each step so that the instantaneous kinetic temperature relaxes exponentially toward target_temperature with time constant \(\tau\) – H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola & J. R. Haak, J. Chem. Phys. 81, 3684 (1984), eq. 11:

\[\lambda^2 = 1 + \frac{\Delta t}{\tau}\left(\frac{T_0}{T} - 1\right), \qquad \frac{dT}{dt} = \frac{T_0 - T}{\tau}.\]

For \(\tau = \Delta t\) this reduces to the immediate, complete rescaling of VelocityRescalingThermostat. Like that scheme it suppresses canonical kinetic-energy fluctuations, so it does not sample the NVT ensemble exactly (see StochasticVelocityRescalingThermostat).

Parameters:

  • target_temperature (float) – \(T_0\); must be positive.

  • tau (float) – Coupling time constant; must be positive.

  • k_b (float)

Examples

With no forces acting, each call moves the temperature a fraction \(\Delta t/\tau\) of the way to the target:

>>> from chemistrykit.md.systems.lj_fluid import LJFluid
>>> fluid = LJFluid.from_lattice(n_per_side=3, cutoff=1.5, density=0.6, temperature=2.0, rng=0)
>>> thermostat = BerendsenThermostat(target_temperature=1.0, tau=0.1)
>>> thermostat.apply(fluid, dt=0.01)
>>> round(fluid.temperature(), 10)  # 2.0 + 0.1 * (1.0 - 2.0)
1.9
apply(system, dt)[source]#

Rescale system.velocities in place by the Berendsen factor.

Parameters:
Return type:

None

class chemistrykit.md.Buckingham(A, B, C)[source]#

Bases: PairPotential

The Buckingham (exp-6) pair potential.

\[U(r) = A e^{-Br} - \frac{C}{r^6}\]

An alternative to Lennard-Jones for the short-range repulsion, replacing the \(r^{-12}\) term with an exponential (closer to the true quantum-mechanical electron-overlap repulsion) at the cost of an unphysical turnover to \(U\to-\infty\) as \(r\to0\) that pure LJ does not have – R. A. Buckingham, Proc. R. Soc. A 168, 264 (1938); used extensively for rare-gas potentials and ionic-crystal lattice energies (see the future chemistrykit.crystal).

Parameters:

  • A (float) – Repulsive prefactor (energy units), must be positive.

  • B (float) – Repulsive exponential decay rate (inverse length), must be positive.

  • C (float) – Dispersion (attractive \(r^{-6}\)) coefficient (energy * length6), must be non-negative.

Examples

At large separation the potential is dominated by the (attractive, slower-decaying) dispersion term, so both the energy and the force are negative:

>>> pot = Buckingham(A=1.0e4, B=3.0, C=1.0)
>>> bool(pot.energy(10.0) < 0.0)
True
>>> bool(pot.force_scalar(10.0) < 0.0)
True
energy(r)[source]#

Return the potential energy U(r) at separation(s) r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

force_scalar(r)[source]#

Return the radial force magnitude \(f(r) = -dU/dr\).

By convention, positive is repulsive (pushes the two particles apart along their separation vector \(\vec r_i - \vec r_j\)) and negative is attractive – so the force on particle i due to particle j is f(r) * (r_i - r_j) / r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

class chemistrykit.md.DiatomicOscillator(potential, m1, m2, r0, v_rel0=0.0, axis=(1.0, 0.0, 0.0))[source]#

Bases: MolecularDynamicsSystem

Two atoms bound by a radial potential: classical bond-vibration MD.

A textbook two-body problem: with no external forces, the center-of-mass moves uniformly and the relative coordinate oscillates as if a single particle of the reduced mass \(\mu=m_1m_2/(m_1+m_2)\) were moving in the bond potential itself. For a HarmonicBond, this gives the exact classical vibration period \(T=2\pi\sqrt{\mu/k}\) (see harmonic_period()), which the simulated trajectory is checked against in the test suite.

Parameters:

  • potential (HarmonicBond or Morse) – The bond potential.

  • m1 (float) – Atomic masses.

  • m2 (float) – Atomic masses.

  • r0 (float) – Initial bond length.

  • v_rel0 (float) – Initial relative (radial) velocity along axis.

  • axis (array-like, shape (3,), default (1, 0, 0)) – Bond direction.

bond_length()[source]#

Current separation between the two atoms.

Return type:

float

Returns:

float

degrees_of_freedom()[source]#

A single relative-coordinate (bond-stretch) vibrational degree of freedom.

Return type:

int

forces_and_potential(positions)[source]#

Return the force on every particle and the total potential energy.

Parameters:

positions (ndarray, shape (n_particles, n_dim))

Returns:

  • forces (ndarray, shape (n_particles, n_dim))

  • potential_energy (float)

static harmonic_period(k, m1, m2)[source]#

Classical vibration period for a harmonic bond, \(T=2\pi\sqrt{\mu/k}\).

Parameters:

  • k (float) – Harmonic force constant.

  • m1 (float) – Atomic masses.

  • m2 (float) – Atomic masses.

Return type:

float

Returns:

float

Examples

>>> round(float(DiatomicOscillator.harmonic_period(k=1.0, m1=1.0, m2=1.0)), 6)
4.442883
class chemistrykit.md.HarmonicAngle(k_theta, theta0)[source]#

Bases: AnglePotential

A harmonic bond-angle-bending potential.

\[U(\theta) = \frac{1}{2}k_\theta(\theta-\theta_0)^2\]
Parameters:

  • k_theta (float) – Force constant (energy per radian squared), must be positive.

  • theta0 (float) – Equilibrium angle, in radians.

Examples

>>> angle = HarmonicAngle(k_theta=100.0, theta0=1.9106)  # ~109.47 deg, tetrahedral
>>> float(angle.energy(angle.theta0))
0.0
>>> float(angle.torque(angle.theta0))
-0.0
energy(theta)[source]#

Return \(U(\theta)\) for angle(s) theta, in radians.

torque(theta)[source]#

Return \(-dU/d\theta\) for angle(s) theta, in radians.

class chemistrykit.md.HarmonicBond(k, r0)[source]#

Bases: PairPotential

A harmonic (Hookean) bond-stretch potential.

\[U(r) = \frac{1}{2}k(r-r_0)^2\]

The simplest bonded term in a classical force field – the small -oscillation limit of any real bond potential, including Morse (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 8, the classical harmonic oscillator).

Parameters:

  • k (float) – Force constant, must be positive.

  • r0 (float) – Equilibrium bond length.

Examples

>>> bond = HarmonicBond(k=50.0, r0=1.0)
>>> float(bond.energy(1.0))
0.0
>>> round(float(bond.force_scalar(1.1)), 10)
-5.0
energy(r)[source]#

Return the potential energy U(r) at separation(s) r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

force_scalar(r)[source]#

Return the radial force magnitude \(f(r) = -dU/dr\).

By convention, positive is repulsive (pushes the two particles apart along their separation vector \(\vec r_i - \vec r_j\)) and negative is attractive – so the force on particle i due to particle j is f(r) * (r_i - r_j) / r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

class chemistrykit.md.HarmonicMolecule(positions, velocities, masses, bonds=(), angles=())[source]#

Bases: MolecularDynamicsSystem

A small, non-periodic cluster of atoms held together by harmonic bonds and angles.

A minimal molecular “force field”: mirrors chemistrykit.kinetics.systems.networks.StoichiometricNetwork’s structural-array-and-factory-closure pattern, but for bonded-force evaluation instead of mass-action reaction rates. The angle-bending force is derived from the chain rule \(\vec F_i=-dU/d\theta\cdot d\theta/d\vec r_i\) applied to \(\theta=\arccos(\hat v_1\cdot\hat v_2)\), a standard force-field result (e.g. Leach, Molecular Modelling: Principles and Applications, 2nd ed., Ch. 4).

Parameters:

  • positions (array-like, shape (n_atoms, 3))

  • velocities (array-like, shape (n_atoms, 3))

  • masses (array-like, shape (n_atoms,))

  • bonds (Sequence) – Harmonic bonds between atom indices i, j with force constant k and equilibrium length r0.

  • angles (Sequence) – Harmonic angle terms at vertex atom j (between bonds i-j and k-j) with force constant k_theta and equilibrium angle theta0 (radians).

forces_and_potential(positions)[source]#

Return the force on every particle and the total potential energy.

Parameters:

positions (ndarray, shape (n_particles, n_dim))

Returns:

  • forces (ndarray, shape (n_particles, n_dim))

  • potential_energy (float)

class chemistrykit.md.LJFluid(positions, velocities, box_length, epsilon=1.0, sigma=1.0, mass=1.0, cutoff=2.5, skin=0.3, rebuild_every=20)[source]#

Bases: MolecularDynamicsSystem

A monatomic Lennard-Jones fluid in a cubic periodic box, in reduced units.

Positions evolve under chemistrykit.integrators.velocity_verlet_step() driven by a @njit pairwise-force loop (_make_lj_accel_njit()) that reads its interacting pairs from a periodically-rebuilt VerletNeighborList, under the minimum-image convention (chemistrykit.md.utils.pbc.minimum_image_displacement()). See the module docstring for the reduced-units convention.

Parameters:

  • positions (array-like, shape (N, 3))

  • velocities (array-like, shape (N, 3))

  • box_length (float) – Cubic periodic box side length; must be at least 2 * cutoff for the minimum-image convention to be unambiguous.

  • epsilon (float) – Lennard-Jones/particle parameters (all particles identical).

  • sigma (float) – Lennard-Jones/particle parameters (all particles identical).

  • mass (float) – Lennard-Jones/particle parameters (all particles identical).

  • cutoff (float) – Interaction cutoff, in units of sigma (2.5 is the conventional standard truncation for a liquid-density LJ fluid; Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., Ch. 1.4).

  • skin (float) – Neighbor-list skin distance, see VerletNeighborList.

  • rebuild_every (int) – Rebuild the neighbor list after this many step() calls.

forces_and_potential(positions)[source]#

Exact (full, non-skin-list) forces, energy, and (cached) virial.

Unlike step()’s force evaluation (which reuses a skin-enlarged neighbor list for performance), this always rebuilds an exact-cutoff pair list from the given positions, so energy, pressure, and g(r) queries are never affected by neighbor-list staleness.

Parameters:

positions (ndarray, shape (N, 3))

Returns:

  • forces (ndarray, shape (N, 3))

  • potential_energy (float)

classmethod from_lattice(n_per_side, density, temperature, epsilon=1.0, sigma=1.0, mass=1.0, cutoff=2.5, rng=None, **kwargs)[source]#

Build an LJFluid on a simple-cubic lattice with Maxwell-Boltzmann velocities.

A standard way to initialize an MD run: a regular lattice avoids any risk of initial particle overlap (which would otherwise blow up the steeply repulsive \(r^{-12}\) term), and velocity components drawn independently from a Gaussian with variance \(k_BT/m\) reproduce the Maxwell-Boltzmann distribution exactly (see chemistrykit.statmech.MaxwellBoltzmannSpeedDistribution); the total momentum is then zeroed and the velocities rescaled so the realized instantaneous temperature matches temperature exactly, rather than only on average.

Parameters:

  • n_per_side (int) – Particles per lattice edge; total particle count is n_per_side**3.

  • density (float) – Number density \(\rho=N/V\), in units of \(\sigma^{-3}\).

  • temperature (float) – Target temperature, in units of \(\epsilon/k_B\).

  • epsilon (float)

  • sigma (float)

  • mass (float)

  • cutoff (float)

  • rng (int, numpy.random.Generator, or None) – Seed or generator for the initial velocities.

  • **kwargs – Forwarded to the LJFluid constructor (e.g. skin, rebuild_every).

Return type:

LJFluid

Returns:

LJFluid

Examples

>>> fluid = LJFluid.from_lattice(n_per_side=4, density=0.5, temperature=1.0, rng=0)
>>> fluid.positions.shape
(64, 3)
>>> round(fluid.temperature(), 6)
1.0
pressure(k_b=1.0)[source]#

Instantaneous pressure via the virial theorem.

\[P = \frac{N k_B T}{V} + \frac{W}{dV}\]

(Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., eq. 2.60, generalized from \(d=3\) to general dimension.)

Parameters:

k_b (float)

Return type:

float

Returns:

float

Examples

A very dilute LJ gas (large box, few particles) has pressure close to the ideal-gas value \(P=\rho k_B T\), since the virial correction from interactions vanishes as the density does:

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> n, L = 20, 40.0
>>> positions = rng.uniform(0.0, L, size=(n, 3))
>>> velocities = rng.normal(0.0, 1.0, size=(n, 3))
>>> velocities -= velocities.mean(axis=0)
>>> fluid = LJFluid(positions, velocities, box_length=L, cutoff=2.5)
>>> fluid.virial()  # no pairs within cutoff at this density
0.0
>>> rho = n / L**3
>>> ideal_pressure = rho * fluid.temperature()
>>> bool(abs(fluid.pressure() - ideal_pressure) < 1e-12)
True
radial_distribution_function(r_max=None, n_bins=100, positions=None)[source]#

Compute the radial distribution function g(r) for a configuration.

\[g(r) = \frac{2\,\langle n_{\text{pairs}}(r, r+dr)\rangle}{N\,\rho\,4\pi r^2 dr}\]

the (3D) local density of particles at distance r from a reference particle, relative to the bulk density \(\rho\) (Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., Ch. 2.6; Frenkel & Smit, Understanding Molecular Simulation, 2nd ed., Ch. 3.6). \(g(r) \to 1\) for an ideal (non-interacting, spatially uniform) gas, and develops structure (a first peak near the LJ minimum, decaying oscillations beyond it) in a liquid or solid.

Parameters:

  • r_max (float | None) – Maximum separation to histogram; defaults to half the box length (beyond which the minimum-image convention is ambiguous).

  • n_bins (int)

  • positions (ndarray, shape (N, 3), optional) – Configuration to analyze; defaults to the system’s current self.positions.

Returns:

  • r (ndarray, shape (n_bins,)) – Bin-center separations.

  • g (ndarray, shape (n_bins,))

step(dt)[source]#

Advance the stored state by one velocity-Verlet step, in place.

Delegates the actual kick-drift-kick arithmetic to the shared chemistrykit.integrators.velocity_verlet_step() – this method exists only to supply that function with this system’s _accel_njit/params, exactly as chemistrykit.kinetics.core.base_system.ReactionNetwork.integrate() wraps chemistrykit.integrators.rk4_integrate().

Parameters:

dt (float) – Time step.

Return type:

None

virial()[source]#

The virial sum \(W=\sum_{i<j}\vec r_{ij}\cdot\vec F_{ij}\) at the current configuration.

Return type:

float

Returns:

float

class chemistrykit.md.LennardJones(epsilon=1.0, sigma=1.0)[source]#

Bases: PairPotential

The 12-6 Lennard-Jones pair potential.

\[U(r) = 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]\]

The \(r^{-12}\) term models short-range Pauli repulsion (steeply, for computational convenience – it has no deeper physical justification), and the \(r^{-6}\) term is the genuine London-dispersion attraction (J. E. Lennard-Jones, Proc. R. Soc. A 106, 463 (1924); Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., Ch. 1.4).

Parameters:

  • epsilon (float) – Well depth.

  • sigma (float) – Finite distance at which U(sigma) == 0.

Examples

The potential has its minimum exactly at \(r=2^{1/6}\sigma\) (r_min), where \(U=-\epsilon\) and the force vanishes – the mechanical-equilibrium separation of an isolated pair:

>>> lj = LennardJones(epsilon=1.0, sigma=1.0)
>>> round(float(lj.energy(lj.r_min)), 10)
-1.0
>>> round(float(lj.force_scalar(lj.r_min)), 10)
-0.0

At r=sigma the potential is exactly zero by construction:

>>> round(float(lj.energy(1.0)), 10)
0.0
energy(r)[source]#

Return the potential energy U(r) at separation(s) r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

force_scalar(r)[source]#

Return the radial force magnitude \(f(r) = -dU/dr\).

By convention, positive is repulsive (pushes the two particles apart along their separation vector \(\vec r_i - \vec r_j\)) and negative is attractive – so the force on particle i due to particle j is f(r) * (r_i - r_j) / r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

property r_min: float#

Separation at the potential minimum, \(2^{1/6}\sigma\).

class chemistrykit.md.MDResult(t, positions, velocities, kinetic_energy, potential_energy, temperature, box_length=None, extra=<factory>)[source]#

Bases: object

Container for the output of a MolecularDynamicsSystem.run() call.

Mirrors chemistrykit.kinetics.core.base_system.KineticsResult (a stable, dataclass return type consumed by visualizers), specialized to a particle trajectory instead of a concentration one.

Parameters:
box_length: float | None = None#

Cubic periodic box side length, or None for a non-periodic system.

Type:

float or None

extra: dict#

Free-form slot for any additional diagnostics a system chooses to attach (e.g. pressure).

Type:

dict

kinetic_energy: ndarray#

Total kinetic energy at each frame.

Type:

ndarray, shape (n_frames,)

positions: ndarray#

Particle positions at each sampled frame.

Type:

ndarray, shape (n_frames, n_particles, n_dim)

potential_energy: ndarray#

Total potential energy at each frame.

Type:

ndarray, shape (n_frames,)

speeds()[source]#

Return per-particle speeds at every sampled frame.

Return type:

ndarray

Returns:

ndarray, shape (n_frames, n_particles)

Examples

>>> import numpy as np
>>> result = MDResult(
...     t=np.array([0.0]),
...     positions=np.zeros((1, 2, 3)),
...     velocities=np.array([[[3.0, 4.0, 0.0], [0.0, 0.0, 0.0]]]),
...     kinetic_energy=np.array([0.0]),
...     potential_energy=np.array([0.0]),
...     temperature=np.array([0.0]),
... )
>>> result.speeds()
array([[5., 0.]])
t: ndarray#

Time at each sampled frame.

Type:

ndarray, shape (n_frames,)

temperature: ndarray#

Instantaneous kinetic temperature (see MolecularDynamicsSystem.temperature()) at each frame.

Type:

ndarray, shape (n_frames,)

property total_energy: ndarray#

kinetic_energy + potential_energy.

For a system with no thermostat attached this should be constant to within integration error – the standard NVE energy-conservation check for a symplectic integrator like velocity-Verlet.

Type:

ndarray, shape (n_frames,)

velocities: ndarray#

Particle velocities at each sampled frame.

Type:

ndarray, shape (n_frames, n_particles, n_dim)

class chemistrykit.md.MolecularDynamicsSystem(positions, velocities, masses)[source]#

Bases: ABC

Common engine for a system of particles integrated by velocity-Verlet.

Concrete subclasses must set, in __init__ (before calling super().__init__(...)):

  • self._accel_njit : an @njit dispatcher with signature (positions, t, params) -> acceleration (the chemistrykit.integrators.RHSFunc convention, specialized so the “state” is a position array and the return value is acceleration, i.e. force already divided by mass – exactly the convention chemistrykit.integrators.velocity_verlet_step() expects).

  • self.params : ndarray – the parameter vector read by self._accel_njit (np.empty(0) if every numeric parameter is instead baked into the njit closure itself).

  • self.box_length : float or None – cubic periodic box side length, for systems that use periodic boundary conditions (None for a finite, non-periodic cluster).

and call super().__init__(positions, velocities, masses).

Parameters:

  • positions (array-like, shape (n_particles, n_dim))

  • velocities (array-like, shape (n_particles, n_dim))

  • masses (array-like, shape (n_particles,))

box_length: float | None = None#
degrees_of_freedom()[source]#

Number of kinetic degrees of freedom, assuming zero net momentum.

n_dim * n_particles - n_dim: total Cartesian degrees of freedom minus the n_dim removed by fixing the center-of-mass velocity – standard practice whenever the total momentum is conserved/zeroed at initialization (Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., Ch. 2.6). Subclasses with different constraints (e.g. a fixed number of internal vibrational modes) should override this.

Return type:

int

Returns:

int

abstractmethod forces_and_potential(positions)[source]#

Return the force on every particle and the total potential energy.

Parameters:

positions (ndarray, shape (n_particles, n_dim))

Returns:

  • forces (ndarray, shape (n_particles, n_dim))

  • potential_energy (float)

kinetic_energy()[source]#

Total kinetic energy, \(K = \sum_i \frac{1}{2} m_i v_i^2\).

Return type:

float

Returns:

float

params: ndarray = array([], dtype=float64)#
potential_energy()[source]#

Total potential energy at the system’s current configuration.

Return type:

float

Returns:

float

run(dt, n_steps, thermostat=None, sample_every=1)[source]#

Integrate forward n_steps velocity-Verlet steps, sampling as it goes.

Parameters:

  • dt (float) – Time step.

  • n_steps (int) – Number of integration steps (must be a multiple of sample_every).

  • thermostat (object, optional) – An object with an apply(system, dt) method, called once after every step (e.g. VelocityRescalingThermostat, NoseHooverThermostat). Omit for a microcanonical (NVE) run.

  • sample_every (int) – Record a frame every this many steps (plus the initial state).

Return type:

MDResult

Returns:

MDResult

step(dt)[source]#

Advance the stored state by one velocity-Verlet step, in place.

Delegates the actual kick-drift-kick arithmetic to the shared chemistrykit.integrators.velocity_verlet_step() – this method exists only to supply that function with this system’s _accel_njit/params, exactly as chemistrykit.kinetics.core.base_system.ReactionNetwork.integrate() wraps chemistrykit.integrators.rk4_integrate().

Parameters:

dt (float) – Time step.

Return type:

None

temperature(k_b=1.0)[source]#

Instantaneous kinetic temperature via the equipartition theorem.

\[T = \frac{2K}{k_B \times \text{dof}}\]

(Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., eq. 2.60.) In the reduced units used by LJFluid, k_b=1.0 (the default) is the correct choice.

Parameters:

k_b (float) – Boltzmann constant in the system’s unit system.

Return type:

float

Returns:

float

class chemistrykit.md.Morse(De, a, re)[source]#

Bases: PairPotential

The Morse interatomic potential.

\[U(r) = D_e\left[1 - e^{-a(r-r_e)}\right]^2 - D_e\]

A three-parameter anharmonic bond potential: dissociation energy \(D_e\), equilibrium separation \(r_e\), and a width parameter a controlling the curvature at the minimum (related to the harmonic force constant there by \(k = 2D_ea^2\), see force_constant) – P. M. Morse, Phys. Rev. 34, 57 (1929); Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12 (used there for vibrational anharmonicity; see also DiatomicOscillator, whose classical vibration this potential drives).

Parameters:

  • De (float) – Dissociation energy (well depth), must be positive.

  • a (float) – Width parameter (inverse length), must be positive.

  • re (float) – Equilibrium separation.

Examples

Exactly reproduces the well depth and zero force at the minimum:

>>> morse = Morse(De=4.0, a=1.5, re=1.0)
>>> round(float(morse.energy(morse.re)), 10)
-4.0
>>> round(float(morse.force_scalar(morse.re)), 10)
0.0

And the dissociation limit \(U(r\to\infty) \to 0\):

>>> round(float(morse.energy(20.0)), 6)
-0.0
energy(r)[source]#

Return the potential energy U(r) at separation(s) r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

property force_constant: float#

Harmonic force constant at the minimum, \(k=2D_ea^2\).

force_scalar(r)[source]#

Return the radial force magnitude \(f(r) = -dU/dr\).

By convention, positive is repulsive (pushes the two particles apart along their separation vector \(\vec r_i - \vec r_j\)) and negative is attractive – so the force on particle i due to particle j is f(r) * (r_i - r_j) / r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

classmethod from_force_constant(De, k, re)[source]#

Build a Morse potential from \(D_e\), a harmonic force constant k, and re.

Inverts \(k=2D_ea^2\) for \(a=\sqrt{k/(2D_e)}\), so the curvature at the minimum matches a given harmonic bond exactly while still capturing anharmonicity/dissociation away from it.

Parameters:
Return type:

Morse

Returns:

Morse

Examples

>>> morse = Morse.from_force_constant(De=4.0, k=10.0, re=1.0)
>>> round(morse.force_constant, 6)
10.0
class chemistrykit.md.NoseHooverThermostat(target_temperature, Q, dof=None, k_b=1.0)[source]#

Bases: object

The Nose-Hoover extended-system thermostat.

Couples the physical system to a fictitious “friction” degree of freedom \(\xi\) whose equation of motion drives the time-averaged kinetic temperature to target_temperature while (in the exact, continuous-time extended Lagrangian) still generating genuine canonical-ensemble fluctuations – S. Nose, J. Chem. Phys. 81, 511 (1984); W. G. Hoover, Phys. Rev. A 31, 1695 (1985).

\[\dot{\vec v}_i = \frac{\vec F_i}{m_i} - \xi \vec v_i, \qquad Q\dot\xi = \left(\sum_i m_i v_i^2\right) - \text{dof}\cdot k_BT\]

This implementation integrates the friction/velocity-scaling part with the explicit, first-order update

\[\xi \mathrel{+}= \frac{2K - \text{dof}\cdot k_BT}{Q}\,dt, \qquad \vec v_i \mathrel{*}= e^{-\xi\,dt}\]

applied once after each ordinary velocity-Verlet step – an approximation to the reference operator-split (“Trotter factorization”) integrator of Martyna, Tuckerman, Tobias & Klein, Mol. Phys. 87, 1117 (1996) (see also Frenkel & Smit, Understanding Molecular Simulation, 2nd ed., Ch. 6.1.2, Algorithm 15), which instead splits the thermostat update into two half-steps symmetrically wrapped around the position/velocity update for improved (time-reversible) long-time behavior. The simplified, single-step coupling used here still relaxes the time-averaged kinetic temperature to target_temperature, but without that scheme’s more favorable long-time energy-conservation properties – flagged explicitly here as the approximation it is.

Parameters:

  • target_temperature (float) – Must be positive.

  • Q (float) – Thermostat “mass” (coupling strength); a larger Q couples more weakly, giving a slower thermostat response. Must be positive.

  • dof (int) – Degrees of freedom used for the target kinetic energy \(\text{dof}\cdot k_BT/2\); defaults to the system’s own degrees_of_freedom() at the first call to apply().

  • k_b (float)

Examples

>>> import numpy as np
>>> from chemistrykit.md.systems.lj_fluid import LJFluid
>>> fluid = LJFluid.from_lattice(n_per_side=3, cutoff=1.5, density=0.6, temperature=2.0, rng=0)
>>> fluid.velocities *= 2.0  # perturb away from the target temperature
>>> thermostat = NoseHooverThermostat(target_temperature=2.0, Q=10.0)
>>> result = fluid.run(dt=0.001, n_steps=200, thermostat=thermostat, sample_every=20)
>>> bool(abs(result.temperature[-1] - 2.0) < abs(result.temperature[0] - 2.0))
True
apply(system, dt)[source]#

Advance the friction coefficient and rescale system.velocities in place.

Parameters:
Return type:

None

class chemistrykit.md.PairPotential[source]#

Bases: ABC

Common interface for a radially symmetric two-body potential U(r).

Concrete subclasses (LennardJones, Morse, Buckingham, HarmonicBond) implement energy() and force_scalar() so they can be swapped in and compared directly, mirroring the role chemistrykit.thermo.core.base_system.EquationOfState plays for equations of state.

abstractmethod energy(r)[source]#

Return the potential energy U(r) at separation(s) r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

abstractmethod force_scalar(r)[source]#

Return the radial force magnitude \(f(r) = -dU/dr\).

By convention, positive is repulsive (pushes the two particles apart along their separation vector \(\vec r_i - \vec r_j\)) and negative is attractive – so the force on particle i due to particle j is f(r) * (r_i - r_j) / r.

Parameters:

r (float or array-like of float)

Returns:

float or ndarray

class chemistrykit.md.ShakeMolecule(positions, velocities, masses, constraints, bonds=(), angles=(), tol=1e-10)[source]#

Bases: HarmonicMolecule

A HarmonicMolecule with rigid (SHAKE) bonds.

Each velocity-Verlet step makes an unconstrained half-kick and drift, corrects the positions with shake(), recovers the half-step velocities from the constrained displacement, completes the second half-kick, and finally removes each constrained bond’s relative velocity component (RATTLE’s velocity stage).

Parameters:

Examples

A rigid rotating diatomic keeps its bond length exactly:

>>> import numpy as np
>>> pos = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
>>> vel = np.array([[0.0, -1.0, 0.0], [0.0, 1.0, 0.0]])
>>> rotor = ShakeMolecule(pos, vel, [1.0, 1.0], constraints=[(0, 1, 1.0)])
>>> for _ in range(500):
...     rotor.step(0.01)
>>> round(float(np.linalg.norm(rotor.positions[0] - rotor.positions[1])), 8)
1.0
bond_lengths()[source]#

Current length of every constrained bond.

Return type:

ndarray

Returns:

ndarray, shape (n_constraints,)

degrees_of_freedom()[source]#

Cartesian degrees of freedom minus center-of-mass motion and one per constraint.

Return type:

int

Returns:

int

step(dt)[source]#

Advance one constrained velocity-Verlet (SHAKE/RATTLE) step, in place.

Parameters:

dt (float)

Return type:

None

class chemistrykit.md.StochasticVelocityRescalingThermostat(target_temperature, tau, rng=None, dof=None, k_b=1.0)[source]#

Bases: object

The Bussi-Donadio-Parrinello stochastic velocity-rescaling thermostat.

Berendsen’s exponential relaxation of the kinetic energy \(K\) plus a correctly scaled stochastic term, chosen so that the kinetic energy samples its exact canonical distribution – G. Bussi, D. Donadio & M. Parrinello, J. Chem. Phys. 126, 014101 (2007), eq. 7 and Appendix (eq. A7):

\[\alpha^2 = c + (1-c)\frac{\bar K}{N_f K}\sum_{i=1}^{N_f} R_i^2 + 2 R_1 \sqrt{c(1-c)\frac{\bar K}{N_f K}}, \qquad c = e^{-\Delta t/\tau},\]

where \(\bar K = N_f k_B T_0/2\), the \(R_i\) are independent standard normal numbers (the sum over \(i \ge 2\) is drawn as a single \(\chi^2_{N_f-1}\) variate), and every velocity is scaled by \(\alpha\). The stationary distribution of \(K\) is the canonical gamma distribution with mean \(\bar K\) and variance \(2\bar K^2/N_f\).

Parameters:

  • target_temperature (float) – \(T_0\); must be positive.

  • tau (float) – Relaxation time; must be positive.

  • rng (int, numpy.random.Generator, or None)

  • dof (int) – \(N_f\); defaults to the system’s degrees_of_freedom().

  • k_b (float)

Examples

>>> from chemistrykit.md.systems.lj_fluid import LJFluid
>>> fluid = LJFluid.from_lattice(n_per_side=3, cutoff=1.5, density=0.6, temperature=3.0, rng=0)
>>> thermostat = StochasticVelocityRescalingThermostat(target_temperature=1.0, tau=0.01, rng=0)
>>> for _ in range(200):
...     thermostat.apply(fluid, dt=0.01)
>>> bool(0.5 < fluid.temperature() < 1.5)
True
apply(system, dt)[source]#

Rescale system.velocities in place by a stochastic factor.

Parameters:
Return type:

None

class chemistrykit.md.VelocityRescalingThermostat(target_temperature, interval=1, k_b=1.0)[source]#

Bases: object

Simple (deterministic) velocity-rescaling thermostat.

Rescales every particle’s velocity by a common factor \(\lambda=\sqrt{T_{\text{target}}/T_{\text{current}}}\) every interval MD steps, so the instantaneous kinetic temperature matches target_temperature exactly at each rescaling (Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., Ch. 2.3; Frenkel & Smit, Understanding Molecular Simulation, 2nd ed., Ch. 6.1.1).

This is the crude, deterministic member of the “velocity rescaling” thermostat family: it does not itself sample the canonical (NVT) ensemble correctly (the kinetic-energy fluctuations of a true canonical ensemble are suppressed to zero at every rescaling step), unlike the stochastic Bussi-Donadio-Parrinello (2007) variant, or (to leading order, with a finite relaxation time) the Berendsen (1984) weak-coupling thermostat. It is nonetheless the standard first introduction to temperature control in MD and is exactly what is meant by “velocity rescaling” here.

Parameters:

  • target_temperature (float) – Must be positive.

  • interval (int) – Rescale every this many calls to apply() (i.e. MD steps).

  • k_b (float) – Boltzmann constant in the system’s unit system (1.0 for the reduced units used by LJFluid).

Examples

>>> import numpy as np
>>> from chemistrykit.md.systems.lj_fluid import LJFluid
>>> fluid = LJFluid.from_lattice(n_per_side=3, cutoff=1.5, density=0.6, temperature=2.0, rng=0)
>>> fluid.velocities *= 2.0  # perturb away from the target temperature
>>> thermostat = VelocityRescalingThermostat(target_temperature=2.0)
>>> result = fluid.run(dt=0.001, n_steps=50, thermostat=thermostat, sample_every=10)
>>> bool(abs(result.temperature[-1] - 2.0) < 0.5)
True
apply(system, dt)[source]#

Rescale system.velocities in place if a rescaling is due.

Parameters:
Return type:

None

class chemistrykit.md.VerletNeighborList(cutoff, skin=0.3, rebuild_every=20)[source]#

Bases: object

A Verlet (“skin”) neighbor list, rebuilt only every few calls.

Building the pair list with a cutoff enlarged by a skin distance means the list stays a superset of the true cutoff-pairs for several integration steps, since a pair can only enter the true cutoff sphere after its separation has changed by up to skin – the standard MD performance trick (Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., Ch. 5.3). This implementation rebuilds unconditionally every rebuild_every calls rather than tracking accumulated displacement, which is simpler but requires the caller to choose rebuild_every/skin conservatively for their timestep and temperature – a stricter implementation would instead track the two largest particle displacements since the last rebuild and force an early rebuild once their sum exceeds skin.

Parameters:

  • cutoff (float) – The physical interaction cutoff.

  • skin (float) – Extra distance added to cutoff when building the list.

  • rebuild_every (int) – Rebuild after this many calls to pairs().

pairs(positions, box_length)[source]#

Return (possibly cached) (pairs_i, pairs_j) for positions.

Parameters:

  • positions (ndarray, shape (N, d))

  • box_length (float | None)

Returns:

pairs_i, pairs_j (ndarray of int64)

reset()[source]#

Force the next call to pairs() to rebuild.

Return type:

None

chemistrykit.md.build_neighbor_list(positions, box_length, cutoff)[source]#

Build the list of all particle pairs within cutoff of each other.

A direct, vectorized \(O(N^2)\) all-pairs distance computation under the minimum-image convention – exact, and fast enough for the particle counts this package targets (up to a few thousand). A production MD code instead uses a cell (linked-list) decomposition to reduce this to \(O(N)\); that further optimization is out of scope here. See VerletNeighborList for the one optimization this module does provide: reusing a list built with a distance “skin” across several integration steps instead of rebuilding every step.

Parameters:

  • positions (ndarray, shape (N, d))

  • box_length (float | None) – Cubic periodic box side length; None or <= 0 for an unbounded (non-periodic) system.

  • cutoff (float) – Pairs farther apart than this are excluded.

Returns:

pairs_i, pairs_j (ndarray of int64, shape (n_pairs,)) – Indices such that pairs_i < pairs_j for every returned pair.

Examples

>>> import numpy as np
>>> positions = np.array([[0.0, 0.0, 0.0], [0.5, 0.0, 0.0], [5.0, 0.0, 0.0]])
>>> pairs_i, pairs_j = build_neighbor_list(positions, box_length=None, cutoff=1.0)
>>> list(zip(pairs_i.tolist(), pairs_j.tolist()))
[(0, 1)]
chemistrykit.md.einstein_diffusion_coefficient(t, msd, n_dim=3, fit_from=0.0)[source]#

Self-diffusion coefficient from the Einstein relation \(\mathrm{MSD}=2dDt\).

Fits a straight line to msd against t for t >= fit_from (skipping the short-time ballistic regime) and returns slope/(2d).

Parameters:
Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> t = np.linspace(0.0, 10.0, 11)
>>> round(einstein_diffusion_coefficient(t, 6.0 * 0.25 * t + 0.1), 12)
0.25
chemistrykit.md.green_kubo_diffusion_coefficient(t, vacf, n_dim=3)[source]#

Self-diffusion coefficient from the Green-Kubo integral of the VACF.

\[D = \frac{1}{d}\int_0^{t_{\max}} \langle\vec v(0)\cdot\vec v(t)\rangle\,dt\]

evaluated by the trapezoidal rule.

Parameters:

  • t (array-like) – Lag times, starting at 0.

  • vacf (array-like) – Unnormalized velocity autocorrelation (see velocity_autocorrelation()).

  • n_dim (int)

Return type:

float

Returns:

float

Examples

An exponentially decaying VACF \(C(t)=3(k_BT/m)e^{-\gamma t}\) (Langevin dynamics) integrates to the Stokes-Einstein-like value \(k_BT/(m\gamma)\):

>>> import numpy as np
>>> t = np.linspace(0.0, 40.0, 40001)
>>> round(green_kubo_diffusion_coefficient(t, 3.0 * 1.5 * np.exp(-2.0 * t)), 6)
0.75
chemistrykit.md.mean_squared_displacement(positions, box_length=None, max_lag=None)[source]#

Time-origin-averaged mean-squared displacement.

\[\mathrm{MSD}(\ell) = \left\langle |\vec r_i(t_0+\ell) - \vec r_i(t_0)|^2 \right\rangle_{i,\,t_0}\]
Parameters:
Returns:

ndarray, shape (max_lag + 1,) – MSD at lags 0, 1, ..., max_lag frames.

Examples

Particles moving ballistically at unit speed have MSD equal to the lag squared:

>>> import numpy as np
>>> t = np.arange(5.0)
>>> positions = np.stack([t, np.zeros(5), np.zeros(5)], axis=-1)[:, None, :]
>>> mean_squared_displacement(positions, max_lag=3)
array([0., 1., 4., 9.])
chemistrykit.md.minimum_image_displacement(r_i, r_j, box_length)[source]#

Displacement \(\vec r_i - \vec r_j\) under the minimum-image convention.

\[\Delta\vec r = (\vec r_i - \vec r_j) - L \, \mathrm{round}\!\left(\frac{\vec r_i - \vec r_j}{L}\right)\]

for a cubic periodic box of side L – the displacement to the closest periodic image of particle j relative to particle i, valid whenever L is at least twice the interaction cutoff (Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., Ch. 1.5.2). If box_length is None or <= 0, plain (non-periodic) displacement is returned unchanged.

Parameters:

  • r_i (array-like) – Broadcastable position array(s).

  • r_j (array-like) – Broadcastable position array(s).

  • box_length (float | None)

Returns:

ndarray

Examples

Two particles near opposite edges of a box of length 10 are actually close together once the box wraps around:

>>> import numpy as np
>>> r_i = np.array([0.5, 0.0, 0.0])
>>> r_j = np.array([9.5, 0.0, 0.0])
>>> minimum_image_displacement(r_i, r_j, box_length=10.0)
array([1., 0., 0.])

With no box, this is just the ordinary displacement:

>>> minimum_image_displacement(r_i, r_j, box_length=None)
array([-9.,  0.,  0.])
chemistrykit.md.shake(positions_new, positions_old, constraints, masses, tol=1e-10, max_iter=1000)[source]#

Correct unconstrained positions so every bond-length constraint holds (SHAKE).

Parameters:

  • positions_new (array-like, shape (n_atoms, n_dim)) – Positions after an unconstrained update.

  • positions_old (array-like, shape (n_atoms, n_dim)) – Positions at the previous step (which satisfy the constraints); their bond vectors set the correction directions.

  • constraints (sequence of (i, j, d)) – Atom pairs and their fixed separations.

  • masses (array-like, shape (n_atoms,))

  • tol (float) – Relative tolerance on \(|s^2 - d^2|/d^2\).

  • max_iter (int)

Returns:

  • positions (ndarray, shape (n_atoms, n_dim)) – Corrected positions.

  • n_iter (int) – Number of sweeps over the constraints needed to converge.

Examples

Two equal masses pulled apart to separation 1.2 are pulled back symmetrically to the constrained length 1.0:

>>> import numpy as np
>>> old = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
>>> new = np.array([[-0.1, 0.0, 0.0], [1.1, 0.0, 0.0]])
>>> fixed, _ = shake(new, old, [(0, 1, 1.0)], masses=[1.0, 1.0])
>>> np.round(fixed[:, 0], 10)
array([0., 1.])
chemistrykit.md.unwrap_trajectory(positions, box_length)[source]#

Undo periodic wrapping of a sampled trajectory.

Accumulates the minimum-image displacement between consecutive frames, so it is exact provided no particle moves more than half a box length between two samples.

Parameters:

  • positions (array-like, shape (n_frames, n_particles, n_dim))

  • box_length (float | None) – Cubic periodic box side length; None or <= 0 returns the positions unchanged.

Returns:

ndarray, shape (n_frames, n_particles, n_dim)

Examples

>>> import numpy as np
>>> wrapped = np.array([[[9.5]], [[0.3]], [[1.1]]])  # crosses the box edge at 10
>>> unwrap_trajectory(wrapped, box_length=10.0).ravel()
array([ 9.5, 10.3, 11.1])
chemistrykit.md.velocity_autocorrelation(velocities, max_lag=None, normalize=False)[source]#

Time-origin-averaged velocity autocorrelation function.

\[C(\ell) = \left\langle \vec v_i(t_0)\cdot\vec v_i(t_0+\ell) \right\rangle_{i,\,t_0}\]
Parameters:

  • velocities (array-like, shape (n_frames, n_particles, n_dim))

  • max_lag (int | None) – Largest lag, in frames; defaults to n_frames // 2.

  • normalize (bool) – Divide by \(C(0)\).

Returns:

ndarray, shape (max_lag + 1,)

Examples

>>> import numpy as np
>>> v = np.array([[[1.0, 0.0]], [[0.0, 1.0]], [[-1.0, 0.0]]])  # a velocity rotating by 90 degrees per frame
>>> velocity_autocorrelation(v, max_lag=2)
array([ 1.,  0., -1.])
chemistrykit.md.wrap_positions(positions, box_length)[source]#

Wrap positions into the primary periodic cell [0, box_length).

Parameters:

  • positions (array-like)

  • box_length (float | None) – If None or <= 0, positions is returned unchanged (as a float array).

Returns:

ndarray

Examples

>>> import numpy as np
>>> wrap_positions(np.array([-0.5, 10.2, 5.0]), box_length=10.0)
array([9.5, 0.2, 5. ])
On this page
Show Source