计算化学公社

标题: xtb输入文件如何写才能同时约束多个键的键长（已解决） [打印本页]
作者
Author: MercuryLamp    时间: 2023-8-12 10:32
标题: xtb输入文件如何写才能同时约束多个键的键长（已解决）
本帖最后由 MercuryLamp 于 2025-6-11 00:27 编辑

老师们好，我在计算一个较大分子（94个原子）与另一个小分子（29个原子）的过渡态时，得到的活化自由能是一个很大的负值（用B3LYP-D3/def2-SVP优化和计算频率，PWPB95-D3/def2-QZVPP计算的单点能），考虑到较大的分子中有一个较长的链取代基，可能柔性较大，我之前优化得到的反应物结构和过渡态结构可能不是最优构象，需要先对过渡态和反应物先做构象搜索。

在论坛中搜索相关帖子时，看到了王老师提到的一篇文献（https://pubs.acs.org/doi/10.1021/jacs.9b13725），我想仿照该文献以及卢老师的使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索(http://bbs.keinsci.com/forum.php?mod=viewthread&tid=16255)教程对过渡态的键长进行约束再用xtb跑动力学产生一些初始构象，但看了xtb的手册，只看到约束两个原子间距离（一个键的键长）的示例
  1. $constrain
  2.    distance: 1, 2, 1.4
  3.    angle: 5, 7, 8, auto
  4.    dihedral: 3, 4, 1, 7, 180
  5. $end
复制代码


现在我想同时约束两个键的键长，尝试了这两种写法
  1. $constrain
  2.    distance: 73, 123, 1.24
  3.    distance: 31, 73, 1.14
  4. $end
复制代码
  1. $constrain
  2.    distance: 73, 123, 1.24; 31, 73, 1.14
  3. $end
复制代码

但跑了几步查看轨迹，发现只能约束一个键的键长，不知应该如何写这个输入文件，还请各位老师指点，十分感谢

已解决：默认的force constant较小，将其设定为更大值后即可解决。
作者
Author: wzkchem5    时间: 2023-8-13 03:36
本帖最后由 wzkchem5 于 2023-8-12 20:40 编辑

如果xtb+molclus确实无法解决的话，有两种可能的替代方法：
（1）用orca的MD模块调用xtb跑动力学，orca的MD模块支持冻结不止一根键
（2）用CREST（https://crest-lab.github.io/crest-docs/）结合xtb做构象搜索
另：一个程序不太可能支持冻结一根键却不支持冻结两根键，因为多冻结几根键从代码实现角度没有任何额外难度。如果确定xtb说明书里没写如何冻结不止一根键，按照冻结一根键的方式简单类推写出的输入文件又不起作用的话，最好跟xtb官方反映一下，让他们更新一下说明书。不然的话将来的用户还会白白地在同一个坑里浪费时间
作者
Author: MercuryLamp    时间: 2023-8-13 10:54
本帖最后由 MercuryLamp 于 2023-8-13 10:58 编辑
wzkchem5 发表于 2023-8-13 03:36
如果xtb+molclus确实无法解决的话，有两种可能的替代方法：
（1）用orca的MD模块调用xtb跑动力学，orca的M ...

好的，谢谢老师您的建议
作者
Author: MercuryLamp    时间: 2023-8-15 21:20
本帖最后由 MercuryLamp 于 2023-8-15 21:33 编辑
wzkchem5 发表于 2023-8-13 03:36
如果xtb+molclus确实无法解决的话，有两种可能的替代方法：
（1）用orca的MD模块调用xtb跑动力学，orca的M ...

不好意思老师，有个问题想再请教一下您，在您之前回答其他坛友提到的这篇构象搜索的文献（https://pubs.acs.org/doi/10.1021/jacs.9b13725）中，我看到在step3时，作者提到对step2得到的过渡态构象使用相同的约束条件，用PBE-D3/def2-SVP级别进一步优化，并提到对能量较低或结构差异比较大的一些构象放开限制进行全优化，发现全优化得到的结构和限制性优化的结构能量差异很小，所以限制性优化得到的过渡态结构可以近似认为是全优化得到的过渡态结构。

我个人理解是指限制性优化得到的结构可以近似认为是最后的过渡态结构，不用放开限制再进行过渡态计算了，即后续计算高精度单点能时，可以直接使用限制性优化得到的结构。不知这么理解是否正确，麻烦老师指点，非常感谢。
作者
Author: wzkchem5    时间: 2023-8-15 22:57
MercuryLamp 发表于 2023-8-15 14:20
不好意思老师，有个问题想再请教一下您，在您之前回答其他坛友提到的这篇构象搜索的文献（https://pubs.a ...

不是的。你可以把限制性优化看作比优化过渡态稍低的一个理论级别，限制性优化和优化过渡态的结果相近，相当于是用高级别对低级别做了一个benchmark，证明以低级别的能量高低筛选出能量较低的过渡态构象是合理的。但是实际算自由能，仍然要用优化过渡态优化收敛的结构，原因有三：（1）如果没有优化到驻点结构，算出来的自由能校正是没有意义的；（2）你在验证限制性优化和优化过渡态的结果是否相近的时候，已经优化了最重要的那些构象的过渡态了，不用白不用；（3）剩下能量足够低但是还没优化过渡态的构象，优化收敛也很快，因为初猜结构离收敛结构很近。
作者
Author: MercuryLamp    时间: 2023-8-15 23:49
wzkchem5 发表于 2023-8-15 22:57
不是的。你可以把限制性优化看作比优化过渡态稍低的一个理论级别，限制性优化和优化过渡态的结果相近，相 ...

哦哦，我明白了，非常感谢老师您的解答
作者
Author: for_you    时间: 2023-9-25 10:53
要不试试这样固定
distance: 2,3,auto
distance: 2,4,auto
作者
Author: zh7729    时间: 2024-1-13 17:24
使用$constrain时把force constant调高一点或者直接用$fix应该可以满足要求
作者
Author: 紫薯糯米小丸子    时间: 2024-6-24 09:55
你好，请问你的问题解决了吗？我想问一下 在进行约束优化时，除了1.xyz坐标信息外，我另建了一个1.inp文件 包含了约束键长 键角和二面角的信息，输入是xtb 1.xyz --input 1.inp --opt tight --gfn 2 > 1.out，为什么在out文件里没有看到约束信息呀？
作者
Author: MercuryLamp    时间: 2024-7-4 21:19
紫薯糯米小丸子 发表于 2024-6-24 09:55
你好，请问你的问题解决了吗？我想问一下 在进行约束优化时，除了1.xyz坐标信息外，我另建了一个1.inp文件  ...

我最后是用orca调用xtb解决的，我不知道您1.inp文件是怎么写的，不太清楚您的问题出在哪，抱歉
作者
Author: Eudaimonia    时间: 2025-6-9 13:34
本帖最后由 Eudaimonia 于 2025-6-9 13:35 编辑

在较新的xtb版本(6.7.1)中，约束的相关信息写在引用文件的后面，可以通过查看该部分确定约束是否正常；按贴主的写法来束缚多根键应该是没有问题的

* started run on 2025/06/09 at 12:57:57.339
   ID    Z sym.   atoms
    1    6 C      1, 3, 4, 6, 9, 12-18, 20, 22, 24, 26, 28, 29, 31, 33, 35,
                  37, 39-41, 48, 49, 56, 58, 59, 61, 64, 67-69, 71, 73, 75,
                  77-80, 82, 84, 86, 90, 97, 98
    2    7 N      2, 5, 57, 60
    3    1 H      7, 8, 10, 11, 19, 21, 23, 25, 27, 30, 32, 34, 36, 38,
                  42-47, 50-55, 62, 63, 65, 66, 70, 72, 74, 76, 81, 83, 85,
                  87-89, 91-96, 99-104, 106
    4    8 O      105
constraining bond 2 9 to    1.4500000 Å, actual value:    1.4514170 Å
constraining bond 5 6 to    1.4500000 Å, actual value:    1.4520452 Å
constraining bond 60 61 to    1.4500000 Å, actual value:    1.4643481 Å
constraining bond 57 64 to    1.4500000 Å, actual value:    1.4648622 Å
spherical wallpotential with radius   11.6941115 Å

           -------------------------------------------------
          |                Calculation Setup                |
           -------------------------------------------------

此外，猜测还是能够观察到键的变动可能是默认的force constant太小导致的，可以适当放大这个数值保证约束够强



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