量子化学与分子模拟之三：势能面、构型优化与构象搜索

势能面PES-Potential energy surface

    根据Born-Oppenheimer近似，分子基电子态的能量可以看做只是核坐标的函数；分子力学中的所有定义的函数均只是核坐标的函数，体系能量的变化可以看成是在一个多维面上的运动。

3.1势能面及独立坐标数

   （1）势能面和独立坐标数

    对于N个原子的体系，势能面是一个超平面，由势能对体系的坐标自由度（独立坐标数）构成。体系的坐标自由度为3N-6（唯一特殊的是双原子分子，其坐标自由度为3N-5）

    对于任何原子数为N的分子：

        其直角坐标数为3N；#每个原子可沿xyz三个方向移动，从而改变能量，势能面

        平动坐标数为3；  #即整体向x，y，或z方向平移，该平移与N无关且不影响原子相对位置，不改变能量

        转动坐标数为3    #即整体绕x轴，y轴或者z轴旋转，与平动一样不改变能量

        独立坐标数为3N-3-3，即3N-6    #独立坐标数是真正改变势能面能量的自由度，双原子分子的独立坐标数是1，因为其转动坐标数少了一个，所以是3N-3-2=3*2-3-2=1。

    PS：关于平动和转动

    平动：在运动过程中，物体上任意一条直线在各时刻的位置都保持平行且相等

    转动：物体绕同一条直线做圆周运动。

（平动和转动示意图）

    （2）体系的内坐标数其实就是体系的独立坐标数（3N-6）

    对于两个原子体系，独立坐标数是1，也就是势能面只一个键长的变化相关；对于三个原子体系，独立坐标数是3，势能面和两个键长和一个键角相关；对于四个原子体系，独立坐标数为6，势能面和三个键长，两个键角还有一个二面角相关；以此类推。。。

    另外，三原子体系即使不用键角去表示，也是有三个自由度，需要三个键长；四原子体系即使不用二面角去表示，也是有六个自由度，需要四个键长，两个键角，如上图所示。

    （3）势能面的描绘

    两个原子：自由度为1，即只有一个键长对应能量的变化，所以可以用二维平面来绘图。

    三个原子：自由度为3N-6=3，即两个键长+1个键角，3个自由度对应能量，难以绘图，已经是个四维的超平面。

    ……

    不过，虽然3个自由度对应1个能量不好描绘，但是我们可以描述某1个自由度或者某2个自由度与能量之间的关系，比如只改变键角，来考察能量的变化，也就是势能面扫描（http://sobereva.com/473）。  改变一个自由度，对应能量做出来二维图；改变两个自由度对应能量做出来三维图。

3.2势能面上的点

    势能面上最令我们感兴趣的点是势能对坐标一阶微商为0的点。势能对坐标的一阶微商对应着力，因此处于势能面这样的点所受到的力为0，这样的点称为不动点（stationary point）。——极小点或者鞍点

（两个自由度对应的势能面示意图）

     PS：微商就是导数，一阶导数为0，证明该点斜率为0，是一个平面的点；二阶微商>0，证明是个凹函数，该点是个极小点，二阶微商<0，证明是个凸函数，该点是个极大点。

1，为什么关注极小点？

    极小点是出现频率最高的结构，是真实体系的代表

    数学意义上极小点是势能对坐标一阶导数为0，二阶导数为正（Hessian矩阵本征值为正），因此可用数学方法搜索（如优化等）。需要注意的是，一般优化方法仅可以找到初始构型附近的极小点，因此优化的初始构型非常重要。

    如果偏离极小点位置则受到相反方向的力，因此可以计算出振动频率，振动频率对应分子振动光谱（红外，拉曼）

    一般优化过程为了节省时间，其Hessian矩阵本征值，采用的是估算，因此要严格确定所优化的结果是否是真正的极小点需要作频率分析，所计算出的频率应均为正（受力与运动方向相反）。如果出现负值，可能由对称性限制引起。

2，鞍点

    势能面上的另一类重要的不动点为鞍点（严格意义讲应该是一阶鞍点），这些鞍点是连接两个极小点中间最低的“山口”，对应于化学反应的过渡态（或构型变化中的中间态）

    从数学意义上，在鞍点处势能对坐标的一阶导数为0，而Hessian矩阵本征值只有一个负值（如果有多个负值就是高阶鞍点，不是我们关注的一阶鞍点）。鞍点是在一个方向上具有极大，而其它方向均为极小。严格的鞍点需要进行频率分析验证，必须有且只有一个虚频率（频率为负，受力与运动方向相同）

3，最小能量途径（minimum energy path）

    最小能量途径（MEP）是连接势能面上两个极小点之间最低的能量途径，MEP也称内禀反映坐标（intrinsic reaction coordinate - IRC）。

AG构象是全局最小点，C，E为局域最小点，BF为过渡态（鞍点）

PS:对体系能量的几点说明

（1）能量的绝对值

    分子力学是以标准的平衡位置（处于平衡键长键角等）做为零点，能量有正有负；从头算能量的零点是所有核和电子相距无穷远，因此所计算出的体系能量都是负值；一般来讲能量的绝对值是没有讨论价值的，比较能量的大小需要使用在同一个水平下计算的结果

（2）能量的比较

    对于不同的体系，更准确地说，对于含有不同原子数的体系，能量的绝对值的比较是毫无意义的。分子模拟中比较的能量值必须是同一体系，在变化前后不能有原子个数，种类的变化

   

    硝基苯的能量值低于苯分子的能量值，但是不能说硝基苯比苯稳定：因为硝基苯比苯分子多了一个硝基少了一个氢原子（两个分子根本不在一个势能面上），按照量子化学能量的计算方法（每增加一个原子，其电子肯定不在距离核无穷远的位置，所以能量肯定是减少的），其能量总值肯定是要小于苯分子的，三硝基甲苯的能量值肯定也是要小于硝基苯的，这是肯定的，不用比。

    苯分子可以和第三个分子进行比较，因为他们是同分异构体，在同一个势能面上。

  PS：如果甲烷变成了甲基自由基和氢原子，其能量变化是：

3.3 能量的优化-01

    3.3.1 能量极小化、算法

    求极小点，我们希望能量对坐标的一阶微商等于0，二阶微商大于0（很难实现完美值）

 

    微分极小化法最常用，能量一阶微分的指向（梯度）表示了极小点的位置所在的方向，二阶微分表示函数的曲率，两者可以用来预测函数变化的方向（？？不太明白---明白了）。微分极小化法包括最速下降法，共轭梯度法，Newton-Raphson法和拟Newton法等。

    3.3.2 收敛的判据

    收敛就是寻找谷底那个点，在Gaussian中，使用四个判据来确定优化的终止：

 最大的力（势能对坐标的偏微分，满足条件说明这个地方是一个斜率接近0的点），力的均方根（这个点所有方向的力的均方根，满足条件的话说明这个地方已经整体比较平了）

最大的位移；位移的均方根（如果位移比较大的话，表示这是一个平坡，并不是一个谷底）

    我们后期验证：可以做频率分析

   3.3.3 初始构型的问题

    不同的初始构型得到的优化结果不同，因为存在局部最小点

    3.3.4 构象搜索

    （1）网格法

    特点：只要网格取的够细，就可以找到能量最低点。缺点：分的很细的话，计算量巨大

     该方法有一个变种，比如：HyperChem使用随机搜索的方式，即将给定的二面角随机产生初始构型，理论上只要搜索足够的数目，可以得到全局极小

    （2）蒙特卡洛（MC）搜索

    特点：随机产生概率。。。不太清楚，总之比网格搜索更容易找到低能构象

    （3）分子动力学模拟搜索-淬火

   给与温度，使原子具有动能，加上动能后，其能量将会超过其势能面的势垒。通过MD后产生多个构象，对这些构象直接进行结构优化（这时候是没有动能的，所以叫淬火），将会优化到对应的坑内。遍历的结构足够多，就可以找到能量最小点。   比蒙特卡洛效率更高

    （4）模拟退火

分为MD模拟退火和MC模拟退火：即将体系升高到高温，逐渐降温，最后体系会落到一个极小点。

完结撒花✿✿ヽ(°▽°)ノ✿