过渡态优化
在本教程中，我们将展示如何使用MLatom 优化 和 检查 过渡状态(TS)几何的性质。此处我们以乙烯和1,3-丁二烯之间的Diels-Alder反应为例。
我们提供了一个演示视频，您可以在 Youtube 上观看。
MLatom支持各种方法，包括QM方法、预定义的或用户定义的机器学习方法以及混合量子力学/机器学习方法，如 概览 所示。在这里，我们选择快速和流行的半经验方法 GFN2-xTB ，您可以 参阅下文 来了解如何使用，例如 AIQM1 (通用的混合QM/ML方法)和B3LYP/6-31G*，以及纯ML模型(预训练模型或用户训练的模型)。
过渡态的几何构型优化
过渡态的几何优化可能相当复杂(此种情况下的处理方法请参见 下文 )。但是，如果对过渡态的结构有一个很好的初始猜测，优化过渡态可能就会很容易。这里我们提供了一个 初始猜测 (以埃为单位):
16
C 0.48462430 -0.55755495 1.43729151
C 0.48462430 -0.55755495 -1.43729151
C -0.27595797 -1.44977527 0.70359025
C -0.27595797 -1.44977527 -0.70359025
C -0.27595797 1.45086377 0.69299925
C -0.27595797 1.45086377 -0.69299925
H 0.37292526 -0.50748993 2.51767690
H 1.44526264 -0.21636383 1.06867438
H 0.37292526 -0.50748993 -2.51767690
H 1.44526264 -0.21636383 -1.06867438
H -1.05536225 -2.01444047 1.21328943
H -1.05536225 -2.01444047 -1.21328943
H 0.51071931 1.96707995 1.23581344
H -1.20625330 1.32768072 1.23591744
H 0.51071931 1.96707995 -1.23581344
H -1.20625330 1.32768072 -1.23591744
用户需要提供输入文件，与最小能量的 几何优化 相似，只需要在输入文件(命令行选项)中将 geomopt 替换为 TS :
GFN2-xTB # 1. define QM method or ML or QM/ML model
TS # 2. request TS geometry optimization
xyzfile=da_ts_init.xyz # 3. initial geometry guess
optxyz=da_ts_opt.xyz # 4. file with optimized geometry
在Python API中，代码也需要做出修改，即添加 ts=True 作为优化例程的选项：
import mlatom as ml
# Get the initial guess for the molecules to optimize
initmol = ml.data.molecule.from_xyz_file('da_ts_init.xyz')
# initmol charge and multiplicity can be changed if required, e.g., by uncommenting the lines below:
#initmol.charge = 1 ; initmol.multiplicity = 2
# Choose one of the predifined (automatically recognized) methods
mymodel = ml.models.methods(method='GFN2-xTB')
# or QM method, e.g., B3LYP with Gaussian
#mymodel = ml.models.methods(method='B3LYP/6-31G*', program='Gaussian')
# or hybrid QM/ML method
#mymodel = ml.models.methods(method='AIQM1')
# or ML model, e.g., MACE:
#mymodel = ml.models.mace(model_file='mace.pt')
# Optimize the geometry with the choosen optimizer:
geomopt = ml.optimize_geometry(model=mymodel, initial_molecule=initmol, ts=True, program='ASE', maximum_number_of_steps=10000)
# Get the final geometry
final_mol = geomopt.optimized_molecule
# and its XYZ coordinates
final_mol.xyz_coordinates
# Check and save the final geometry
print('Optimized coordinates:')
print(final_mol.get_xyz_string())
final_mol.write_file_with_xyz_coordinates(filename='da_ts_opt.xyz')
# You can also check the optimization trajectory
print('Number of optimization steps:', len(geomopt.optimization_trajectory.steps))
# and access any geometry during optimization as
intermediate_step_geom = geomopt.optimization_trajectory.steps[2].molecule
print('Intermediate coordinates:')
print(intermediate_step_geom.get_xyz_string())
在命令行模式下使用MLatom 输出文件 中的代码片段如下(使用Gaussian优化器，其他选项和注意事项请参阅 此文 ):
==============================================================================
Optimization of molecule 1
==============================================================================
Iteration Energy (Hartree)
1 -17.8096315123480
2 -17.8106949632610
3 -17.8116884258690
4 -17.8120367568860
5 -17.8122673179990
6 -17.8122596563090
7 -17.8122592324260
Final energy of molecule 1: -17.8122592324260 Hartree
==============================================================================
The final optimized geometry (in da_ts_opt.xyz file obtained in either way) is:
16
C 0.4813024374103 -1.4349227612885 0.5137386536117
C 0.4859702737149 1.4344889056846 0.5129530071391
C 1.2680466082266 -0.7096089630187 -0.3161249847856
C 1.2698482503147 0.7063075171378 -0.3168590194909
C -1.5743749532249 -0.6736952829099 -0.2267301172749
C -1.5741402429035 0.6770058850694 -0.2253736447211
H 0.4107820982713 -2.5048539561628 0.4024004256024
H 0.1545008650554 -1.0465393829088 1.4604758694958
H 0.4182444345553 2.5045006285090 0.4008820574544
H 0.1575374270938 1.0474822961018 1.4596556145266
H 1.7590868456084 -1.2036512202944 -1.1442175065287
H 1.7624326531584 1.1982575357438 -1.1452660986280
H -1.9999454494089 -1.2291294376367 0.5911010005786
H -1.4029477841077 -1.2229351940318 -1.1358103575475
H -1.9978808718237 1.2312093010793 0.5942018379777
H -1.4017244596312 1.2282076255495 -1.1330462098024
检查过渡态
当找到了过渡态结构以后，需要检查其正确性。最简单的方式是对优化得到的结构进行频率计算(请参阅我们之前关于 频率计算 的教程)，然后分析虚频的正常模式。MLatom会保存频率计算的 gaussian.log 文件。用户可以通过多种可视化软件（如GaussView）进行查看。以下是GaussView中虚频(在XACS云平台上使用GFN2-xTB计算得到在正态模式下的振动：
另一种方法是进行内禀反应坐标(IRC)计算(详见 IRC教程 )，这通过Gaussian接口运行，输入文件如下:
IRC
GFN2-xTB
xyzfile=da_ts_opt.xyz
在Python API中，我们通过 mlatom.simulations.irc() 来执行IRC计算。用户需要提供 model 和 ts_molecule ，即待检查的过渡态构型。
import mlatom as ml
# Get optimized molecules
optmol = ml.data.molecule.from_xyz_file('da_ts_opt.xyz')
# Choose one of the predifined (automatically recognized) methods
mymodel = ml.models.methods(method='GFN2-xTB')
# or QM method, e.g., B3LYP with Gaussian
#mymodel = ml.models.methods(method='B3LYP/6-31G*', program='Gaussian')
# or hybrid QM/ML method
#mymodel = ml.models.methods(method='AIQM1')
# or ML model, e.g., MACE:
#mymodel = ml.models.mace(model_file='mace.pt')
IRC = ml.simulations.irc(model=mymodel,ts_molecule=optmol)
IRC计算结束以后，用户可以通过 gaussian.log 文件查看IRC曲线（由GFN2-xTB优化的几何构型生成）
复杂情况的处理
过渡态优化并不能总是找到正确的过渡态。因此，用户可能需要对计算的设置和使用做出进一步修改，例如，使用将反应物和生成物的几何构型考虑在内的算法。目前，如果您熟悉Gaussian程序，我们建议先使用MLatom进行初步优化，然后修改其生成的 gaussian.com ，以此作为输入文件，例如使用QST2或QST3方法。以下是Diels-Alder反应的例子:
%nproc=1
#p opt(qst3) external='/share/apps/anaconda3/bin/python3 /export/home/xacscloud/mlatom_versions/mlatom/interfaces/gaussian_interface.py'
reactants
0 1
C 0.484624300000 -0.557554950000 1.437291510000
C 0.484624300000 -0.557554950000 -1.437291510000
C -0.275957970000 -1.449775270000 0.703590250000
C -0.275957970000 -1.449775270000 -0.703590250000
C -0.311976759298 2.284843109406 0.344043544386
C -0.311976759298 2.284843109406 -1.041954955614
H 0.372925260000 -0.507489930000 2.517676900000
H 1.445262640000 -0.216363830000 1.068674380000
H 0.372925260000 -0.507489930000 -2.517676900000
H 1.445262640000 -0.216363830000 -1.068674380000
H -1.055362250000 -2.014440470000 1.213289430000
H -1.055362250000 -2.014440470000 -1.213289430000
H -0.921357965784 1.567908780096 0.886857734386
H 0.545932937341 2.665125132700 0.886961734386
H -0.921357965784 1.567908780096 -1.584769145614
H 0.545932937341 2.665125132700 -1.584873145614
products
0 1
C -0.661894740000 -0.059709790000 -1.393910490000
C -1.486209330000 -0.113548400000 -0.144719260000
C -0.691627710000 0.317536480000 1.081726060000
C 0.691627710000 -0.317536480000 1.081726060000
C 1.486209330000 0.113548400000 -0.144719260000
C 0.661894740000 0.059709790000 -1.393910490000
H -1.190276150000 -0.118797780000 -2.341761460000
H -2.373254720000 0.518104980000 -0.263582030000
H -1.871283790000 -1.133268850000 -0.007315520000
H -1.235540190000 0.060874970000 1.994969250000
H -0.583769070000 1.408342230000 1.074593440000
H 1.235540190000 -0.060874970000 1.994969250000
H 0.583769070000 -1.408342230000 1.074593440000
H 2.373254720000 -0.518104980000 -0.263582030000
H 1.871283790000 1.133268850000 -0.007315520000
H 1.190276150000 0.118797780000 -2.341761460000
TS
0 1
C 0.484624300000 -0.557554950000 1.437291510000
C 0.484624300000 -0.557554950000 -1.437291510000
C -0.275957970000 -1.449775270000 0.703590250000
C -0.275957970000 -1.449775270000 -0.703590250000
C -0.275957970000 1.450863770000 0.692999250000
C -0.275957970000 1.450863770000 -0.692999250000
H 0.372925260000 -0.507489930000 2.517676900000
H 1.445262640000 -0.216363830000 1.068674380000
H 0.372925260000 -0.507489930000 -2.517676900000
H 1.445262640000 -0.216363830000 -1.068674380000
H -1.055362250000 -2.014440470000 1.213289430000
H -1.055362250000 -2.014440470000 -1.213289430000
H 0.510719310000 1.967079950000 1.235813440000
H -1.206253300000 1.327680720000 1.235917440000
H 0.510719310000 1.967079950000 -1.235813440000
H -1.206253300000 1.327680720000 -1.235917440000
备注
由于Python和MLatom的安装路径应该不同，故以上输入文件仅作为说明的例子。您还需要在执行初步过渡态优化的文件夹(或其副本)中运行该文件，因为它生成了额外的必需文件，如 model.json 。当然，专业用户也可以自行创建 model.json 文件。
选择方法或模型
MLatom提供了多种用于优化的方法或模型。用户可以查看 单点能计算教程 来了解如何定义这些方法。这里我们给出一些简单的例子来快速入门。MLatom可以自动识别的方法有很多，如ANI-1ccx或AIQM1，如 手册 所示。
定义QM方法的一般步骤是使用 method 和相应的QM程序(MLatom使用这些程序的接口)。例如，可以选择从头算方法(如HF或MP2)或DFT方法(如B3LYP/6-31G*):
method=B3LYPG/6-31G* # 1. request running DFT optimization with B3LYP
qmprog=PySCF # 2. request using PySCF for B3LYP calculations; qmprog=Gaussian can be also used
TS # 3. request transition state geometry optimization
xyzfile=da_ts_init.xyz # 4. initial geometry guess
optxyz=da_ts_opt.xyz # 5. file with optimized geometry
这里我们提供了与过渡态参考值(在B3LYP/6-31G*水平)的优化时间、迭代次数和RMSD的对比。所有的计算都在 XACS云平台 上执行。
方法 |
时间（min） |
迭代次数 |
RMSD (\(\unicode{x212B}\)) |
|---|---|---|---|
GFN2-xTB |
0.50 |
7 |
0.03 |
AIQM1 |
1.67 |
18 |
0.07 |
对于这种特殊情况，GFN2-xTB具有出色的性能。但AIQM1也不错，我们的 大规模基准测试 表明，它通常比GFN2-xTB更鲁棒、更快。因此，这两种方法都可用于过渡态的初步优化，然后再使用DFT优化。
备注
AIQM1目前仅适用于含有CHNO元素的化合物。另外，为了释放其全部潜力，我们建议为其QM部分安装MNDO程序。XACS云平台使用的是Sparrow程序，该程序目前只提供数值梯度，限制了其对几何优化的适用性。
使用用户训练的机器学习模型进行优化的一般方法是使用 MLmodelType 和 MLmodelIn 关键字指定机器学习模型类型，例如:
MLmodelType=MACE # 1. request optimization with the MACE-type of machine learning potential
MLmodelIn=mace.pt # 2. the file with the trained model should be provided, here it is mace.pt file
TS # 3. request transition state geometry optimization
xyzfile=init.xyz # 4. initial geometry guess
optxyz=opt.xyz # 5. file with optimized geometry
关于如何训练此类模型，请参阅相关的手册和教程，例如 MACE 。
高级的用户定义混合模型，如基于delta-learning概念的模型，只能通过Python API进行优化(参见 几何优化中的示例 )。
选择优化程序
For transition state geometry optimization, MLatom supports the use of Gaussian and ASE optimizer. They can be chosen by adding optporg=[ASE, geometric or Gaussian] keyword as :
GFN2-xTB # 1. define the method
TS # 2. requests transition state geometry optimization
xyzfile=init.xyz # 3. initial geometry guess
optxyz=opt.xyz # 4. file with optimized geometry
optprog=gaussian # 5. Gaussian is choosen as the optimizer
If no optimizer is requested, MLatom will check the availability of the programs in the order Gaussian -> geometric -> ASE.
各个优化器都有各自的优缺点：
ASE是开源的ASE，并已经经过测试。它可以在XACS云平台中使用，并且可以很容易地安装。
GeomeTRIC is open-source and features its various choices for coordinate systems. Compared with Gaussian, it may take much more iterations to converge.
Gaussian是一个非常通用的程序，MLatom的Gaussian接口允许用户使用Gaussian的所有功能，并且还能访问ML模型。Gaussian应由用户单独获取并安装。由于目前的I/O效率并不高，因此使用非常快的ML模型进行优化将比使用其他优化器花费更长的时间。这再将来可能会得到改进。
备注
MLatom正在快速发展，未来将具备更多功能和更好的几何优化能力。我们的教程也会做出相应修改，请关注我们的更新。
无论您选择哪种优化器，以下注意事项和提示信息都值得了解。
ASE优化器
目前，ASE的过渡态优化只支持dimer方法，后续我们也可能实现微扰弹性带方法(NEB)。
备注
当使用dimer方法时，它在当前的MLatom实现中会给出不同的结果(因为ASE中的默认设置是使用random seed)。因此，如果多次进行优化，每次都会获得不同的结果(其中一个可能是用户需要的结果)。与Gaussian优化器相比，dimer方法也需要大量的迭代。
geomeTRIC optimizer
GeomeTRIC is an open-source software which introduce a new translation-rotation-internal coordinate (TRIC) system to describe molecule structure for efficient geometry optimization. It accepts external software for the energy and gradient, which is provided by various methods in MLatom, through wrapper functions. By default, geomeTRIC uses BFGS for optimization task.
On XACS cloud, geomeTRIC optimizer is the default option if optprog is not specified in input file and program is not set for mlatom.optimize_geometry() with python.
For example, you can try the input below using our newly updated UAIQM method to perform previous task on XACS cloud:
uaiqm # 1. define QM method or ML or QM/ML model
TS # 2. request TS geometry optimization
xyzfile=da_ts_init.xyz # 3. initial geometry guess
optxyz=da_ts_opt_uaiqm.xyz # 4. file with optimized geometry
And we provide the output file and optimized molecule as reference, as well as the results from frequency calcuclation as validation.
Gaussian优化器
对于Gaussian优化器，计算等效于在Gaussian输入文件中使用 opt(ts,calcfc,noeigen,nomicro) (参见工作目录中的 gaussian.com 文件)。
当使用Gaussian优化器时，您将看到与普通Gaussian输出文件相似的 gaussian.log 文件，不同之处在于，它将使用MLatom的模型来预测能量和梯度以执行几何优化。用户可以自行选择工具来分析这个输出，这对Gaussian用户来说是一个很大的优势。
完整的 Gaussian输出文件 可供下载。
另一个重要的优点是用户可以修改 gaussian.com 文件，该文件包含使用MLatom模型进行过渡态优化所需的输入。例如，如 前一节 所示，如果您想在Gaussian中使用各种优化算法，可以将 opt(ts) 替换为 opt(qst3) 或其他算法。此外，通过这种方式，用户可以使用Gaussian的所有功能进行几何优化，如冻结坐标或在内部坐标中定义输入，等等。
在修改完 gaussian.com 以后，用户可以照常执行Gaussian任务，例如，在用户的系统中使用 g16 gaussian.com 命令行。
问题或建议
如果您有进一步的问题、批评和建议，我们很乐意通过 电子邮件 、 Slack (首选)或微信(请发送电子邮件请求将您添加到XACS用户支持组)以英文或中文接收您的意见。