2005
Computational Chemistry laboratory
Beijing Normal university
计算化学理论和应用
－第九讲分子几何结构优化 (3)
优化收敛问题 (2)
问题 引起原因 解决方案力过大 (a)初始构型太差
(b)分子坐标选择得不好
(a)改变分子构型重新输入
(b)重新构造分子坐标优化结构的 Hessian
矩阵出现负本征值
(a)结构没有优化到极小点
(b)Hessian矩阵更新中数值计算错误
(a)沿负本征值方向优化得到能量更低的结构
(b)改变方法初猜 Hessian矩阵重新优化
(c)如果 (b)仍不行，冻结负本征值所在的坐标矢量优化收敛后再放开优化
Hessian矩阵的本征值太大 (少见 )
(a)Hessian矩阵初猜不好
(b)Hessian矩阵更新不好
(c)不同坐标之间偶合过强
(a)(b) 重新初猜 Hessian矩阵
(c) 重新构造分子坐标优化步数太多 (a)优化比较困难，需要的步数确实很多
(b)输入内坐标中存在冗余
(c)存在很松弛的坐标
(d)不同坐标之间偶合过强
(e)初猜不好
(a)增加优化次数 maxcyc=n
(b)重新输入内坐标
(c)冻结松弛的坐标优化收敛后再放开优化
(a)重新构造分子坐标优化中的对称性问题优化到高对称性一般可以做到，但是收敛非常慢，而且受算法的限制，
只能得到近似的结果
*解决办法：提高到相应对称性直接优化然后作比较优化到低对称性优化中计算得到的梯度矢量属于分子点群的对称表示，
优化过程中分子对称性不会降低
*解决办法：检查 Hessian矩阵，依负本征值所在矢量方向降低对称性过渡状态理论与反应速率计算过渡状态理论简介
R e a cti o n co o r d i n a te s
A A B
K
a c t
K
d e a c t
K
E
传统过渡态理论的反应速率公式
,0,0( ) /( ) / ABAB U U k TG G k T
A
AdB
k A k A k k
d t A
Q
k k e k e
Q
势能面的特性
1,二次区域 (qudratic region)
2,过渡矢量 (transition vector)
二次区域例,函数在点 (0,0)处的梯度为零，
Hessian矩阵为点 (0,0)是一个鞍点,Hessian矩阵的本征值和本征矢分别是
4 (0,1); 和 -4 (1,0)
22
2
401 2 8 4 1 6
041 6 8 4
x y x y
x y x
4 2 2 2 2(,) 4 2 2f x y x x y x y
-1 - 0,8 - 0,6 - 0,4 - 0,2 0 0,2 0,4 0,6 0,8 1
0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
1
-2 - 1,5 -1 - 0,5 0 0,5 1 1,5 2
-1
0
1
2
3
4
5
6
7
8
过渡结构的优化
1,一般问题沿过渡矢量方向优化到极大，而其他方向优化到极小
2,考虑对称性后的简化
H
C
H
H
XX
过渡结构的优化方法
势能面扫描
梯度模法
Gaussian中的 Berny优化
线性或二次准同步过渡法
坐标驱动法 (本征矢跟踪法 )
线性和二次准同步过渡法
R,P,反应物和产物 TS：过渡结构点 1：线性准同步过渡法寻找到的过渡结构点 3：二次准同步过渡法寻找到的过渡结构
2
( ) /
[ ( ) / ( ) /
/[ 2 ( ) ( ) ]
22
2 2 2 2
||
| | | |
| | | | | | + | |
T P R P R
T a P X P X R X R X
a R X P X R X P X R X P X
坐标驱动法和本征矢跟踪法存在的问题梯度模法
0Tgg?
梯度模法存在的问题梯度模极小值点不一定是驻点，有可能是势能函数的拐点优化过渡结构所使用关键词和选项
*J,Comp,Chem,Vol,3,No,2,214-218(1982)
Berny优化默认使用 RFO法控制优化步长，此时可以用
iop(1/8)=n (0.01*n Bohr)来控制 TSR
如果在使用 Berny优化时选用 Newton法，此时需通过
NrScale来控制步长算法 关键词 选项
Berny优化 * Opt=TS NoEigentest,iop(1/8=6)
STQN Opt=QST2,QST3
本征矢跟踪 Opt=EF CalcAll,CalcFc,EnOnly
不同方法的比较
Re ac t i o n Z - m at r i x i n t e r n al s Re d u n d an t i n t e r n al s
r e g u l ar Ca l cF C Q S T 3 Ca l cF C Q S T 2 Q S T 3
CH F CH HF
4 3
6 4 6 5 8 5
CH O CH OH
3 2
12 9 9 8 8 9
S i H H S i H
2 2 4
11 7 11 7 8 8
C H F C H HF
2 5 2 4
16 12 15 13 17 11
D i e l s- A l d e r r e ac t i o n 56 11 23 8 13 14
Cl ai se n r e ac t i o n 38 8 15 7 15 15
En e r e ac t i o n f ai l 15 28 13 18 18
优化过渡态时可能遇到的问题及解决问题 引起原因 解决方案优化过程中 Hessian
矩阵负本征值过多
(a)优化得到二级鞍点结构
(b)Hessian矩阵有数值错误
(a)沿非反应坐标方向寻找极小点
(b)改变方法初猜 Hessian矩阵优化过程中没有负本征值
(a)结构离鞍点的二次区域太远
(b)Hessian矩阵有数值错误
(a)使用更好的初始结构
(b)改变方法初猜 Hessian矩阵重新优化主反应内坐标变量后加 3计算相应 Hessian矩阵元例,0 1C
N 1 r1
H1 r2 2 a1
r1=1.2
r2=1.4 3
a1=59.0 3
D,K,Malick,G,A,Petersson,and J,A,Montgomery Jr.,J,Chem,Phys,108,5704 (1998).
IRCMax (MP2/6-31G(d):HF/3-21G*,Zero,Stepsize=12)
P,Y,Ayala and H,B,Schlegel,J,Chem,Phys,107,375 (1997).
通过反应途径寻找过渡结构
Opt(QST2,Path=M)
使用两种方法分别优化反应途径并计算反应途径上的能量最高点输入：过渡结构同时优化反应途径上的 M个点和过渡结构输入：同 QST2或 QST3
过渡结构算例
#T UHF/6-31G(d) Opt=QST2
H3CO --> H2COH Reactants
0,2
C
O 1 1.48
H 1 R 2 A
H 1 1.08 2 110,3 120.
H 1 1.08 2 110,3 -120.
R=1.08
A=110.
H3CO --> H2COH Reactants
0,2
C
O 1 1.48
H 1 R 2 A
H 1 1.08 2 110,3 120.
H 1 1.08 2 110,3 -120.
R=1.9
A=30.
E
-1 1 4,4 1 4 6 0
-1 1 4,3 3 0 5 9
-1 1 4,3 8 8 5 1
过渡结构的确认
计算完整的 Hessian矩阵，确保有且仅有一个负本征值通过振动分析计算完成，可以使用连续作业
# B3LYP /6-31g opt freq
确认本征矢量的方向正确连接到反应物和产物通过计算反应途径进行确认
IRC(内禀反应坐标 )原理和计算
IRC的定义和计算方法反应途径：
势能面上，经过过渡结构，连接反应物和产物的结构变化途径
特性：
1,满足 最小作用 原理
2,原子核的运动速度无限小
推论：
1,反应途径是势能面上连接反应物,过渡结构和产物的最小能量途径
2,反应途径上每点的梯度方向与运动方向一致
寻找反应途径的方式：
从过渡态出发寻找连接反应物和产物的最陡下降途径最陡下降途径的定义式：
运动方向,梯度：
反应途径：
选择不同的坐标系，可以有不同的最陡下降途径
i i i
ii i
jik
i j k
V
x d x x
x
dxd x d x
c
V V V
x x x
Fukui的 IRC方程
jik
i j k
dXd X d X
c
V V V
X X X
质权坐标系下的最陡下降途径，就是最小能量途径
K.Fukui,Acc,Chem,Res,1981,Vol,14,12,363
1,确认过渡结构
2,寻找反应机理
IRC计算的意义
Gaussian程序中 IRC的算法
C,Gonzalez and H,B,Schlegel,J,Chem,Phys,90,2154 (1989).
C,Gonzalez and H,B,Schlegel,J,Phys,Chem,94,5523 (1990).
IRC计算的步骤和选项
1,优化出过渡结构
2,对优化结果进行振动分析计算确认过渡结构，得到零点能，获得力常数矩阵
3,进行 IRC计算力常数选项：
从振动分析计算的 chk文件中读取 IRC=(RcFc)
使用当前方法计算 IRC=(CalcFc)
反应坐标数的选取 IRC(MaxPoints=n) 默认,n=6
点数不够时可以重起作业：
IRC=(Restart,MaxPoints=n)
步骤：
选项：
IRC算例
%chk=freqhessian.chk
#p UHF/6-31G(d) IRC=(maxpoints=10,rcfc,forward)
H3CO --> H2COH IRC
0 2
6 0.03104 0.63055 0,
8 0.03104 -0.73696 0,
1 -0.99138 -0.13555 0,
1 0.27839 1.12398 0.92614
1 0.27839 1.12398 -0.92614
振动分析计算 (2)
Hessian矩阵本征值与振动模式的关系
1,Hessian矩阵的本征值对应于简正振动的力常数
2，简正振动的频率等于力常数的平方根
3，简正振动的振动模式是相应本征值的本征矢热力学函数的统计表达式
*
*
/
l n l n ( ! )
!
i
i
i
n
i
i
kT
ii
g
S k k N
n
N
n g e
f
2
,
,
ln
ln
ln
B
NV
BB
NV
f
E Nk T
T
H E P V
f
S Nk T f Nk T
T
G H TS
配分函数和热力学函数各个分量的计算平动配分函数，平动能和平动熵转动配分函数，转动能和转动熵
e t rf f f f f
析因子性质：
32
2
2
32
2
2
8 3
()
2
8 3
l n ( )
2
a b c B
rr
a b c B
r
I I I kT
f E RT
h
I I I kT
SR
h
3 / 2 3 / 2
22
2 3 2 3( ) l n ( )
22t t t
m k T m k Tf V E R T S R V
hh
振动配分函数，振动能和振动熵热容的计算
V
V
UC
T
36
1 / 2
/
1
36
/
1
36
/
/
1
1
( 1 )
1
( 1 )
l n ( 1 )
( 1 )
i iB
iB
iB
iB
n
i i i i v h k T
i
n
i
v h k T
i B
n
h k Ti
v h k T
i B
f
e
h
ER
ke
h
S R e
k T e
振动分析算例
%chk=freqhessian.chk
#p UHF/6-31G(d) freq=(noraman,readisotopes)
H3CO --> H2COH TRANSITION STATE
0 2
6 0.03104 0.63055 0.0
8 0.03104 -0.73696 0.0
1 -0.99138 -0.13555 0.0
1 0.27839 1.12398 0.92614
1 0.27839 1.12398 -0.92614
400 1.5 0.8929
12
16
1
1
1
ReadIsotopes选项：
设置同位素和其它环境参量
M.W.Wong,K.B.Wiberg,and M.J.Frisch,
J.Chem.Phys.,95,8991,(1991)
Computational Chemistry laboratory
Beijing Normal university
计算化学理论和应用
－第九讲分子几何结构优化 (3)
优化收敛问题 (2)
问题 引起原因 解决方案力过大 (a)初始构型太差
(b)分子坐标选择得不好
(a)改变分子构型重新输入
(b)重新构造分子坐标优化结构的 Hessian
矩阵出现负本征值
(a)结构没有优化到极小点
(b)Hessian矩阵更新中数值计算错误
(a)沿负本征值方向优化得到能量更低的结构
(b)改变方法初猜 Hessian矩阵重新优化
(c)如果 (b)仍不行，冻结负本征值所在的坐标矢量优化收敛后再放开优化
Hessian矩阵的本征值太大 (少见 )
(a)Hessian矩阵初猜不好
(b)Hessian矩阵更新不好
(c)不同坐标之间偶合过强
(a)(b) 重新初猜 Hessian矩阵
(c) 重新构造分子坐标优化步数太多 (a)优化比较困难，需要的步数确实很多
(b)输入内坐标中存在冗余
(c)存在很松弛的坐标
(d)不同坐标之间偶合过强
(e)初猜不好
(a)增加优化次数 maxcyc=n
(b)重新输入内坐标
(c)冻结松弛的坐标优化收敛后再放开优化
(a)重新构造分子坐标优化中的对称性问题优化到高对称性一般可以做到，但是收敛非常慢，而且受算法的限制，
只能得到近似的结果
*解决办法：提高到相应对称性直接优化然后作比较优化到低对称性优化中计算得到的梯度矢量属于分子点群的对称表示，
优化过程中分子对称性不会降低
*解决办法：检查 Hessian矩阵，依负本征值所在矢量方向降低对称性过渡状态理论与反应速率计算过渡状态理论简介
R e a cti o n co o r d i n a te s
A A B
K
a c t
K
d e a c t
K
E
传统过渡态理论的反应速率公式
,0,0( ) /( ) / ABAB U U k TG G k T
A
AdB
k A k A k k
d t A
Q
k k e k e
Q
势能面的特性
1,二次区域 (qudratic region)
2,过渡矢量 (transition vector)
二次区域例,函数在点 (0,0)处的梯度为零，
Hessian矩阵为点 (0,0)是一个鞍点,Hessian矩阵的本征值和本征矢分别是
4 (0,1); 和 -4 (1,0)
22
2
401 2 8 4 1 6
041 6 8 4
x y x y
x y x
4 2 2 2 2(,) 4 2 2f x y x x y x y
-1 - 0,8 - 0,6 - 0,4 - 0,2 0 0,2 0,4 0,6 0,8 1
0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
1
-2 - 1,5 -1 - 0,5 0 0,5 1 1,5 2
-1
0
1
2
3
4
5
6
7
8
过渡结构的优化
1,一般问题沿过渡矢量方向优化到极大，而其他方向优化到极小
2,考虑对称性后的简化
H
C
H
H
XX
过渡结构的优化方法
势能面扫描
梯度模法
Gaussian中的 Berny优化
线性或二次准同步过渡法
坐标驱动法 (本征矢跟踪法 )
线性和二次准同步过渡法
R,P,反应物和产物 TS：过渡结构点 1：线性准同步过渡法寻找到的过渡结构点 3：二次准同步过渡法寻找到的过渡结构
2
( ) /
[ ( ) / ( ) /
/[ 2 ( ) ( ) ]
22
2 2 2 2
||
| | | |
| | | | | | + | |
T P R P R
T a P X P X R X R X
a R X P X R X P X R X P X
坐标驱动法和本征矢跟踪法存在的问题梯度模法
0Tgg?
梯度模法存在的问题梯度模极小值点不一定是驻点，有可能是势能函数的拐点优化过渡结构所使用关键词和选项
*J,Comp,Chem,Vol,3,No,2,214-218(1982)
Berny优化默认使用 RFO法控制优化步长，此时可以用
iop(1/8)=n (0.01*n Bohr)来控制 TSR
如果在使用 Berny优化时选用 Newton法，此时需通过
NrScale来控制步长算法 关键词 选项
Berny优化 * Opt=TS NoEigentest,iop(1/8=6)
STQN Opt=QST2,QST3
本征矢跟踪 Opt=EF CalcAll,CalcFc,EnOnly
不同方法的比较
Re ac t i o n Z - m at r i x i n t e r n al s Re d u n d an t i n t e r n al s
r e g u l ar Ca l cF C Q S T 3 Ca l cF C Q S T 2 Q S T 3
CH F CH HF
4 3
6 4 6 5 8 5
CH O CH OH
3 2
12 9 9 8 8 9
S i H H S i H
2 2 4
11 7 11 7 8 8
C H F C H HF
2 5 2 4
16 12 15 13 17 11
D i e l s- A l d e r r e ac t i o n 56 11 23 8 13 14
Cl ai se n r e ac t i o n 38 8 15 7 15 15
En e r e ac t i o n f ai l 15 28 13 18 18
优化过渡态时可能遇到的问题及解决问题 引起原因 解决方案优化过程中 Hessian
矩阵负本征值过多
(a)优化得到二级鞍点结构
(b)Hessian矩阵有数值错误
(a)沿非反应坐标方向寻找极小点
(b)改变方法初猜 Hessian矩阵优化过程中没有负本征值
(a)结构离鞍点的二次区域太远
(b)Hessian矩阵有数值错误
(a)使用更好的初始结构
(b)改变方法初猜 Hessian矩阵重新优化主反应内坐标变量后加 3计算相应 Hessian矩阵元例,0 1C
N 1 r1
H1 r2 2 a1
r1=1.2
r2=1.4 3
a1=59.0 3
D,K,Malick,G,A,Petersson,and J,A,Montgomery Jr.,J,Chem,Phys,108,5704 (1998).
IRCMax (MP2/6-31G(d):HF/3-21G*,Zero,Stepsize=12)
P,Y,Ayala and H,B,Schlegel,J,Chem,Phys,107,375 (1997).
通过反应途径寻找过渡结构
Opt(QST2,Path=M)
使用两种方法分别优化反应途径并计算反应途径上的能量最高点输入：过渡结构同时优化反应途径上的 M个点和过渡结构输入：同 QST2或 QST3
过渡结构算例
#T UHF/6-31G(d) Opt=QST2
H3CO --> H2COH Reactants
0,2
C
O 1 1.48
H 1 R 2 A
H 1 1.08 2 110,3 120.
H 1 1.08 2 110,3 -120.
R=1.08
A=110.
H3CO --> H2COH Reactants
0,2
C
O 1 1.48
H 1 R 2 A
H 1 1.08 2 110,3 120.
H 1 1.08 2 110,3 -120.
R=1.9
A=30.
E
-1 1 4,4 1 4 6 0
-1 1 4,3 3 0 5 9
-1 1 4,3 8 8 5 1
过渡结构的确认
计算完整的 Hessian矩阵，确保有且仅有一个负本征值通过振动分析计算完成，可以使用连续作业
# B3LYP /6-31g opt freq
确认本征矢量的方向正确连接到反应物和产物通过计算反应途径进行确认
IRC(内禀反应坐标 )原理和计算
IRC的定义和计算方法反应途径：
势能面上，经过过渡结构，连接反应物和产物的结构变化途径
特性：
1,满足 最小作用 原理
2,原子核的运动速度无限小
推论：
1,反应途径是势能面上连接反应物,过渡结构和产物的最小能量途径
2,反应途径上每点的梯度方向与运动方向一致
寻找反应途径的方式：
从过渡态出发寻找连接反应物和产物的最陡下降途径最陡下降途径的定义式：
运动方向,梯度：
反应途径：
选择不同的坐标系，可以有不同的最陡下降途径
i i i
ii i
jik
i j k
V
x d x x
x
dxd x d x
c
V V V
x x x
Fukui的 IRC方程
jik
i j k
dXd X d X
c
V V V
X X X
质权坐标系下的最陡下降途径，就是最小能量途径
K.Fukui,Acc,Chem,Res,1981,Vol,14,12,363
1,确认过渡结构
2,寻找反应机理
IRC计算的意义
Gaussian程序中 IRC的算法
C,Gonzalez and H,B,Schlegel,J,Chem,Phys,90,2154 (1989).
C,Gonzalez and H,B,Schlegel,J,Phys,Chem,94,5523 (1990).
IRC计算的步骤和选项
1,优化出过渡结构
2,对优化结果进行振动分析计算确认过渡结构，得到零点能，获得力常数矩阵
3,进行 IRC计算力常数选项：
从振动分析计算的 chk文件中读取 IRC=(RcFc)
使用当前方法计算 IRC=(CalcFc)
反应坐标数的选取 IRC(MaxPoints=n) 默认,n=6
点数不够时可以重起作业：
IRC=(Restart,MaxPoints=n)
步骤：
选项：
IRC算例
%chk=freqhessian.chk
#p UHF/6-31G(d) IRC=(maxpoints=10,rcfc,forward)
H3CO --> H2COH IRC
0 2
6 0.03104 0.63055 0,
8 0.03104 -0.73696 0,
1 -0.99138 -0.13555 0,
1 0.27839 1.12398 0.92614
1 0.27839 1.12398 -0.92614
振动分析计算 (2)
Hessian矩阵本征值与振动模式的关系
1,Hessian矩阵的本征值对应于简正振动的力常数
2，简正振动的频率等于力常数的平方根
3，简正振动的振动模式是相应本征值的本征矢热力学函数的统计表达式
*
*
/
l n l n ( ! )
!
i
i
i
n
i
i
kT
ii
g
S k k N
n
N
n g e
f
2
,
,
ln
ln
ln
B
NV
BB
NV
f
E Nk T
T
H E P V
f
S Nk T f Nk T
T
G H TS
配分函数和热力学函数各个分量的计算平动配分函数，平动能和平动熵转动配分函数，转动能和转动熵
e t rf f f f f
析因子性质：
32
2
2
32
2
2
8 3
()
2
8 3
l n ( )
2
a b c B
rr
a b c B
r
I I I kT
f E RT
h
I I I kT
SR
h
3 / 2 3 / 2
22
2 3 2 3( ) l n ( )
22t t t
m k T m k Tf V E R T S R V
hh
振动配分函数，振动能和振动熵热容的计算
V
V
UC
T
36
1 / 2
/
1
36
/
1
36
/
/
1
1
( 1 )
1
( 1 )
l n ( 1 )
( 1 )
i iB
iB
iB
iB
n
i i i i v h k T
i
n
i
v h k T
i B
n
h k Ti
v h k T
i B
f
e
h
ER
ke
h
S R e
k T e
振动分析算例
%chk=freqhessian.chk
#p UHF/6-31G(d) freq=(noraman,readisotopes)
H3CO --> H2COH TRANSITION STATE
0 2
6 0.03104 0.63055 0.0
8 0.03104 -0.73696 0.0
1 -0.99138 -0.13555 0.0
1 0.27839 1.12398 0.92614
1 0.27839 1.12398 -0.92614
400 1.5 0.8929
12
16
1
1
1
ReadIsotopes选项：
设置同位素和其它环境参量
M.W.Wong,K.B.Wiberg,and M.J.Frisch,
J.Chem.Phys.,95,8991,(1991)
