.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated_examples/0-beginner/05-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_generated_examples_0-beginner_05-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:`arch-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-ase.yaml :language: yaml We first train the model same model but and before import the necessary libraries and run the training process and the integration of ASE. .. GENERATED FROM PYTHON SOURCE LINES 23-36 .. code-block:: Python import subprocess 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 metatomic.torch.ase_calculator import MetatomicCalculator .. GENERATED FROM PYTHON SOURCE LINES 38-45 .. code-block:: Python # Here, we run training as a subprocess, in reality you would run this from the command # line as ``mtt train options-ase.yaml --output model-md.pt``. subprocess.run( ["mtt", "train", "options-ase.yaml", "--output", "model-md.pt"], check=True ) .. rst-class:: sphx-glr-script-out .. code-block:: none CompletedProcess(args=['mtt', 'train', 'options-ase.yaml', '--output', 'model-md.pt'], returncode=0) .. GENERATED FROM PYTHON SOURCE LINES 46-54 A detailed step-by-step introduction on how to train a model is provided in the :ref:`label_basic_usage` tutorial. 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 55-59 .. code-block:: Python train_frames = ase.io.read("ethanol_reduced_100.xyz", ":") atoms = train_frames[0].copy() .. GENERATED FROM PYTHON SOURCE LINES 60-61 Below we show the initial configuration of a single ethanol molecule in vacuum. .. GENERATED FROM PYTHON SOURCE LINES 62-71 .. code-block:: Python ase.visualize.plot.plot_atoms(atoms) plt.xlabel("Å") plt.ylabel("Å") plt.show() .. image-sg:: /generated_examples/0-beginner/images/sphx_glr_05-run_ase_001.png :alt: 05 run ase :srcset: /generated_examples/0-beginner/images/sphx_glr_05-run_ase_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 72-74 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 75-78 .. code-block:: Python ase.md.velocitydistribution.MaxwellBoltzmannDistribution(atoms, temperature_K=300) .. GENERATED FROM PYTHON SOURCE LINES 79-81 We now register our exported model as the energy calculator to obtain energies and forces. .. GENERATED FROM PYTHON SOURCE LINES 82-85 .. code-block:: Python atoms.calc = MetatomicCalculator("model-md.pt", extensions_directory="extensions/") .. GENERATED FROM PYTHON SOURCE LINES 86-88 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 89-93 .. code-block:: Python integrator = ase.md.VelocityVerlet(atoms, timestep=0.5 * ase.units.fs) .. GENERATED FROM PYTHON SOURCE LINES 94-103 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 104-122 .. 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 123-134 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 135-147 .. 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:: /generated_examples/0-beginner/images/sphx_glr_05-run_ase_002.png :alt: 05 run ase :srcset: /generated_examples/0-beginner/images/sphx_glr_05-run_ase_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 148-153 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 154-158 .. 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 159-167 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 168-177 .. 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 178-181 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 182-189 .. 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 190-192 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 193-198 .. 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 199-201 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 202-211 .. 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:: /generated_examples/0-beginner/images/sphx_glr_05-run_ase_004.png :alt: 05 run ase :srcset: /generated_examples/0-beginner/images/sphx_glr_05-run_ase_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 31.577 seconds) .. _sphx_glr_download_generated_examples_0-beginner_05-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: 05-run_ase.ipynb <05-run_ase.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 05-run_ase.py <05-run_ase.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 05-run_ase.zip <05-run_ase.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_