技術文章 » MindSponge分子動力學模擬——計算單點能（2023.08）

MindSponge分子動力學模擬——計算單點能（2023.08）

2023-08-31 18:03:24

技術背景

在前面的幾篇文章中，我們分別介紹了MindSponge的軟體架構、MindSponge的安裝與使用和如何在MindSponge中定義一個分子系統。那麼就像深度學習中的損失函數，或者目標函數，這裡分子力學的主要目標函數就是勢能（也有動能項，但動能項更多的來源於分子動力學模擬的過程，而不是實驗中的引數）。這篇文章我們主要探討一下，在MindSponge中如何計算給定體系的單點能。

淺析ForceField類

ForceField是用於構建分子力場的基礎類，首先我們還是來看一下ForceField這個類的內部實現，主要關注下需要哪些輸入，以及這些輸入的含義是什麼。

class ForceField(ForceFieldBase):
    r"""Potential energy of classical force field. It is a subclass of `ForceFieldBase`.

        The `ForceField` class use force field parameter files to build the potential energy.

    Args:

        system (Molecule):      Simulation system.

        parameters (Union[dict, str, List[Union[dict, str]]]):
                                Force field parameters. It can be a `dict` of force field parameters,
                                a `str` of filename of a force field file in MindSPONGE YAML format,
                                or a `list` or `tuple` containing multiple `dict` or `str`.
                                If a filename 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 MindSPONGE's built-in force field.
                                If multiple sets of parameters are given and the same atom type
                                is present in different parameters, then the atom type in the parameter
                                at the back of the array will replace the one at the front.

        cutoff (float):         Cutoff distance. Default: None

        rebuild_system (bool):  Whether to rebuild the atom types and bond connection of the system
                                based on the template in parameters.
                                Default: True

        length_unit (str):      Length unit. If None is given, it will be assigned with the global length unit.
                                Default: None

        energy_unit (str):      Energy unit. If None is given, it will be assigned with the global energy unit.
                                Default: None

    Returns:

        energy (Tensor):    Tensor of shape `(B, E)`. Data type is float.

    Supported Platforms:

        ``Ascend`` ``GPU``

    Symbols:

        B:  Batchsize, i.e. number of walkers in simulation.
        E:  Number of energy terms.

    """

在這裡面重點關注下這三個輸入資訊：

  1. 分子系統。這裡麵包含了原子型別、成鍵關係等構建力場所必須的引數，因此我們可以直接把構建好的Molecule物件直接傳入ForceField。
  2. 力場引數。我們已經在MindSponge內部轉化了一些常用的Amber力場引數檔案，包含amber.ff03，amber.ff14sb，amber.ff96和amber.gaff等，還有兩個常用的水分子力場spce和tip3p。我們在構建ForceField物件的時候，只需要用字串給parameters這一項賦值即可呼叫相關的力場檔案。如果是自己做了一個力場檔案，那麼傳入該力場檔案的絕對路徑即可。能被成功呼叫的前提條件是，完全按照MindSponge的力場檔案格式要求做的力場。
  3. 截斷半徑。在力場裡面的長程相互作用力的計算中會使用到近鄰表的計算，因為對大體系而言我們不太可能考慮全連線圖，因此需要用一個截斷半徑對長程相互作用的範圍進行一定的截斷。

除了ForceField類之外，我們還需要關注一個相關的類是MindSponge中的WithEnergyCell，這也是在計算單點能中一個非常重要的類。

class WithEnergyCell(Cell):
    r"""
    Cell that wraps the simulation system with the potential energy function.
    This Cell calculates the value of the potential energy of the system at the current coordinates and returns it.

    Args:
        system(Molecule):               Simulation system.
        potential(PotentialCell):       Potential energy function cell.
        bias(Union[Bias, List[Bias]]):  Bias potential function cell. Default: None
        cutoff(float):                  Cut-off distance for neighbour list. If None is given, it will be assigned
                                        as the cutoff value of the of potential energy.
                                        Defulat: None
        neighbour_list(NeighbourList):  Neighbour list. Default: None
        wrapper(EnergyWrapper):         Network to wrap and process potential and bias.
                                        Default: None

    Inputs:
        - **\*inputs** (Tuple(Tensor)) - Tuple of input tensors of 'WithEnergyCell'.

    Outputs:
        energy, Tensor of shape `(B, 1)`. Data type is float. Total potential energy.

    Supported Platforms:
        ``Ascend`` ``GPU``

    Symbols:
        B:  Batchsize, i.e. number of walkers of the simulation.
        A:  Number of the atoms in the simulation system.
        N:  Number of the maximum neighbouring atoms.
        U:  Number of potential energy terms.
        V:  Number of bias potential terms.
    """

在前面一部我們構建好的ForceField裡面，實際上更多的是建模過程，在分子構象空間和拓撲連線性上面查詢相應的物件，並匹配到相應的力場引數，最終作為一個能量計算環節的引數輸入。而WithEnergyCell則可以根據system和forcefield兩者提供的資訊，直接計算能量。關於WithEnergyCell類的輸入，我們重點關注下這幾個引數：

  1. 分子系統。就是前面構造好的system引數，或者Molecule物件。
  2. 力場模型。根據上一個介紹的ForceField類中構建的力場模型，有了這個力場模型，就可以直接計算各項力場的能量和力的大小。
  3. 偏置勢。這一項雖然看起來只是一個引數，但其實最能夠體現使用AI框架來做分子動力學模擬的優勢，後面會專門寫一篇文章來介紹。其主要作用是新增一些偏置勢能項，用於約束分子系統，向期望的方向演化。比如我們可以新增一個球形諧振子勢，把體系約束在一個給定體積大小的球體中。也可以通過bias新增Meta Dynamics，用於快速模擬化學反應的過程。
  4. 近鄰表。相比於給定一個cutoff，這裡給定一個近鄰表會更加直接。由於近鄰表在分子模擬的過程中有可能直接決定了分子模擬的速度，或者是可模擬體系的大小，因為近鄰表所產生的最大Tensor維度有可能為(B,A,N,D)，是所有Tensor中最大的。熟悉GPU計算的童鞋可能都知道，在GPU上佔用視訊記憶體過大有可能導致不可計算的問題。
  5. 修飾器。wrapper這個引數，為力場的使用增添了許多的靈活性。跟bias的功能是很相似的，但是bias相當於是新增了一個新的勢能函數，而wrapper則是直接對原始力場的能量進行修飾。

以上，就是ForceField和WithEnergyCell的基本介紹，有了這兩個工具，我們就可以匯入一些已有的力場，來完成單點能的計算。

計算單點能

給定一個初始的分子構象，求這個構象的分子勢能。一般情況，這個問題可能需要通過計算化學的方式來進行求解，比如CCSD(T)之類的方法。但是我們做分子模擬，需要快速的演化和迭代，如果使用計算化學的方法，速度是無法滿足計算的需求的。因此將實驗資料或者是計算化學的資料擬合成一系列的力場引數，目前來說還是一個相對比較快速而且不損失太多精度的方案。接下來我們展示一下如何使用MindSponge來計算單點能。

from sponge import Protein, ForceField, WithEnergyCell
system = Protein('case1.pdb', rebuild_hydrogen=True)
energy = ForceField(system, parameters='amber.ff14sb')
with_energy = WithEnergyCell(system, energy)
print (with_energy.energy_names)
print (with_energy.calc_energies())

在上面這個案例中，先是用MindSponge內建的Protein類載入了case1.pdb這個案例作為範例。這個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
TER

執行結束之後會在case1.pdb同路徑下生成一個補氫的檔案case1_addH.pdb，內容如下所示：

MODEL     1
ATOM      1  N   ALA A   1      -0.095 -11.436  -0.780   1.0   0.0           N
ATOM      2  CA  ALA A   1      -0.171 -10.015  -0.507   1.0   0.0           C
ATOM      3  CB  ALA A   1       1.201  -9.359  -0.628   1.0   0.0           C
ATOM      4  C   ALA A   1      -1.107  -9.319  -1.485   1.0   0.0           C
ATOM      5  O   ALA A   1      -1.682  -9.960  -2.362   1.0   0.0           O
ATOM      6  H   ALA A   1      -0.994 -11.866  -0.701   1.0   0.0           H
ATOM      7  HA  ALA A   1      -0.767  -9.928   0.292   1.0   0.0           H
ATOM      8  HB1 ALA A   1       1.816  -9.816   0.015   1.0   0.0           H
ATOM      9  HB2 ALA A   1       1.407  -8.427  -0.329   1.0   0.0           H
ATOM     10  HB3 ALA A   1       1.797  -9.446  -1.427   1.0   0.0           H
TER
ENDMDL
END

上述單點能的計算結果如下所示：

[MindSPONGE] Adding 10 hydrogen atoms for the protein molecule in 0.003 seconds.
Warrning! The head_atom of residue 0 is not None but the tail_atom of residue -1 is None.
Warrning! The head_atom of residue 0 is not None but the tail_atom of residue -1 is None.
['bond_energy', 'angle_energy', 'dihedral_energy', 'coulomb_energy', 'lj_energy', 'nb_pair_energy']
[[  45.856487     76.87538       2.9642215  -100.20307      -0.50251234
   206.13483   ]]

而這裡如果我們還是依舊使用這個程式碼，只是替換一個輸入的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] Adding 57 hydrogen atoms for the protein molecule in 0.008 seconds.
['bond_energy', 'angle_energy', 'dihedral_energy', 'improper_energy', 'coulomb_energy', 'lj_energy', 'nb_pair_energy']
[[ 231.55423  1049.3224    194.43379    20.708637 -689.10315  4051.5076
   163.09871 ]]

在MindSponge的單點能的計算過程中，我們保留了每一項能量的name和value，便於結果的比對。

需要注意的是，不同的分子模擬軟體所使用的預設輸出單位也是不一樣的。所以我們用MindSponge進行分子模擬的時候，最好在程式碼的最前頭把單位設定加上：

from sponge import set_global_units
set_global_units('nm', 'kj/mol')

比如像上面這樣操作一下，就把全域性的能量單位設定成了kj/mol，全域性的長度單位設定成了nm。常用的長度單位還有A，常用的能量單位還有kcal/mol，主要就是這裡兩組單位之間進行切換。

總結概要

本文主要銜接前面的文章，繼「MindSponge的安裝與使用」、「MindSponge軟體架構」以及「MindSponge中定義一個分子系統」系列文章之後，再講解一下如何根據一個定義好的分子系統進行力場建模，使用力場來計算單點能，就是一個比較簡單的案例。

版權宣告

本文首發連結為：https://www.cnblogs.com/dechinphy/p/single-point-energy.html

作者ID：DechinPhy

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

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