.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/dynamics/02_velocity_verlet_integration.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_examples_dynamics_02_velocity_verlet_integration.py: Velocity Verlet Dynamics with Lennard-Jones Potential ===================================================== This example demonstrates GPU-accelerated molecular dynamics using the Velocity Verlet integrator with a Lennard-Jones (LJ) potential. We simulate argon in the **NVE** ensemble (no thermostat) to demonstrate: 1. Building an FCC lattice system 2. Computing neighbor lists with periodic boundaries 3. Using the LJ potential for energy and forces 4. Running Velocity Verlet dynamics 5. Monitoring energy conservation and basic stability diagnostics Notes ----- Velocity Verlet is stable and time-reversible, but for LJ systems it is still sensitive to: - timestep size, - neighbor list rebuild logic, - force discontinuities at the cutoff (consider using `switch_width > 0`). .. GENERATED FROM PYTHON SOURCE LINES 41-43 Imports ------- .. GENERATED FROM PYTHON SOURCE LINES 43-59 .. code-block:: Python from __future__ import annotations import matplotlib.pyplot as plt import numpy as np import warp as wp from _dynamics_utils import ( DEFAULT_CUTOFF, DEFAULT_SKIN, EPSILON_AR, SIGMA_AR, MDSystem, create_fcc_argon, run_velocity_verlet, ) .. GENERATED FROM PYTHON SOURCE LINES 60-62 System setup ------------ .. GENERATED FROM PYTHON SOURCE LINES 62-92 .. code-block:: Python print("=" * 95) print("VELOCITY VERLET DYNAMICS WITH LENNARD-JONES POTENTIAL (NVE)") print("=" * 95) print() device = "cuda:0" if wp.is_cuda_available() else "cpu" print(f"Using device: {device}") print("\n--- Creating FCC Argon System ---") positions, cell = create_fcc_argon(num_unit_cells=4, a=5.26) # 256 atoms print(f"Created {len(positions)} atoms in {cell[0, 0]:.2f} ų box") print("\n--- Initializing MD System ---") system = MDSystem( positions=positions, cell=cell, epsilon=EPSILON_AR, sigma=SIGMA_AR, cutoff=DEFAULT_CUTOFF, skin=DEFAULT_SKIN, # NVE can benefit from smooth cutoffs; keep hard cutoff by default for now. switch_width=0.0, device=device, dtype=np.float64, ) print("\n--- Setting Initial Temperature ---") system.initialize_temperature(temperature=94.4, seed=42) .. rst-class:: sphx-glr-script-out .. code-block:: none =============================================================================================== VELOCITY VERLET DYNAMICS WITH LENNARD-JONES POTENTIAL (NVE) =============================================================================================== Using device: cuda:0 --- Creating FCC Argon System --- Created 256 atoms in 21.04 ų box --- Initializing MD System --- Initialized MD system with 256 atoms Cell: 21.04 x 21.04 x 21.04 Å Cutoff: 8.50 Å (+ 0.50 Å skin) LJ: ε = 0.0104 eV, σ = 3.40 Å Device: cuda:0, dtype: Units: x [Å], t [fs], E [eV], m [eV·fs²/Ų] (from amu), v [Å/fs] --- Setting Initial Temperature --- Initialized velocities: target=94.4 K, actual=90.4 K .. GENERATED FROM PYTHON SOURCE LINES 93-95 Velocity Verlet (NVE) --------------------- .. GENERATED FROM PYTHON SOURCE LINES 95-104 .. code-block:: Python print("\n--- NVE Run (3000 steps) ---") stats = run_velocity_verlet( system=system, num_steps=3000, dt_fs=1.0, log_interval=200, ) .. rst-class:: sphx-glr-script-out .. code-block:: none --- NVE Run (3000 steps) --- Running 3000 Velocity Verlet steps (NVE) dt = 1.000 fs =============================================================================================== Step KE (eV) PE (eV) Total (eV) T (K) Neighbors min r (Å) max|F| =============================================================================================== 0 2.9787 -21.5621 -18.5833 90.37 9984 3.713 2.313e-03 200 1.2443 -19.7214 -18.4771 37.75 10304 3.324 1.422e-01 400 1.5604 -20.0374 -18.4771 47.34 10256 3.313 1.771e-01 600 1.3336 -19.8084 -18.4749 40.46 10279 3.344 1.640e-01 800 1.5552 -20.0421 -18.4869 47.18 10253 3.239 1.964e-01 1000 1.5486 -20.0457 -18.4971 46.98 10265 3.316 1.433e-01 1200 1.5071 -19.9964 -18.4893 45.72 10266 3.322 1.814e-01 1400 1.4764 -19.9599 -18.4835 44.79 10275 3.253 1.971e-01 1600 1.5812 -20.0671 -18.4859 47.97 10245 3.318 1.719e-01 1800 1.5937 -20.0755 -18.4818 48.35 10269 3.317 1.582e-01 2000 1.5902 -20.0884 -18.4983 48.24 10268 3.351 1.498e-01 2200 1.4885 -19.9817 -18.4932 45.16 10260 3.302 1.611e-01 2400 1.4969 -19.9773 -18.4805 45.41 10257 3.331 1.949e-01 2600 1.4802 -19.9617 -18.4815 44.91 10267 3.296 1.753e-01 2800 1.4200 -19.9046 -18.4845 43.08 10257 3.292 1.663e-01 2999 1.5927 -20.0865 -18.4939 48.32 10252 3.272 1.715e-01 .. GENERATED FROM PYTHON SOURCE LINES 105-107 Analysis -------- .. GENERATED FROM PYTHON SOURCE LINES 107-135 .. code-block:: Python print("\n--- Analysis ---") temps = np.array([s.temperature for s in stats]) total_energies = np.array([s.total_energy for s in stats]) steps = np.array([s.step for s in stats]) e0 = total_energies[0] drift = total_energies - e0 print(f" Mean Temperature: {temps.mean():.2f} ± {temps.std():.2f} K") print(f" Mean Total Energy: {total_energies.mean():.4f} eV") print(f" Energy Fluctuation: {total_energies.std():.4f} eV") print( f" Energy Drift (last): {drift[-1]:.6f} eV (relative: {drift[-1] / e0 * 100:.3e}%)" ) fig, ax = plt.subplots(2, 1, figsize=(7.0, 5.0), sharex=True, constrained_layout=True) ax[0].plot(steps, temps, lw=1.5) ax[0].set_ylabel("Temperature (K)") ax[1].plot(steps, drift, lw=1.5) ax[1].axhline(0.0, color="k", ls="--", lw=1.0) ax[1].set_xlabel("Step") ax[1].set_ylabel(r"$E(t) - E(0)$ (eV)") fig.suptitle("Velocity Verlet (NVE): Temperature and Energy Drift") print("\n" + "=" * 95) print("SIMULATION COMPLETE") print("=" * 95) .. image-sg:: /examples/dynamics/images/sphx_glr_02_velocity_verlet_integration_001.png :alt: Velocity Verlet (NVE): Temperature and Energy Drift :srcset: /examples/dynamics/images/sphx_glr_02_velocity_verlet_integration_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none --- Analysis --- Mean Temperature: 48.25 ± 11.24 K Mean Total Energy: -18.4918 eV Energy Fluctuation: 0.0246 eV Energy Drift (last): 0.089460 eV (relative: -4.814e-01%) =============================================================================================== SIMULATION COMPLETE =============================================================================================== .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.670 seconds) .. _sphx_glr_download_examples_dynamics_02_velocity_verlet_integration.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 02_velocity_verlet_integration.ipynb <02_velocity_verlet_integration.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 02_velocity_verlet_integration.py <02_velocity_verlet_integration.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 02_velocity_verlet_integration.zip <02_velocity_verlet_integration.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_