求助迭代自洽电荷加ONIOM计算重组能的流程细节

杲杲出日

我想按照J. Am. Chem. Soc., 2008, 130 (37), pp 12377–12384中的流程去考虑中心分子和环境分子的极化作用，以计算重组能，关于电荷自洽部分，我的做法如下：
1. 通过QM方法计算单分子的ESP电荷，并将其作为体系中所有分子的初始电荷。
2. 采用QEq方法对各MM区域分子的电荷逐个进行更新，其中其余分子的电荷作为背景点电荷引入计算。Q1：在这一步中，我发现g16计算QEq时无法考虑背景电荷，那按照博文sobereva.com/333中的做法：具体来说，是先让QM描述的中心分子的初始电荷设为单独计算时得到的CHELPG电荷，然后在Gaussian算QEq电荷时固定QM原子的电荷不变去确定MM原子的电荷以表现环境分子对中心分子电荷分布的响应。这该怎么去实现呢？后面步骤中低层区域的电荷我是用QM方法去做的
3. 接着利用ONIOM方法计算体系能量，并判断电荷迭代是否收敛。若未收敛，则以更新后的MM区域电荷作为背景点电荷，重新计算QM区域分子的电荷，并继续对MM区域电荷进行迭代更新。关键词为ONIOM(B3LYP/6-31+G(d):UFF)=EmbedCharge geom=connectivity
4. 若计算的是s0或cation最优结构的能量，我在第3步的电荷自洽后，固定电荷对其使用ONIOM方法进行结构优化。然后我考虑到其结构发生变化，那么其esp电荷也会发生变化，所以我比较初始构型与优化后构型的能量差：若二者差异较大，则采用新的优化构型重新进行上述电荷迭代过程；若能量差较小，则认为体系的结构与电荷均已达到收敛，最终输出体系能量。关键词为ONIOM(B3LYP/6-31+G(d):UFF)=EmbedCharge geom=connectivity nosymm
Q2：按照上述流程计算得到的S0结构下Cation态的能量Es0-cation反而比cation态最优结构的能量Ecation还要更低一些，这是不是因为我每次优化时改变了MM区域的电荷分布，相当于改变了QM的哈密顿量，导致在不同的势能面上优化？
Q3：如果说不需要上述过程的第4步，那么我是应该先优化结构，在最优结构下平衡电荷计算能量呢，还是先平衡电荷后固定电荷去优化结构呢？