.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/statphys/molecular_dynamics/plot_alder_wainwright_molecular_dynamics.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_statphys_molecular_dynamics_plot_alder_wainwright_molecular_dynamics.py: Alder and Wainwright: molecular dynamics simulation ====================================================== In 1957 Alder and Wainwright took a different route from Monte Carlo: integrate Newton's equations for a few hundred interacting particles and read equilibrium properties directly off the trajectory. Their simulations found that hard spheres with no attraction at all freeze into a crystal when packed densely enough. Order came from geometry alone. This example runs the same kind of experiment with :class:`~physicskit.statphys.chapters.molecular_dynamics.LennardJonesGas`, which integrates :func:`~physicskit.statphys.core.md_engine.lj_forces` with the time-reversible Velocity Verlet scheme (:func:`~physicskit.statphys.core.md_engine.velocity_verlet_step`). It checks that the total energy is conserved, as a microcanonical simulation requires, and far better than with a naive Euler step. It then compares the pair distribution function :math:`g(r)` of a dilute gas with that of a dense system, whose sharp shells reveal a solid. .. GENERATED FROM PYTHON SOURCE LINES 25-49 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from physicskit.statphys.chapters.molecular_dynamics import LennardJonesGas from physicskit.statphys.core.md_engine import lj_forces N = 144 def make_gas(density, T0, seed=1957): box = np.sqrt(N / density) return LennardJonesGas(n_particles=N, box_size=box, temperature_init=T0, dt=0.004, initial_velocity_distribution="maxwell_boltzmann", seed=seed) def pair_distribution(positions, box, r_max, bins=120): d = positions[:, None, :] - positions[None, :, :] d -= box * np.round(d / box) r = np.sqrt((d**2).sum(-1))[np.triu_indices(len(positions), 1)] hist, edges = np.histogram(r, bins=bins, range=(0, r_max)) shell = np.pi * (edges[1:] ** 2 - edges[:-1] ** 2) ideal = shell * len(positions) * (len(positions) - 1) / 2 / box**2 return 0.5 * (edges[1:] + edges[:-1]), hist / ideal .. GENERATED FROM PYTHON SOURCE LINES 50-56 Energy conservation: Verlet vs Euler ---------------------------------------- Same dense system, same time step. The time-reversible Verlet energy stays within about a percent (the small steps come from pairs crossing the interaction cutoff); forward Euler pumps energy in until particles are thrown into each other's cores and the simulation explodes. .. GENERATED FROM PYTHON SOURCE LINES 56-76 .. code-block:: Python gas = make_gas(0.7, 0.6) E_verlet = [] for _ in range(400): gas.step(5) E_verlet.append(gas.total_energy()) E_verlet = np.array(E_verlet) euler = make_gas(0.7, 0.6) x, v = euler.positions.copy(), euler.velocities.copy() E_euler = [] for k in range(2000): F, U = lj_forces(x, euler.box_size) if k % 5 == 4: E_euler.append(0.5 * np.sum(v**2) + U) x = np.mod(x + v * euler.dt, euler.box_size) v = v + F * euler.dt E_euler = np.array(E_euler) print(f"Verlet: relative energy drift over 2000 steps {abs(E_verlet[-1] - E_verlet[0]) / abs(E_verlet[0]):.1e}") print(f"Euler: relative energy drift over 2000 steps {abs(E_euler[-1] - E_euler[0]) / abs(E_euler[0]):.1e}") .. rst-class:: sphx-glr-script-out .. code-block:: none Verlet: relative energy drift over 2000 steps 8.0e-03 Euler: relative energy drift over 2000 steps 4.8e+57 .. GENERATED FROM PYTHON SOURCE LINES 77-87 Structure from dynamics: gas vs solid ----------------------------------------- Equilibrate a dilute gas and a dense system, then average :math:`g(r)` over the trajectory. The gas shows one weak neighbour peak and then nothing. The dense system shows sharp shells at the spacings of a triangular crystal, :math:`r_1`, :math:`\sqrt3 r_1` and :math:`2r_1`. The dense system starts on a square lattice, which is not its lowest energy state; as it rearranges into a triangular one it releases energy and warms up, so its temperature is reported as measured. .. GENERATED FROM PYTHON SOURCE LINES 87-128 .. code-block:: Python results = {} for label, density, T0 in [("dilute gas, rho = 0.05", 0.05, 1.0), ("dense, rho = 0.9", 0.9, 0.3)]: sim = make_gas(density, T0) sim.step(1500) g_acc, temps = 0, [] for _ in range(40): sim.step(25) r, g = pair_distribution(sim.positions, sim.box_size, r_max=4.0) g_acc = g_acc + g / 40 temps.append(sim.temperature()) results[f"{label}, T = {np.mean(temps):.2f}"] = (r, g_acc, sim.positions.copy(), sim.box_size) print(f"{label}: kinetic temperature {np.mean(temps):.3f} +/- {np.std(temps):.3f}") dense_key = [key for key in results if key.startswith("dense")][0] r1 = results[dense_key][0][np.argmax(results[dense_key][1])] print(f"dense system: first peak r1 = {r1:.3f}; triangular-lattice shells expected at {r1 * np.sqrt(3):.3f} and {2 * r1:.3f}") fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 4)) ax1.semilogy(np.arange(len(E_verlet)) * 5 * gas.dt, np.abs(E_verlet / E_verlet[0] - 1) + 1e-16, color="steelblue", label="Velocity Verlet") ax1.semilogy(np.arange(len(E_euler)) * 5 * euler.dt, np.abs(E_euler / E_euler[0] - 1) + 1e-16, color="firebrick", label="forward Euler") ax1.set_xlabel("time") ax1.set_ylabel(r"$|E(t)/E(0) - 1|$") ax1.set_title("Microcanonical energy conservation") ax1.legend(fontsize=8) for label, (r, g, _, _) in results.items(): ax2.plot(r, g, label=label) for shell in (r1, np.sqrt(3) * r1, 2 * r1): ax2.axvline(shell, color="0.6", ls=":", lw=0.8) ax2.set_xlabel("r / sigma") ax2.set_ylabel("g(r)") ax2.set_title("Pair distribution function") ax2.legend(fontsize=8) _, _, pos, box = results[dense_key] ax3.plot(pos[:, 0], pos[:, 1], "o", ms=6, color="darkorchid") ax3.set_xlim(0, box) ax3.set_ylim(0, box) ax3.set_aspect("equal") ax3.set_title("Dense system: a crystal") fig.tight_layout() plt.show() .. image-sg:: /api/gallery/statphys/molecular_dynamics/images/sphx_glr_plot_alder_wainwright_molecular_dynamics_001.png :alt: Microcanonical energy conservation, Pair distribution function, Dense system: a crystal :srcset: /api/gallery/statphys/molecular_dynamics/images/sphx_glr_plot_alder_wainwright_molecular_dynamics_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none dilute gas, rho = 0.05: kinetic temperature 0.989 +/- 0.018 dense, rho = 0.9: kinetic temperature 0.795 +/- 0.041 dense system: first peak r1 = 1.083; triangular-lattice shells expected at 1.876 and 2.167 .. GENERATED FROM PYTHON SOURCE LINES 129-134 Check ----- Velocity Verlet conserves energy to < 2%; forward Euler explodes. The dilute gas keeps its temperature; the dense system crystallizes with its first neighbour at the triangular-lattice spacing sqrt(2 / (sqrt(3) rho)). .. GENERATED FROM PYTHON SOURCE LINES 134-137 .. code-block:: Python assert abs(E_verlet[-1] - E_verlet[0]) / abs(E_verlet[0]) < 0.02 assert abs(E_euler[-1] - E_euler[0]) / abs(E_euler[0]) > 1 assert abs(r1 - np.sqrt(2 / (np.sqrt(3) * 0.9))) < 0.07 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.433 seconds) .. _sphx_glr_download_api_gallery_statphys_molecular_dynamics_plot_alder_wainwright_molecular_dynamics.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_alder_wainwright_molecular_dynamics.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_alder_wainwright_molecular_dynamics.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_alder_wainwright_molecular_dynamics.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_