layout: post title: 水溶液体系成笼分析程序cagen categories:
在对水溶液体系, 特别是气体水合物体系的模拟结果进行分析时, 常需要分析水分子借助氢键网络所形成的笼子. 甲烷水合物的常见构型sI, sII等就是以笼子形状区分的. 常见的相关模拟论文中通常会统计各种笼子数目与形态的变化. 目前我能查到的成笼分析程序有
GRADE: Mahmoudinobar F, Dias C L. GRADE: A code to determine clathrate hydrate structures. Computer Physics Communications, 2019, 244: 385-391.FSICA: Guo G-J, Zhang Y-G, Liu C-J, et al. Using the face-saturated incomplete cage analysis to quantify the cage compositions and cage linking structures of amorphous phase hydrates. Physical Chemistry Chemical Physics, 2011, 13(25): 12048.ICO: Hao Y, Xu Z, Du S, et al. Iterative Cup Overlapping: An Efficient Identification Algorithm for Cage Structures of Amorphous Phase Hydrates. The Journal of Physical Chemistry B, 2021, 125(4): 1282-1292.HTR: Liu Y, Xu K, Xu Y, et al. HTR: An ultra-high speed algorithm for cage recognition of clathrate hydrates. Nanotechnology Reviews, 2022, 11(1): 699-711. 测试版本见HTRgenice-cage: M. Matsumoto, A. Baba, I. Ohmine; Topological building blocks of hydrogen bond network in water; J. Chem. Phys. 127(13):134504, 2007; 10.1063/1.2772627. 可基于其开源的python代码进行二次开发, 但可能慢且繁琐NetworkX: 虽然没有直接的笼型识别功能, 但具有很多作为其基础的图论算法其中GRADE是开源的, 其他程序要么不开源, 要么过于久远不易获得. 可惜的是GRADE只能识别几种常规的笼形(常见sI水合物为5¹²和5¹²-6², sII为5¹²和5¹²-6⁴). 实际应用中, 仅识别5, 6元环并不够, 总有其他各种奇形怪状的多面体值得研究, 一般需要识别4, 5, 6, 8元环及其组合. 然而要将GRADE代码扩展到能够识别任意笼形并不容易.
鉴于很多做水合物模拟的人需要进行成笼分析, 我就思考了一下这个问题, 得到了一些想法, 并初步编写了一个C++程序以验证想法的可行性.
汝果欲学诗，工夫在诗外
我们先来考虑一下这个问题的定义. 处理任何问题的起点都是明确它的定义. 有些问题不易解决, 就在于没有良好定义. 良好地定义问题等于解决了问题的一半.
要进行成笼分析, 必须首先明确笼子的定义. 广义地说, 所谓笼子就是一个多面体区域, 由顶点, 边, 面等几何元素构成区域的内外分界面, 区域内部也可以存在一些东西. 这样定义的广义笼子最普遍通用, 但很难处理, 因为有太多可能性. 基于我们的目的, 可以将笼子定义得更狭义具体一些. 首先, 我们只考虑近似凸多面体, 因为非凸多面体更难处理, 数学上除了将其划分为更小的凸多面体外也没有太多处理办法. 其次, 我们不考虑复杂的笼子, 只考虑单笼, 即笼子只能相邻而不能互相嵌套穿插重叠. 一个笼子内部也不会包含其他笼子或者片段. 这样问题就简单了不少. 我们要找的笼子是凸多面体, 所以只要考虑由几个多边形封闭形成的空间即可. 为此, 我们可以首先找到各种多边形. 在此基础上, 根据边的相连接性拼合起来形成笼子.
接下来的问题是找多边形. 多边形是由几条边连接成的环, 所以要首先找到这些边, 然后根据边找到环, 形成多边形. 要找到边, 那首先需要定义顶点之间存在边的条件. 这样我们就将问题倒推到了要处理的最小元素--顶点. 解决问题的时候将这个过程反过来, 从顶点(原子)开始, 连点成边(键), 连边成面(多边形, 环), 连面成体(笼).
问题定义明确了, 下面就要看看能否将其化归为数学上已知的问题. 很多时候一些特定领域的问题都是数学上早已研究过并解决了的问题, 不需要重复发明更拙劣的轮子. 我们所考虑的这个要寻找顶点, 边, 面, 体的问题, 实际上是图论中的经典问题. 在图论中, 环也称为圈(cycle, loop), 无交叉的圈称为无弦圈(chordless cycle). 根据顶点连接关系找到圈的方法, 属于遍历问题, 采用深度优先搜索即可. 虽然也可以像GRADE那样采用多重循环的方法, 但代码写起来很麻烦, 也不好处理任意的多元环. 可以参考的资料罗列如下:
如何根据得到的多边形(环)组合得到多面体(笼子)是另一个难点. 目前文献上有两类方法. 一类是GRADE采用的cup方法或者其变体, 就是从一个环开始, 将与其有邻边的其他环拼接起来形成杯, 再由杯拼接成笼子. 这种方法容易理解, 处理起来却麻烦, 因为要根据环的连接信息和取向不断调整. 另一类方法是FSICA采用的遮掩方法, 简单说来就是将空间分格子, 得到很多格点, 以每个格点为起始点, 测试其周围是否存在笼子. 这类方法思路可取, 但做法过于简陋笨重, 找到的笼子也未必符合我们的要求. 我觉得除上述两类方法外还可以试试以下两类方法:
Qhull程序可以获得半平面的交集, 所以可以基于其结果获得笼子的构型. 我目前暂时不打算涉及Qhull过多, 所以只是在此提出这种我觉得可行的方法, 实际验证其可行性, 有待来者了.FSICA的思路, 但不是基于格点, 而是计算每个面的中心, 法向, 从每个面A的中心向其他面B的中心(使用所有顶点更可靠)发出射线, 检查射线是否被其他面遮挡. 如果被遮挡, 说明面B和面A不属于同一个笼子. 对每个面, 最多可以同时处于两个笼子之中, 所以取其法线正反两个方向检查是否与其他面相交即可. 计算遮掩属于计算机图形学或计算几何中的经典问题. 计算射线与平面的交点, 或更具体的, 计算射线与三角形的交点, 有很多高效的算法, 可以参考的代码/资料很多. 其中最常用的是MT方法, 见Moller T, Trumbore B. Fast, Minimum Storage Ray-Triangle Intersection. Journal of Graphics Tools, 1997, 2(1): 21-28, 近年来还有BW方法, 见Baldwin D, Weber M. Fast ray-triangle intersections by coordinate transformation. 2016, 5(3): 39-49等. 对各种方法的系统比较, 可以参考Fast CPU Ray-Triangle Intersection Method. 计算交点是任何渲染器, 光线追踪都要用到的一步处理. 所以即便到今天还是有很多人在研究这个问题. 我一直对计算机图形学感兴趣, 所以优先选择实现这个方法.关于射线, 三角形相交的参考:
可以想见的
其流程如下
Windows下可执行程序及其测试文件见gmxtools.
命令行选项:
-f : 输入gro文件, 目前只处理文件中的第一帧构型-wat : 水分子的残基名称-oxy : 水分子氧的原子名称, 与前一个选项一起用于确认水分子氧原子-or : 输出环构型的pdb文件-oc : 输出笼构型的pdb文件-od : 输出每种笼的数目-hb : 水分子氢键的标准(gmx), 距离(Å), 角度(°)-ring: 成环标准, 最少原子数, 最多原子数, 环中原子偏离平面的最大距离(Å)-cage: 成笼标准, 遮掩射线起点与环中心距离(Å), 起点与笼中各环中心的连线偏离法线的最大角度(°)conf1_512-64.gro: 常规笼型conf2_512-64.groconf3_512-62.groconf4_512-62.groconf5_512.groconf6_512.groconf7_MH1.gro: 甲烷水溶液conf8_MH2.gro: 甲烷水合物第1个测试文件 只含一个笼子5¹²-6⁴,
<div class="highlight"><pre style="line-height:125%"><span style="color:#A2F">cagen</span> <span style="color:#666">-f</span> conf1_512-64.gro <span style="color:#A2F">>>>>>>>>>>>>>>>>>>>>>>>>>></span> cagen <<<<<<<<<<<<<<<<<<<<<<<<<< <span style="color:#A2F">>>>>>>>>>>>>>>>>>>>>>>>>>></span> Jicun Li <<<<<<<<<<<<<<<<<<<<<<<<<< <span style="color:#A2F">>>>>>>>>>>>>>>>>>>>>>>>>>></span> 2022-09-02 12:16:48 <<<<<<<<<<<<<<<<<<<<<<<<<< <span style="color:#A2F">>></span> Usage: cagen [-f FILE.gro] [-wat WATER] [-oxy WATEROXYGEN] <span style="color:#A2F">[-or</span> RING.pdb] [-oc CAGE.pdb] [-od CAGE.dat] <span style="color:#A2F">[-hb</span> DISTACE ANGLE] <span style="color:#A2F">[-ring</span> MIN MAX DEVIATION] <span style="color:#A2F">[-cage</span> DISTACE ANGLE] <span style="color:#A2F">>></span> Default: cagen <span style="color:#666">-f</span> conf.gro <span style="color:#666">-wat</span> SOL <span style="color:#666">-oxy</span> OW <span style="color:#666">-or</span> conf~ring.pdb <span style="color:#666">-oc</span> conf~cage.pdb <span style="color:#666">-od</span> conf~cage.dat <span style="color:#666">-hb</span> 3.5 30 <span style="color:#666">-ring</span> 4 7 100 <span style="color:#666">-cage</span> 1 85 <span style="color:#A2F">-------------------------------------------------------------------------------</span> <span style="color:#A2F">>></span> Option: <span style="color:#666">-f</span> : input configuration file <span style="color:#666">-wat</span> : residue name of water <span style="color:#666">-oxy</span> : atom name of water oxygen <span style="color:#666">-or</span> : output ring file <span style="color:#666">-oc</span> : output cage file <span style="color:#666">-od</span> : output cage data <span style="color:#666">-hb</span> : criteria of hydrogen bonding(A, degree) <span style="color:#666">-ring:</span> min, max number of members, max deviation of atoms from plane(A) <span style="color:#666">-cage:</span> distance from center along normal to define ray starting point(A) <span style="color:#A2F">max</span> angle between normal and ray in the same cage(degree) <span style="color:#A2F">-------------------------------------------------------------------------------</span> <span style="color:#A2F">>></span> Warning: <span style="color:#A2F">!!!</span> 1. Eigen is needed to compile the code <span style="color:#A2F">!!!</span> 2. Use JMOL to view output pdbs <span style="color:#A2F">-------------------------------------------------------------------------------</span> <span style="color:#A2F">>>>></span> Settings <span style="color:#A2F">conf</span> File: conf1_512-64.gro <span style="color:#A2F">ring</span> File: conf~ring.pdb <span style="color:#A2F">cage</span> File: conf~cage.pdb <span style="color:#A2F">data</span> File: conf~cage.dat <span style="color:#A2F">water</span> Name: SOL <span style="color:#A2F">OW</span> Name: OW <span style="color:#A2F">hb</span> Criteria: 3.5(A) 30(degree) <span style="color:#A2F">ring</span> Criteria: 4 7 100(A) <span style="color:#A2F">cage</span> Criteria: 1(A) 85(degree) <span style="color:#A2F">>>>></span> System <span style="color:#A2F">numAtom:</span> 112 <span style="color:#A2F">degree:</span> size= 112 capacity= 112 <span style="color:#A2F">resName:</span> size= 112 capacity= 112 <span style="color:#A2F">atmName:</span> size= 112 capacity= 112 <span style="color:#A2F">atmXYZ:</span> size= 112 3 capacity= 112 3 <span style="color:#A2F">isOW:</span> size= 112 capacity= 128 <span style="color:#A2F">isEdge:</span> size= 112 112 capacity= 112 128 <span style="color:#A2F">>>>></span> Water <span style="color:#A2F">OW:</span> 28 <span style="color:#A2F">numNode:</span> size= 28 capacity= 28 <span style="color:#A2F">distRing:</span> size= 28 capacity= 28 <span style="color:#A2F">centRing:</span> size= 28 3 capacity= 28 3 <span style="color:#A2F">normRing:</span> size= 28 3 capacity= 28 3 <span style="color:#A2F">tangRing:</span> size= 28 7 12 capacity= 28 7 12 <span style="color:#A2F">nearRing:</span> size= 28 28 capacity= 28 28 <span style="color:#A2F">edgeList:</span> size= 112 capacity= 112 <span style="color:#A2F">ringList:</span> size= 28 28 capacity= 28 28 <span style="color:#A2F">cageList:</span> size= 56 28 capacity= 56 28 <span style="color:#A2F">cageNorm:</span> size= 56 capacity= 56 <span style="color:#A2F">>>>></span> Node degree: <span style="color:#A2F">>>>></span> Edge Size: 84 <span style="color:#A2F">>>>></span> Edge Size (degree>1): 84 <span style="color:#A2F">>>>></span> Find all Rings Iter: 1601 Time(s): 0.002 <span style="color:#A2F">Ring:</span> 16 <span style="color:#A2F">>>>></span> Ring Info. <span style="color:#A2F">numNode:</span> size= 16 capacity= 28 <span style="color:#A2F">distRing:</span> size= 16 capacity= 28 <span style="color:#A2F">centRing:</span> size= 16 3 capacity= 28 3 <span style="color:#A2F">normRing:</span> size= 16 3 capacity= 28 3 <span style="color:#A2F">tangRing:</span> size= 16 7 12 capacity= 28 7 12 <span style="color:#A2F">nearRing:</span> size= 16 16 capacity= 28 16 <span style="color:#A2F">ringList:</span> size= 16 6 capacity= 28 6 <span style="color:#A2F">cageList:</span> size= 32 16 capacity= 56 16 <span style="color:#A2F">>>>></span> Fit Ring to Plane Time(s): 0.002 <span style="color:#A2F">Ring:</span> 16 <span style="color:#A2F">>>>></span> Neighbor all Rings Time(s): 0.001 <span style="color:#A2F">>>>></span> Write Rings to File <span style="color:#A2F">>>>></span> Determin Ring Occulasion Time(s): 0 <span style="color:#A2F">>>>></span> Conect Rings to Cage Time(s): 0 <span style="color:#A2F">Cages:</span> 32 <span style="color:#A2F">>>>></span> Deduplicate Cages Time(s): 0 <span style="color:#A2F">Cages:</span> 1 <span style="color:#A2F">>>>></span> Write Cages to File</pre></div>识别笼型时先根据氢键判据确定两个水分子之间是否存在氢键, 再由此确定出所有可能的环(默认为4元环到7元环), 并拟合每个环平面的方程, 计算其法向, 以及环中原子偏离环平面的最大距离, 如果此最大距离小于设定值(默认100 Å), 则认为此环存在.
环构型输出文件中, 紫红原子对应环的中心, 白色原子对应环的法向, 由环中心出发的法线长度为1 Å(默认值).
以每个环的法线点为起始点, 向其他环中心连线, 如果此连线与法线的夹角小于85°(默认值), 则认为其他环可能处于同一笼中, 因此进行遮掩计算. 如果连线不会与任何一个环三角形相交, 且其他环与起始环可相连, 则它们处于同一笼中.
第7个测试文件 为甲烷与水的混合溶液, 基本相当于甲烷的水溶液, 使用默认条件识别的笼子很多, 且有些奇怪的笼子. 这说明单纯使用欧拉多面体公式作为封闭笼子的的判据有时并不够, 因为有些非凸的笼子也是符合这个公式的. 初步的处理方法是调整成环判据的偏离值, 剔除那些偏离平面过大的环.
第8个测试文件 为甲烷水合物, 使用默认条件识别的笼子并不多, 大部分都是比较常规的笼子.
