""".. _introduction:

ASE Introduction: Nitrogen on copper
====================================

This section gives a quick (and incomplete) overview of what ASE can do.
You can download the code shown in this tutorial (and others in the same style)
as python scripts or jupyter notebooks at the bottom of this page.

We will calculate the adsorption energy of a nitrogen
molecule on a copper surface.
This is done by calculating the total
energy for the isolated slab and for the isolated molecule. The
adsorbate is then added to the slab and relaxed, and the total energy
for this composite system is calculated. The adsorption energy is
obtained as the sum of the isolated energies minus the energy of the
composite system.

You will be able to see an image of the system after relaxation,
later in the "Visualization" section.

Please have a look at the following script:
"""

from ase import Atoms
from ase.build import add_adsorbate, fcc111
from ase.calculators.emt import EMT
from ase.constraints import FixAtoms
from ase.optimize import QuasiNewton

h = 1.85
d = 1.10

slab = fcc111('Cu', size=(4, 4, 2), vacuum=10.0)

slab.calc = EMT()
e_slab = slab.get_potential_energy()

molecule = Atoms('2N', positions=[(0.0, 0.0, 0.0), (0.0, 0.0, d)])
molecule.calc = EMT()
e_N2 = molecule.get_potential_energy()

add_adsorbate(slab, molecule, h, 'ontop')
constraint = FixAtoms(mask=[a.symbol != 'N' for a in slab])
slab.set_constraint(constraint)
dyn = QuasiNewton(slab, trajectory='N2Cu.traj')
dyn.run(fmax=0.05)

print('Adsorption energy:', e_slab + e_N2 - slab.get_potential_energy())

# %%
# Assuming you have ASE setup correctly (:ref:`download_and_install`)
# you can copy it into a python file (e.g. N2Cu.py) and run the script::
#
#  python N2Cu.py
#
# Please read below what the script does.
#
# Atoms
# -----
#
# The :class:`~ase.Atoms` object is a collection of atoms.  Here
# is how to define a N2 molecule by directly specifying the position of
# two nitrogen atoms
#

d = 1.10
molecule = Atoms('2N', positions=[(0.0, 0.0, 0.0), (0.0, 0.0, d)])

# %%
# You can also build crystals using, for example, the lattice module
# which returns :class:`~ase.Atoms` objects corresponding to
# common crystal structures. Let us make a Cu (111) surface
#

from ase.build import fcc111

slab = fcc111('Cu', size=(4, 4, 2), vacuum=10.0)

# %%
# Adding calculator
# -----------------
#
# In this overview we use the effective medium theory (EMT) calculator,
# as it is very fast and hence useful for getting started.
#
# We can attach a calculator to the previously created
# :class:`~ase.Atoms` objects
#

from ase.calculators.emt import EMT

slab.calc = EMT()
molecule.calc = EMT()

# %%
# and use it to calculate the total energies for the systems by using
# the :meth:`~ase.Atoms.get_potential_energy` method from the
# :class:`~ase.Atoms` class
#

e_slab = slab.get_potential_energy()
e_N2 = molecule.get_potential_energy()

# %%
# Structure relaxation
# --------------------
#
# Let's use the :class:`~ase.optimize.QuasiNewton` minimizer to optimize the
# structure of the N2 molecule adsorbed on the Cu surface. First add the
# adsorbate to the Cu slab, for example in the on-top position
#

h = 1.85
add_adsorbate(slab, molecule, h, 'ontop')

# %%
# In order to speed up the relaxation, let us keep the Cu atoms fixed in
# the slab by using :class:`~ase.constraints.FixAtoms` from the
# :mod:`~ase.constraints` module. Only the N2 molecule is then allowed
# to relax to the equilibrium structure
#

from ase.constraints import FixAtoms

constraint = FixAtoms(mask=[a.symbol != 'N' for a in slab])
slab.set_constraint(constraint)

# %%
# Now attach the :class:`~ase.optimize.QuasiNewton` minimizer to the
# system and save the trajectory file. Run the minimizer with the
# convergence criteria that the force on all atoms should be less than
# some ``fmax``
#

from ase.optimize import QuasiNewton

dyn = QuasiNewton(slab, trajectory='N2Cu.traj')
dyn.run(fmax=0.05)

# %%
# .. note::
#
#   The general documentation on
#   :ref:`structure optimizations <structure_optimizations>` contains
#   information about different algorithms, saving the state of an optimizer
#   and other functionality which should be considered when performing
#   expensive relaxations.
#
# Input-output
# ------------
#
# Writing the atomic positions to a file is done with the
# :func:`~ase.io.write` function
#

from ase.io import write

write('slab.xyz', slab)

# %%
# This will write a file in the xyz-format.  Possible formats are:
#
# ========  ===========================
# format    description
# ========  ===========================
# ``xyz``   Simple xyz-format
# ``cube``  Gaussian cube file
# ``pdb``   Protein data bank file
# ``traj``  ASE's own trajectory format
# ``py``    Python script
# ========  ===========================
#
# Reading from a file is done like this
#

from ase.io import read

slab_from_file = read('slab.xyz')

# %%
# If the file contains several configurations, the default behavior of
# the :func:`~ase.io.write` function is to return the last
# configuration. However, we can load a specific configuration by
# doing
#

read('N2Cu.traj')  # last configuration
read('N2Cu.traj', -1)  # same as above
read('N2Cu.traj', 0)  # first configuration

# %%
# Visualization
# -------------
#
# The simplest way to visualize the atoms is the :func:`~ase.visualize.view`
# function
#

from ase.visualize import view

view(slab)

# %%
#
# This will pop up a :mod:`ase.gui` window.  Alternative viewers can be used
# by specifying the optional keyword ``viewer=...`` - use one of
# 'ase.gui', 'gopenmol', 'vmd', or 'rasmol'. (Note that these alternative
# viewers are not a part of ASE and will need to be installed by the user
# separately.) The VMD viewer can take an optional ``data`` argument to
# show 3D data
#
#
# .. code-block:: python
#
#    view(slab, viewer='VMD', data=array)
#
#
#
# If you do not want a gui to open and plot this directly, you can do
# this with plot_atoms in Matplotlib


import matplotlib.pyplot as plt

from ase.visualize.plot import plot_atoms

fig, ax = plt.subplots()
plot_atoms(slab_from_file, ax)
ax.set_axis_off()

# %%
#
# Molecular dynamics
# ------------------
#
# Let us look at the nitrogen molecule as an example of molecular
# dynamics with the :class:`VelocityVerlet <ase.md.verlet.VelocityVerlet>`
# algorithm. We first create the :class:`VelocityVerlet
# <ase.md.verlet.VelocityVerlet>` object giving it the molecule and the time
# step for the integration of Newton's law. We then perform the dynamics
# by calling its :meth:`~ase.md.verlet.VelocityVerlet.run` method and
# giving it the number of steps to take:
#

from ase import units
from ase.md.verlet import VelocityVerlet

dyn = VelocityVerlet(molecule, timestep=1.0 * units.fs)
for i in range(10):
    pot = molecule.get_potential_energy()
    kin = molecule.get_kinetic_energy()
    print('%2d: %.5f eV, %.5f eV, %.5f eV' % (i, pot + kin, pot, kin))
    dyn.run(steps=20)