{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Atoms and calculators\n\nASE allows atomistic calculations to be scripted with different\ncomputational codes. In this introductory exercise, we go through the\nbasic concepts and workflow of ASE and will eventually\ncalculate the binding curve of :mol:`N_2`.\n\n\n## Python\n\nIn ASE, calculations are performed by writing and running Python\nscripts.  A very short primer on Python can be found in the\n`ASE documentation <what is python>`.\nIf you are new to Python it would be wise to look through\nthis to understand the basic syntax, datatypes, and\nthings like imports.  Or you can just wing it --- we won't judge.\n\n\n## Atoms\n\nLet's set up a molecule and run a calculation.\nWe can create\nsimple molecules by manually typing the chemical symbols and a\nguess for the atomic positions in \u00c5ngstr\u00f6m.  For example\n:mol:`N_2`:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase import Atoms\n\natoms = Atoms('N2', positions=[[0, 0, -1], [0, 0, 1]])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Just in case we made a mistake, we should visualize our molecule\nusing the :mod:`ASE GUI <ase.gui>`:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.visualize import view"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "```python\nview(atoms)\n```\nEquivalently we can save the atoms in some format, often ASE's own\n:mod:`~ase.io.trajectory` format:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.io import write\n\nwrite('myatoms.traj', atoms)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then run the GUI from a terminal::\n\n $ ase gui myatoms.traj\n\nASE supports quite a few different formats.   For the full list, run::\n\n $ ase info --formats\n\nAlthough we won't be using all the ASE commands any time soon,\nfeel free to get an overview::\n\n $ ase --help\n\n\n## Calculators\n\nNext, let us perform a calculation.  ASE uses\n:mod:`~ase.calculators` to perform calculations. Calculators are\nabstract interfaces to different backends which do the actual computation.\nNormally, calculators work by calling an external electronic structure\ncode or force field code.  To run a calculation, we must first create a\ncalculator and then attach it to the :class:`~ase.Atoms` object.\nFor demonstration purposes, we use the :class:`~ase.calculators.emt.EMT`\ncalculator which is implemented in ASE.\nHowever, there are many other internal and external\ncalculators to choose from (see :mod:`~ase.calculators`).\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.calculators.emt import EMT\n\ncalc = EMT()\natoms.calc = calc"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Once the :class:`~ase.Atoms` object have a calculator with appropriate\nparameters, we can do things like calculating energies and forces:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "e = atoms.get_potential_energy()\nprint('Energy', e)\nf = atoms.get_forces()\nprint('Forces', f)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This will give us the energy in eV and the forces in eV/\u00c5\n(see :mod:`~ase.units` for the standard units ASE uses).\n\nDepending on the calculator, other properties are also available to calculate.\nFor this check the documentation of the respective calculator\nor print the implemented properties the following way:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(EMT.implemented_properties)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Binding curve of :mol:`N_2`\n\nThe strong point of ASE is that things are scriptable.\n``atoms.positions`` is a numpy array containing the atomic positions:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(atoms.positions)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can move the nitrogen atoms by adding or assigning other values into some\nof the array elements.  ASE understands that the state of the atoms object has\nchanged and therefore we can trigger a new calculation by calling\n:meth:`~ase.Atoms.get_potential_energy` or :meth:`~ase.Atoms.get_forces`\nagain, without reattatching a calculator.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "atoms.positions[0, 2] += 0.1  # z-coordinate change of atom 0\ne = atoms.get_potential_energy()\nprint('Energy', e)\nf = atoms.get_forces()\nprint('Forces', f)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This way we can implement any series of calculations by changing the atoms\nobject and subsequently calculating a property.  When running\nmultiple calculations, we often want to write them into a file.\nWe can use the standard :mod:`~ase.io.trajectory` format to write multiple\ncalculations in which the atoms objects and their respective properties\nsuch as energy and forces are contained. Here for a single\ncalculation:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.io.trajectory import Trajectory\n\nwith Trajectory('mytrajectory.traj', 'w') as traj:\n    traj.write(atoms)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we can displace one of the atoms in small steps to trace out a binding\nenergy curve $E(d)$ around the equilibrium\ndistance.  We safe each step to a single trajectory file so that we can\nevaluate the results later on separately.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "atoms = Atoms('N2', positions=[[0, 0, -1], [0, 0, 1]])\n\ncalc = EMT()\natoms.calc = calc\n\nstep = 0.1\nnsteps = int(6 / step)\n\nwith Trajectory('binding_curve.traj', 'w') as traj:\n    for i in range(nsteps):\n        d = 0.5 + i * step\n        atoms.positions[1, 2] = atoms.positions[0, 2] + d\n\n        e = atoms.get_potential_energy()\n        f = atoms.get_forces()\n        print('distance, energy', d, e)\n        print('force', f)\n        traj.write(atoms)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As before, you can use the command line interface to visualize\nthe dissociation process::\n\n $ ase gui binding_curve.traj\n\nAlthough the GUI will plot the energy curve for us, publication\nquality plots usually require some manual tinkering.\nASE provides two functions to read trajectories or other files:\n\n * :func:`ase.io.read` reads and returns the last image, or possibly a list of images if the ``index`` keyword is also specified.\n\n * :func:`ase.io.iread` reads multiple images, one at a time.\n\nUse :func:`ase.io.iread` to read the images back in, e.g.:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.io import iread\n\nfor atoms in iread('binding_curve.traj'):\n    print(atoms.get_potential_energy())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we can plot the binding curve (energy as a function of distance)\nwith matplotlib and calculate the dissociation energy.\nWe first collect the energies and the distances when looping\nover the trajectory.  The atoms already have the energy.  Hence, calling\n``atoms.get_potential_energy()`` will simply retrieve the energy\nwithout calculating anything.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n\nenergies = []\ndistances = []\n\nfor atoms in iread('binding_curve.traj'):\n    energies.append(atoms.get_potential_energy())\n    distances.append(atoms.positions[1, 2] - atoms.positions[0, 2])\n\nax = plt.gca()\nax.plot(distances, energies)\nax.set_xlabel('Distance [\u00c5]')\nax.set_ylabel('Total energy [eV]')\nplt.show()\n\nprint('Dissociation energy [eV]: ', energies[-1] - min(energies))"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.13.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}