技術(shù)背景

在前面的幾篇博客中，我們分別介紹過 Histogram算法的使用 、 Plumed安裝與簡單使用 。Plumed一般就是兩種用法：要么在運行分子動力學(xué)模擬的過程中實時的對接，要么就是把分子模擬的相關(guān)軌跡保存下來，然后再使用Plumed進行后分析，本文介紹的是后面這種使用方法。

軌跡準備

做后分析，我們要先準備一手軌跡。比如我們做Histogram，那么就需要保留一條CV的軌跡，或者說反應(yīng)坐標的軌跡。一般為了歸一化的需求，我們可能還需要保留反應(yīng)坐標所對應(yīng)的單點能，或者稱之為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定義一個歸一化的系數(shù)rw，對應(yīng)的溫度為300K
4~13. 定義Histogram參數(shù)，使用高斯核，區(qū)間最小值為0.3，區(qū)間最大值為1.65，區(qū)間內(nèi)分50個格子，波包帶寬為0.05，使用rw進行歸一化
14. 將輸出數(shù)據(jù)保存到名為histo的文件內(nèi)

運行Plumed

安裝好plumed以后，之間在對應(yīng)文件的目錄下運行：

$ 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 文件，內(nèi)容為：

#! 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

這個結(jié)果的三列數(shù)據(jù)分別為：cv值、密度值和標準差，對于我們而言，主要關(guān)注下前兩列的結(jié)果即可，可以寫一個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()

輸出結(jié)果為：

總結(jié)概要

Plumed是一個強大的分子模擬數(shù)據(jù)處理工具，可以在模擬的過程中逐步分析，也可以保存模擬的軌跡做后分析。本文緊接前面的“增強采樣軟件PLUMED的安裝與使用”文章，還有“直方圖與核密度估計”文章。介紹了如何使用Plumed后分析工具，對輸出的反應(yīng)坐標的軌跡進行核密度估計。

版權(quán)聲明

本文首發(fā)鏈接為： https://www.cnblogs.com/dechinphy/p/plumed-post.html

作者ID：DechinPhy

更多原著文章： https://www.cnblogs.com/dechinphy/

請博主喝咖啡： https://www.cnblogs.com/dechinphy/gallery/image/379634.html