图论角度简化交联建模问题

Graphite

最近接到个需求，简单来说建模时情况如下：
1、有三接头的交联剂X和二接头的线性分子Y，随机交联，任意指定交联度。X只会和Y交联，不存在X-X和Y-Y。
2、这个交联反应基团转化比较复杂，不是简单加成或者脱去一份子，但是X-Y交联后linker结构（比如酰胺基团之类）是固定的。

“传统”做法是用fix bond/create或者交联脚本抽卡，但是那对于加成反应或者简单消去一分子的还好，我这个体系估计要rubbish in rubbish out，也没法控制交联度。

正好最近在搞reax_tools，学了一点图论，发现就是个二分图问题。
以4X+6Y为例。

图一是只显示X的情况，图二三加了linker，图二三同构。从图三实际上我们已经发现这就是个有边数约束的列表匹配，很容易用代码写出来。如果是没有完全反应的，其实就是在X上面加一个封端Xt，Y同理。

所以其实要写出完整的拓扑文件，只需要以下信息：

• 每种组件（X核，Y核，LINKER，X封端，Y封端）的拓扑，对应图中的点
• 因为连接而增加的拓扑，对应图中的边
• 一个相当于图三的表格
  

这样我们对于lammps很容易写出基于XY数量、反应概率（边的生成概率）来生成moltemplate文件的脚本，对于gromacs就是整合top/itp文件的脚本。（moltemplate其实就相当于是lammps版的top/itp处理器）

对于初始坐标，用packmol就可以，毕竟原子数肯定是保守的，当然，位置必然是紊乱的，这个时候就要祭出pair_style soft揉搓大法了。（http://bbs.keinsci.com/thread-42434-1-1.html）

一个含有20个X、30个Y、60个linker的体系

揉搓前、揉搓后、结构展开之后
（在处理这种网络的展开时，只能在边缘断掉，比如地球仪->地图，只是显示问题）

然后我们就可以正常做NPT、NVT、力学性能等模拟了。

对于不完全反应的情况，也就是需要多定义封端和封端相关的拓扑而已。这个是输入5X+9Y设定95%反应概率的情况。解不唯一，不过不满意再抽也就2秒钟。


当然了，这点图论在数学和计算机科学中大概就是小学生水平...
题外话，这个问题其实和岗位匹配问题等价：X个程序员，每个程序员能干3个活，Y个程序，每个程序需要2个程序员，程序员工作量不饱和的，去做一份外包（X封端），程序人手不足的，喊一个外包（Y封端）。