```python
%load_ext autoreload
%autoreload 2
```


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

# Double Well Potential

Our double well potential is described by the following expression:

\begin{equation}
V(x,y,z)=E_{0}\left[ \left(\frac{x}{a}\right)^4 -2\left(\frac{x}{a}\right)^2 \right]-\frac{b}{a}x + \frac{1}{2}k\left( y^2 + z^2 \right)
\end{equation}

This potential can be split in two summands, the first one with the potential for the $X$ axis, and a second term with the potential for the $Y$ and $Z$ axes:

\begin{equation}
V(x,y,z)=V_{x}(x)+V_{y,z}(y,z)
\end{equation}

## A double well potential along $X$

The one dimensional potential $V_{x}(x)$ is a function of $x$ as well as three parameters: $E_{0}$, $a$, $b$.

\begin{equation}
V_{x}(x)=E_{0}\left[ \left(\frac{x}{a}\right)^4 -2\left(\frac{x}{a}\right)^2 \right]-\frac{b}{a}x + \frac{1}{2}k\left( y^2 + z^2 \right)
\end{equation}


Let's see the meaning of these constants. For this purpose the following method is defined to evaluate the contribution of the potential energy terms corresponding to the double well in $X$.


```python
def double_well_potential_1D(x,Eo,a,b):
    return Eo*((x/a)**4-2*(x/a)**2)-(b/a)*x
```

## Symmetric case ($b=0$)

A symmetric double well potential can be represented by the former mathematical expression when $b=0$. In this situation, the double well minima are placed in $x=-a$ and $x=a$. And $E_{0}$ is the value of the barrier height (equal for both basins). Lets verify these statements with an example:


```python
Eo=3.0 * unit.kilocalories_per_mole # barrier height when the double well is symmetric.
a=1.0 * unit.nanometers # Absolute value of the coordinates of minima (the potential is an even function).
b=0.0 * unit.kilocalories_per_mole # No need to explanation at this moment

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

plt.plot(x_serie, double_well_potential_1D(x_serie,Eo,a,b), 'r-')
plt.ylim(-4,1)
plt.xlim(-2,2)
plt.grid()
plt.xlabel("X ({})".format(unit.nanometers))
plt.ylabel("Energy ({})".format(unit.kilocalories_per_mole))
plt.title("Symmetric Double Well")

plt.hlines(-3, -1.7, -1, color='gray', linestyle='--')
plt.hlines(0, -1.7, 0, color='gray', linestyle='--')
plt.vlines(-1.6, -3, 0, color='gray', linestyle='--')
plt.vlines(1, -4, -3, color='gray', linestyle='--')
plt.text(-1.85, -1.5, 'Eo', fontsize=12)
plt.text(1.1, -3.6, 'a', fontsize=12)

plt.show()
```

You can play with the last cell, changing the $E_{0}$ and $a$ values, to check the consistency of their description. Or you can instead work with the first derivative of the potential to do the same analytically:


```python
x, Eo, a = sy.symbols('x Eo a')
f = Eo*((x/a)**4-2*(x/a)**2)

g=sy.diff(f,x) # Primera derivada con respecto a x de la funcion f
```


```python
g
```




$\displaystyle Eo \left(- \frac{4 x}{a^{2}} + \frac{4 x^{3}}{a^{4}}\right)$



The first derivative can be factorized to unveil the value of its three roots: $x=0$, the barrier position, and $x=a$ y $x=-a$, the minima positions.


```python
sy.factor(g)
```




$\displaystyle \frac{4 Eo x \left(- a + x\right) \left(a + x\right)}{a^{4}}$



The height of the barrier, from the bottom of the basins to its top, can then be calculated:

\begin{equation}
V_{x}(0)-V_{x}(c)=0-E_{0}\left[ 1-2 \right]=E_{0}
\end{equation}

### Frequency of small oscillations around the minima

Now that the role of $E_{0}$ and $a$ is clear for the symmetric double well, we are interested in the frequency of the small oscillations at the energy minima. This frequency will be our time of reference to choose an appropriate integration time step of at the time of simulate the dynamics of a particle in this potential. This can be done attending to the second derivative of the potential.

The value of any mathematical function close enough to a minimum can be approximated by the value of the function at the minimum plus the resulting contribution of an harmonic potential. The stiffness of this harmonic potential is equal the to value of the second derivative of the function at the minimum. This is what is known as the Taylor expansion of any function (truncated in the third grade):

\begin{equation}
f(x) \approx f(x_{0}) + f'(x_{0})(x-x_{0}) + \frac{1}{2} f''(x_{0})(x-x_{0})^{2}
\end{equation}

And by definition of minimum:

\begin{equation}
f'(x_{0})=0
\end{equation}

So:

\begin{equation}
f(x) \approx f(x_{0}) + \frac{1}{2} f''(x_{0})(x-x_{0})^{2}
\end{equation}

At this point lets make a break to talk about the harmonic potential. We all know the hooks law to describe the force  suffered by a mass attached to an ideal spring around the equilibrium position:

\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. The potential energy $V(x)$ is now deduced given that:

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

So, 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 frequency of oscillation of the 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}

Now, going back to our Taylor expansion of any mathematical function $f(x)$. If the shape of the $f(x)$ around a minimum $x_{0}$ can be approximated with an harmonic potential, this means that the characteristic frequency $\omega$ of the oscillations around the near surroundings of the minimum is -by comparison with the potential behind the Hooke's law-:

\begin{equation}
\omega = \sqrt{\frac{f''(x_{0})}{m}}
\end{equation}

This way the frequency of the small oscillations of a particle with mass $m$ around a minimum can then by obtained from the value of the second derivative at the minimum:


```python
x, Eo, a = sy.symbols('x Eo a')
f = Eo*((x/a)**4-2*(x/a)**2)

gg=sy.diff(f,x,x) # Second derivative of f with respect x

gg
```




$\displaystyle - \frac{4 Eo \left(1 - \frac{3 x^{2}}{a^{2}}\right)}{a^{2}}$



Let's look at the minimum at $x=a$:


```python
gg.subs({x:a})
```




$\displaystyle \frac{8 Eo}{a^{2}}$



In this case the frequency of the oscillations of a particle with mass $m$ is:

\begin{equation}
\omega = \sqrt{\frac{8E_{0}}{ma^{2}}}
\end{equation}

And the period is:

\begin{equation}
T = \frac{2\pi}{\omega} = 2\pi \sqrt{\frac{ma^{2}}{8E_{0}}}
\end{equation}

We can see graphically how the Taylor expansion is a good aproximation when close enough to a minimum:


```python
def harmonic_well_potential_1D(x,k,a,Eo):
    return 0.5*k*(x-a)**2-Eo

Eo=3.0 * unit.kilocalories_per_mole # barrier height when the double well is symmetric.
a=1.0 * unit.nanometers # Absolute value of the coordinates of minima (the potential is an even function).
b=0.0 * unit.kilocalories_per_mole # No need to explanation at this moment
k=(8*Eo)/a**2 # harmonic stiffness

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

plt.plot(x_serie, double_well_potential_1D(x_serie,Eo,a,b), 'r-')
plt.plot(x_serie, harmonic_well_potential_1D(x_serie,k,a,Eo), color='k', linestyle='--')
plt.ylim(-4,1)
plt.xlim(-2,2)
plt.grid()
plt.xlabel("X ({})".format(unit.nanometers))
plt.ylabel("Energy ({})".format(unit.kilocalories_per_mole))
plt.title("Symmetric Double Well and Harmonic Potential")
plt.show()

x_serie = np.arange(-0.5, 1.5, 0.005) * unit.nanometers

plt.plot(x_serie, double_well_potential_1D(x_serie,Eo,a,b), 'r-')
plt.plot(x_serie, harmonic_well_potential_1D(x_serie,k,a,Eo), color='k', linestyle='--')
plt.ylim(-3.05,-2.8)
plt.xlim(0.9,1.1)
plt.grid()
plt.xlabel("X ({})".format(unit.nanometers))
plt.ylabel("Energy ({})".format(unit.kilocalories_per_mole))
plt.title("Harmonic approximation nearby a minimum")
plt.show()
```

## Assymetric case ($b\neq0$)

In the case of $b\neq 0$, our double well turns in to an assymetric potential. In this situation $E_{0}$ and $a$ have **approximately** the same interpretation, and $b$ can be **approximately** understood as the amount energy that basins shift up or down, depending on the relative position to $x=0$. In our double well, the left basin raises and the right basin drops when $b>0$. Lets see this in a plot:


```python
Eo=3.0 * unit.kilocalories_per_mole # barrier height when the double well is symmetric.
a=1.0 * unit.nanometers # Absolute value of the coordinates of minima (the potential is an even function).
b=1.0 * unit.kilocalories_per_mole # vertical shift of basins (approx.)

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

plt.plot(x_serie, double_well_potential_1D(x_serie,Eo,a,b), 'r-')
plt.ylim(-4,1)
plt.xlim(-2,2)
plt.grid()
plt.xlabel("X ({})".format(unit.nanometers))
plt.ylabel("Energy ({})".format(unit.kilocalories_per_mole))
plt.title("Asymmetric Double Well")

plt.hlines(-2, -1.45, -1, color='gray', linestyle='--')
plt.hlines(0.05, -1.45, 0, color='gray', linestyle='--')
plt.vlines(-1.4, -2, 0.05, color='gray', linestyle='--')
plt.vlines(-0.95, -4, -2, color='gray', linestyle='--')
plt.text(-1.75, -2.4, r'$\approx Eo-b$', fontsize=12)
plt.text(-0.9, -3.6, r'$\approx a$', fontsize=12)

plt.show()
```

The value of the energy barrier from the left minimum is $\approx E_{0}+b$, while $\approx E_{0}-b$ accounts for the barrier from the right minimum.

### Frequency of small oscillations around the minima.

Given that $b$ is included as a linear factor in the potential, its effect vanishes in the second derivative. As such, the harmonic approximation described for the symmetric double well is exact also for this case. But notice that the positions of the minima, and the position of the barrier, are slightly shifted with respect to those for the symmetric double well.

## An harmonic potential along $Y$ and $Z$

The behaviour of a particle in the double well potential along $X$ is independent of what happens along $Y$ and $Z$. We could keep the three dimensional potential equal to $V_{x}(x)$, but in this case the particle would diffuse freely in the subspace $Y-Z$. To avoid these, for aesthetic purposes only, let's add an harmonic well in those two axes:

\begin{equation}
V_{yz}(y,z)= \frac{1}{2}k\left( y^2 + z^2 \right)
\end{equation}

The value of the elastic constant $k$ should be lower, or at least equal, to the harmonic approximation of $V_{x}(x)$ around the minima. Otherwise we should be aware of $V_{yz}$ limiting the integration time step. It's then suggested that:

\begin{equation}
k\le \frac{8E_{0}}{a^2}
\end{equation}

## Working with this test system

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


```python
from openmolecularsystems import DoubleWell

double_well_symmetric = DoubleWell(n_particles = 1, mass = 64 * unit.amu,
                                   Eo=3.0 * unit.kilocalories_per_mole,
                                   a=5.0 * unit.angstroms, b=0.0 * unit.kilocalories_per_mole,
                                   k=1.0 * unit.kilocalories_per_mole/unit.angstrom**2)
```

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


```python
double_well_symmetric.potential_expression
```




$\displaystyle Eo \left(- \frac{2.0 x^{2}}{a^{2}} + \frac{x^{4}}{a^{4}}\right) + 0.5 k \left(y^{2} + z^{2}\right) - \frac{b x}{a}$




```python
double_well_symmetric.parameters
```




    {'n_particles': 1,
     'mass': Quantity(value=64, unit=dalton),
     'Eo': Quantity(value=3.0, unit=kilocalorie/mole),
     'a': Quantity(value=5.0, unit=angstrom),
     'b': Quantity(value=0.0, unit=kilocalorie/mole),
     'k': Quantity(value=1.0, unit=kilocalorie/(angstrom**2*mole))}



And there is a method to evaluate the potential at a given positions:


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




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




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

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

And the system to be used with OpenMM is accesible as:


```python
double_well_symmetric.system
```




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




```python
double_well_symmetric.get_coordinates_minima()
```




    [Quantity(value=array([-0.5,  0. ,  0. ]), unit=nanometer),
     Quantity(value=array([0.5, 0. , 0. ]), unit=nanometer)]




```python
double_well_symmetric.get_coordinates_maximum()
```




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




```python
double_well_symmetric.get_small_oscillations_time_periods_around_minima()
```




    ([Quantity(value=array([-0.5,  0. ,  0. ]), unit=nanometer),
      Quantity(value=array([0.5, 0. , 0. ]), unit=nanometer)],
     [Quantity(value=array([2.50806277, 2.45738961, 2.45738961]), unit=picosecond),
      Quantity(value=array([2.50806277, 2.45738961, 2.45738961]), unit=picosecond)])




```python
double_well_symmetric.get_harmonic_standard_deviations_around_minima(temperature=300*unit.kelvin)
```




    ([Quantity(value=array([-0.5,  0. ,  0. ]), unit=nanometer),
      Quantity(value=array([0.5, 0. , 0. ]), unit=nanometer)],
     [Quantity(value=array([0.07880364, 0.07721148, 0.07721148]), unit=nanometer),
      Quantity(value=array([0.07880364, 0.07721148, 0.07721148]), unit=nanometer)])



### Newtonian dynamics

The newtonian dynamics of the system can be integrated with the method 'tools.md.newtonian()':


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

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

traj_dict = newtonian(double_well_symmetric, time=20.0*unit.picoseconds,
                      saving_timestep=0.1*unit.picoseconds,
                      integration_timestep=0.02*unit.picoseconds,
                      initial_coordinates=initial_positions)
```

    1000it [00:00, 8839.92it/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()
```

Now we can wonder, is the period of the oscillations is in agreement with the value calculated before? A shorter run with more points could be enough to say, by inspection, that the calculus is correct.


```python
mass = 64 * unit.amu
Eo=3.0 * unit.kilocalories_per_mole
a=5.0 * unit.angstroms
T = 2*np.pi*np.sqrt((mass*a**2)/(8.0*Eo))

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

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


As we can also check with the method:


```python
double_well_symmetric.get_small_oscillations_time_periods_around_minima()
```




    ([Quantity(value=array([-0.5,  0. ,  0. ]), unit=nanometer),
      Quantity(value=array([0.5, 0. , 0. ]), unit=nanometer)],
     [Quantity(value=array([2.50806277, 2.45738961, 2.45738961]), unit=picosecond),
      Quantity(value=array([2.50806277, 2.45738961, 2.45738961]), 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 approximation
plt.xlabel('time ({})'.format(traj_dict['time'].unit))
plt.ylabel('X ({})'.format(traj_dict['coordinates'].unit))
plt.show()
```

It seems a good approximation to what is observed graphically in the integrated trajectory. Keep in mind that this is an approximation valid when the particle moves very close the bottom of the basins. In an real harmonic oscillator this period does not change with amplitude, but in this case does it. If the amplitude is large enough the oscillations are far from being well approximated by the cuadratic term of the Taylor expansion. Try to play with larger and shorter initial distances to the minimum to see this effect:


```python
amplitudes = 0.5*unit.nanometers + np.array([0.02, 0.05, 0.07, 0.1, 0.15]) * unit.nanometers
traj_dicts = []

for amplitude in amplitudes:
    initial_positions[0,0] = amplitude
    double_well_symmetric.set_coordinates(initial_positions)
    traj_dict = newtonian (double_well_symmetric, time = 5 * unit.picoseconds,
                           integration_timestep = 0.02 * unit.picoseconds,
                           saving_timestep = 0.1 * unit.picoseconds,
                           tqdm=False)
    traj_dicts.append(traj_dict)
    
for traj_dict in traj_dicts:
    plt.plot(traj_dict['time'], traj_dict['coordinates'][:,0,0])

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

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


```python
traj_dict = newtonian(double_well_symmetric, time=20.0*unit.picoseconds,
                      saving_timestep=0.1*unit.picoseconds,
                      integration_timestep=0.02*unit.picoseconds,
                      friction = 0.5 / unit.picoseconds,
                      initial_coordinates=initial_positions,
                      tqdm=False)

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(double_well_symmetric, time=20.0*unit.picoseconds,
                      saving_timestep=0.1*unit.picoseconds,
                      integration_timestep=0.02*unit.picoseconds,
                      friction = 5.0 / unit.picoseconds,
                      initial_coordinates=initial_positions,
                      tqdm=False)

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()
```

And the same can be checked for the asymmetric double well potential. In the case of the period of short oscillations around a minimum, keep in mind that there the minimum is not longer $a$, and $T$ was here calculated with the value of the second derivative (equal, no matter the value of $b$) in $a$. But still is a good approximation.

### Stochastic Dynamics

Thanks to the `tools.md.langevin_NVT` method in this library, a simple newtonian dynamics can be run with few efforce:


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

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

traj_dict = langevin_NVT(double_well_symmetric, time=10.0*unit.nanoseconds,
                      saving_timestep=1.0*unit.picoseconds,
                      integration_timestep=0.1*unit.picoseconds,
                      friction=1.0/unit.picoseconds, temperature=300.0*unit.kelvin,
                      initial_coordinates=initial_positions)
```

    100000it [00:06, 16466.46it/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()
```

The trajectory is long enough to observe the system hopping on and off from the initial basin.

14ef9c4d8f09d0eaa3f7d18d7c1163602504a014
ipynb
Jupyter Notebook
docs/contents/double_well/Double_Well.ipynb
dprada/OpenMolecularSystems
5787fc159f87091ec498cf23abd07c1c2aec6138
MIT
2021-07-02T14:42:08.000Z
2021-07-02T14:42:08.000Z
docs/contents/double_well/Double_Well.ipynb
dprada/OpenMolecularSystems
5787fc159f87091ec498cf23abd07c1c2aec6138
MIT
2021-07-25T02:28:07.000Z
2021-07-25T02:28:07.000Z
docs/contents/double_well/Double_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