Plumed分子模拟后分析

技术背景

在前面的几篇博客中，我们分别介绍过Histogram算法的使用、Plumed安装与简单使用。Plumed一般就是两种用法：要么在运行分子动力学模拟的过程中实时的对接，要么就是把分子模拟的相关轨迹保存下来，然后再使用Plumed进行后分析，本文介绍的是后面这种使用方法。

轨迹准备

做后分析，我们要先准备一手轨迹。比如我们做Histogram，那么就需要保留一条CV的轨迹，或者说反应坐标的轨迹。一般为了归一化的需求，我们可能还需要保留反应坐标所对应的单点能，或者称之为Bias偏置势。如下是一条轨迹的示例record_cv_bias.txt，含有100个点：

#! FIELDS time rbias dcomb
0 1.0 0.722307
1 1.0 0.929853
2 1.0 1.046353
3 1.0 0.849326
4 1.0 0.635665
5 1.0 1.133585
6 1.0 1.602735
7 1.0 1.11138
8 1.0 0.815045
9 1.0 0.744088
10 1.0 0.631533
11 1.0 1.006049
12 1.0 0.507855
13 1.0 0.620414
14 1.0 1.43557
15 1.0 0.798832
16 1.0 1.007135
17 1.0 0.684447
18 1.0 1.073844
19 1.0 0.952023
20 1.0 0.754293
21 1.0 0.530406
22 1.0 1.078594
23 1.0 1.325044
24 1.0 1.418314
25 1.0 1.110357
26 1.0 0.964378
27 1.0 0.893131
28 1.0 1.473515
29 1.0 1.103729
30 1.0 1.019812
31 1.0 1.037889
32 1.0 1.092906
33 1.0 0.602312
34 1.0 1.016394
35 1.0 0.845226
36 1.0 1.210281
37 1.0 0.744589
38 1.0 0.666467
39 1.0 1.276913
40 1.0 0.976801
41 1.0 0.92263
42 1.0 1.386575
43 1.0 1.243241
44 1.0 1.442755
45 1.0 1.284636
46 1.0 0.756184
47 1.0 1.162758
48 1.0 1.177665
49 1.0 0.717487
50 1.0 1.193379
51 1.0 0.798996
52 1.0 1.045093
53 1.0 1.489541
54 1.0 1.067574
55 1.0 1.10109
56 1.0 1.135074
57 1.0 1.049557
58 1.0 0.362635
59 1.0 1.030856
60 1.0 1.323538
61 1.0 1.405822
62 1.0 0.935292
63 1.0 0.868032
64 1.0 0.946401
65 1.0 1.468123
66 1.0 1.062565
67 1.0 1.05637
68 1.0 0.962974
69 1.0 1.50403
70 1.0 0.95933
71 1.0 1.218142
72 1.0 1.303102
73 1.0 0.876645
74 1.0 1.188313
75 1.0 1.31572
76 1.0 0.693456
77 1.0 0.54377
78 1.0 1.026873
79 1.0 1.3925
80 1.0 1.317733
81 1.0 0.972162
82 1.0 1.373311
83 1.0 1.244547
84 1.0 1.00565
85 1.0 1.180062
86 1.0 1.221193
87 1.0 1.046702
88 1.0 1.409161
89 1.0 1.132955
90 1.0 0.428334
91 1.0 0.890674
92 1.0 0.63586
93 1.0 0.997099
94 1.0 0.969676
95 1.0 1.280118
96 1.0 1.19793
97 1.0 1.112535
98 1.0 1.388056
99 1.0 0.946911

有了轨迹之后，写一个简单的Plumed配置文件plumed.dat：

cv: READ FILE=record_cv_bias.txt VALUES=dcomb IGNORE_TIME
w: READ FILE=record_cv_bias.txt VALUES=rbias IGNORE_TIME
rw: REWEIGHT_BIAS ARG=w TEMP=300
hh: HISTOGRAM ...
   ARG=cv
   KERNEL=gaussian
   GRID_MIN=0.3
   GRID_MAX=1.65
   GRID_BIN=50
   BANDWIDTH=0.05
   NORMALIZATION=true
   LOGWEIGHTS=rw
...
DUMPGRID GRID=hh FILE=histo

这个文件的逐行释义为：

1. 读取record_cv_bias.txt文件中标签为dcomb的一列，作为cv
2. 读取record_cv_bias.txt文件中标签为rbias的一列，作为w
3. 使用w定义一个归一化的系数rw，对应的温度为300K
4~13. 定义Histogram参数，使用高斯核，区间最小值为0.3，区间最大值为1.65，区间内分50个格子，波包带宽为0.05，使用rw进行归一化
14. 将输出数据保存到名为histo的文件内

运行Plumed

安装好plumed以后，之间在对应文件的目录下运行：

$ plumed driver --plumed plumed.dat --noatoms
PLUMED: PLUMED is starting
PLUMED: Version: 2.7.1 (git: Unknown) compiled on Jul 12 2021 at 09:24:30
PLUMED: Please cite these papers when using PLUMED [1][2]
PLUMED: For further information see the PLUMED web page at http://www.plumed.org
PLUMED: Root: /usr/local/lib/plumed
PLUMED: For installed feature, see /usr/local/lib/plumed/src/config/config.txt
PLUMED: Molecular dynamics engine: driver
PLUMED: Precision of reals: 8
PLUMED: Running over 1 node
PLUMED: Number of threads: 1
PLUMED: Cache line size: 512
PLUMED: Number of atoms: 0
PLUMED: File suffix: 
PLUMED: FILE: plumed.dat
PLUMED: Action READ
PLUMED:   with label cv
PLUMED:   with stride 1
PLUMED:   reading data from file record_cv_bias.txt
PLUMED:   reading value dcomb and storing as cv
PLUMED: Action READ
PLUMED:   with label w
PLUMED:   with stride 1
PLUMED:   reading data from file record_cv_bias.txt
PLUMED:   reading value rbias and storing as w
PLUMED: Action REWEIGHT_BIAS
PLUMED:   with label rw
PLUMED:   with arguments w
PLUMED: Action HISTOGRAM
PLUMED:   with label hh
PLUMED:   with stride 1
PLUMED:   with arguments cv
PLUMED:   reweighting using weights from rw 
PLUMED:   grid of 51 equally spaced points between (0.3) and (1.65)
PLUMED: Action DUMPGRID
PLUMED:   with label @4
PLUMED:   with stride 0
PLUMED:   outputting grid calculated by action hh to file named histo with format %f 
PLUMED:   outputting data for replicas 0 END FILE: plumed.dat
PLUMED: Timestep: 1.000000
PLUMED: KbT has not been set by the MD engine
PLUMED: It should be set by hand where needed
PLUMED: Relevant bibliography:
PLUMED:   [1] The PLUMED consortium, Nat. Methods 16, 670 (2019)
PLUMED:   [2] Tribello, Bonomi, Branduardi, Camilloni, and Bussi, Comput. Phys. Commun. 185, 604 (2014)
PLUMED: Please read and cite where appropriate!
PLUMED: Finished setup
PLUMED:                                               Cycles        Total      Average      Minimum      Maximum
PLUMED:                                                    1     0.007538     0.007538     0.007538     0.007538
PLUMED: 1 Prepare dependencies                            99     0.000267     0.000003     0.000002     0.000008
PLUMED: 2 Sharing data                                    99     0.000007     0.000000     0.000000     0.000000
PLUMED: 3 Waiting for data                                99     0.000010     0.000000     0.000000     0.000000
PLUMED: 4 Calculating (forward loop)                      99     0.001478     0.000015     0.000014     0.000057
PLUMED: 5 Applying (backward loop)                        99     0.000130     0.000001     0.000001     0.000003
PLUMED: 6 Update                                          99     0.003019     0.000030     0.000014     0.000093

运行完成后，如果没有报错，会在当前目录下生成一个histo文件，内容为：

#! FIELDS cv hh dhh_cv
#! SET normalisation  146.331611
#! SET min_cv 0.3
#! SET max_cv 1.65
#! SET nbins_cv  50
#! SET periodic_cv false
 0.300000 0.040185 1.087032
 0.327000 0.073737 1.333674
 0.354000 0.108112 1.138788
 0.381000 0.132720 0.673529
 0.408000 0.146143 0.387485
 0.435000 0.158725 0.647624
 0.462000 0.185775 1.399849
 0.489000 0.233253 2.034814
 0.516000 0.290241 2.109999
 0.543000 0.347304 2.200385
 0.570000 0.415370 2.931229
 0.597000 0.504292 3.503145
 0.624000 0.592054 2.763309
 0.651000 0.645747 1.208607
 0.678000 0.663637 0.297321
 0.705000 0.670123 0.268624
 0.732000 0.677724 0.233552
 0.759000 0.679871 -0.077180
 0.786000 0.677381 0.020735
 0.813000 0.688778 0.956408
 0.840000 0.733653 2.431984
 0.867000 0.822586 4.209499
 0.894000 0.963284 6.230239
 0.921000 1.155414 7.841243
 0.948000 1.373556 8.056159
 0.975000 1.575612 6.691082
 1.002000 1.725112 4.227284
 1.029000 1.795526 0.856452
 1.056000 1.767797 -2.871105
 1.083000 1.648824 -5.650900
 1.110000 1.480563 -6.442901
 1.137000 1.318483 -5.309715
 1.164000 1.198476 -3.600376
 1.191000 1.115220 -2.744803
 1.218000 1.043134 -2.600741
 1.245000 0.977289 -2.155749
 1.272000 0.928752 -1.443136
 1.299000 0.894999 -1.109955
 1.326000 0.868484 -0.743408
 1.353000 0.859332 0.120417
 1.380000 0.869809 0.444883
 1.407000 0.867177 -0.906627
 1.434000 0.811246 -3.257711
 1.461000 0.694147 -5.274589
 1.488000 0.536089 -6.215685
 1.515000 0.370839 -5.747150
 1.542000 0.236681 -4.040289
 1.569000 0.154469 -2.128090
 1.596000 0.112702 -1.142470
 1.623000 0.083908 -1.071380
 1.650000 0.053970 -1.101246

这个结果的三列数据分别为：cv值、密度值和标准差，对于我们而言，主要关注下前两列的结果即可，可以写一个Python脚本做一下简单的可视化：

import numpy as np
with open('histo', 'r') as hfile:
    hlines = hfile.readlines()
hcv = []
hbar = []
for hline in hlines[6:]:
    hcv.append(float(hline.strip().split()[0]))
    hbar.append(float(hline.strip().split()[1]))
hcv = np.array(hcv)
hbar = np.array(hbar)
from matplotlib import pyplot as plt
plt.figure()
plt.plot(hcv, hbar, color='black')
plt.show()

输出结果为：