MindSponge分子動力學模擬——定義一個分子系統（2023.08）
2023-08-30 18:00:48
技術背景
在前面兩篇文章中，我們分別介紹了分子動力學模擬軟體MindSponge的軟體架構和安裝與使用教學。這裡我們進入到實用化階段，假定大家都已經在本地部署好了基於MindSpore的MindSponge的程式設計環境，開始用MindSponge去做一些真正的分子模擬的工作。那麼分子模擬的第一步，我們就需要在MindSponge中去定義一個分子系統Molecule()。
基礎類Molecule的解析
我們先來看一下原始碼中的Molecule這個類的自我介紹：
class Molecule(Cell):
r"""
Base class for molecular system, used as the "system module" in MindSPONGE.
The `Molecule` Cell can represent a molecule or a system consisting of multiple molecules.
The major components of the `Molecule` Cell is the `Residue` Cell. A `Molecule` Cell can
contain multiple `Residue` Cells.
Args:
atoms(Union[List[Union[str, int]], ndarray]): Array of atoms. The data in array can be str of atom
name or int of atomic number. Defulat: None
atom_name(Union[List[str], ndarray]): Array of atom name with data type `str`. Defulat: None
atom_type(Union[List[str], ndarray]): Array of atom type with data type `str`. Defulat: None
atom_mass(Union[Tensor, ndarray, List[float]]): Array of atom mass of shape `(B, A)` with data type
`float`. Defulat: None
atom_charge(Union[Tensor, ndarray, List[float]]): Array of atom charge of shape `(B, A)` with data type
`float`. Defulat: None
atomic_number(Union[Tensor, ndarray, List[float]]): Array of atomic number of shape `(B, A)` with data type
`int`. Defulat: None
bond(Union[Tensor, ndarray, List[int]]): Array of bond connection of shape `(B, b, 2)` with data
type `int`. Defulat: None
coordinate(Union[Tensor, ndarray, List[float]]): Tensor of atomic coordinates :math:`R` of shape
`(B, A, D)` with data type `float`. Default: None
pbc_box(Union[Tensor, ndarray, List[float]]): Tensor of box size :math:`\vec{L}` of periodic boundary
condition (PBC). The shape of tensor is `(B, D)`,
and the data type is `float`. Default: None
template(Union[dict, str, List[Union[dict, str]]]): Template for molecule. It can be a `dict` in MindSPONGE
template format or a `str` for the filename of a
MindSPONGE template file. If a `str` is given,
it will first look for a file with the same name in the
current directory. If the file does not exist, it will
search in the built-in template directory of
MindSPONGE (`mindsponge.data.template`).
Default: None.
residue(Union[Residue, List[Residue]]): Residue or a list of residues. If template is not None,
only the residues in the template will be used.
Default: None.
length_unit(str): Length unit. If `None` is given, the global length
units will be used. Default: None
Outputs:
- coordinate, Tensor of shape `(B, A, D)`. Data type is float.
- pbc_box, Tensor of shape `(B, D)`. Data type is float.
Supported Platforms:
``Ascend`` ``GPU``
Symbols:
B: Batchsize, i.e. number of walkers in simulation
A: Number of atoms.
b: Number of bonds.
D: Spatial dimension of the simulation system. Usually is 3.
"""
可以先看一下Molecule所接收的資訊，其實可以主要分為以下幾大類別：
- 原子特徵資訊。用於區分不同原子之間的差異性，比如atom_name原子名稱、atom_type原子型別、atomic_number原子序數等。
- 拓撲資訊。在構建Molecule的時候需要傳入鍵連資訊bond，否則不帶鍵連關係的Molecule計算出來的力場能量是錯誤的。
- 構象資訊。主要是原子座標coordinate和週期性邊界條件pbc_box，作為近鄰表計算和力場能量計算的輸入，但不作為拓撲連線資訊的輸入。
- 模組化資訊。除了逐個原子的去構建一個Molecule，還可以定義好一系列完整的殘基Residue再輸入給Molecule進行構建，或者通過模板template來進行構建。
- 單位資訊。主要包含長度單位length_unit和能量單位energy_unit。
上述主要是給Molecule的輸入資訊，輸入給Molecule之後在內部構建build一次，才能得到一個最終的分子系統物件。接下來看看構建之後的Molecule的一些重要內建屬性（self.xxx）：
- 原子特徵資訊。除了上述傳入的那些資訊之外，還有原子數num_atoms，batch數量num_walker以及靈活的維度數量dimension。除了每個原子的基本型別外，還儲存了一個heavy_atom_mask重原子的資訊，便於快速區分重原子和氫原子。
- 拓撲資訊。除了鍵連關係bonds資訊，還有h_bonds氫原子成鍵的資訊。
- 構象資訊。主要就是coordinate原子座標，因為需要在Updater中更新迭代，因此這裡的coordinate需要是一個Parameter的型別，而不是普通的Tensor。
- 模組化資訊。在構建的過程中，對傳入的Residue也都進行了extend，因此最終Residue內部的這些資訊，都會被合併到前面提到的Molecule的原子特徵資訊和拓撲資訊、構象資訊中，同時會保留一個atom_resid用於追溯原子所在的residue。如果在template模板中有設定一些約束限制，比如settle約束演演算法相關的引數settle_index和settle_length，也會儲存在Molecule的屬性中，用於後續約束演演算法的計算。
- 單位資訊。units把相關的單位資訊都儲存在一個Units物件中，支援從global units中呼叫，可以隨時呼叫。
除了內建屬性，Molecule還有一些內建函數可以關注一下：
- 單位轉化。主要是convert_length_from和convert_length_to兩個函數，用於執行長度單位的變換。
- 系統擴充套件函數。比如copy函數，可以用於將本系統拷貝一份，但是該拷貝的過程會生成一個新的物件，而不是原有的Molecule物件。但如果是多個的Molecule物件，可以用內建函數append進行合併。如果需要節省一些麻煩，想對系統進行擴充套件，可以直接使用內建函數reduplicate，在系統內部複製一份。類似於append的功能，可以使用內建函數add_residue來新增新的residue。上述幾種方法主要針對於非週期性的體系，如果是帶有周期性邊界條件的體系，直接使用repeat_box函數即可完成對體系的快速複製。
- 構建函數。一般情況下對於只是想做MD的童鞋而言，沒有必要使用到build_system構建系統和build_space構建構象這些函數，但是如果有需要自行調整Molecule的內容時，就需要重新build一次。
- 補媒介函數。一般給定的pdb檔案會丟失一些氫原子和溶劑分子的資訊，這些都可以在做模擬之前手動補上。目前MindSponge支援的是對體系加水分子fill_water，可以指定溶劑層的厚度，或者指定一個盒子的大小。
- 回撥函數。在深度學習或者MindSponge分子動力學模擬的過程中，我們會使用到回撥函數CallBack來對輸出內容進行追蹤。但是CallBack本身是不儲存任何體系相關的資訊的，因此追蹤的內容其實也是從Molecule和ForceField內部進行回撥。比如在Molecule中，可以呼叫get_atoms，get_coordinate，get_pbc_box等等函數，而如果直接使用MindSpore的Cell中所特有的construct函數，這裡也會返回coordinate和pbc_box兩個資訊，這些都可以認為是Molecule類的「回撥函數」。
從模板定義一個分子
關於MindSponge的安裝和使用，可以參考這篇文章，在這裡我們就不重複贅述了，假設大家已經完成了MindSponge的安裝。但是需要提一句的是，在開始MindSponge模擬前，最好在python指令碼的最前面加上這樣一些環境變數的設定，否則容易報錯：
import os
os.environ['GLOG_v']='4'
os.environ['MS_JIT_MODULES']='sponge'
接下來我們就可以簡單的使用模板檔案去建立一個新的分子：
from sponge import Molecule
system = Molecule(template='water.spce.yaml')
print ('The number of atoms in the system is: ', system.num_atoms)
print ('All the atom names in the system are: ', system.atom_name)
print ('The coordinates of atoms are: \n{}'.format(system.coordinate.asnumpy()))
輸出的結果如下所示：
The number of atoms in the system is: 3
All the atom names in the system are: [['O' 'H1' 'H2']]
The coordinates of atoms are:
[[[ 0. 0. 0. ]
[ 0.08164904 0.0577359 0. ]
[-0.08164904 0.0577359 0. ]]]
這裡因為water.spce.yaml是一個預置的模板，類似的還有water.tip3p.yaml。這種預置的模板我們可以直接當做template來建立，但如果是使用者自行定義的模板檔案，最好在這裡寫清楚yaml檔案的絕對路徑，否則會導致報錯。相關yaml檔案的內容如下所示：
template:
base: water_3p.yaml
WAT:
atom_mass: [15.9994, 1.008, 1.008]
atom_charge: [-0.8476, 0.4238, 0.4238]
settle:
mandatory: false
length_unit: nm
distance:
OW-HW: 0.1
HW-HW: 0.16330
molecule:
residue:
- WAT
length_unit: nm
coordinate:
- [0.0, 0.0, 0.0]
- [0.081649043, 0.057735897, 0.0]
- [-0.081649043, 0.057735897, 0.0]
這裡的base是指向了另外一個較為基礎的yaml引數檔案：
template:
WAT:
atom_name: [O, H1, H2]
atom_type: [OW, HW, HW]
atom_mass: [16.00, 1.008, 1.008]
atomic_number: [8, 1, 1]
bond:
- [0, 1]
- [0, 2]
head_atom: null
tail_atom: null
有了這些參考，使用者就可以自行定義一些模板，用於計算。
從檔案定義一個分子
MindSponge也支援一些特定格式的分子匯入，比如mol2格式的分子和pdb格式的蛋白質分子，這個章節介紹一下如何將檔案匯入為一個MindSponge的Molecule。比如我這裡有一個非常簡單的pdb格式的多肽鏈：
REMARK Generated By Xponge (Molecule)
ATOM 1 N ALA 1 -0.095 -11.436 -0.780
ATOM 2 CA ALA 1 -0.171 -10.015 -0.507
ATOM 3 CB ALA 1 1.201 -9.359 -0.628
ATOM 4 C ALA 1 -1.107 -9.319 -1.485
ATOM 5 O ALA 1 -1.682 -9.960 -2.362
ATOM 6 N ARG 2 -1.303 -8.037 -1.397
ATOM 7 CA ARG 2 -2.194 -7.375 -2.328
ATOM 8 CB ARG 2 -3.606 -7.943 -2.235
ATOM 9 CG ARG 2 -4.510 -7.221 -3.228
ATOM 10 CD ARG 2 -5.923 -7.789 -3.136
ATOM 11 NE ARG 2 -6.831 -7.111 -4.087
ATOM 12 CZ ARG 2 -8.119 -7.421 -4.205
ATOM 13 NH1 ARG 2 -8.686 -8.371 -3.468
ATOM 14 NH2 ARG 2 -8.844 -6.747 -5.093
ATOM 15 C ARG 2 -2.273 -5.882 -2.042
ATOM 16 O ARG 2 -1.630 -5.388 -1.119
ATOM 17 N ALA 3 -3.027 -5.119 -2.777
ATOM 18 CA ALA 3 -3.103 -3.697 -2.505
ATOM 19 CB ALA 3 -1.731 -3.041 -2.625
ATOM 20 C ALA 3 -4.039 -3.001 -3.483
ATOM 21 O ALA 3 -4.614 -3.643 -4.359
ATOM 22 N ALA 4 -4.235 -1.719 -3.394
ATOM 23 CA ALA 4 -5.126 -1.057 -4.325
ATOM 24 CB ALA 4 -6.538 -1.625 -4.233
ATOM 25 C ALA 4 -5.205 0.436 -4.039
ATOM 26 O ALA 4 -4.561 0.930 -3.116
ATOM 27 OXT ALA 4 -5.915 1.166 -4.728
TER
使用MindSponge來讀取該pdb檔案的方法為[*注：由於一般pdb檔案中會忽略氫原子，因此載入的時候需要使用rebuild_hydrogen將其重構成一個完整的分子]：
from sponge import Protein
system = Protein('case1.pdb', rebuild_hydrogen=True)
print ('The number of atoms in the system is: ', system.num_atoms)
print ('All the atom names in the system are: ', system.atom_name)
相應的輸出結果為：
[MindSPONGE] Adding 57 hydrogen atoms for the protein molecule in 0.007 seconds.
The number of atoms in the system is: 57
All the atom names in the system are: [['N' 'CA' 'CB' 'C' 'O' 'H1' 'H2' 'H3' 'HA' 'HB1' 'HB2' 'HB3' 'N' 'CA'
'CB' 'CG' 'CD' 'NE' 'CZ' 'NH1' 'NH2' 'C' 'O' 'H' 'HA' 'HB2' 'HB3' 'HG2'
'HG3' 'HD2' 'HD3' 'HE' 'HH11' 'HH12' 'HH21' 'HH22' 'N' 'CA' 'CB' 'C'
'O' 'H' 'HA' 'HB1' 'HB2' 'HB3' 'N' 'CA' 'CB' 'C' 'O' 'OXT' 'H' 'HA'
'HB1' 'HB2' 'HB3']]
可以看到的是，在對應的位置上，我們將氫原子補在了一個相對合適的位置。一般情況下，重構完氫原子之後，需要對系統進行一個能量極小化，否則會導致初始系統的能量過於不穩定。具體的加氫效果可以看一下這個執行的結果：
其實加氫是很難做到一步到位的，但是我們可以儘可能的將氫原子擺放在一個相對合理的位置，便於後續的能量計算和優化。
自定義分子
由於python這一程式語言的靈活性，使得我們不僅支援從檔案和模板檔案中去定義一個分子系統，還可以直接用指令碼的形式傳一系列的python列表給Molecule來構建一個分子系統。比如我們只傳原子型別和座標還有鍵連關係，就能構建一個簡單的水分子：
from sponge import Molecule
system = Molecule(atoms=['O', 'H', 'H'],
coordinate=[[0, 0, 0], [0.1, 0, 0], [-0.0333, 0.0943, 0]],
bonds=[[[0, 1], [0, 2]]])
print ('The number of atoms in the system is: ', system.num_atoms)
print ('All the atom names in the system are: ', system.atom_name)
print ('The coordinates of atoms are: \n{}'.format(system.coordinate.asnumpy()))
相應的輸出結果如下所示：
The number of atoms in the system is: 3
All the atom names in the system are: [['O' 'H' 'H']]
The coordinates of atoms are:
[[[ 0. 0. 0. ]
[ 0.1 0. 0. ]
[-0.0333 0.0943 0. ]]]
總結概要
本文通過解析MindSponge的原始碼實現，詳細介紹了在MindSponge中Molecule基礎分子類的內建屬性和內建函數，以及三種相應的分子系統定義方法：我們既可以使用yaml模板檔案來定義一個分子系統，也可以從mol2和pdb檔案格式中直接載入一個Molecule，還可以直接使用python列表的形式傳入一些手動定義的內容，直接構建一個Molecule。有了最基礎的分子系統之後，後面就可以開始定義一些能量項和迭代器，開始分子動力學模擬。
版權宣告
本文首發連結為：https://www.cnblogs.com/dechinphy/p/mol-system.html
作者ID：DechinPhy
更多原著文章請參考：https://www.cnblogs.com/dechinphy/
打賞專用連結：https://www.cnblogs.com/dechinphy/gallery/image/379634.html
騰訊雲專欄同步：https://cloud.tencent.com/developer/column/91958
CSDN同步連結：https://blog.csdn.net/baidu_37157624?spm=1008.2028.3001.5343
51CTO同步連結：https://blog.51cto.com/u_15561675