注:该Python代码可以实现所有高斯(Gauss)计算输入文件的处理(特别是针对超分子体系的路径建模非常有用,也可以作为处理其它体系的参考)。具体来说,实现客体分子围绕主体分子的任意平面旋转任意角度或者实现客体分子从主体分子的任意平面的垂直方向穿过,并批量产生该路径上的高斯计算的输入文件。本实例以客体分子穿过主体分子为例。
。
正文如下:
超分子建模过程中,如果需要模拟客体分子穿过主体分子孔洞的过程(如小分子穿过大环分子)的能量变化,或者相互作用力的变化,具体实例如图一所示,球状富勒烯C60穿过一个分子环的过程。
图一
假设大环是由10个苯环通过σ键连接而成,形成[10]CPP,中文名称10环对苯撑,富勒烯C60分子作为客体分子,与大环会存在相互作用力,显而易见,在尺寸上这两者非常匹配,类似于在富勒烯球外面套上一根腰带,而这类分子是典型的超分子结构,理论计算对于研究这类体系尤为重要,特别是研究主体分子和客体分子间的相互作用力,扯远了,主客体相互作用不是本帖子讨论的重点。。。。。
回归正题,实现计算客体分子穿过主体分子这一过程的能量变化,设计这个模型的逻辑如下:
1. 总体思路是,富勒烯从一端到另一端的路径上设置多个点,可密可稀,对路径上的每一个点对应的超分子结构进行平面内的优化,求单点能做成二维图即可(也可以是刚性,默认中心的每一个点就是在该平面上的能量最低点,这样不用再优化结构,大大降低计算耗时)。
2. 这个路径是在[10]CPP大环的环平面对应垂直方向上,也就是说需要确定一个[10]CPP的环平面,可以默认为苯环之间的碳碳σ键均处在大平面上,也就是说,只要在这个位置找三个点就可以确定大环的化学平面(如C13-C4-C43三个碳原子对应的坐标点),再求三个坐标形成平面的法向量,再将得到的法向量转化为单位向量,以方便后续对坐标进行操作。
3. 在给定坐标基础上,只需要将富勒烯分子所有的坐标值在这个法向量上进行数学加减,即可代表在这个方向上的移动,移动的位置通过实际距离长度乘以单位法向量,即可得实际移动距离的具体坐标。
4. 通过距离批量命名输入文件名,得到对应高斯计算的输入文件。
5. 通过Gauss软件(主流是g09或者g16版本)计算得到每一个点能量,提取能量作图(非本文讨论内容)。
具体实现如下:
球形富勒烯C60穿过大环的过程的模拟,第一步是先得到优化后的超分子结构,如图二所示:
图二
用Gaussview打开该分子,存储为高斯输入文件(后缀为.gjf文件),这里需要注意,坐标顺序必须保持主体分子在前面,客体分子在后面,如图二中的数字编号,[10]CPP在前,C60在后。具体坐标如下所示:
%nprocshared=12
%mem=2500mb
%rwf=1,2000mb,2,2000mb,3,2000mb,4,2000mb,5,2000mb,6,2000mb,7,2000mb,8,2000mb,9,-1
%chk=xxxx.chk
# b3lyp/3-21g*
c60-cycloparaphenylene
0 1
C 1.73732300 -6.65945200 0.02809900
C 1.02755200 -7.09867900 -1.09886500
C -0.35923100 -7.15166600 -1.09301500
C -1.09214200 -6.76784500 0.03963400
C -0.37754700 -6.49784200 1.21181700
C 1.01002500 -6.44545800 1.20553100
C -2.51785300 -6.37048700 -0.04431700
C -2.97016100 -5.74036500 -1.21000000
C -4.13845500 -4.99392400 -1.21062100
C -4.90312800 -4.84814500 -0.04645400
C -4.53591800 -5.60946900 1.07118000
C -3.36411000 -6.35602500 1.07264900
C -5.85139600 -3.71116200 0.03786900
C -5.90440600 -2.95677700 1.21537900
C -6.40004400 -1.65956700 1.21155600
C -6.85703100 -1.06756500 0.02933800
C -6.97558800 -1.88698300 -1.10270300
C -6.48354100 -3.18455800 -1.09801100
C -6.90248000 0.41140400 -0.06167900
C -6.40817300 1.02325000 -1.21997500
C -6.03985300 2.36013100 -1.21989300
C -6.14496300 3.13798100 -0.06066500
C -6.78453300 2.56780500 1.04839600
C -7.15754000 1.22949300 1.04716500
C