.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/ase/run_ase.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_ase_run_ase.py: Running molecular dynamics with ASE =================================== This tutorial demonstrates how to use an already trained and exported model to run an ASE simulation of a single ethanol molecule in vacuum. We use a model that was trained using the :ref:`architecture-soap-bpnn` architecture on 100 ethanol systems containing energies and forces. You can obtain the :download:`dataset file ` used in this example from our website. The dataset is a subset of the `rMD17 dataset `_. The model was trained using the following training options. .. literalinclude:: options.yaml :language: yaml You can train the same model yourself with .. literalinclude:: train.sh :language: bash A detailed step-by-step introduction on how to train a model is provided in the :ref:`label_basic_usage` tutorial. .. GENERATED FROM PYTHON SOURCE LINES 28-30 First, we start by importing the necessary libraries, including the integration of ASE calculators for metatensor atomistic models. .. GENERATED FROM PYTHON SOURCE LINES 31-42 .. code-block:: Python import ase.md import ase.md.velocitydistribution import ase.units import ase.visualize.plot import matplotlib.pyplot as plt import numpy as np from ase.geometry.analysis import Analysis from metatensor.torch.atomistic.ase_calculator import MetatensorCalculator .. GENERATED FROM PYTHON SOURCE LINES 43-48 Setting up the simulation ------------------------- Next, we initialize the simulation by extracting the initial positions from the dataset file which we initially trained the model on. .. GENERATED FROM PYTHON SOURCE LINES 49-53 .. code-block:: Python train_frames = ase.io.read("ethanol_reduced_100.xyz", ":") atoms = train_frames[0].copy() .. GENERATED FROM PYTHON SOURCE LINES 54-55 Below we show the initial configuration of a single ethanol molecule in vacuum. .. GENERATED FROM PYTHON SOURCE LINES 56-65 .. code-block:: Python ase.visualize.plot.plot_atoms(atoms) plt.xlabel("Å") plt.ylabel("Å") plt.show() .. image-sg:: /examples/ase/images/sphx_glr_run_ase_001.png :alt: run ase :srcset: /examples/ase/images/sphx_glr_run_ase_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 66-68 Our initial coordinates do not include velocities. We initialize the velocities according to a Maxwell-Boltzmann Distribution at 300 K. .. GENERATED FROM PYTHON SOURCE LINES 69-72 .. code-block:: Python ase.md.velocitydistribution.MaxwellBoltzmannDistribution(atoms, temperature_K=300) .. GENERATED FROM PYTHON SOURCE LINES 73-75 We now register our exported model as the energy calculator to obtain energies and forces. .. GENERATED FROM PYTHON SOURCE LINES 76-79 .. code-block:: Python atoms.calc = MetatensorCalculator("model.pt", extensions_directory="extensions/") .. rst-class:: sphx-glr-script-out .. code-block:: none the model suggested to use CUDA devices before CPU, but we are unable to find it .. GENERATED FROM PYTHON SOURCE LINES 80-82 Finally, we define the integrator which we use to obtain new positions and velocities based on our energy calculator. We use a common timestep of 0.5 fs. .. GENERATED FROM PYTHON SOURCE LINES 83-87 .. code-block:: Python integrator = ase.md.VelocityVerlet(atoms, timestep=0.5 * ase.units.fs) .. GENERATED FROM PYTHON SOURCE LINES 88-97 Run the simulation ------------------ We now have everything ready to run the MD simulation at constant energy (NVE). To keep the execution time of this tutorial small we run the simulations only for 100 steps. If you want to run a longer simulation you can increase the ``n_steps`` variable. During the simulation loop we collect data about the simulation for later analysis. .. GENERATED FROM PYTHON SOURCE LINES 98-116 .. code-block:: Python n_steps = 100 potential_energy = np.zeros(n_steps) kinetic_energy = np.zeros(n_steps) total_energy = np.zeros(n_steps) trajectory = [] for step in range(n_steps): # run a single simulation step integrator.run(1) trajectory.append(atoms.copy()) potential_energy[step] = atoms.get_potential_energy() kinetic_energy[step] = atoms.get_kinetic_energy() total_energy[step] = atoms.get_total_energy() .. GENERATED FROM PYTHON SOURCE LINES 117-128 Analyse the results ------------------- Energy conservation ################### For a first analysis, we plot the evolution of the mean of the kinetic, potential, and total energy which is an important measure for the stability of a simulation. As shown below we see that both the kinetic, potential, and total energy fluctuate but the total energy is conserved over the length of the simulation. .. GENERATED FROM PYTHON SOURCE LINES 129-141 .. code-block:: Python plt.plot(potential_energy - potential_energy.mean(), label="potential energy") plt.plot(kinetic_energy - kinetic_energy.mean(), label="kinetic energy") plt.plot(total_energy - total_energy.mean(), label="total energy") plt.xlabel("step") plt.ylabel("energy / kcal/mol") plt.legend() plt.show() .. image-sg:: /examples/ase/images/sphx_glr_run_ase_002.png :alt: run ase :srcset: /examples/ase/images/sphx_glr_run_ase_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 142-147 Inspect the systems ################### Even though the total energy is conserved, we also have to verify that the ethanol molecule is stable and the bonds did not break. .. GENERATED FROM PYTHON SOURCE LINES 148-152 .. code-block:: Python animation = ase.visualize.plot.animate(trajectory, interval=100, save_count=None) plt.show() .. container:: sphx-glr-animation .. raw:: html
.. GENERATED FROM PYTHON SOURCE LINES 153-161 Carbon-hydrogen radial distribution function ############################################ As a final analysis we also calculate and plot the carbon-hydrogen radial distribution function (RDF) from the trajectory and compare this to the RDF from the training set. To use the RDF code from ase we first have to define a unit cell for our systems. We choose a cubic one with a side length of 10 Å. .. GENERATED FROM PYTHON SOURCE LINES 162-171 .. code-block:: Python for atoms in train_frames: atoms.cell = 10 * np.ones(3) atoms.pbc = True for atoms in trajectory: atoms.cell = 10 * np.ones(3) atoms.pbc = True .. GENERATED FROM PYTHON SOURCE LINES 172-175 We now can initilize the :py:class:`ase.geometry.analysis.Analysis` objects and compute the the RDF using the :py:meth:`ase.geometry.analysis.Analysis.get_rdf` method. .. GENERATED FROM PYTHON SOURCE LINES 176-183 .. code-block:: Python ana_traj = Analysis(trajectory) ana_train = Analysis(train_frames) rdf_traj = ana_traj.get_rdf(rmax=5, nbins=50, elements=["C", "H"], return_dists=True) rdf_train = ana_train.get_rdf(rmax=5, nbins=50, elements=["C", "H"], return_dists=True) .. GENERATED FROM PYTHON SOURCE LINES 184-186 We extract the bin positions from the returned values and and averege the RDF over the whole trajectory and dataset, respectively. .. GENERATED FROM PYTHON SOURCE LINES 187-192 .. code-block:: Python bins = rdf_traj[0][1] rdf_traj_mean = np.mean([rdf_traj[i][0] for i in range(n_steps)], axis=0) rdf_train_mean = np.mean([rdf_train[i][0] for i in range(n_steps)], axis=0) .. GENERATED FROM PYTHON SOURCE LINES 193-195 Plotting the RDF verifies that the hydrogen bonds are stable, confirming that we performed an energy-conserving and stable simulation. .. GENERATED FROM PYTHON SOURCE LINES 196-205 .. code-block:: Python plt.plot(bins, rdf_traj_mean, label="trajectory") plt.plot(bins, rdf_train_mean, label="training set") plt.legend() plt.xlabel("r / Å") plt.ylabel("radial distribution function") plt.show() .. image-sg:: /examples/ase/images/sphx_glr_run_ase_004.png :alt: run ase :srcset: /examples/ase/images/sphx_glr_run_ase_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 19.019 seconds) .. _sphx_glr_download_examples_ase_run_ase.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: run_ase.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: run_ase.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: run_ase.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_