让Gaussian的优化过程看起来更流畅

1、让Gaussian的优化过程看起来更流畅（2007-12-30 18:13:45）转载Q标签：分类：计算化学计算化学教育在许多时候，我们很想浏览或演示一下Gaussian程序对几何构型优化的过程， 例如，起始构型是怎样一步步地变化成最终平衡构型的，反应的IRC是怎么走 的。一些软件（例如大家最常用的GaussView）能够读取Gaussian优化输出文件中各个优化步骤的分子几何构型，然后在界面上显示出其3D模型来。在 GaussView中，你可以点击鼠标一步步地观察整个优化的过程。如下图所示：然而，往往我们会发现，我们所看到的Gaussian优化过程并不那么“流畅”， 例如上图动画演示中，出现

2、了好几次“剧烈的”构型翻转。其实各个优化步骤的 构型之间并不存在太大的构型差异，出现上述“剧烈的”构型翻转的原因是这些 构型的坐标轴取向出现了很大变化。原来，Gaussian在优化分子构型时，所选 取的坐标轴的取向往往并非保持不变（不过，对于z-matrix即内坐标优化，情 况是否也如此我还没考察过，大家有兴趣可以留意一下）。一般来说，Gaussian 会选取分子的三个主轴作为直角坐标的三个轴，不过，有时也不总是这样。总之， 尽管相邻优化步骤的两个构型很相近，但如果它们的坐标轴的取向差别很大的 话，那么它们的直角坐标数值也会相差悬殊，于是在软件中所显示出来的3D模 型就会呈现“剧烈的”构型翻转

3、，从而使我们无法“流畅”地观看优化全过程。下面我就来重点介绍一下，如何校正Gaussian构型优化的输出文件，使我们可 以“流畅”地观看优化全过程。总体思路很简单，那就是将所有优化步骤的坐标轴取向都取成一致即可。具体可 以这样实施：以优化的第N步分子构型作为参照，考察第N+1步分子构型。如果 我们不断旋转第N+1步分子，总可以找到一个最佳的分子取向使得旋转后的第 N+1步分子坐标与第N步分子坐标最为接近。下面我们就来解析这个问题。设分子的第N个几何构型的各个原子的位置坐标矩阵为X = : x . x ，第 N+1个几何构型的各个原子的位置坐标矩阵为X = x . / 。现在，我 们假定X和X的

4、几何中心都被事先平移在原点。我们需要寻找一个适当的旋转 操作入，将X变换为AX 2 X，即我们需要求解下面不相容方程组的最小二乘 解：X = AX根据最小二乘法，上述方程的求解步骤如下:X Xt = AX, X, T因此，A = X Xt(X Xt) -1然而接下来的问题是：上述满足最小二乘解的变换矩阵A并不一定是正交规范 阵。而对称操作矩阵一定是正交规范阵(证明见附录1)。因此，我们需要寻找 一个最接近于最小二乘解A的正交规范阵B。求这个最优的正交规范阵B的问题等价于这样的矢量极值问题：已知三维空间中的三个彼此非正交的单位矢量x, y, z，求三个彼此正交的单位矢 量x,y,z，使得这组新矢

5、量最接近于旧矢量，即要求满足：f (x,y,z)= |x - x|2 + |y - y|2 + |z - z|2最小，同时束缚条件是：X |2 = |y |2 = |z |2 = 1x y = y z = z x = 0(注意：可以将A的三个列向量先归一化，这样A = xyz，B = x y z)一个自然的思路是利用Lagrange条件极值法求解，但得到的方程组很复杂，直 接寻求解析解比较困难。因此，我们不妨寻求数值解法。我们首先考虑二维的情况，即：已知二维空间中的两个彼此非正交的单位矢量x,y，求两个彼此正交的单位矢 量x,y，使得这组新矢量最接近于旧矢量。0 0对于这个二维问题我们很容易利

6、用图解法求解，见下图：如上图所示，最接近于旧矢量x, y的一组正交矢量x, y是将矢量x, y分别旋 转q和-q角度而得到的矢量，其中q = p/4 - a/2，而a则是矢量x, y的夹角。现在，我们基于二维的结果采用迭代法来解决三维问题。首先我们由旧矢量x, y求得二维最优新矢量x,y，然后取新矢量z = x Xy； 接着由新矢量y, z求得二维最优新矢量y, z，再取新矢量x=yXz； 最后由新矢量z, x求得二维最优新矢量z, x，再取新矢量y= z Xx。至此，我们从旧矢量x, y, z得到一组更接近于正交化的一组新矢量 x, y, z,从而完成了一次迭代。将旧矢量x, y, z更新为

7、新矢量x, y, z去做下一次迭代，这样反复迭代直到得到一组正交性相当好的矢量 x*, y*, z*。这就是我们所求的最优正交规范矢量。上述算法很容易通过编程实现，从而得到最优的旋转操作矩阵B。下面是我按照上述算法将Gaussian构型优化的输出文件进行了重新取向，得到 的优化过程演示图如下：是不是流畅了很多？附录1：点对称操作矩阵一定是正交规范阵的证明。设分子的各个原子的位置坐标矩阵为X = : x . x (不妨事先将坐标原点 取在分子的中心位置)，其中3X1列向量xi是每个原子的位置坐标矢量。我们注意看乘积x t x = wij :,得到的矩阵具有如下特点：对角元Wii = r：，是位置坐标矢量长度的平方，也就是各个原子到分子中心距离的平方。1非对角元W = rr，是两个不同位置坐标矢量的点积，它只与这两位置坐 标矢量的夹角有关；由于点对称操作不改变各个原子之间的相对距离和夹角，因此点对称操作不改变 矩阵X T X。设点对称操作的变换矩阵为A (3X3矩阵)，则：(AX) t (