.. 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:`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-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 24-37 .. 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 39-46 .. 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 47-55 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 56-60 .. code-block:: Python train_frames = ase.io.read("ethanol_reduced_100.xyz", ":") atoms = train_frames[0].copy() .. GENERATED FROM PYTHON SOURCE LINES 61-62 Below we show the initial configuration of a single ethanol molecule in vacuum. .. GENERATED FROM PYTHON SOURCE LINES 63-72 .. 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 73-75 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 76-79 .. code-block:: Python ase.md.velocitydistribution.MaxwellBoltzmannDistribution(atoms, temperature_K=300) .. GENERATED FROM PYTHON SOURCE LINES 80-82 We now register our exported model as the energy calculator to obtain energies and forces. .. GENERATED FROM PYTHON SOURCE LINES 83-86 .. code-block:: Python atoms.calc = MetatomicCalculator("model-md.pt", extensions_directory="extensions/") .. GENERATED FROM PYTHON SOURCE LINES 87-89 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 90-94 .. code-block:: Python integrator = ase.md.VelocityVerlet(atoms, timestep=0.5 * ase.units.fs) .. GENERATED FROM PYTHON SOURCE LINES 95-104 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 105-123 .. 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 124-135 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 136-148 .. 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 149-154 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 155-159 .. 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 160-168 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 169-178 .. 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 179-182 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 183-190 .. 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 191-193 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 194-199 .. 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 200-202 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 203-212 .. 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 22.711 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 `_