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

.. only:: html

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

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

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

.. _sphx_glr_examples_generated_tutorials_equil_tipnp.py:

.. _tipnp water box equilibration:

Equilibrating a TIPnP Water Box
===============================

This tutorial shows how to use the TIP3P and TIP4P force fields in
ASE.

.. GENERATED FROM PYTHON SOURCE LINES 9-15

.. note::

  Due to GitLab limit we have to cut the simulation steps in the tutorial.
  Please adjust the number of steps to your own system.
  We recommend multiply the number of steps we provided by 1000 for
  a realistic use.

.. GENERATED FROM PYTHON SOURCE LINES 17-24

.. code-block:: Python

    import ase.units as units
    from ase.build import molecule
    from ase.calculators.tip3p import TIP3P
    from ase.constraints import FixBondLengths
    from ase.io.trajectory import Trajectory
    from ase.md import Langevin








.. GENERATED FROM PYTHON SOURCE LINES 25-27

Let's first create a starting point of the simulaiton.
We will create a water box at 20 °C density.

.. GENERATED FROM PYTHON SOURCE LINES 27-38

.. code-block:: Python

    atoms = molecule('H2O')
    density = 0.9982  # denisty of water in g/cm^3 at 20°C
    box_length = ((atoms.get_masses().sum() / units.mol) / (density * 1e-24)) ** (
        1 / 3
    )  # box length in Å
    atoms.set_cell((box_length, box_length, box_length))
    atoms.center()
    # Repeat the water molecule we just created to end up with a PBC cell
    atoms *= (3, 3, 3)
    atoms.set_pbc(True)








.. GENERATED FROM PYTHON SOURCE LINES 39-40

We can visualise the starting box

.. GENERATED FROM PYTHON SOURCE LINES 40-47

.. code-block:: Python

    import matplotlib.pyplot as plt

    from ase.visualize.plot import plot_atoms

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



.. image-sg:: /examples_generated/tutorials/images/sphx_glr_equil_tipnp_001.png
   :alt: equil tipnp
   :srcset: /examples_generated/tutorials/images/sphx_glr_equil_tipnp_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 48-54

Since the TIPnP type water interpotentials are for rigid
molecules, there are no intramolecular force terms, and we need to
constrain all internal degrees of freedom. For this, we're
using the RATTLE-type constraints of the :ref:`FixBondLengths` class to
constrain all internal atomic distances (O-H1, O-H2, and H1-H2) for
each molecule.

.. GENERATED FROM PYTHON SOURCE LINES 54-61

.. code-block:: Python

    atoms.constraints = FixBondLengths(
        [(3 * i + j, 3 * i + (j + 1) % 3) for i in range(3**3) for j in [0, 1, 2]]
    )
    # RATTLE-type constraints on O-H1, O-H2, H1-H2.
    tag = 'tip3p_27mol_equil'
    atoms.calc = TIP3P(rc=4.5)  # set the calculator to be the TIP3P force field








.. GENERATED FROM PYTHON SOURCE LINES 62-68

For efficiency, we first equillibrate a smaller box, and then repeat that
once more for the final equillibration. However, the potentials are not
parallelized, and are mainly included for testing and for use with QM/MM
tasks, so expect to let it run for some time.
For illustration, we will equillibrate our system with the Langevin
thermostat.

.. GENERATED FROM PYTHON SOURCE LINES 68-105

.. code-block:: Python


    # Equillibrate in a small box first.
    md = Langevin(
        atoms,
        1 * units.fs,
        temperature_K=293.15,  # 20 °C
        friction=0.01,
        logfile=tag + '.log',
    )

    traj = Trajectory(tag + '.traj', 'w', atoms)
    md.attach(traj.write, interval=1)
    md.run(4)  # please use 4000 to better equilibrate

    # Repeat box and equilibrate further.
    tag = 'tip3p_216mol_equil'
    atoms.set_constraint()  # repeat not compatible with FixBondLengths currently.
    atoms *= (2, 2, 2)
    atoms.constraints = FixBondLengths(
        [
            (3 * i + j, 3 * i + (j + 1) % 3)
            for i in range(len(atoms) // 3)
            for j in [0, 1, 2]
        ]
    )
    atoms.calc = TIP3P(rc=7.0)
    md = Langevin(
        atoms,
        2 * units.fs,
        temperature_K=293.15,  # 20 °C
        friction=0.01,
        logfile=tag + '.log',
    )

    traj = Trajectory(tag + '.traj', 'w', atoms)
    md.attach(traj.write, interval=1)
    md.run(2)  # please use 2000 to better equilibrate




.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    /build/python-ase-iXAd1o/python-ase-3.28.0/ase/md/langevin.py:110: FutureWarning: The implementation of `fixcm=True` in `Langevin` does not strictly sample the correct NVT distributions. The deviations are typically small for large systems but can be more pronounced for small systems. Use `fixcm=False` together with `ase.constraints.FixCom`. `fixcm` is deprecated since ASE 3.28.0 and will be removed in a future release.
      warnings.warn(msg, FutureWarning)

    True



.. GENERATED FROM PYTHON SOURCE LINES 106-112

.. note::

 The temperature calculated by ASE is assuming all degrees of freedom
 are available to the system. Since the constraints have removed the 3
 vibrational modes from each water, the shown temperature will be 2/3
 of the actual value.

.. GENERATED FROM PYTHON SOURCE LINES 114-121

The procedure for the TIP4P force field is the same, with the following
exception: the atomic sequence **must** be OHH, OHH, ... .
So to perform the same task using TIP4P, you simply have to import
that calculator instead:
``from ase.calculators.tip4p import TIP4P, rOH, angleHOH``

More info about the TIP4P potential: :mod:`ase.calculators.tip4p`


.. _sphx_glr_download_examples_generated_tutorials_equil_tipnp.py:

.. only:: html

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

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

      :download:`Download Jupyter notebook: equil_tipnp.ipynb <equil_tipnp.ipynb>`

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

      :download:`Download Python source code: equil_tipnp.py <equil_tipnp.py>`

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

      :download:`Download zipped: equil_tipnp.zip <equil_tipnp.zip>`


.. only:: html

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

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