diff --git a/.nojekyll b/.nojekyll new file mode 100644 index 0000000..e69de29 diff --git a/404.html b/404.html new file mode 100644 index 0000000..eee5687 --- /dev/null +++ b/404.html @@ -0,0 +1,1699 @@ + + + +
+ + + + + + + + + + + + + + +首先将需要分析的团簇轨迹文件读入到程序中,用户可以指定 path
参数为文件路径,
+其他的参数请以可选参数形式传入。
from catflow.structure.cluster import Cluster
+
+trajfile = "./dump.lammpstrj"
+c = Cluster(trajfile, topology_format="LAMMPSDUMP", dt=0.0005)
+
提示
+注意默认的格式为"XYZ",如需要导入其他格式,请参考 +MDAnalysis文档。
+或将 MDAnalysis 的一个 Universe
实例导入:
from MDAnalysis import Universe
+u = Universe(trajfile, topology_format="LAMMPSDUMP", dt=0.0005))
+d = Cluster.convert_universe(u)
+
为分析体系的相变性质,Lindemann等人[1]提出了对结构随时间均方键长的演变进行分析来确定相变性质的方法, +一般计算公式如下:
+这里参照 Welford 算法[2]计算方差,简单实现了 Lindemann Index的计算方法,调用如下:
+lpf = c.lindemann_per_frames(u, select_lang="name Pt")
+
即可得到关于整条轨迹Lindemann index的变化趋势,对 lpf
作图,可以帮助判断相变情况以及确定MD是否收敛。
注意这里 lindemann_per_frames
读入的是MDAnalysis中的Universe对象,通常用 select_lang
来指定需要对哪些原子进行分析,
+即给定对应的Atom selection language来选取,
+请参考官方文档的说明。
为对团簇相变行为的温度依赖进行分析,常对得到的曲线进行拟合,可以采用如下函数,以dataframe形式输出拟合的曲线:
+import numpy as np
+
+temperature = np.array([300., 400., 500., 600., 700.])
+lindemann = np.array([0.05, 0.10, 0.20, 0.25, 0.32])
+bounds = ([-np.inf, -np.inf, -np.inf, -np.inf, 400, 15.],
+ [np.inf, np.inf, np.inf, np.inf, 700., 100.])
+df = fitting_lindemann_curve(temperature, lindemann, bounds, function='func2')
+
其中,
+func1
函数形式为:
+ $$ f(x) = b + (a - b) x + \frac{d}{1 + \exp({\frac{x - x_0}{\mathrm{d}x})}} + cx $$func2
函数形式为
+ $$ f(x) = \frac{ax+b}{1 + \exp({\frac{x - x_0}{\mathrm{d}x})}} + \frac{cx+d}{1 + \exp({\frac{x_0 - x}{\mathrm{d}x})}}$$默认采用 func2
进行拟合,bounds
中上下界的参数对应即为a, b, c, d, x0, dx
的取值范围。
[1] F. A. Lindemann, Zeitschrift für, Phys. 1910, 11, 609–612.
+[2] Donald E. Knuth, The art of computer programming, volume 2 (3rd ed.): seminumerical algorithms, Addison-Wesley Longman Publishing Co, 1997, 232.
+ + + + + + +首先,导入环境:
+from catflow.tesla.dpgen import DPTask
+
加载DP-GEN工作目录:
+t = DPTask(
+ path='/path/to/dpgen/',
+ param_file='param.json',
+ machine_file='machine.json',
+ record_file='record.dpgen'
+)
+
便可根据所需分析的部分,对训练情况进行分析。
+导入分析器(DPAnalyzer
),这里我们选择训练,即:
from catflow.tesla.dpgen.training import DPTrainingAnalyzer
+
从任务初始化分析器实例:
+ana = DPTrainingAnalyzer(t)
+
即可利用ana
的内置函数进行作图:
fig = ana.plot_lcurve(
+ iteration=28, test=False, style='ticks', context='talk'
+)
+fig.set_size_inches((12,12))
+
类似地,我们也可以对模型的model deviation分布进行分析:
+from catflow.tesla.dpgen.exploration import DPExplorationAnalyzer
+ana = DPExplorationAnalyzer(t)
+
利用分析器自带的方法进行作图:
+fig = ana.plot_single_iteration(
+ iteration=41,
+ temps=[400, 600, 800, 1000, 1200],
+ xlimit=1000000,
+ f_trust_lo=t.param_data['model_devi_f_trust_lo'],
+ f_trust_hi=t.param_data['model_devi_f_trust_hi'],
+ style='ticks',
+ group_by='temps',
+ label_unit='K',
+ context='talk'
+)
+
其中:
+iteration
对应为所需分析的轮数,默认为最新进行过Exploartion的轮数。
f_trust_lo
和f_trust_hi
即对应的最大力偏差上下限设置。
通过 group_by
指定所需作图的参数,对应到 param.json
中 model_devi_jobs
中该轮数需要迭代的List,例如:
{
+ "template": {
+ "lmp": "lmp/input-meta.lammps",
+ "plm": "lmp/input-meta.plumed"
+ },
+ "sys_idx": [
+ 53
+ ],
+ "traj_freq": 1000,
+ "rev_mat": {
+ "lmp": {
+ "V_NSTEPS": [
+ 1000000
+ ],
+ "V_TEMP": [
+ 400,
+ 600,
+ 800,
+ 1000,
+ 1200
+ ]
+ }
+ },
+ "model_devi_f_trust_lo": 0.23,
+ "model_devi_f_trust_hi": 0.75
+ }
+
若指定group_by
参数为 V_TEMP
,则根据该轮的热浴温度分组作图,若指定V_TEMP=[400, 600, 800, 1000, 1200]
,则可由400、600、800、1000、1200K分别对model deviation作图。
label_unit
即 group_by
参数的单位,例如这里是温度,故为"K"。效果如下:
+ + + + + + + +本软件包含了一些可能有用的分析模块,具体使用说明如下:
+First, load necessary modules to be load in the following steps.
+import numpy as np
+
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.
from catflow.metad.hills import Hills
+
#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]])
+
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).
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.
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')
+
(<Figure size 960x672 with 2 Axes>, + <Axes: xlabel='CV1 - phi', ylabel='CV2 - psi'>)+
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.
# do sum_hills
+fes2 = fes.remove_cvs([1])
+
2023-07-13 21:58:53,381 - INFO : Removing CV 1. ++
# plot the FES
+fes2.plot(cmap="RdYlBu", levels=20, dpi=96, style='ticks', context='talk')
+
(<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.
# do sum_hills
+fes_1d = FreeEnergySurface.from_hills(h1, resolution=256)
+
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')
+
(<Figure size 960x672 with 1 Axes>, + <Axes: xlabel='CV1 - phi', ylabel='Free Energy (kJ/mol)'>)+
We could then use fes.find_minima
to analyze the FES acquired, and to label them in the plot.
fes.find_minima()
+
+# plot the minimas
+fes.plot_minima(mark_color="white", png_name=None, style='ticks', context='notebook')
+
(<Figure size 960x672 with 2 Axes>, + <Axes: xlabel='CV1 - phi', ylabel='CV2 - psi'>)+
Local minima are stored in pandas.DataFrame
, as shown in below.
fes.minima
+
+ | 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 could draw the time-dependent profile of free energies of local minima from FEProfile
.
from catflow.metad.profile import FreeEnergyProfile
+fe_profile = FreeEnergyProfile(fes, h2)
+fe_profile.plot(energy_unit="kJ/mol", cmap="viridis", style='ticks', context='notebook')
+
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
+
#from catflow.metad.string import StringMethod
+s = StringMethod(fes)
+
Load the minima from the upper block.
+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.
+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.
+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:
+s.plot_string_evolution(cmap="RdYlBu_r", levels=20, dpi=96)
+
(<Figure size 960x672 with 2 Axes>, + <Axes: xlabel='CV1 - phi', ylabel='CV2 - psi'>)+
And draw the final MEP on the plot:
+s.plot_mep(cmap="RdYlBu_r", levels=20, dpi=96)
+
(<Figure size 960x672 with 2 Axes>, + <Axes: xlabel='CV1 - phi', ylabel='CV2 - psi'>)+
Last but not least, the free energy profile of the MEP:
+s.plot_mep_energy_profile(dpi=96)
+
(<Figure size 614.4x460.8 with 1 Axes>, + <Axes: xlabel='Reaction coordinate', ylabel='Free Energy (kJ/mol)'>)+
{"use strict";/*!
+ * escape-html
+ * Copyright(c) 2012-2013 TJ Holowaychuk
+ * Copyright(c) 2015 Andreas Lubbe
+ * Copyright(c) 2015 Tiancheng "Timothy" Gu
+ * MIT Licensed
+ */var Ha=/["'&<>]/;Nn.exports=$a;function $a(e){var t=""+e,r=Ha.exec(t);if(!r)return t;var o,n="",i=0,s=0;for(i=r.index;i