技術背景

增強取樣（Enhanced Sampling）是一種在分子動力學模擬中常用的技術，其作用是幫助我們更加快速的在時間軸上找到儘可能多的體系結構及其對應的能量。比如一個氫氣的燃燒反應，在中間過程中會產生眾多的反應產物，但是我們光從結果來看的話，就是從\(H_2\)加\(O_2\)，經過燃燒以後變成了\(H_2O\)，那麼中間的過程就被我們全部忽略了。實際上化學反應是一個整體的過程，是非常複雜的，如果要控制一個化學反應的過程，或者製備某個中間態的產物，或者需要提升最終預期產物的佔比，我們不得不考慮整個化學反應的路徑。

由於分子動力學模擬是基於牛頓力學的，因此中間發生的位移和碰撞我們都認為是連續的過程，但是為了在時間軸上更快的讀取資料，我們不能無限制的對時間進行分割，只能通過對時間進行取樣來獲得一個近似的過程。最典型的就是均勻取樣，比如每隔10ps採一個樣本點。但是這種方法在保障計算效率的前提下，很容易忽略了中間過程中出現時間極其短暫的一個體系狀態。因此就需要使用到增強抽樣的方法，對於高簡併度的狀態，我們降低其被取樣的概率，而對於低簡併度的狀態，我們提升其被採到的概率。常見的方法有：Meta Dynamics、VES和ITS等。這裡我們探索一下分子動力學模擬軟體PLUMED的安裝，該軟體已經集成了很大一部分的CV和增強取樣的方法。

PLUMED的下載與安裝

首先訪問其官方下載網站找到一個合適的版本進行下載，比如這裡我是直接下載的最新的版本：



下載完成後可以用tar -xvf plumed-2.7.1.tgz指令來進行解壓縮，可以看到解壓後的目錄如下所示：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed$ ll

總用量 103736

drwxrwxr-x  3 dechin dechin      4096 7月  12 09:21 ./

drwxrwxr-x 10 dechin dechin      4096 7月  12 09:20 ../

drwxrwxr-x 18 dechin dechin      4096 4月  16 14:48 plumed-2.7.1/

-rw-rw-r--  1 dechin dechin 106210796 7月  12 09:21 plumed-2.7.1.tgz

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed$ cd plumed-2.7.1/

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/plumed-2.7.1$ ll

總用量 504

drwxrwxr-x 18 dechin dechin   4096 4月  16 14:48 ./

drwxrwxr-x  3 dechin dechin   4096 7月  12 09:21 ../

drwxrwxr-x  6 dechin dechin   4096 4月  16 14:48 astyle/

-rw-rw-r--  1 dechin dechin     71 4月  16 14:48 .astyle.options

drwxrwxr-x  2 dechin dechin   4096 4月  16 14:48 CHANGES/

-rw-rw-r--  1 dechin dechin     94 4月  16 14:48 .codecov.yml

drwxrwxr-x  4 dechin dechin   4096 4月  16 14:48 conda/

-rwxrwxr-x  1 dechin dechin 323087 4月  16 14:48 configure*

-rw-rw-r--  1 dechin dechin  43692 4月  16 14:48 configure.ac

-rw-rw-r--  1 dechin dechin   7639 4月  16 14:48 COPYING.LESSER

drwxrwxr-x  2 dechin dechin   4096 4月  16 14:48 developer-doc/

drwxrwxr-x  2 dechin dechin   4096 4月  16 14:48 docker/

drwxrwxr-x  3 dechin dechin   4096 4月  16 14:48 .github/

-rw-rw-r--  1 dechin dechin    179 4月  16 14:48 .gitignore

-rw-rw-r--  1 dechin dechin    245 4月  16 14:48 .lgtm.yml

drwxrwxr-x  2 dechin dechin   4096 4月  16 14:48 macports/

-rw-rw-r--  1 dechin dechin   2485 4月  16 14:48 Makefile

-rw-rw-r--  1 dechin dechin   1200 4月  16 14:48 Makefile.conf.in

drwxrwxr-x  7 dechin dechin   4096 4月  16 14:48 patches/

-rw-rw-r--  1 dechin dechin    527 4月  16 14:48 PEOPLE

drwxrwxr-x  3 dechin dechin   4096 4月  16 14:48 python/

-rw-rw-r--  1 dechin dechin   8146 4月  16 14:48 README.md

drwxrwxr-x 26 dechin dechin   4096 4月  16 14:48 regtest/

-rwxrwxr-x  1 dechin dechin   5397 4月  16 14:48 release.sh*

drwxrwxr-x  2 dechin dechin   4096 4月  16 14:48 scripts/

-rw-rw-r--  1 dechin dechin    319 4月  16 14:48 sourceme.sh.in

drwxrwxr-x 45 dechin dechin   4096 4月  16 14:48 src/

-rw-rw-r--  1 dechin dechin      0 4月  16 14:48 stamp-h.in

drwxrwxr-x  6 dechin dechin   4096 4月  16 14:48 test/

drwxrwxr-x  2 dechin dechin   4096 4月  16 14:48 .travis/

-rw-rw-r--  1 dechin dechin  10558 4月  16 14:48 .travis.yml

drwxrwxr-x  4 dechin dechin   4096 4月  16 14:48 user-doc/

-rw-rw-r--  1 dechin dechin    322 4月  16 14:48 VERSION

drwxrwxr-x  2 dechin dechin   4096 4月  16 14:48 vim/

直接在這個plumed-2.7.1（注意對應自己安裝的版本）目錄下執行如下的指令進行安裝：

$ ./configure --prefix=/usr/local

$ make -j 4

$ sudo make install

安裝驗證

安裝完成後，可以在命令列中直接使用plumed指令，可以用如下幫助指令，確認軟體是否被安裝成功：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples$ plumed help

Usage: plumed [options] [command] [command options]

  plumed [command] -h|--help: to print help for a specific command

Options:

  [help|-h|--help]          : to print this help

  [--is-installed]          : fails if plumed is not installed

  [--has-mpi]               : fails if plumed is running without MPI

  [--has-dlopen]            : fails if plumed is compiled without dlopen

  [--load LIB]              : loads a shared object (typically a plugin library)

  [--standalone-executable] : tells plumed not to look for commands implemented as scripts

Commands:

  plumed completion : dump a function usable for programmable completion

  plumed driver : analyze trajectories with plumed

  plumed driver-float : analyze trajectories with plumed (single precision version)

  plumed gen_example : construct an example for the manual that users can interact with

  plumed gentemplate : print out a template input for a particular action

  plumed info : provide informations about plumed

  plumed kt : print out the value of kT at a particular temperature

  plumed manual : print out a description of the keywords for an action in html

  plumed pathtools : print out a description of the keywords for an action in html

  plumed pdbrenumber : Modify atom numbers in a PDB, possibly using hybrid-36 coding

  plumed pesmd : Langevin dynamics on PLUMED energy landscape

  plumed simplemd : run lj code

  plumed sum_hills : sum the hills with  plumed

  plumed patch : patch an MD engine

  plumed vim2html : convert plumed input file to colored html using vim syntax

  plumed selector : create lists of serial atom numbers

  plumed config : inquire plumed about how it was configure

  plumed newcv : create a new collective variable from a template

  plumed mklib : compile a .cpp file into a shared library

  plumed partial_tempering : scale parameters in a gromacs topology to implement solute or partial tempering

簡單示例

使用plumed，需要有一個完整的反應路徑檔案，這裡我們可以用wget直接下載plumed倉庫中的示例軌跡檔案：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples$ wget https://plumed.github.io/doc-v2.4/user-doc/html/tutorial-resources/trieste-1.tar.gz

--2021-07-12 10:18:18--  https://plumed.github.io/doc-v2.4/user-doc/html/tutorial-resources/trieste-1.tar.gz

正在解析主機 plumed.github.io (plumed.github.io)... ::1, 185.199.109.153, 185.199.108.153, ...

正在連線 plumed.github.io (plumed.github.io)|::1|:443... 失敗：拒絕連線。

正在連線 plumed.github.io (plumed.github.io)|185.199.109.153|:443... 已連線。

已發出 HTTP 請求，正在等待迴應... 301 Moved Permanently

位置：https://www.plumed.org/doc-v2.4/user-doc/html/tutorial-resources/trieste-1.tar.gz [跟隨至新的 URL]

--2021-07-12 10:18:19--  https://www.plumed.org/doc-v2.4/user-doc/html/tutorial-resources/trieste-1.tar.gz

正在解析主機 www.plumed.org (www.plumed.org)... ::1, 185.199.109.153, 185.199.108.153, ...

正在連線 www.plumed.org (www.plumed.org)|::1|:443... 失敗：拒絕連線。

正在連線 www.plumed.org (www.plumed.org)|185.199.109.153|:443... 已連線。

已發出 HTTP 請求，正在等待迴應... 200 OK

長度： 335374 (328K) [application/gzip]

正在儲存至: “trieste-1.tar.gz”

trieste-1.tar.gz    100%[==================>] 327.51K  --.-KB/s    用時 0.04s 

2021-07-12 10:18:20 (8.02 MB/s) - 已儲存 “trieste-1.tar.gz” [335374/335374])

同樣的執行解壓：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples$ ll

總用量 336

drwxrwxr-x 2 dechin dechin   4096 7月  12 10:18 ./

drwxrwxr-x 4 dechin dechin   4096 7月  12 09:38 ../

-rw-rw-r-- 1 dechin dechin 335374 11月 18  2020 trieste-1.tar.gz

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples$ tar -xvf trieste-1.tar.gz

trieste-1/

trieste-1/.solutions/

trieste-1/.solutions/ref-rna.pdb

trieste-1/.solutions/plumed-ex2c.dat

trieste-1/.solutions/plumed-ex2.dat

trieste-1/.solutions/plumed-ex2b.dat

trieste-1/.solutions/plumed-ex1b.dat

trieste-1/.solutions/plumed-ex1.dat

trieste-1/traj-broken.xtc

trieste-1/ref.pdb

trieste-1/traj-whole.xtc

解壓後可以在目錄下看到一個pdb的配置檔案和兩個xtc的軌跡檔案，這個是用gromacs生成的資料檔案格式：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples/trieste-1$ ll

總用量 752

drwxr-xr-x 3 dechin dechin   4096 11月 18  2020 ./

drwxrwxr-x 3 dechin dechin   4096 7月  12 10:19 ../

-rw-r--r-- 1 dechin dechin 519979 11月 18  2020 ref.pdb

drwxr-xr-x 2 dechin dechin   4096 11月 18  2020 .solutions/

-rw-r--r-- 1 dechin dechin 117620 11月 18  2020 traj-broken.xtc

-rw-r--r-- 1 dechin dechin 116036 11月 18  2020 traj-whole.xtc

在具備了軌跡檔案之後，就可以按照官方的文件示例來編寫plumed.dat配置輸入檔案，用於定義需要計算的內容，比如這裡定義的是計算1號原子和2號原子之間的距離，並每10個step將其寫入到名為colvar的檔案下：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples/trieste-1$ cat plumed.dat

d1: DISTANCE ATOMS=1,2

PRINT ARG=d1 FILE=colvar STRIDE=10

執行plumed指令得到的結果如下：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples/trieste-1$ plumed driver --mf_xtc traj-whole.xtc --plumed plumed.dat 

DRIVER: Found molfile format trajectory xtc with name traj-whole.xtc

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: 6580

PLUMED: File suffix:

PLUMED: FILE: plumed.dat

PLUMED: Action DISTANCE

PLUMED:   with label d1

PLUMED:   between atoms 1 2

PLUMED:   using periodic boundary conditions

PLUMED: Action PRINT

PLUMED:   with label @1

PLUMED:   with stride 10

PLUMED:   with arguments d1

PLUMED:   on file colvar

PLUMED:   with format  %f

PLUMED: 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.003111     0.003111     0.003111     0.003111

PLUMED: 1 Prepare dependencies                             5     0.000005     0.000001     0.000000     0.000003

PLUMED: 2 Sharing data                                     1     0.001410     0.001410     0.001410     0.001410

PLUMED: 3 Waiting for data                                 1     0.000022     0.000022     0.000022     0.000022

PLUMED: 4 Calculating (forward loop)                       1     0.000017     0.000017     0.000017     0.000017

PLUMED: 5 Applying (backward loop)                         1     0.000047     0.000047     0.000047     0.000047

PLUMED: 6 Update                                           1     0.000052     0.000052     0.000052     0.000052

在這個給定的軌跡檔案中其實包含了5個step，但是由於我們設定了每10個step列印一次，因此最終儲存到colvar中的資料只有1個值：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples/trieste-1$ ll

總用量 760

drwxr-xr-x 3 dechin dechin   4096 7月  12 10:22 ./

drwxrwxr-x 3 dechin dechin   4096 7月  12 10:19 ../

-rw-rw-r-- 1 dechin dechin     37 7月  12 10:22 colvar

-rw-rw-r-- 1 dechin dechin     60 7月  12 10:21 plumed.dat

-rw-r--r-- 1 dechin dechin 519979 11月 18  2020 ref.pdb

drwxr-xr-x 2 dechin dechin   4096 11月 18  2020 .solutions/

-rw-r--r-- 1 dechin dechin 117620 11月 18  2020 traj-broken.xtc

-rw-r--r-- 1 dechin dechin 116036 11月 18  2020 traj-whole.xtc

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples/trieste-1$ cat colvar

#! FIELDS time d1

 0.000000 0.095760

但是如果我們把上面的儲存步長STRIDE修改為1的話，得到的colvar就是如下所示的結果：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples/trieste-1$ cat colvar

#! FIELDS time d1

 0.000000 0.095760

 1.000000 0.096845

 2.000000 0.095885

 3.000000 0.096010

 4.000000 0.095321

這裡是將5個step中的所有對應的原子距離都返回了出來。

軌跡檔案的格式轉換

上面用到的xtc檔案是GROMACS生成的軌跡副檔名，由於是二進位制檔案並不方便讀取，這裡我們可以將其轉換成hdf5的格式，然後就可以用python直接來讀取其中的資料。這裡我們需要藉助於mdtraj這個工具。

用pip安裝mdtraj

經過嘗試用conda來裝失敗了，因此建議直接使用pip來進行安裝：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples/trieste-1$ python3.9 -m pip install mdtraj

Collecting mdtraj

  Downloading mdtraj-1.9.6-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.whl (6.0 MB)

     |████████████████████████████████| 6.0 MB 1.1 MB/s

Collecting astunparse

  Using cached astunparse-1.6.3-py2.py3-none-any.whl (12 kB)

Requirement already satisfied: numpy>=1.6 in /home/dechin/anaconda3/envs/AmberTools21/lib/python3.9/site-packages (from mdtraj) (1.21.0)

Requirement already satisfied: scipy in /home/dechin/anaconda3/envs/AmberTools21/lib/python3.9/site-packages (from mdtraj) (1.7.0)

Requirement already satisfied: pyparsing in /home/dechin/anaconda3/envs/AmberTools21/lib/python3.9/site-packages (from mdtraj) (2.4.7)

Requirement already satisfied: six<2.0,>=1.6.1 in /home/dechin/anaconda3/envs/AmberTools21/lib/python3.9/site-packages (from astunparse->mdtraj) (1.16.0)

Requirement already satisfied: wheel<1.0,>=0.23.0 in /home/dechin/anaconda3/envs/AmberTools21/lib/python3.9/site-packages (from astunparse->mdtraj) (0.36.2)

Installing collected packages: astunparse, mdtraj

Successfully installed astunparse-1.6.3 mdtraj-1.9.6

安裝完成後會生成一個名為mdconvert的可執行檔案，一般情況是在系統路徑下可以直接指令呼叫的，下述是mdconvert的幫助文件：

$ mdconvert -h

usage: mdconvert [-h]

                 -o OUTPUT         # 必需，指定輸出檔案

                 [-c CHUNK]        # 可選，指定一次讀入記憶體的幀數，預設1000.

                 [-f]              # 可選，如果輸出檔案已存在，則強制覆蓋

                 [-s STRIDE]       # 可選，只加載 every stride-th frame

                 [-i INDEX]        # 可選，更加靈活地指定載入特定的幀(Python格式)，eg： -i N; -i -1; -i N:M

                 [-a ATOM_INDICES] # 可選，通過檔案指定只加載特定的某些原子，檔案內容為原子序號(從0開始計數)

                 [-t TOPOLOGY]     # 可選，指定拓撲檔案(PDB/prmtop 格式)，一旦指定即可將軌跡輸出為pdb格式。

                 input [input ...] # 必需，指定輸入檔案(可以是多條軌跡)

通過mdconvert，可以將上面的案例中所提到的xtc檔案配合pdb檔案轉化成hdf5格式的檔案：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples/trieste-1$ mdconvert -o test.h5 -t ref.pdb traj-whole.xtc

Warning: 'step' data from input file(s) will be discarded. output format only supports fields: 'xyz', 'time', 'cell_lengths', 'cell_angles', 'velocities', 'kineticEnergy', 'potentialEnergy', 'temperature', 'lambda', 'topology'

converted 5 frames, 6580 atoms

而讀取hdf5的檔案，可以用python中的h5py來實現，沒有安裝h5py的可以通過pip簡單的安裝一下：

(base) dechin@ubuntu2004:~/projects/gitlab/dechin/src/plumed/examples/trieste-1$ python3 -m pip install h5py==2.9

Collecting h5py==2.9

  Downloading h5py-2.9.0-cp38-cp38-manylinux1_x86_64.whl (2.8 MB)

     |████████████████████████████████| 2.8 MB 997 kB/s

Requirement already satisfied: numpy>=1.7 in /home/dechin/anaconda3/lib/python3.8/site-packages (from h5py==2.9) (1.20.2)

Requirement already satisfied: six in /home/dechin/.local/lib/python3.8/site-packages (from h5py==2.9) (1.16.0)

Installing collected packages: h5py

  Attempting uninstall: h5py

    Found existing installation: h5py 3.2.1

    Uninstalling h5py-3.2.1:

      Successfully uninstalled h5py-3.2.1

Successfully installed h5py-2.10.0

具體的h5py的使用方法這裡就不進行展開了，可以參考官方的使用文件。

總結概要

本文作為一個入門級的文章，主要介紹了分子動力學模擬中增強取樣的基本概念與相應工具的安裝和使用。PLUMED是業界比較出名的一款增強取樣開源軟體，能夠對接多個分子動力學模擬軟體，如GROMACS等，並利用這些軟體生成的路徑資訊來進行取樣。並且為了可以在python上也能看到路徑資訊等重要資料，可以考慮使用mdconvert將路徑資料轉化成python上常用的hdf5格式並用h5py進行讀寫。

版權宣告

本文首發連結為：https://www.cnblogs.com/dechinphy/p/plumed.html

作者ID：DechinPhy

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

打賞專用連結：https://www.cnblogs.com/dechinphy/gallery/image/379634.html

騰訊雲專欄同步：https://cloud.tencent.com/developer/column/91958

參考連結

1. https://www.jianshu.com/p/edbcd57ffa03
2. https://blog.chembiosim.com/alter-md-trajectory-file-format-with-mdconvert/
3. https://blog.chembiosim.com/plumed-basics-01/