Metadynamics 分析¶

Load Environment¶

First, load necessary modules to be load in the following steps.

In [2]:

import numpy as np

import numpy as np

Load Hills¶

We could load hills from HILLS produced by Plumed, to do more analysis. Here, we just use the examples provided by V. Spiwok, which is trajectories of Alanine Dipeptide in water with 1, 2 or 3 Ramachandran angles, respectively.

In [3]:

from catflow.metad.hills import Hills

from catflow.metad.hills import Hills

In [4]:

#load hills
h1 = Hills(name="../../tests/metad/data/hills/acealanme1d", periodic=[True], cv_per=[[-np.pi, np.pi]])
h2 = Hills(name="../../tests/metad/data/hills/acealanme", periodic=[True,True], cv_per=[[-np.pi, np.pi], [-np.pi, np.pi]])
h3 = Hills(name="../../tests/metad/data/hills/acealanme3d", periodic=[True,True,True], cv_per=[[-np.pi, np.pi], [-np.pi, np.pi], [-np.pi, np.pi]])

#load hills h1 = Hills(name="../../tests/metad/data/hills/acealanme1d", periodic=[True], cv_per=[[-np.pi, np.pi]]) h2 = Hills(name="../../tests/metad/data/hills/acealanme", periodic=[True,True], cv_per=[[-np.pi, np.pi], [-np.pi, np.pi]]) h3 = Hills(name="../../tests/metad/data/hills/acealanme3d", periodic=[True,True,True], cv_per=[[-np.pi, np.pi], [-np.pi, np.pi], [-np.pi, np.pi]])

Sum Hills using Fes¶

We could just use metadynminer.fes to sum the hills to get the Free Energy Surface (FES).

For example, here we just use the fast approach to draw the FES of acealanme (with 2 CVs).

In [5]:

from catflow.metad.fes import FreeEnergySurface

# do sum_hills
fes = FreeEnergySurface.from_hills(h2, resolution=256)

from catflow.metad.fes import FreeEnergySurface # do sum_hills fes = FreeEnergySurface.from_hills(h2, resolution=256)

For here, there are two approaches. One is the fast approach, quickly do the sum in lack of accuracy. The other is the original approach, just do the same as plumed sum_hills command. Here, we use the original method with 4 workers to accelerate the build of FES.

In [6]:

fes.make_fes_original(resolution=256, n_workers=4)
# plot the FES
fes.plot(cmap="RdYlBu_r", levels=20, dpi=96, style='ticks', context='talk')

fes.make_fes_original(resolution=256, n_workers=4) # plot the FES fes.plot(cmap="RdYlBu_r", levels=20, dpi=96, style='ticks', context='talk')

Out[6]:

(<Figure size 960x672 with 2 Axes>,
 <Axes: xlabel='CV1 - phi', ylabel='CV2 - psi'>)

Remove CVs¶

You could use fes.remove_cv(cv_index) or fes.remove_cvs([cv_index_0, cv_index_1]) to remove CV(s) from the built FES, to plot the relationship between one or two CVs and free energy surface.

In [6]:

# do sum_hills
fes2 = fes.remove_cvs([1])

# do sum_hills fes2 = fes.remove_cvs([1])

2023-07-13 21:58:53,381 - INFO : Removing CV 1.

In [7]:

# plot the FES
fes2.plot(cmap="RdYlBu", levels=20, dpi=96, style='ticks', context='talk')

# plot the FES fes2.plot(cmap="RdYlBu", levels=20, dpi=96, style='ticks', context='talk')

Out[7]:

(<Figure size 960x672 with 1 Axes>,
 <Axes: xlabel='CV1 - phi', ylabel='Free Energy (kJ/mol)'>)

Compared with FES from h1, we could see similar trend.

In [8]:

# do sum_hills
fes_1d = FreeEnergySurface.from_hills(h1, resolution=256)

# do sum_hills fes_1d = FreeEnergySurface.from_hills(h1, resolution=256)

In [9]:

fes_1d.make_fes_original(resolution=256, n_workers=4)
# plot the FES
fes_1d.plot(cmap="RdYlBu", levels=20, dpi=96, style='ticks', context='talk')

fes_1d.make_fes_original(resolution=256, n_workers=4) # plot the FES fes_1d.plot(cmap="RdYlBu", levels=20, dpi=96, style='ticks', context='talk')

Out[9]:

(<Figure size 960x672 with 1 Axes>,
 <Axes: xlabel='CV1 - phi', ylabel='Free Energy (kJ/mol)'>)

Find local minima¶

We could then use fes.find_minima to analyze the FES acquired, and to label them in the plot.

In [6]:

fes.find_minima()

# plot the minimas
fes.plot_minima(mark_color="white", png_name=None, style='ticks', context='notebook')

fes.find_minima() # plot the minimas fes.plot_minima(mark_color="white", png_name=None, style='ticks', context='notebook')

Out[6]:

(<Figure size 960x672 with 2 Axes>,
 <Axes: xlabel='CV1 - phi', ylabel='CV2 - psi'>)

Local minima are stored in pandas.DataFrame, as shown in below.

In [11]:

fes.minima

fes.minima

Out[11]:

	Minimum	free energy	CV1bin	CV2bin	CV1 - phi	CV2 - psi
0	0	0.000000	78.0	236.0	-1.214913	2.662991
1	1	1.635963	27.0	240.0	-2.466641	2.761165
2	2	2.670061	74.0	117.0	-1.313088	-0.257709
3	3	5.255746	166.0	150.0	0.944932	0.552233
4	4	12.537359	170.0	251.0	1.043107	3.031146

Free Energy Time-dependent Profile¶

We could draw the time-dependent profile of free energies of local minima from FEProfile.

In [12]:

from catflow.metad.profile import FreeEnergyProfile
fe_profile = FreeEnergyProfile(fes, h2)
fe_profile.plot(energy_unit="kJ/mol", cmap="viridis", style='ticks', context='notebook')

from catflow.metad.profile import FreeEnergyProfile fe_profile = FreeEnergyProfile(fes, h2) fe_profile.plot(energy_unit="kJ/mol", cmap="viridis", style='ticks', context='notebook')

Using String Method to Find Minimal Energy Path (MEP)¶

We could use string method to find MEP on the plotted FEP.

Reference: Weinan E, et al. "Simplified and improved string method for computing the minimum energy paths in barrier-crossing events", J. Chem. Phys. 126, 164103 (2007), https://doi.org/10.1063/1.2720838

In [46]:

#from catflow.metad.string import StringMethod
s = StringMethod(fes)

#from catflow.metad.string import StringMethod s = StringMethod(fes)

Load the minima from the upper block.

In [47]:

s.load_minima()

s.load_minima()

2023-07-14 12:17:02,456 - INFO :    Minimum  free energy  CV1bin  CV2bin  CV1 - phi  CV2 - psi
0        0     0.000000    78.0   236.0  -1.214913   2.662991
1        1     1.635963    27.0   240.0  -2.466641   2.761165
2        2     2.670061    74.0   117.0  -1.313088  -0.257709
3        3     5.255746   166.0   150.0   0.944932   0.552233
4        4    12.537359   170.0   251.0   1.043107   3.031146

We just try the path from Minima 1 to 2, crossing 0.

In [48]:

s.mep_from_minima(begin_index=1, end_index=2, mid_indices=[0], maxsteps=200000)

s.mep_from_minima(begin_index=1, end_index=2, mid_indices=[0], maxsteps=200000)

2023-07-14 12:17:15,186 - INFO : Change in string: 0.0069036067
2023-07-14 12:17:27,184 - INFO : Change in string: 0.0122908930
2023-07-14 12:17:38,706 - INFO : Change in string: 0.0361591569
2023-07-14 12:17:50,272 - INFO : Change in string: 0.0056718707
2023-07-14 12:18:02,275 - INFO : Change in string: 0.0002050943
2023-07-14 12:18:14,258 - INFO : Change in string: 0.0000045245
2023-07-14 12:18:27,275 - INFO : Change in string: 0.0000006603
2023-07-14 12:18:39,405 - INFO : Change in string: 0.0000002474
2023-07-14 12:18:51,680 - INFO : Change in string: 0.0000000967
2023-07-14 12:19:03,603 - INFO : Change in string: 0.0000000379
2023-07-14 12:19:15,659 - INFO : Change in string: 0.0000000148
2023-07-14 12:19:21,624 - INFO : Change in string lower than tolerance.
2023-07-14 12:19:21,625 - INFO : Converged in 116 steps.

It took 116 iterations for the string to be converged. We could also use the fourth-order Runge-Kutta method to solve the evolution of string.

In [51]:

s.mep_from_minima(begin_index=1, end_index=2, mid_indices=[0], maxsteps=200000, integrator="rk4")

s.mep_from_minima(begin_index=1, end_index=2, mid_indices=[0], maxsteps=200000, integrator="rk4")

2023-07-14 12:20:56,777 - INFO : Change in string: 0.0072241848
2023-07-14 12:21:45,994 - INFO : Change in string: 0.0103739428
2023-07-14 12:22:35,167 - INFO : Change in string: 0.0353083649
2023-07-14 12:23:22,341 - INFO : Change in string: 0.0084581203
2023-07-14 12:24:09,113 - INFO : Change in string: 0.0001819289
2023-07-14 12:24:55,353 - INFO : Change in string: 0.0000021306
2023-07-14 12:25:41,926 - INFO : Change in string: 0.0000004627
2023-07-14 12:26:28,478 - INFO : Change in string: 0.0000001716
2023-07-14 12:27:15,062 - INFO : Change in string: 0.0000000642
2023-07-14 12:28:01,155 - INFO : Change in string: 0.0000000240
2023-07-14 12:28:43,129 - INFO : Change in string lower than tolerance.
2023-07-14 12:28:43,130 - INFO : Converged in 110 steps.

After 110 iteration, the string converged. It took more time than forward Euler method, while the former is more stable. We could plot its evolution:

In [49]:

s.plot_string_evolution(cmap="RdYlBu_r", levels=20, dpi=96)

s.plot_string_evolution(cmap="RdYlBu_r", levels=20, dpi=96)

Out[49]:

(<Figure size 960x672 with 2 Axes>,
 <Axes: xlabel='CV1 - phi', ylabel='CV2 - psi'>)

And draw the final MEP on the plot:

In [50]:

s.plot_mep(cmap="RdYlBu_r", levels=20, dpi=96)

s.plot_mep(cmap="RdYlBu_r", levels=20, dpi=96)

Out[50]:

(<Figure size 960x672 with 2 Axes>,
 <Axes: xlabel='CV1 - phi', ylabel='CV2 - psi'>)

Last but not least, the free energy profile of the MEP:

In [11]:

s.plot_mep_energy_profile(dpi=96)

s.plot_mep_energy_profile(dpi=96)

Out[11]:

(<Figure size 614.4x460.8 with 1 Axes>,
 <Axes: xlabel='Reaction coordinate', ylabel='Free Energy (kJ/mol)'>)