```python
import numpy as np
import sympy as sy
import matplotlib.pyplot as plt
from simtk import unit
```

# Harmonic Well Potential

The harmonic well potential is described by the following expression:

\begin{equation}
f(x)=\frac{1}{2} k \left(x^2 + y^2 + z^2 \right)
\end{equation}

Where $k$ is the only parameter and represents the stiffness of the harmonic potential -or the stiffness of the harmonic spring described by the hookes' law-. Notice that the potential for potential dimensions $Y$ and $Z$, has the same shape. In this way we have a three dimensional harmonic well. But let's see here only the proyection over a single dimension $X$ since $Y$ and $Z$ are decorrelated, and there by will behave as $X$.


```python
def harmonic_well(x,k):
    return 0.5*k*x**2
```


```python
k=5.0 * unit.kilocalories_per_mole/ unit.nanometers**2 # stiffness of the harmonic potential

x_serie = np.arange(-5., 5., 0.05) * unit.nanometers

plt.plot(x_serie, harmonic_well(x_serie, k), 'r-')
plt.ylim(-1,5)
plt.xlim(-2,2)
plt.grid()
plt.xlabel("X ({})".format(unit.nanometers))
plt.ylabel("Energy ({})".format(unit.kilocalories_per_mole))
plt.title("Harmonic Well")
plt.show()
```

Different values of $k$ can be tested to graphically see how this parameter accounts for the openness of the well's arms. Or, as it is described below, the period of oscillations of a particle of mass $m$ in a newtonian dynamics.

The hooks' law describes the force suffered by a mass attached to an ideal spring as:

\begin{equation}
F(x) = -k(x-x_{0})
\end{equation}

Where $k$ is the stiffness of the spring and $x_{0}$ is the equilibrium position. Now, since the force is minus the gradient of the potential energy $V(x)$,

\begin{equation}
F(x) = -\frac{d V(x)}{dx},
\end{equation}

we can proof that the spring force is the result of the first harmonic potential derivative:

\begin{equation}
V(x) = \frac{1}{2} k (x-x_{0})^{2}
\end{equation}

And the angular frequency of oscillations of a spring, or a particle goberned by the former potential, is:

\begin{equation}
\omega = \sqrt{\frac{k}{m}}
\end{equation}

Where $m$ is the mass of the particle. This way the potential can also be written as:

\begin{equation}
V(x) = \frac{1}{2} k (x-x_{0})^{2} = \frac{1}{2} m \omega^{2} (x-x_{0})^{2}
\end{equation}

Finnally, the time period of these oscillations are immediately computed from the mass of the particle, $m$, and the stiffness parameter $k$. Given that by definition:

\begin{equation}
T = 2\pi / \omega
\end{equation}

Then:

\begin{equation}
T = 2\pi \sqrt{\frac{m}{k}}
\end{equation}

## Working with this test system

This test system is fully documented in [HarmonicWell class API](../api/_autosummary/uibcdf_test_systems.HarmonicWell.html). Let's see an example of how to interact with it:


```python
from openmolecularsystems import HarmonicWell

open_molecular_system = HarmonicWell(n_particles = 1, mass = 32 * unit.amu,
                                     k=5.0 * unit.kilocalories_per_mole/unit.nanometers**2)
```

The potential expression and the value of the parameters are stored in `potential`:


```python
open_molecular_system.potential_expression
```




$\displaystyle 0.5 k \left(x^{2} + y^{2} + z^{2}\right)$




```python
open_molecular_system.parameters
```




    {'n_particles': 1,
     'mass': Quantity(value=32, unit=dalton),
     'k': Quantity(value=5.0, unit=kilocalorie/(nanometer**2*mole))}




```python
open_molecular_system.coordinates
```




    Quantity(value=array([[0., 0., 0.]], dtype=float32), unit=nanometer)




```python
open_molecular_system.topology
```


```python
open_molecular_system.system
```




    <simtk.openmm.openmm.System; proxy of <Swig Object of type 'OpenMM::System *' at 0x7f513ab0f300> >



There is a method to evaluate the potential at a given positions:


```python
open_molecular_system.evaluate_potential([-1.5, 0.0, 0.0] * unit.nanometers)
```




    Quantity(value=5.625, unit=kilocalorie/mole)




```python
position = np.zeros((200,3), dtype=float) * unit.nanometers
position[:,0] = np.linspace(-5., 5., 200) * unit.nanometers

plt.plot(position[:,0], open_molecular_system.evaluate_potential(position) , 'r-')
plt.ylim(-1,5)
plt.xlim(-2,2)
plt.grid()
plt.xlabel("X ({})".format(unit.nanometers))
plt.ylabel("Energy ({})".format(unit.kilocalories_per_mole))
plt.title("Harmonic Well")
plt.show()
```


```python
open_molecular_system.get_oscillations_time_period()
```




    Quantity(value=7.770948260727904, unit=picosecond)




```python
open_molecular_system.get_standard_deviation(300.0*unit.kelvin)
```




    Quantity(value=0.3453002396733165, unit=nanometer)



### Newtonian dynamics


```python
initial_positions = np.zeros([1, 3], np.float32) * unit.nanometers
initial_positions[0,0] = 1.0 * unit.nanometers

open_molecular_system.set_coordinates(initial_positions)
```


```python
from openmolecularsystems.tools.md import newtonian
traj_dict = newtonian(open_molecular_system, time=50.0*unit.picoseconds, saving_timestep=0.1*unit.picoseconds,
                      integration_timestep=0.02*unit.picoseconds)
```

    100%|██████████| 2500/2500 [00:00<00:00, 8947.70it/s]


We can now plot the trajectory of the x coordinate:


```python
plt.plot(traj_dict['time'], traj_dict['coordinates'][:,0,0])
plt.xlabel('time ({})'.format(traj_dict['time'].unit))
plt.ylabel('X ({})'.format(traj_dict['coordinates'].unit))
plt.show()
```

We can wonder now if the period of these oscillations is in agreement with the value calculated above.


```python
mass = 32 * unit.amu
T = 2*np.pi*np.sqrt(mass/k)

print('The period of the small oscillations around the minimum is',T)
```

    The period of the small oscillations around the minimum is 7.770948260727904 ps



```python
open_molecular_system.get_oscillations_time_period()
```




    Quantity(value=7.770948260727904, unit=picosecond)




```python
plt.plot(traj_dict['time'], traj_dict['coordinates'][:,0,0])
plt.axvline(T._value, color='gray', linestyle='--') # Period of the harmonic oscillations
plt.xlabel('time ({})'.format(traj_dict['time'].unit))
plt.ylabel('X ({})'.format(traj_dict['coordinates'].unit))
plt.show()
```

Remember that the integration timestep must be smaller than $\sim T/10.0$ to guarantee that no artifacts are introduced by the timestep size.

The newtonian dynamics can also include damping. This way we can simulate damped oscillations around the minimum.


```python
traj_dict = newtonian(open_molecular_system, time=50.0*unit.picoseconds, saving_timestep=0.1*unit.picoseconds,
                      integration_timestep=0.02*unit.picoseconds, friction=0.25/unit.picoseconds)
```

    100%|██████████| 2500/2500 [00:00<00:00, 8990.40it/s]



```python
plt.plot(traj_dict['time'], traj_dict['coordinates'][:,0,0])
plt.xlabel('time ({})'.format(traj_dict['time'].unit))
plt.ylabel('X ({})'.format(traj_dict['coordinates'].unit))
plt.show()
```

What would be the friction value needed to enter in the overdamped regime?


```python
traj_dict = newtonian(open_molecular_system, time=50.0*unit.picoseconds, saving_timestep=0.1*unit.picoseconds,
                      integration_timestep=0.02*unit.picoseconds, friction=5.0/unit.picoseconds)
```

    100%|██████████| 2500/2500 [00:00<00:00, 8904.40it/s]



```python
plt.plot(traj_dict['time'], traj_dict['coordinates'][:,0,0])
plt.xlabel('time ({})'.format(traj_dict['time'].unit))
plt.ylabel('X ({})'.format(traj_dict['coordinates'].unit))
plt.show()
```

### Stochastic Dynamics


```python
from openmolecularsystems.tools.md import langevin_NVT

traj_dict = langevin_NVT(open_molecular_system, time=0.5*unit.nanoseconds, saving_timestep=0.5*unit.picoseconds,
                         integration_timestep=0.02*unit.picoseconds,
                         friction=1.0/unit.picoseconds, temperature=300.0*unit.kelvin,
                         initial_velocities='boltzmann')
```

    100%|██████████| 25000/25000 [00:00<00:00, 30010.47it/s]


Let us see the time evolution of the coordinate $x$ of our single particle:


```python
plt.plot(traj_dict['time'], traj_dict['coordinates'][:,0,0])
plt.xlabel('time ({})'.format(traj_dict['time'].unit))
plt.ylabel('X ({})'.format(traj_dict['coordinates'].unit))
plt.show()
```


```python
open_molecular_system.get_standard_deviation(300.0*unit.kelvin)
```




    Quantity(value=0.3453002396733165, unit=nanometer)




```python
np.std(traj_dict['coordinates'][:,0,0])
```




    Quantity(value=0.36140549642467434, unit=nanometer)



fc09801bd70109b203f6196a63f0fb585fef464f
ipynb
Jupyter Notebook
docs/contents/harmonic_well/Harmonic_Well.ipynb
dprada/OpenMolecularSystems
5787fc159f87091ec498cf23abd07c1c2aec6138
MIT
2021-07-02T14:42:08.000Z
2021-07-02T14:42:08.000Z
docs/contents/harmonic_well/Harmonic_Well.ipynb
dprada/OpenMolecularSystems
5787fc159f87091ec498cf23abd07c1c2aec6138
MIT
2021-07-25T02:28:07.000Z
2021-07-25T02:28:07.000Z
docs/contents/harmonic_well/Harmonic_Well.ipynb
dprada/OpenMolecularSystems
5787fc159f87091ec498cf23abd07c1c2aec6138
MIT
2021-06-17T18:56:55.000Z
2021-06-17T18:56:55.000Z
Qwen/Qwen-72B
1. YES
2. YES