.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/md/pair_potentials/plot_01_morse_vs_harmonic.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

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

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

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

.. _sphx_glr_api_gallery_md_pair_potentials_plot_01_morse_vs_harmonic.py:


Morse vs. harmonic bond potentials, and diatomic bond vibration
==================================================================

The harmonic bond potential (:class:`~chemistrykit.md.systems.pair_potentials.HarmonicBond`)
is the small-oscillation limit of the anharmonic
:class:`~chemistrykit.md.systems.pair_potentials.Morse` potential -- they
agree near the equilibrium separation but diverge away from it (Morse
correctly flattens out toward a finite dissociation energy; harmonic
does not). :class:`~chemistrykit.md.systems.pair_potentials.DiatomicOscillator`
integrates the classical two-body bond-vibration problem via the shared
velocity-Verlet integrator, for either potential.

.. GENERATED FROM PYTHON SOURCE LINES 16-38

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from chemistrykit.md.systems.pair_potentials import DiatomicOscillator, HarmonicBond, Morse

    De, re = 4.0, 1.0
    harmonic_k = 2.0 * De * 1.8**2  # match harmonic force constant to Morse's curvature at re (a=1.8)
    morse = Morse(De=De, a=1.8, re=re)
    harmonic = HarmonicBond(k=harmonic_k, r0=re)

    r = np.linspace(0.7, 3.0, 400)

    fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
    axes[0].plot(r, morse.energy(r), label="Morse")
    axes[0].plot(r, harmonic.energy(r), label="harmonic")
    axes[0].axhline(-De, color="gray", linestyle=":", linewidth=0.8, label="dissociation limit (Morse)")
    axes[0].set_xlabel("r")
    axes[0].set_ylabel("U(r)")
    axes[0].set_ylim(-De * 1.2, De * 1.5)
    axes[0].set_title("Potential energy curves")
    axes[0].legend()




.. image-sg:: /api/gallery/md/pair_potentials/images/sphx_glr_plot_01_morse_vs_harmonic_001.png
   :alt: Potential energy curves
   :srcset: /api/gallery/md/pair_potentials/images/sphx_glr_plot_01_morse_vs_harmonic_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x11b1530e0>



.. GENERATED FROM PYTHON SOURCE LINES 39-44

Simulating the classical bond vibration for each potential, started
stretched away from equilibrium: the harmonic bond oscillates forever
with a fixed period; the (slightly softer, anharmonic) Morse bond
oscillates at very nearly the same period for this small a displacement,
since both potentials share the same curvature at re by construction.

.. GENERATED FROM PYTHON SOURCE LINES 44-69

.. code-block:: Python


    m1 = m2 = 1.0
    r0_displaced = 1.3
    dt = 0.002
    n_steps = 2000

    for potential, label in [(harmonic, "harmonic"), (morse, "Morse")]:
        oscillator = DiatomicOscillator(potential, m1=m1, m2=m2, r0=r0_displaced)
        lengths = [oscillator.bond_length()]
        for _ in range(n_steps):
            oscillator.step(dt)
            lengths.append(oscillator.bond_length())
        axes[1].plot(np.arange(n_steps + 1) * dt, lengths, label=label)

    period = DiatomicOscillator.harmonic_period(harmonic_k, m1, m2)
    axes[1].axvline(period, color="gray", linestyle=":", linewidth=0.8, label="harmonic period")
    axes[1].set_xlabel("t")
    axes[1].set_ylabel("bond length r(t)")
    axes[1].set_title("Classical bond-vibration trajectories")
    axes[1].legend()
    fig.tight_layout()

    print(f"Harmonic vibration period: {period:.4f}")

    plt.show()




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

 .. code-block:: none

    Harmonic vibration period: 0.8727





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.374 seconds)


.. _sphx_glr_download_api_gallery_md_pair_potentials_plot_01_morse_vs_harmonic.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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