.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/2body_potential.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_2body_potential.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_2body_potential.py:

.. _2_body_potential:

Calculating the 2 body potential of Argon using ASE and CP2K
============================================================

In this tutorial you will learn how to calculate the 2 body potential between two Argon atoms with DFT. 
We will use ASE (Atomic Simulation Environment) to setup our simulation. 
To do our energy calculation, ASE will call CP2K, which implements DFT.

First, we import our Python modules

.. GENERATED FROM PYTHON SOURCE LINES 14-28

.. code-block:: Python


    # uncomment the following line if running a jupyter notebook
    # %matplotlib inline

    import numpy as np
    import matplotlib.pyplot as plt
    from ase import Atoms
    from ase.calculators.cp2k import CP2K
    import matplotlib
    import scipy.constants as c

    # uncomment the following line if running on a cluster
    # matplotlib.use("agg")


.. GENERATED FROM PYTHON SOURCE LINES 29-30

We also import pint, which will help us to keep track of units.

.. GENERATED FROM PYTHON SOURCE LINES 30-35

.. code-block:: Python


    import pint

    ureg = pint.UnitRegistry()


.. GENERATED FROM PYTHON SOURCE LINES 36-37

Next, we create an Atoms object with two argon atoms placed a distance apart.

.. GENERATED FROM PYTHON SOURCE LINES 37-41

.. code-block:: Python


    distance = 3.3 * ureg.angstrom
    two_argon_atoms = Atoms("Ar2", [[0, 0, 0], [0, 0, distance.magnitude]])


.. GENERATED FROM PYTHON SOURCE LINES 42-44

We set the cubic simulation box size to something much larger than Ar-Ar distance
and apply periodic boundary conditions in all directions.

.. GENERATED FROM PYTHON SOURCE LINES 44-48

.. code-block:: Python


    two_argon_atoms.center(vacuum=3)
    two_argon_atoms.pbc = [1, 1, 1]


.. GENERATED FROM PYTHON SOURCE LINES 49-51

We will use CP2K as a calculator for our Atoms object.
First, we specify calculation settings, you will not have to worry about these

.. GENERATED FROM PYTHON SOURCE LINES 51-73

.. code-block:: Python


    inp = """
    # Parameters for force calculation.
    &FORCE_EVAL
        # Define the DFT parameters

        &DFT
            &SCF
                &OT
                    MINIMIZER DIIS
                    PRECONDITIONER FULL_SINGLE_INVERSE
                &END OT
                &OUTER_SCF
                    MAX_SCF 100
                    EPS_SCF 1.0E-6
                &END OUTER_SCF
            &END SCF

        &END DFT
    &END FORCE_EVAL
    """


.. GENERATED FROM PYTHON SOURCE LINES 74-75

Next, we specify our basis set, pseudo potential and exchange correlation functional

.. GENERATED FROM PYTHON SOURCE LINES 75-85

.. code-block:: Python


    two_argon_atoms.calc = CP2K(inp=inp,
                                basis_set="DZVP-MOLOPT-SR-GTH",
                                pseudo_potential="GTH-PBE-q8",
                                potential_file="GTH_POTENTIALS",
                                xc="PBE",
                                command="cp2k_shell.psmp",
                                )



.. GENERATED FROM PYTHON SOURCE LINES 86-87

We can now run our first DFT calculation with

.. GENERATED FROM PYTHON SOURCE LINES 87-91

.. code-block:: Python


    E = two_argon_atoms.get_potential_energy()
    print(E)


.. GENERATED FROM PYTHON SOURCE LINES 92-94

This will return us a single potential energy for the specified interatomic distance. 
Next, we want to sample a region between 3.3 and 6.0 Angstrom and get the potential energy for each distance

.. GENERATED FROM PYTHON SOURCE LINES 94-108

.. code-block:: Python


    distances = np.linspace(3.3, 6.0, 20)
    energies = np.zeros(distances.shape)

    for i in range(len(distances)):
        two_argon_atoms.set_positions([[0, 0, 0], [0, 0, distances[i]]])
        two_argon_atoms.center(vacuum=3)
        print(distances[i])
        print(two_argon_atoms.get_positions())

        E = two_argon_atoms.get_potential_energy()
        print(E)
        energies[i] = E


.. GENERATED FROM PYTHON SOURCE LINES 109-110

Finally, we plot the energy as a function of the interatomic distance. 

.. GENERATED FROM PYTHON SOURCE LINES 110-119

.. code-block:: Python


    plt.plot(distances, energies, "x", label="DFT")
    plt.xlabel("$R$ in Angstrom")
    plt.ylabel("$E$ in eV")

    # TODO: Plot the LJ potential
    plt.legend()
    plt.savefig("2body_potential.png", dpi=300)
    #plt.show()


.. _sphx_glr_download_auto_examples_2body_potential.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: 2body_potential.ipynb <2body_potential.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: 2body_potential.py <2body_potential.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: 2body_potential.zip <2body_potential.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_