计算化学公社

标题: 使用机器学习力场UMA优化分子结构与搜寻过渡态初猜 [打印本页]
作者
Author: LPOwenChan    时间: 2025-10-4 22:26
标题: 使用机器学习力场UMA优化分子结构与搜寻过渡态初猜
本帖最后由 LPOwenChan 于 2025-10-8 15:26 编辑

使用机器学习力场UMA优化分子结构与搜寻过渡态初猜
Optimizing Molecular Structures and Searching for Transition State Initial Guesses Using Machine Learning Force Field UMA


0.前言
1.UMA简介与工作流程构建
2.在Windows中的PyCharm IDE中创建分子结构探索环境
3.使用UMA与ASE对分子结构进行优化及代码思路
4.使用UMA与ASE创建NEB对分子内C-C偶联反应进行过渡态初猜搜索及代码思路
5.在探索中发现的问题与讨论
6.总结与展望

**********
2025-10-08 更新：
有关5.4节中的问题在Gau_UMA: 联用Gaussian和UMA机器学习势一贴中进行了更为详细的讨论，代码修改思路也可以根据贴中回复进行。
**********

0.前言
0.0 灵感来源               
        本文受David_R的《Using the Universal Model for Atoms for single point calculations》一文启发，并继续参考了微信公众号"计算化学简讯"的《UMA能预测金属有机催化过渡态》一文，结合了自身研究中需要将其加入工作流程的实际需求，从而写下了本文。本文的附件可以从这里下载。

0.1 注意事项                        
        由于本人习惯将任务的预处理(建模、预优化)与后处理(波函数分析、数据处理与绘图)在个人PC上完成，仅将较繁重的计算任务提交至集群，因此本UMA使用例将全程在没有NVIDIA GPU的笔记本上完成，一方面可以为论坛各位展示在孱弱性能的笔记本上运行模型是什么体验，另一方面也可以为没有高性能GPU服务器的各位提供一种可能性，即近几年的笔记本电脑完全有可能可以仅仅依靠内存而非GPU显存来运行小参数的模型。

0.2 硬件与软件配置
0.3 警告        
        不要将UMA优化的结构作为最终结构！不要将UMA优化的结构作为最终结构！不要将UMA优化的结构作为最终结构！ 很明显，机器学习势的精度在大部分情况下依然无法达到DFT精度，本文仅提供UMA可能的结构预处理应用场景，在有明确的Benchmark文章确认其精度或对该模型有明确的可解释性之前都不要将UMA优化的结构作为最终结构！不要将UMA优化的结构作为最终结构！不要将UMA优化的结构作为最终结构！

1. UMA简介与工作流程构建1.0 UMA是什么                              
       UMA(Universal Models for Atoms)是Meta的FAIRCHEM团队提出的一系列通用模型，旨在替代传统的密度泛函理论计算，以更快、更准确地预测原子系统的能量、力、应力等物理性质，适用于分子、材料、催化剂、分子晶体等多种化学领域。他的训练集超过了5亿个体系，个人比较在意的是UMA中的omol任务类型，它可以用于优化纯粹有机分子、金属络合物等分子体系，其训练集是OMol25数据集在ORCA中以ωB97M-V/def2-TZVPD进行优化得到的结构训练集。

1.1 为什么我尝试将UMA加入工作流当中                                 
        在目前的工作当中，时常需要对金属有机催化物种及其反应路径进行分析，这些物种通常是由有机配体-金属中心-反应底物三者共同构成的物种。对于牵涉到这类物种的反应路径研究，完整的工作流程通常是: 建模(ChemDraw/Chem3D/GaussView) → ”预“预优化(Avogadro) → 构象搜索或预优化(CREST/xTB) → 结构优化与相关性质计算(Gaussian/ORCA) → 数据处理与波函数分析(Origin/Multiwfn)。
        在发现UMA之前，本人经常使用xTB来完成构象搜索或预优化任务。具体地说，使用CREST结合GFN-FF力场进行构象搜索，使用xTB结合GFN2(1/0)-xTB来进行粗糙的结构预优化以缩短DFT的结构优化时间。而目前的问题在于，对金属有机催化物种而言，无论是GFN-FF力场还是GFN2(1/0)-xTB，都在某些情况下无法对体系做到定性正确。在我的实际工作中遇到过GFN2(1)-xTB无法很好维持Pd-O/Pd-C/Pd-P/Pd-N/Pt-X等金属络合物体系的配位键，优化过程中体系倾向于配体从金属中心解离而无法继续优化配体；也遇到过GFN-FF力场在搜寻金属配合物构象、或是进行极为粗糙的优化时，金属中心跟随烷基配体直接脱离合理的配位几何结构的情况。尽管上述的情况可以通过冻结原子位置或者定义键连关系来避免，但是作为“预”优化流程如果依然需要较多的步骤来对体系进行限制，可能并非是一个性价比较高的方法。
        UMA力场似乎可以在粗略优化金属配合物的层面上更稳健。无论是论文中所展示训练集，甚至是官方给出的Demo页面中直接映入眼帘的Simulation inputs，都无不强调该力场在处理这类物种时的出色能力，因此如果有望能将UMA力场加入工作流或者能替换一些不那么稳健的方法，那在金属络合物预处理方面是十分有吸引力的。

2. 在Windows中的PyCharm IDE中创建分子结构探索环境2.1 在PyCharm中创建新的Python解释器        
        为了能让模型在一个干净的环境下运行并解决各种包的兼容性问题，我们可以在PyCharm中创建新的Python解释器。具体的操作可以是：设置(Ctrl + Alt + S) → Python → 解释器 → 添加解释器 → 添加本地解释器。
        在当前的"添加Python解释器"页面，”环境“选择”生成新的“；”类型“就用最简单的Virtualenv；”基础Python“就可以使用当前环境中的Python，例如我这里就是默认的3.12.3；接着选择合适的环境路径，这个就因人而异；然后无需勾选下面的两个选项，点击”确定“创建新环境。例如下图，我在PyCharm默认的项目目录下新创建了一个以Python 3.12.3为基础，名字是.fairchem的新环境：
(, 下载次数 Times of downloads: 1)

        接下来我们需要将当前的工作环境切换到新的虚拟环境，具体做法是：唤起终端(Alt + F12/Alt + Fn + F12)，输入：
  1. .\your_env_name\Scripts\activate
复制代码
       以我的.fairchem环境名为例，应该就是输入：
  1. .\.fairchem\Scripts\activate
复制代码
       这样，你的终端样式就应该会从之前的
  1. (.venv) D:\PyCharm Community\PyCharm Projects\pythonProject
复制代码
       变成
  1. (.fairchem) PS D:\PyCharm Community\PyCharm Projects\pythonProject>
复制代码
       接着可以尝试重启一下你的PyCharm，应该可以看到终端会显示
  1. (.fairchem) D:\PyCharm Community\PyCharm Projects\pythonProject
复制代码
       这样就说明已经成功进入了新的、干净的虚拟环境。如果你此时在终端中输入命令
  1. pip list
复制代码
       你应该会发现，目前只安装了pip这一个库。

2.2 安装Fairchem-Core
        这一步相当简单，唤起终端(Alt + F12/Alt + Fn + F12)后，在终端中输入
  1. pip install fairchem-core
复制代码
       所有相关的包都会自动被安装，并且由于我们是在全新的虚拟环境中进行这一步，理论上极少会有兼容性问题，只需要耐心等待即可！在安装过程中可能PyCharm可能会询问你是否要安装以下的包，此时需要输入‘y’来进行确认。安装完成后，可以继续使用pip list命令在终端中查看各个包的版本信息，在此列举一些以供参考：
  1. ...
  2. ase                     3.26.0
  3. fairchem-core           2.5.0
  4. networkx                3.5
  5. numpy                   2.2.6
  6. scipy                   1.16.2
  7. torch                   2.6.0
  8. ...
复制代码

2.3 UMA模型权重文件的获取
        在LICENSE · facebook/UMA at main中，Meta明确地通过直接法律限制(要求遵守美国的Trade Control Laws，明确禁止受美国制裁的国家、地区和个人使用该模型)与间接程序限制(美国加州法院对该许可证下的产品具有专属管辖权)的方式对特定地区进行限制。如果你认为你或者你所在的组织受到美国的制裁，或者你担心在美国司法体系下解决潜在纠纷的可行性，那么该许可证会很大程度上限制你使用这个模型。尽管在License部分有“额外说明“:UMA is available via HuggingFace globally, except in...，但是较为幽默的是，如果你把这部分“额外内容”与许可证同时喂给Gemini或者Grok并询问两个条款是否存在潜在冲突，你会得到很有意思的答案...
        抛开小扎手下并不FAIR的FAIRCHEM团队不谈，我默认大家都有进入互联网的能力，前往facebook/UMA · Hugging Face，在页面中找到权重申请的入口，进入后会要求你填写姓名与学校。以我个人的经验而言，非国内的IP、常见的美欧人名结合IP所在地的高等学府，基本上申请是秒过的。申请成功之后，就可以下载权重文件了，个人认为当前最优价值的模型应该是uma-s-1.1。其体积小巧，仅有1.1GB大小，可以非常顺畅地在近几年配置的个人电脑上部署并运行。结合0.0节中David_R的方法，可以通过修改fairchem-core的源代码来实现直接运行调用模型。同时，高技术力的读者也可以尝试手搓PyTorch代码来调用该模型权重以实现相同的功能。
         在Windows中结合PyCharm，我们也可以做到只需模型下载权限，而无需手动下载权重文件的方式来调用该模型，这可能也是较为推荐的一种做法，因为即使出现了问题也相比于修改源代码或自己手搓torch代码有着更快的订正速度: 首先来到Hugging Face → 点击右上角自己的用户头像 → Access Tokens → 右上方 Create New token，便可以来到创建新的Access Token的页面。我们可以在Token Type中选择Fine-Grained，Token Name可以填写自己喜欢的，例如我在此处填写UMA，下方Repositories选项中请务必把三个选项都勾选上(如下图所示)，其余选项没有必要打勾，拖动页面至最下方，点击Create token后复制新的token到自己电脑上的笔记本中妥善保存。
(, 下载次数 Times of downloads: 1)
        接下来来到PyCharm中，打开设置(Ctrl + Alt + S) → 外观与行为 → 系统设置 → HTTP，点击手动配置，选择HTTP，主机名与端口号的填写因软件而异。以小猫咪为例，默认的地址是127.0.0.1:7890，那就将地址与端口分别填入主机名与端口之中，下方”检查链接“可以使用常用谷歌学术来测试，如果显示”连接成功“，那么配置就是正确的。
        最后唤起终端(Alt + F12/Alt + Fn + F12)，输入
  1. hf auth login
复制代码
       如果你上述的配置都是正确的且网络环境通畅，那么会显示如下内容：
  1.     _|    _|  _|    _|    _|_|_|    _|_|_|  _|_|_|  _|      _|    _|_|_|      _|_|_|_|    _|_|      _|_|_|  _|_|_|_|
  2.     _|    _|  _|    _|  _|        _|          _|    _|_|    _|  _|            _|        _|    _|  _|        _|
  3.     _|_|_|_|  _|    _|  _|  _|_|  _|  _|_|    _|    _|  _|  _|  _|  _|_|      _|_|_|    _|_|_|_|  _|        _|_|_|
  4.     _|    _|  _|    _|  _|    _|  _|    _|    _|    _|    _|_|  _|    _|      _|        _|    _|  _|        _|
  5.     _|    _|    _|_|      _|_|_|    _|_|_|  _|_|_|  _|      _|    _|_|_|      _|        _|    _|    _|_|_|  _|_|_|_|

  7.     To log in, `huggingface_hub` requires a token generated from https://huggingface.co/settings/tokens.
  8. Token can be pasted using 'Right-Click'.
  9. Enter your token (input will not be visible):
复制代码
       复制你的token，可以直接使用右键或者 Ctrl + V 粘贴进去，然后回车，等待几秒后你会看到
  1. Token is valid (permission: fineGrained).
  2. The token `UMA` has been saved to C:\Users\your_user_name\.cache\huggingface\stored_tokens
  3. Your token has been saved in your configured git credential helpers (manager).
  4. Your token has been saved to C:\Users\your_user_name\.cache\huggingface\token
  5. Login successful.
  6. The current active token is: `UMA`
复制代码
       此时Hugging Face已经成功验证了你的身份，PyCharm上所涉及模型的权重文件现在会在代码运行过程中自动从Hugging Face拉取，权重文件一般会保存在以下路径：
  1. C:\Users\your_user_name\.cache\fairchem\models--facebook--UMA\snapshots\random_string_of_40_digits\checkpoints
复制代码

3. 使用UMA与ASE对分子结构进行优化及代码思路
3.1 代码构建
3.1.1 定义FAIRChemCalculator
        可调用UMA权重的Fairchem包代码与ASE的兼容性相当好，其实只要设置一个FAIRChemCalculator就可以被ASE所调用了。一个FAIRChemCalculator本身的参数与需要传递的参数较多，但是我们在优化特定的物种时是需要一致的Calculator的，所以我们可以先来定义一个函数，确定好在优化物种结构时的参数并返回一个封装好各类参数的FAIRChemCalculator，之后进需要调用即可。我们可以这样来写：
  1. from typing import Literal
  2. from fairchem.core import FAIRChemCalculator, pretrained_mlip

  4. def calculator(device: Literal['cpu', 'cuda'] = 'cpu', task_name: str = 'omol') -> FAIRChemCalculator:
  5.         valid_devices = {'cpu', 'cuda'}
  6.         valid_tasks = {'oc20', 'omat', 'omol', 'odac', 'moc'}

  8.         if device not in valid_devices:
  9.                 raise ValueError(f'Choose one device in {valid_devices}.')

  11.         if task_name not in valid_tasks:
  12.                 raise ValueError(f'Choose one task in {valid_tasks}.')

  14.         predictor = pretrained_mlip.get_predict_unit(model_name='uma-s-1p1', device=device)
  15.         return FAIRChemCalculator(predictor, task_name=task_name)
复制代码
       上述的函数返回一个FAIRChemCalculator，其中确定了运行模型的设备、任务类型、模型传递方法以及模型名称。设备类型会被严格限制在CPU与CUDA之间，由于我们现在在没有NVIDIA GPU的笔记本电脑上运行，因此默认设备类型为CPU；由于我们现在针对的物种是金属络合物，其对应的任务类型应该是omol，官方给出的所有可用的任务类型是：
  1. '''
  2. Set the task for your application and calculation
  3.         oc20 - use this for catalysis
  4.         omat - use this for inorganic materials
  5.         omol - use this for molecules
  6.         odac - use this for MOFs
  7.         omc - use this for molecular crystals
  8. '''
复制代码
       为了提升程序健壮性，在输入除去有效设备类型与有效任务类型之外的设备或任务都会抛出ValueError，请务必选择合适的选项。

3.1.2 结合ASE构建结构优化函数        
        现在我们需要定义如何使用ASE和FAIRChemCalculator来对结构进行优化。很明显，我们需要传入函数的数据有原始结构、Calculator本身，而函数会传出一个优化好的结构。同时，我们也可以设定一些ASE结构优化中的所需要的默认参数，如结构的电荷、自旋多重度、力收敛阈值、最大运行步数、优化的轨迹文件等。整体的参数设置与参数传递都很好理解：
  1. from ase import Atoms
  2. from ase.optimize import LBFGS
  3. from ase.io import write, Trajectory

  5. def geo_opt(reactant: Atoms, calc: FAIRChemCalculator, traj_name: str = '',
  6.              charge: int = 0, spin: int = 1, fmax: float = 0.025, maxsteps: int = 1000, logfile: str = '-') -> Atoms:

  8.         print('==========START REACTANT OPTIMIZATION==========')

  10.         reactant.calc = calc
  11.         reactant.info = {'charge': charge, 'spin': spin}
  12.         traj_filename = 'opt.traj'
  13.         traj = Trajectory(traj_filename, 'w', reactant)

  15.         opt_react = LBFGS(reactant, alpha=50, logfile=logfile, trajectory=traj)
  16.         opt_react.run(fmax=fmax, steps=maxsteps)

  18.         print('==========REACTANT OPTIMIZATION DONE==========')

  20.         print('-----WRITING OPTIMIZATION TRAJECTORY-----')

  22.         frames = read(traj_filename, index=':')
  23.         write(f'{traj_name}_traj.xyz', frames)

  25.         print('-----TRAJECTORY DONE-----')

  27.         return reactant.copy()
复制代码
       对于ASE的结构优化方法而言，LBFGS/BFGS都是可以考虑的，读者也可以试试看ASE提供的其他方法对于自己的体系是否更加高效；LBFGS/BFGS中的alpha确定了Hessian初猜，其默认的值是70，读者可以根据体系自行调整；当logfile被设置为logfile=‘-’时，ASE会在运行中默认输出LBFGS优化过程的各项参数，包括优化方法、步数、时间、体系能量以及结构受力。代码还将默认输出的traj轨迹格式转化成xyz轨迹格式，以方便观察优化轨迹。

3.1.3 定义环境变量        
        在David_R的博文中指出，使用GPU运行UMA时ASE仍倾向于使用掉每一个能用的CPU核心，因此他把CPU线程数设定为了1。我个人认为以纯CPU运行模型时，可以调整的CPU使用率是一个更好的选择，因此可以在代码的开头这样写以实现自定义CPU线程数：
  1. import os

  3. num_cpu_threads = '12'
  4. os.environ['OMP_NUM_THREADS'] = num_cpu_threads
  5. os.environ['MKL_NUM_THREADS'] = num_cpu_threads
  6. os.environ['OPENBLAS_NUM_THREADS'] = num_cpu_threads
  7. os.environ['NUMEXPR_NUM_THREADS'] = num_cpu_threads
复制代码

3.2 结构优化实例
        在上述代码完成之后，我们就差不多可以开始真正优化结构了。这里只需要四行代码就可以进入优化阶段，代码构建也符合直觉：确定要优化结构的路径 → ASE读取结构 → 结构传入Calculator → Calculator传入ASE进行优化。
  1. from ase.io import read
  2. from pathlib import Path

  4. react = r'absolute_path\your_structure'
  5. react = read(react, format=Path(react).suffix[1:])
  6. calculation = calculator()
  7. opted_react = geo_opt(reactant=react, traj_name='react', calc=calculation)
复制代码
        这里使用suffix方法来自动获取传入结构的格式，可以是xyz、mol等这类ASE可以直接读取的结构文件格式。在运行上述代码之后，我们会得到类似的输出：
  1. W1003 12:04:43.975000 13540
  2. .uma\Lib\site-packages\torch\distributed\
  3. elastic\multiprocessing\redirects.py:29]
  4. NOTE: Redirects are currently not supported in Windows or MacOs.
  5.        Step     Time          Energy          fmax
  6. LBFGS:    0 12:04:49   -65865.782109        3.082992
  7. LBFGS:    1 12:04:51   -65866.459520        2.014846
  8. LBFGS:    2 12:04:52   -65866.913820        1.120183
  9. LBFGS:    3 12:04:53   -65867.337472        1.540557
  10. LBFGS:    4 12:04:54   -65867.517059        1.157642
  11. LBFGS:    5 12:04:55   -65867.694608        0.802547
  12. LBFGS:    6 12:04:56   -65867.895905        0.586043
  13. LBFGS:    7 12:04:57   -65867.979146        0.612534
  14. LBFGS:    8 12:04:58   -65868.075387        0.627200
  15. LBFGS:    9 12:04:59   -65868.150023        0.360363
  16. ...
  17. LBFGS:  120 12:07:07   -65869.415791        0.029772
  18. LBFGS:  121 12:07:08   -65869.416312        0.029220
  19. LBFGS:  122 12:07:09   -65869.416791        0.029055
  20. LBFGS:  123 12:07:10   -65869.417148        0.035234
  21. LBFGS:  124 12:07:11   -65869.417477        0.024485
复制代码
程序在到达收敛限后会自动停止并输出优化轨迹文件。
        我们以Adv. Synth. Catal. 2025, 00, e70103中一个较为复杂的Pd催化中间体IM4R为例，其存在Pd-C/Pd-P/Pd-N配位键，并且存在弱相互作用。我们分别使用PBE0-D3(BJ)/6-31G*/def2-SVP(For Pd)以及UMA力场对原始输入结构进行优化，Gaussian的输入文件为：
  1. %nproc=96
  2. %mem=360GB
  3. #p opt PBE1PBE/genecp em=gd3bj

  5. Title

  7. 0 1
  8. P                  1.32242366    1.50787567    0.74995307
  9. O                  2.78177566    3.73559167    0.28834107
  10. O                  4.05917017   -1.48600229    1.64674044
  11. O                  2.29832463   -0.73980731   -2.70395439
  12. C                  1.58136666    1.32546267    2.63393507
  13. C                  2.98763166    1.39018267   -0.01752193
  14. C                  1.40387666    3.35999767    0.46490407
  15. C                  2.76126166    2.16624367    3.15243207
  16. C                  1.83381866   -0.16097933    2.94674407
  17. C                  0.26980066    1.77172567    3.31109507
  18. C                  3.54171066    2.67987767   -0.11895793
  19. C                  3.72807466    0.28263067   -0.47042593
  20. C                  3.17623266   -1.10732033   -0.49928493
  21. C                  4.83159266    2.89549267   -0.60289593
  22. C                  5.02728766    0.50108067   -0.95469193
  23. C                  5.57061366    1.78706467   -1.01037093
  24. C                  3.38126775   -2.00119988    0.57625495
  25. C                  2.50373882   -1.58763325   -1.63372149
  26. C                  2.90047328   -3.31349572    0.52066650
  27. C                  2.01715285   -2.89814168   -1.69747921
  28. C                  2.21657515   -3.74980427   -0.61717293
  29. C                  5.47441053   -1.59477971    1.66881515
  30. C                  3.34378337   -0.74990985   -3.67680312
  31. H                  0.99873066    3.94170767    1.29477207
  32. H                  2.85493966    2.02353567    4.23726107
  33. H                  2.62553166    3.23721467    2.97121607
  34. H                  3.70733066    1.86180567    2.69386907
  35. H                  1.89309666   -0.29467833    4.03542507
  36. H                  2.77696966   -0.50864533    2.51563907
  37. H                  1.02036766   -0.79570033    2.57908407
  38. H                 -0.57632334    1.15228867    2.99736407
  39. H                  0.37425666    1.67253367    4.39942607
  40. H                  0.02452966    2.81907267    3.10019507
  41. H                  5.22553466    3.90482667   -0.66040593
  42. H                  5.61057366   -0.34874333   -1.29548493
  43. H                  6.57773666    1.92766267   -1.39371793
  44. H                  3.04652857   -3.99359913    1.35120781
  45. H                  1.48211295   -3.22145553   -2.58410936
  46. H                  1.83667076   -4.76700351   -0.65199913
  47. Pd                -0.57244323    0.33446891   -0.14032557
  48. H                 -3.41675243    5.43032725   -1.03473335
  49. C                 -2.88911583    4.49594407   -0.86149479
  50. C                 -2.06950861    3.95479797   -1.85354462
  51. H                 -1.95558832    4.46598243   -2.80722446
  52. C                 -1.39000807    2.75072270   -1.63134292
  53. H                 -0.76089047    2.34261472   -2.41883397
  54. C                 -1.52279570    2.07351192   -0.40969447
  55. C                 -2.35153077    2.61773461    0.58131707
  56. H                 -2.48738771    2.10320882    1.52776037
  57. C                 -3.02909012    3.82155814    0.35285045
  58. H                 -3.67206195    4.22799569    1.13083934
  59. H                 -3.48397582   -1.37386042   -4.15827632
  60. C                 -3.34925694   -1.70528917   -3.12063550
  61. H                 -2.72104903   -2.60579780   -3.14234439
  62. H                 -4.99030070   -3.06855181   -2.68522117
  63. H                 -1.32202290   -1.39276653   -2.95357095
  64. H                 -5.48661017   -1.39744104   -2.88467123
  65. C                 -4.69398900   -2.03466600   -2.45831400
  66. C                 -2.29096943   -1.32258954   -2.44800610
  67. C                 -4.57499500   -1.84542900   -0.96913400
  68. N                 -2.46731012   -0.83893784   -1.03060420
  69. C                 -3.42105100   -1.44183500   -0.36813300
  70. H                 -5.34777802   -2.17777331    1.39309783
  71. C                 -3.34981200   -1.43306000    1.12429100
  72. C                 -4.43488091   -1.87407323    1.89736338
  73. C                 -2.17834984   -1.03146914    1.79028724
  74. H                 -1.32266932   -0.71067246    1.20344604
  75. C                 -4.35598894   -1.91214766    3.29153683
  76. H                 -5.21156086   -2.25218295    3.87008819
  77. C                 -2.09637519   -1.07556982    3.18261224
  78. C                 -3.18669495   -1.51252864    3.94078147
  79. H                 -1.17463909   -0.77822841    3.67736879
  80. H                 -3.12295851   -1.54423253    5.02551332
  81. H                 -2.60563866   -0.26722007   -2.49539031
  82. H                 -5.40392600   -2.15630400   -0.33975900
  83. H                  4.28572253   -0.38725253   -3.24570877
  84. H                  5.79867720   -2.64325267    1.69008799
  85. C                  0.61800319    3.70729508   -0.81313638
  86. H                 -0.06720711    2.91640850   -1.03648646
  87. C                  1.59904198    3.88336446   -1.98708803
  88. H                  2.50662635    4.32382853   -1.63046215
  89. H                  1.81215950    2.92838603   -2.42009741
  90. H                  1.15997062    4.52021271   -2.72637190
  91. C                 -0.16271818    5.01680833   -0.59576757
  92. H                  0.15107380    5.47120080    0.32072688
  93. H                  0.02853298    5.68527674   -1.40907746
  94. H                 -1.21021038    4.80380364   -0.54792322
  95. C                  5.92760588   -0.89525092    2.96375663
  96. H                  6.98496944   -0.73458253    2.93107698
  97. H                  5.42650814    0.04575006    3.05496080
  98. H                  5.68622239   -1.51088430    3.80496383
  99. C                  6.11267635   -0.89022261    0.45728056
  100. H                  7.00823536   -0.39223889    0.76531675
  101. H                  6.34796086   -1.61516822   -0.29371549
  102. H                  5.42530302   -0.17373097    0.05846026
  103. C                  3.55885616   -2.16238016   -4.25149549
  104. H                  4.05356346   -2.09098958   -5.19757586
  105. H                  4.16077084   -2.73422513   -3.57651990
  106. H                  2.61166774   -2.64331709   -4.37968980
  107. C                  2.89628063    0.19921057   -4.80397289
  108. H                  1.92637780   -0.09085816   -5.15045323
  109. H                  2.85697231    1.20156903   -4.43163919
  110. H                  3.59456473    0.14637437   -5.61298937

  112. C H O P N 0
  113. 6-31G*
  114. ****
  115. Pd 0
  116. def2svp
  117. ****

  119. Pd 0
  120. def2svp
复制代码
      将gjf文件转换为xyz文件作为ASE输入，UMA的收敛标准设定为：
  1. fmax = 0.0025
  2. maxstep = 2000
复制代码
       Gaussian用时76m41s 73步收敛，UMA用时9m21s 450步收敛。对于这个103个原子的体系，在纯CPU上运行，UMA平均一步耗时1.25s，占用内存约1.7GB。观察一下用两种不同方法优化得到的结构，以及其中的细节：
(, 下载次数 Times of downloads: 2)
        可以看到，作为一类机器学习力场，UMA对于配位键键长的描述已经近乎达到了DFT级别，在不做构象搜索的情况下，配体上的若干个甲基的朝向几乎也是相同的。并且对于Pd-O之间的静电吸引也作出了与DFT相似的表述。相关的输入输出文件在附件中，有兴趣的读者也可以度量一下各个配位键以及Pd-O的键长，UMA相比DFT的误差甚至可以做到0.01埃左右，在这个结构中毫无疑问地有着十分惊人的表现。

4.使用UMA与ASE创建NEB对分子内C-C偶联反应进行过渡态初猜搜索及代码思路
4.1使用UMA与ASE构建NEB
4.1.1使用ASE构建NEB       
        UMA作为一种力场在金属配合物结构预优化中的表现堪称顶级，如果该方法有着相当的精度，使用它来帮助我们完成金属配合物反应当中的过渡态优化是否可行？NEB方法通常被用于寻找能量最小路径，也被我们用来寻找一些反应的过渡态结构初猜。在仅牵涉到分子反应的过渡态搜寻过程中，NEB可能是opt=TS/QST或者刚性/柔性扫描都失效后的最终手段。由于当NEB images的数量增加时，优化耗时的增加会非常恐怖，因此NEB并非是首选。而有了UMA之后，能量计算和优化的成本是大幅降低的，所以对于某些过渡态搜寻，我们可以直接大大方方地上NEB。
        首先定义一个函数来专门处理NEB问题，这个函数接受原子序号对应的反应物结构与生成物结构、NEB image的数量，并通过NEB优化返回一条NEB路径，返回类型是由结构构成的列表；同时，该函数也要接受FAIRChemCalculator的传入，以及电荷、自旋多重度、力收敛阈值和最大优化步数，因此代码可以这样写：
  1. # noinspection PyTypeChecker
  2. def neb_opt(head: Atoms, tail: Atoms, calc_factory,charge: int = 0, spin: int = 1,
  3.             fmax: float = 0.025, maxsteps: int = 2000, num_images: int = 5, logfile: str = '-') -> list[Atoms]:

  5.         print('==========START NEB PATH OPTIMIZATION==========')
  6.         print(f'=========={num_images + 2} IMAGES IN TOTAL==========')

  8.         images = [head]
  9.         for i in range(num_images):
  10.                 image = head.copy()
  11.                 image.info.update({'charge': charge, 'spin': spin})
  12.                 image.calc = calc_factory()
  13.                 images.append(image)
  14.         images.append(tail)

  16.         neb = NEB(images)
  17.         neb.interpolate()

  19.         neb_traj_filename = 'neb.traj'
  20.         neb_opt = LBFGS(neb, logfile=logfile, trajectory=neb_traj_filename)
  21.         neb_opt.run(fmax=fmax, steps=maxsteps)

  23.         print('==========NEB PATH DONE==========')

  25.         print('-----WRITING OPTIMIZATION TRAJECTORY-----')

  27.         frames = read(neb_traj_filename, index=':')
  28.         write(f'neb_traj.xyz', frames)
  29.         write('neb_final_path.xyz', images)

  31.         print('-----TRAJECTORY DONE-----')

  33.         return images
复制代码
       我们将所有需要优化的结构存在一个列表当中，通过复制若干的反应物结构(head)来填充这个列表，并在最后将生成物结构(tail)也加入这个列表。列表中结构的个数取决于你想要的结构个数与NEB路径的平滑程度二者的平衡：想要算的快，就减少NEB images的数量，但是得到的路径肯定相当不平滑；想要清晰光滑的反应路径，就需要很多NEB images，但是不可避免运算量与耗时会增加。整个列表中一个会有{num_images + 2}个结构，多出来的2个结构分别是一头一尾的反应物与生成物。
        接着使用这一个结构列表来创建NEB并对其中的结构进行interpolate，ASE默认的方法是linear interpolate，如有需要可以自行修改，同时也可以通过在NEB()中添加climb=True参数来实现CI-NEB。
        然后确定一下NEB优化轨迹的文件路径，并使用LBFGS来对NEB路径进行优化。这里的代码存在一个小问题是，直接将结构列表传入LBFGS优化器中会引发PyCharm的类型警告：LBFGS只能接受Atoms类型而不接受list[Atoms]类型。但是在我的实际操作中发现除了logfile中的一些小问题以外(在下文5.2节中叙述)，传入NEB的list[Atoms]也可以得到正确的NEB路径，因此为了消除这个无用的弱警告，需要在该函数前一行注释
  1. # noinspection PyTypeChecker
复制代码
以跳过该函数内部的类型检查。

4.1.2 NEB的能量分析
        在完成优化函数、NEB函数之后，我们还需要一个函数来对NEB的结果进行分析，具体而言是分析NEB路径上每个image的相对能量大小以确定鞍点的位置。结合numpy、matplotlib和ASE中的单位转换功能也可以很方便地实现这一点：传入最后一次优化得到的NEB路径，返回一张相对能量 - 反应路径图。代码参考了AI4Chem仓库中的内容：
  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. from ase.units import kcal, mol

  5. def neb_plot(images: list[Atoms], filename: str ='neb_energy_profile.png') -> None:
  6.         energies = np.array([atoms.get_total_energy() / (kcal / mol) for atoms in images])
  7.         energies = energies - energies[0]
  8.         num_images = len(images)

  10.         distances = []
  11.         cumulative_distance = 0.0
  12.         if num_images > 1:
  13.                 for i in range(num_images - 1):
  14.                         cumulative_distance =  (cumulative_distance + np.linalg.norm(images[i + 1].get_positions() - images[i].get_positions()))
  15.                         distances.append(cumulative_distance)
  16.                 distances.insert(0, 0.0)
  17.         else:
  18.                 distances = [0.0]

  20.         plt.figure(figsize=(10, 6))
  21.         plt.plot(distances, energies, marker='o', linestyle='-')
  22.         plt.xlabel('Distance along Reaction Path (Å)')
  23.         plt.ylabel('Relative Energy (kcal/mol)')
  24.         plt.title('NEB Energy Profile')
  25.         plt.grid(True)
  26.         plt.tight_layout()
  27.         plt.savefig(filename, dpi=300)
  28.         plt.close()
复制代码
       NEB路径中所有的images都会与反应物进行能量的比较，从而得到整个路径上每一个image的相对能量，最后输出一张能量曲线图，方便我们观察能量变化趋势。

4.2 结合UMA的NEB示例                        
        现在我们的路线也是明确的，我们有一个反应物结构与一个生成物结构，便可以创建一个NEB并优化得到可能的过渡态结构。我们仍然使用Adv. Synth. Catal. 2025, 00, e70103中的例子，用NEB来搜寻文中TS6R过渡态结构。反应物IM4R已在上文中做过叙述，而生成物FS2在原文中已经是从Pd物种上脱离下来的状态，我们可以将FS2的N原子配位至与IM4R相同的位置上来表示其未从催化物种上脱离的状态。为了保证反应路径中所有原子标号一致，可以在IM4R的结构基础上进行修改。现在，我们就有了一个已经优化好的反应物IM4R和一个还未进行优化的生成物FS2X，在这两个物种上构建NEB。
        在上文中我们已经拥有了Calculator函数、结构优化函数、NEB优化函数以及NEB能量分析函数，只要结合反应物与生成物结构的具体路径，我们就可以有完整的基于UMA的NEB代码了，下面将展示项目的完整代码结构，其中的各个函数在上文已有展示：
  1. import os
  2. from pathlib import Path
  3. from typing import Literal

  5. num_cpu_threads = '12'
  6. os.environ['OMP_NUM_THREADS'] = num_cpu_threads
  7. os.environ['MKL_NUM_THREADS'] = num_cpu_threads
  8. os.environ['OPENBLAS_NUM_THREADS'] = num_cpu_threads
  9. os.environ['NUMEXPR_NUM_THREADS'] = num_cpu_threads

  11. from ase import Atoms
  12. from ase.mep import NEB
  13. from ase.units import kcal, mol
  14. from ase.optimize import LBFGS
  15. from ase.io import read, write, Trajectory
  16. from fairchem.core import FAIRChemCalculator, pretrained_mlip

  18. import numpy as np
  19. import matplotlib.pyplot as plt

  21. def calculator(device: Literal['cpu', 'cuda'] = 'cpu', task_name: str = 'omol') -> FAIRChemCalculator:
  22.     ...

  24. def geo_opt(reactant: Atoms, calc: FAIRChemCalculator, traj_name: str = '',
  25.              charge: int = 0, spin: int = 1, fmax: float = 0.0025, maxsteps: int = 1000, logfile: str = '-') -> Atoms:   
  26.     ...

  28. # noinspection PyTypeChecker
  29. def neb_opt(head: Atoms, tail: Atoms, calc_factory,charge: int = 0, spin: int = 1,
  30.             fmax: float = 0.2 , maxsteps: int = 2000, num_images: int = 5, logfile: str = '-') -> list[Atoms]:
  31.     ...

  33. def neb_plot(images: list[Atoms], filename: str ='neb_energy_profile.png') -> None:
  34.     ...

  36. if __name__ == '__main__':
  37.         react = r'absolute_path\your_reactant'
  38.         react = read(react, format=Path(react).suffix[1:])

  40.         prod = r'absolute_path\your_product'  
  41.         prod = read(prod, format=Path(prod).suffix[1:])

  43.         calculation_1 = calculator()
  44.         calculation_2 = calculator()

  46.         opted_react = geo_opt(reactant=react, traj_name='react', calc=calculation_1)
  47.         opted_prod = geo_opt(reactant=prod, traj_name='prod', calc=calculation_2)

  49.         neb_images = neb_opt(react, prod, calculator, num_images=18)
  50.         neb_plot(neb_images)
复制代码
        在主程序中可以调用上述的四个函数，结构较为简单易懂。首先载入反应物与生成物的具体路径并让ASE读取分子结构；接着创建两个FAIRChemCalculator在之后的计算中传递给反应物与生成物；然后分别对二者进行几何优化得到合适的结构；最后针对二者进行NEB的优化并得到NEB的能量分析结果，寻找到潜在的过渡态结构。我在此例中插入18个NEB image，设置了很宽松的NEB优化力收敛阈值以节省时间，且目前的反应物、生成物结构已经提前优化。运行上述代码，可以得到类似如下内容：
  1. ==========START REACTANT OPTIMIZATION==========
  2.        Step     Time          Energy          fmax
  3. LBFGS:    0 22:33:45   -65869.475626        0.002484
  4. ==========REACTANT OPTIMIZATION DONE==========
  5. -----WRITING OPTIMIZATION TRAJECTORY-----
  6. -----TRAJECTORY DONE-----
  7. ==========START REACTANT OPTIMIZATION==========
  8.        Step     Time          Energy          fmax
  9. LBFGS:    0 22:33:47   -65869.911009        0.002391
  10. ==========REACTANT OPTIMIZATION DONE==========
  11. -----WRITING OPTIMIZATION TRAJECTORY-----
  12. -----TRAJECTORY DONE-----
  13. ==========START NEB PATH OPTIMIZATION==========
  14. ==========20 IMAGES IN TOTAL==========
  15.        Step     Time          Energy          fmax
  16. LBFGS:    0 22:46:37   -65863.632958       16.867927
  17. LBFGS:    1 22:46:59   -65866.423111       11.438179
  18. LBFGS:    2 22:47:21   -65867.640214        4.062935
  19. LBFGS:    3 22:47:43   -65868.055238        1.875868
  20. LBFGS:    4 22:48:05   -65868.179259        1.412737
  21. LBFGS:    5 22:48:28   -65868.261376        1.433932
  22. LBFGS:    6 22:48:50   -65868.318152        1.038011
  23. LBFGS:    7 22:49:12   -65868.352802        0.816594
  24. LBFGS:    8 22:49:34   -65868.377574        1.009610
  25. LBFGS:    9 22:49:55   -65868.402276        0.724988
  26. LBFGS:   10 22:50:17   -65868.419051        0.525514
  27. LBFGS:   11 22:50:39   -65868.428215        0.477047
  28. LBFGS:   12 22:51:01   -65868.436577        0.599866
  29. LBFGS:   13 22:51:23   -65868.445362        0.426750
  30. LBFGS:   14 22:51:45   -65868.455088        0.372402
  31. LBFGS:   15 22:52:07   -65868.465119        0.352490
  32. LBFGS:   16 22:52:29   -65868.472281        0.374115
  33. LBFGS:   17 22:52:50   -65868.477368        0.351974
  34. LBFGS:   18 22:53:12   -65868.482401        0.344732
  35. LBFGS:   19 22:53:34   -65868.488276        0.367053
  36. LBFGS:   20 22:53:56   -65868.494220        0.352513
  37. LBFGS:   21 22:54:18   -65868.498228        0.328189
  38. LBFGS:   22 22:54:39   -65868.501957        0.327379
  39. LBFGS:   23 22:55:01   -65868.507301        0.318098
  40. LBFGS:   24 22:55:22   -65868.514008        0.316225
  41. LBFGS:   25 22:55:44   -65868.519416        0.221761
  42. LBFGS:   26 22:56:05   -65868.523299        0.269194
  43. LBFGS:   27 22:56:27   -65868.526838        0.221552
  44. LBFGS:   28 22:56:48   -65868.529811        0.169127
  45. ==========NEB PATH DONE==========
  46. -----WRITING OPTIMIZATION TRAJECTORY-----
  47. -----TRAJECTORY DONE-----
复制代码
       对于较大的收敛阈值，NEB进需要28步就完成了优化，加载20个结构(插入18个image)共使用了约14.2GB内存，平均每个image需要约800MB内存。运行完毕后得到结果分析：
(, 下载次数 Times of downloads: 1)
        打开随NEB运行输出的‘neb_final_path.xyz’轨迹文件，寻找最高点所在的第11帧结构，将其与手动摆放过渡态结构、并用PBE0-D3(BJ)/6-31G*/def2-SVP水平优化出的过渡态结构进行一下比较：
(, 下载次数 Times of downloads: 1)
        可以看出仅在分子内C-C键偶联的过渡态搜寻上，UMA结合NEB在fmax=0.2这样极为粗糙的收敛条件下，只需要约10分钟就可以给出一个相当合理的过渡态初猜结构了。尽管对于这类反应而言可能根本无需使用DFT级别的NEB来寻找过渡态，但是就耗时而言，其仅需DFT级别的NEB耗时零头里的零头了，更何况目前的模型根本没有运行在GPU+VRAM上，而是在慢一个数量级的CPU+RAM当中。可以想象，假如模型在合适的设备上以更低的收敛阈值运行NEB等方法时，极有可能以更短的时间、更高的精度给出过渡态结构初猜，非常有益于提高后续DFT过渡态优化的效率。

5.在探索中发现的问题与讨论
5.1 UMA是在“背题”吗？               
        在omol任务中调用UMA，其对应的训练集是OMol25训练集。OMol25数据集与UMA模型一样，都存在着针对地域的获取限制。况且OMol25数据集相当庞大，我目前还没有研究过这个数据集中是否包含着本文当中的结构优化示例与NEB示例数据。如果在之后的深入探索中发现了这个数据集中存有示例数据，那么本文当中的示例可能是不可靠的，这也从侧面反映我们需要更多“纯净”的数据来测试UMA的稳健性。
        上述的情况我在选择示例时其实也有考虑到，因此特意选择了一篇非常新鲜的文章中的结构来尝试规避这个问题，这也如各位所见，如果模型不是在“背题”，那么说明这个模型具有相当强大的泛化能力，UMA力场的精度也可以说是达到一个新的台阶。

5.2 NEB中的LBFGS logfile
        对于结构优化而言，LBFGS给出的logfile是非常清晰易懂的，其中的Energy和fmax项分别表示了分子的能量与受力。但是当使用LBFGS优化NEB时，其输出就会有一些令人疑惑了，如4.2节中的logfile：
  1.        Step     Time          Energy          fmax
  2. LBFGS:    0 22:46:37   -65863.632958       16.867927
  3. LBFGS:    1 22:46:59   -65866.423111       11.438179
  4. LBFGS:    2 22:47:21   -65867.640214        4.062935
  5. LBFGS:    3 22:47:43   -65868.055238        1.875868
  6. LBFGS:    4 22:48:05   -65868.179259        1.412737
  7. LBFGS:    5 22:48:28   -65868.261376        1.433932
  8. ...
  9. LBFGS:   23 22:55:01   -65868.507301        0.318098
  10. LBFGS:   24 22:55:22   -65868.514008        0.316225
  11. LBFGS:   25 22:55:44   -65868.519416        0.221761
  12. LBFGS:   26 22:56:05   -65868.523299        0.269194
  13. LBFGS:   27 22:56:27   -65868.526838        0.221552
  14. LBFGS:   28 22:56:48   -65868.529811        0.169127
复制代码
从每一轮的优化耗时来看，不难看出此时优化器正在同时优化band上的所有结构。而此处的Energy和fmax项就会有一些理解上的疑问，这两项分别代表了什么结构的能量和受力？
        众所周知，NEB方法通过优化一条由imgaes构成的band，使得每个结构所受的真实原子力在反应路径切线方向上的分量为零。当NEB收敛时，整条反应路径就近似为最小能量路径，其上每个点都只在垂直于路径的方向上受力，而沿着反应路径方向的净受力为零。因此，这里的fmax所代表的受力可能是：某个image在优化过程中的受力、所有images在优化过程中的平均受力、优化过程中band的总受力。但是结合Energy项的数量级，显然不会是整个band上所有image的能量加和。所以这里的两项究竟是代表了某一个结构的数据，还是所有结构平均之后的数据，希望有精通ASE的老师同学帮助解惑，并针对NEB中fmax项的合理设置给出见解！

5.3 模型对于计算资源的消耗情况
        从我的测试结果中可以看出，尽管都使用LBFGS方法，结构优化与NEB优化中独立结构所消耗的内存并不是一样的。对于32GB内存的笔记本电脑而言，做一次约100个原子的NEB，假如不清空后台，极限约是25个结构即23个image，最多不超过28个image。其计算资源消耗量可供各位参考。

5.4 UMA在双分子反应过渡态初猜搜寻中的表现               
        虽然上文中的分子内C-C偶联过渡态初猜搜寻中，UMA表现迅猛，但是经过我的一例测试发现，其在双分子过渡态的表现存在明显的问题，我使用的测试例子是最简单的1,3-丁二烯与乙烯的Diels-Alder加成反应，对其进行反应物、生成物优化以及过渡态初猜搜索，可以得到如下的数据：
(, 下载次数 Times of downloads: 1)
        观察能量最高点过渡态初猜结构：
(, 下载次数 Times of downloads: 1)
        初看这个过渡态初猜感觉有一丝合理性，但是回想一下基础的有机知识就可以发现[4+2]的Diels-Alder加成反应并非是在同一个平面内发生的(附件中对应的结构，读者可以体系侧面观察扭转的丁二烯与乙烯的位置，二者在主轴线位置几乎被拉到一个平面内了)，而是非共平面的，如下图在B3LYP-D3(BJ)/6-311G*水平下优化的[4+2]加成过渡态：
(, 下载次数 Times of downloads: 0)
        以二者自身的主轴线而言，两个分子是不在同一平面上的。UMA在这一例简单的双分子过渡态初猜优化中表现不佳，出现了一些定性错误。
        对于此情况，我个人分析的原因是，OMol25训练集当中的数据有一个隐含前提，就是优化好的结构是“一个分子”，UMA可能没有办法从都是“一个分子”的数据中泛化出“一个丁二烯与一个乙烯在靠近时二者的位置关系”，因为后者完完全全是“两个分子”，不包含任何的键连关系，所以UMA在这个体系中失效了而DFT的描述是准确的。在这个体系中，UMA似乎在尽全力将其视为“一个分子”，因此NEB路径中的最高点结构也已经完全倾向于生成环己烯的结构了。

6.总结与展望
        总的来说，这次探索让我个人觉得UMA模型在金属配合物反应当中的结构预优化和过渡态初猜方面具有相当大的开发潜力。即使在普通的笔记本电脑上，也能高效地对复杂的金属有机分子进行结构优化，给出的键长甚至能和DFT结果媲美，而且表现得比一些传统半经验方法更稳健。UMA结合NEB搜索一个分子内反应的过渡态，在十分钟内就能搜寻到相当靠谱的初猜结构，其效率令人印象深刻。当然，该模型可能由于训练集的限制，会出现一些定性的错误，例如在处理像双分子Diels-Alder反应这种需要精确把握两个分子空间取向的过渡态时，就显得有些力不从心。
        目前来看，该模型可以通过精确的预优化极大地节省DFT优化的时间，这是其最大的优势。从代码方面来说，目前迫切的需求是结合ASE的Trajectory或者某些Checkpoint机制，开发出健壮的“断点续算”功能，可以我们更加放心地去进行较长时间的NEB初筛计算；或者是使用该机器学习势，开发出类似于CREST一样的构象搜索系统，在某些极端的体系中可能会比GFN-FF表现更好。

(, 下载次数 Times of downloads: 62)




作者
Author: Daniel_Arndt    时间: 2025-10-20 10:27
楼主找出来的Diels–Alder反应的第一个过渡态对应于J. Org. Chem. 2000, 65(21), 7134–7138的Scheme 1中的TS3，楼主找出来的Diels–Alder反应的第二个过渡态对应于J. Org. Chem. 2000, 65(21), 7134–7138的Scheme 1中的TS6。
作者
Author: RAL    时间: 2025-10-20 10:38
D-A反应恐怕是初猜或者NEB的问题，我们组用ADCR结合UMA的“前身”，即eSEN架构在OMol25上训练的baseline模型，对181个教科书中性闭壳层/中性单自由基反应测试，都有定性+定量非常精确的描述。
贴图中找到的是典型的反式二烯体+亲二烯体的环加成过渡态。。
作者
Author: Whitedwarf    时间: 2025-10-23 11:08
请问PBE0-D3(BJ/6-31G*/def2-SVP和UMA(uma-s-1.1)/NEB(ASE-LBFGS)比较的可视化是用的什么程序，感觉标记非常清晰
作者
Author: slxc920113    时间: 2025-10-23 12:31
Whitedwarf 发表于 2025-10-23 11:08
请问PBE0-D3(BJ/6-31G*/def2-SVP和UMA(uma-s-1.1)/NEB(ASE-LBFGS)比较的可视化是用的什么程序，感觉标记非 ...

应该是VMD，可以参考这篇帖子：http://bbs.keinsci.com/thread-1080-1-1.html
作者
Author: adong    时间: 2025-10-23 16:44
Whitedwarf 发表于 2025-10-23 11:08
请问PBE0-D3(BJ/6-31G*/def2-SVP和UMA(uma-s-1.1)/NEB(ASE-LBFGS)比较的可视化是用的什么程序，感觉标记非 ...

CYLView，https://www.cylview.org/
作者
Author: Whitedwarf    时间: 2025-10-23 20:54
adong 发表于 2025-10-23 16:44
CYLView，https://www.cylview.org/

感谢
作者
Author: archer    时间: 2025-11-5 10:05
UMA的准确度和xtb比起来如何
作者
Author: LPOwenChan    时间: 2025-11-8 11:10
archer 发表于 2025-11-5 10:05
UMA的准确度和xtb比起来如何

       目前还没有Benchmark文章来比较最新的高精度机器学习力场(UMA等)、半经验(xTB等)与DFT三者之间的精度差异。不过就我自己研究的金属有机催化方面而言(主要牵涉Pd、Fe、Cu、Co与Ni)，UMA的预优化精度可以与DFT媲美了，远远高于GFN-FF与GFN2-xTB，主要的差异来源于配位键的长度。在我个人研究的一些体系中，GFN-FF与GFN2-xTB对于配位键长度的误判带来的最大问题就是配体解离与金属中心的逃逸(尤其是在一些平面四配位的体系中)，而UMA目前是很好地解决了这一痛点的。
作者
Author: maoxinxina    时间: 2025-11-15 19:46
可以把这个势文件下载好放在服务器上面用么。
作者
Author: LPOwenChan    时间: 2025-11-15 23:30
maoxinxina 发表于 2025-11-15 19:46
可以把这个势文件下载好放在服务器上面用么。

可以的，修改一下源代码就可以
作者
Author: maoxinxina    时间: 2025-11-17 09:47
LPOwenChan 发表于 2025-11-15 23:30
可以的，修改一下源代码就可以

能说一下大概怎么改么，谢谢。
作者
Author: LPOwenChan    时间: 2025-11-17 13:08
maoxinxina 发表于 2025-11-17 09:47
能说一下大概怎么改么，谢谢。

需要两个文件才能让模型在本地跑起来，一个是uma-s-1p1.pt，即模型权重文件本身；另一个是iso_atom_elem_refs.yaml，是原子索引文件（附于文末）。假设把他们放在以下两个路径：
  1. /public/home/username/work_path/model/iso_atom_elem_refs.yaml
  2. /public/home/username/work_path/model/uma-s-1p1.pt
复制代码
我的fairchem版本是2.7.1，假设fairchem-core核心代码文件夹的安装位置是：
  1. /public/home/username/miniconda3/envs/fairchem/lib/python3.12/site-packages/fairchem/core/calculate
复制代码
现将上述路径下的pretrained_mlip.py文件备份，接着将我改写好的pretrained_mlip.py覆盖掉原来的文件（附于文末），然后继续修改pretrained_mlip.py文件；
首先寻找到第52行，内容是：
  1. checkpoint_path = r""
复制代码
上下文是：
  1. def pretrained_checkpoint_path_from_name(model_name: str):
复制代码
进入编辑模式，将模型权重文件绝对路径填入双引号中：
  1. checkpoint_path = r"/public/home/username/work_path/model/uma-s-1p1.pt"
复制代码
接下来找到第91行，内容是：
  1. checkpoint_path = r""
复制代码
上下文是：
  1. try:
  2.         model_checkpoint = _MODEL_CKPTS.checkpoints[model_name]
  3. except KeyError as err:
  4.         raise KeyError(
  5.                 f"Model '{model_name}' not found. Available models: {available_models}"
  6.         ) from err
  7. checkpoint_path = r""
  8. atom_refs = get_isolated_atomic_energies(model_name, cache_dir)
  9. return load_predict_unit(
  10.         checkpoint_path, inference_settings, overrides, device, atom_refs
  11. )
复制代码
将模型权重文件绝对路径填入双引号中：
  1. checkpoint_path = r"/public/home/username/work_path/model/uma-s-1p1.pt"
复制代码
最后找到第112行：
  1. atomic_refs_path = r""
复制代码
上下文为：
  1. model_checkpoint = _MODEL_CKPTS.checkpoints[model_name]
  2. atomic_refs_path = r""
  3. return OmegaConf.load(atomic_refs_path)
复制代码
将原子索引文件绝对路径填入双引号中，注意这里填写的是原子索引文件路径而非模型权重路径
  1. atomic_refs_path = r"/public/home/username/work_path/model/iso_atom_elem_refs.yaml"
复制代码
保存并退出，即完成源代码的修改!
(, 下载次数 Times of downloads: 14) (, 下载次数 Times of downloads: 12)




作者
Author: maoxinxina    时间: 2025-11-17 15:57
LPOwenChan 发表于 2025-11-17 13:08
需要两个文件才能让模型在本地跑起来，一个是uma-s-1p1.pt，即模型权重文件本身；另一个是iso_atom_elem_ ...

试过了，真的可行，太厉害了。想问下要是那个uma-m-1p1.pt势，也是在这里面改吗。
作者
Author: LPOwenChan    时间: 2025-11-17 21:49
maoxinxina 发表于 2025-11-17 15:57
试过了，真的可行，太厉害了。想问下要是那个uma-m-1p1.pt势，也是在这里面改吗。

是的，要用什么别的模型应该就是写入对应模型权重文件的绝对路径应该就可以了。
作者
Author: maoxinxina    时间: 2025-11-18 09:42
LPOwenChan 发表于 2025-11-17 21:49
是的，要用什么别的模型应该就是写入对应模型权重文件的绝对路径应该就可以了。

谢谢。
作者
Author: quanta    时间: 2026-1-14 15:57
博主好，我重复了一下您文中的结果，但是得到的两个分子的差别很大，RMSD有1.0+Å了，不知道是怎么回事儿
(, 下载次数 Times of downloads: 0)
作者
Author: LPOwenChan    时间: 2026-1-14 18:31
quanta 发表于 2026-1-14 15:57
博主好，我重复了一下您文中的结果，但是得到的两个分子的差别很大，RMSD有1.0+Å了，不知道是怎么回事 ...

可能有很多种原因值得考虑的：
1. 初始的输入结构是否合理(以DFT最终优化结果为参考)？
2. 多构象（我个人认为你的结果可能就是有构象干扰，力场可能陷入了一个局部最小值，这种问题可能是优化器导致的）
3. 数值精度问题（CPU推理和GPU推理数值精度相差很大，但是这种可能性比较小）
作者
Author: quanta    时间: 2026-1-15 08:50
LPOwenChan 发表于 2026-1-14 18:31
可能有很多种原因值得考虑的：
1. 初始的输入结构是否合理(以DFT最终优化结果为参考)？
2. 多构象（我 ...

博主您好，我用UMA起始优化的结构就是您帖子里的那个结构，您用uma优化的时候，是拿DFT的优化结果开始的嘛？
作者
Author: LPOwenChan    时间: 2026-1-15 10:39
quanta 发表于 2026-1-15 08:50
博主您好，我用UMA起始优化的结构就是您帖子里的那个结构，您用uma优化的时候，是拿DFT的优化结果开始的 ...

初始结构基于Avogadro中UFF力场粗糙优化的结果，然后分别用UMA和DFT优化
作者
Author: Huschein    时间: 2026-1-15 23:46
推荐一下MAPLE程序，一款专门为各种MLFF(包含UMA)开发的计算化学平台，功能涵盖大部分Gaussian/ORCA的基础任务，用户可以不需要会任何代码也可以使用MLFF进行各种计算

https://chemrxiv.org/engage/chem ... 9bbafb3b11118e9807e
http://bbs.keinsci.com/thread-57731-1-1.html



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3