Gau_UMA: 联用Gaussian和UMA机器学习势

zhouoh

0 前言
之前一直在关注通用机器学习力场的发展，但是受制于模型普适性的种种限制（如只适用于C, H, N, O元素，只能用于中性分子等）一直没有上手尝试。最近看到meta公布了他们创建的83M OMol25数据集以及基于其训练的UMA机器学习势能面，据称对全周期表和各种电荷，多重的的体系都有良好的表现，同时收到了论坛里之前帖子的启发（使用机器学习力场UMA优化分子结构与搜寻过渡态初猜 - 量子化学 (Quantum Chemistry) - 计算化学公社），萌生了使用UMA代替xTB进行体系预优化的想法。然而UMA 原生只提供了基于ASE的计算借口，直接使用的话每次都要写脚本不是很方便，而且相较于Gaussian，ASE对于过渡态优化，频率分析，IRC等的支持也有限。于是我写了一个接口使得Gaussian可以调用UMA模型进行几何优化，振动分析和IRC计算。

1 接口设计
程序具体的设计参考了sob老师的gau_xtb程序（将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析 - 量子化学 (Quantum Chemistry) - 计算化学公社），其中大致的思路是在Gaussian每次调用接口程序时将Gaussian传递的信息转化成xTB的输入，之后在命令行中调用xTB执行计算并提取结果。然而，和xTB这类的cpu计算程序不同，UMA在冷启动（即从命令行直接开始运行计算脚本）的时候涉及初始化GPU计算环境，将参数文件读入显存等等操作，而这些会导致几秒钟的延迟。考虑到在GPU上UMA一次推理的时间一般只有十几到几十毫秒，在每轮计算中重新启动UMA环境会大大拖慢整体的速度。所以我把整个接口程序分成了两个部分：在后台常驻运行的UMA服务器端和用于被Gaussian调用的客户端脚本。具体的结构如下：

其中GauUMA_server.py为常驻后台的服务器，而GauUMA_client.py则是可以直接被Gaussian调用，两者之间通过socket套接字进行沟通。这样就保证了UMA模型只需要在server里被加载一次，同时附带的好处是执行Gaussian计算的机器可以和执行GPU推理的机器相互分离，两者之间只通过局域网相互通信。在写server端代码的时候我尽量只使用了ASE包的一些标准操作，这样保证了理论上其他提供ASE接口的模型，比如orbmol，DPA3等只需要稍微修改代码就可以实现相同的接口。

2 使用方法和示例

2.1 机器配置与环境

• CPU: Ryzen 7950X3D
• RAM: 64GB DDR5 6400MHz
• GPU: RTX 4090
• OS: Ubuntu 22.04
• Python Version: 3.10
  

2.2  Pd催化中间体IM4R的优化


这里我使用了和@LPOweChan在他的帖子使用机器学习力场UMA优化分子结构与搜寻过渡态初猜 - 量子化学 (Quantum Chemistry) - 计算化学公社中一样的IM4R体系。在计算之前需要先获取UMA的权重文件以及安装meta提供的Fairchem库，具体的安装方法可以参考上面的帖子。安装完成后，执行

1. python3 GauUMA_server.py (或者 nohup python3 GauUMA_server.py &! 放入后台运行)

复制代码

便会启动接口服务端，程序默认使用50423端口进行通信，也可以在代码中自行修改。
优化时使用的Gaussian输入文件为：

1. %nproc=1
   
2. #P opt=(nomicro)  external='python path/to/GauUMA_client.py'
   
3. 
   
4. t
   
5. 
   
6. 0 1
   
7. 坐标略...

复制代码

使用Gaussian运行后经历82步收敛，用时仅为19.6s。平均下来对于这个104原子的体系，每步仅耗时0.24s。观察优化期间的GPU占用发现GPU占用率仅为10%，功耗也只有70w左右。Gaussian输出文件报告Gaussian本身占用了14.2s的cpu时间，说明决速步应该是在Gaussian内部的计算了。计算过程的显存占用约为4.3 GB，应该大部分近几年的游戏显卡都可以无压力运行上百原子体系的几何优化。

几何优化之后，也同样可以通过柔性扫描来初步定位C-N还原消除的过渡态：

1. %nproc=1
   
2. #P opt=(nomicro,modredundant,loose)  external='python ../GauUMA_client.py'
   
3. 
   
4. t
   
5. 
   
6. 0 1
   
7. 坐标略...
   
8. 
   
9. 47 61 S 15 -0.1

复制代码

15步扫描用时51.6s。扫描成功定位了一个最高点和过渡态的初始结构：

接下来同样可以使用通用的关键词进行过渡态优化，振动分析和IRC计算：

1. %nproc=1
   
2. #P opt=(nomicro,calcfc,ts,noeigen)  external='python ../GauUMA_client.py'
   
3. 
   
4. t
   
5. 
   
6. 0 1

复制代码

1. #P freq  external='python ../GauUMA_client.py'

复制代码

1. #P irc=(lqa,calcfc)  external='python ../GauUMA_client.py'

复制代码

需要注意的是UMA本身只能预测分子的受力，因此calcfc和freq所需的Hessian矩阵是由程序自动进行有限差分得到的。不过由于有限差分是在内部实现不涉及和Gaussian的通信，所以整体效率较高，对于这个过渡态体系一次频率计算仅耗时17.8s。这里理论上还可以通过批量推理或者使用@wzkchem5老师的O1NumHess进一步优化，不过考虑到程序的兼容性以及本身计算的绝对时间也不长，就没有进一步实现。以下是过渡态结构中的一些关键键长和PBE0结果的对比，可以说作为初猜已经非常可用了：

  

得到的IRC曲线也非常平滑：

3 更多示例和测试

上一节中测试的体系虽然比较大，但是涉及的还原消除过渡态总体比较简单，然后Pd金属本身的电子结构也比较友善。于是我又测试了一个Cr催化的乙烯插入过渡态，其中涉及多个C-C键，Cr-C键的形成和断裂，同时Cr也是三重态的开壳层结构。测试体系来源于Chem. Sci., 2020,11, 9665-9674, 可以看到整体还是成功维持了过渡态的形状。

接着又几个对称破缺体系的表现，首先测试了乙烷中C-C键的解理曲线，UMA成功描述了C-C键在接近断裂时候的行为。给出的解离能和实验值90 kcal/mol也大致相符（粗略估计，实验上的键能还需要分别优化并且考虑热力学矫正）。对于乙烯扭转的结构，UMA成功优化出了两个CH2垂直的过渡态结构，但是其给出的能垒高达92 kcal/mol。倒是和典型的RKS结果类似。这表明在这一类电子结构比较复杂的体系上使用此类机器学习势还是要谨慎对待结果。

最后对于@LPOwenChan提出的UMA不适用于分子间Diels-Alder反应的观察，我这里的测试发现UMA完全可以成功的定位乙烯和丁二烯[4+2]的过渡态，其中的键长为2.23 A，与DFT计算值2.25 A符合极好。我猜测问题问题可能出现在NEB的参数或者初始结构设置上，贴中提供的过渡态结构经过calcall优化下得到的是丁二烯s-cis 到 s-trans构型变化的虚频。具体输出文件放在附件中可自行查看。

4 代码和所有实例输入输出文件
见附件。欢迎大家提出意见和建议。

LPOwenChan

      感谢老师的倾情开发与关注，针对[4+2]环加成的问题我在之后也进行了进一步的分析，确实如老师所言是“初始结构设置”的问题，更确切地说是我的博文当中的“先优化反应物与生成物结构，然后优化NEB”这一思路可能并不适用于双分子加成这一类型反应。我被限制在了之前算单分子内C-C偶联的思路里，因为反应物和生成物都是“单一分子”，用UMA进行优化理论上不会有任何问题，NEB也能很好地给出过渡态初猜。但是[4+2]加成的“反应体系”显然不是单分子体系，两分子接近的方向与相距的距离依赖于手动摆放，并且我们也希望反应体系的初始结构就是“手动摆放的位置”，因此对于手动摆放好的体系再去使用UMA优化在逻辑上是不合理的。在我观察被错误优化的1,3-丁二烯+乙烯的初始结构后发现，该结构已经非常倾向于1,3-丁二烯构象变化的过渡态。对这样的结构进行NEB也不可能得到正确的加成过渡态，因为反应体系的原始摆放位置已经被完全破坏了，因此不应该对初始体系进行优化。
      仅对于我的代码而言，解决方案也是很简单的，就是对于“双分子反应体系”的反应物仅计算体系能量而不去优化体系结构，但是依然优化生成物分子结构（将原来的优化反应物与生成物变为仅优化生成物），然后将二者分别设置为NEB的一头一尾，就可以得到合理的能量高点与合理的过渡态初猜。

Daniel_Arndt

LPOwenChan 发表于 2025-10-8 15:20
感谢老师的倾情开发与关注，针对[4+2]环加成的问题我在之后也进行了进一步的分析，确实如老师所言是 ...

我之前用GFN-xTB和一些软件进行NEB找过渡态时，有个trick就是观察GFN-xTB优化后的反应物和产物，对于其中结构较“松散”的一个，我会用 https://xtb-docs.readthedocs.io/ ... training-potentials 中的constrain功能重新优化出一个较“紧密”的结构，这样的话，NEB的两端，一端是势能面上的极小点，另一端不是势能面上的极小点，但比较靠近过渡态。

Marucs

请问 GauUMA_server.py和GauUMA_client.py应该放在哪个文件下，我执行的时候显示ModuleNotFoundError: No module named 'ase'
谢谢老师！

zhouoh

Marucs 发表于 2025-10-23 05:41
请问 GauUMA_server.py和GauUMA_client.py应该放在哪个文件下，我执行的时候显示ModuleNotFoundError: No m ...

理论上放在任意文件夹都可以。报错应该是没有安装ASE包。执行pip install ase就好了

Marucs

zhouoh 发表于 2025-10-24 11:29
理论上放在任意文件夹都可以。报错应该是没有安装ASE包。执行pip install ase就好了

感谢解答

Huschein

感觉倒数第二张图蛮好看的 想请问一下用什么工具做的 就是两个分子重叠的那个图

zhouoh

Huschein 发表于 2025-11-21 06:58
感觉倒数第二张图蛮好看的 想请问一下用什么工具做的 就是两个分子重叠的那个图

是用Chimera X结合里面的soft风格画的

Huschein

来推广一下MAPLE 现在不再需要手写脚本接入ORCA或者Gaussian 在MAPLE中即可直接使用各种机器学习势（包含UMA），进行结构优化 势能面扫描 过渡态搜索等诸多任务
http://bbs.keinsci.com/thread-57731-1-1.html