计算化学公社

标题: 计算分子内质子转移的轨迹，该用什么软件呢？大家来给我的想法提提意见吧 [打印本页]

作者
Author: 算一算 时间: 2018-8-8 07:43
标题: 计算分子内质子转移的轨迹，该用什么软件呢？大家来给我的想法提提意见吧
苯环上有一个氨基（-NH3）、一个羧基（-COOH），质子在这两个基团之间转移，如果要计算这个质子的运动轨迹，有没有哪款免费软件能算全过程呢？
感觉已有的各种分子模拟软件都是针对周期性体系的，那我这个不是周期性的，就是单个分子（目前不考虑溶剂），不知能不能用那些软件，如果能用，用哪款最适合？并且，该选哪种算法，velocity verlet? Leapfrog? 还是别的？
如果真的没有一款软件能算这种东西，需要我手动一步一步的计算的话，大家看看我这个方案可行不可行：假设每一步的步长是1飞秒（这个是不是太长了呢？），我每一步用Gaussian里的DFT方法（DFT算这个是不是足够了？）算质子受力，然后根据这个力把质子移动到下一飞秒它该在的位置，然后再进行DFT计算，这样一步一步的算，直到质子最终和另一基团成键。
但是我对这个一步一步手算的方案有一个很大的疑问：每一步DFT计算时，Gaussian算出的都是能量最低的电子结构，但是在现实中，这个质子转移的过程中，不一定总是能量最低的结构（在几个飞秒内它不能迅速转换电子结构以达到能量最低），那我这样算出来的质子运动轨迹，是不是就和现实有很大偏差了？这种情况有解决办法吗？如果用分子模拟软件，会不会也出现这种情况？
谢谢大家花时间阅读我这么多问题！

作者
Author: alonewolfyang 时间: 2018-8-8 08:30
个人感觉，既然说到分子受力，为何还要用DFT？分子力学不行吗？

作者
Author: sobereva 时间: 2018-8-8 09:13
直接在Gaussian中通过ADMP或BOMD关键词做从头算动力学模拟即可

作者
Author: liyuanhe211 时间: 2018-8-8 15:56
你想通过“质子转移的轨迹”说明什么问题？

要用“分子模拟”的软件做，分子模拟即使用周期性盒子也可以做“孤立”的分子（稀溶液），盒子足够大就可以了，这不是问题。问题是大部分力场不会断键的。

简单的算受力、按受力走的方法显然有问题，因为如果从极小点出发、受力是0，哪儿都走不到。简单的从头算动力学的问题也不是“电子结构是能量最低的”，因为电子的弛豫快，往往认为电子结构处在基态造成的问题不太大。

真去做、你会发现你得到的不是质子转移的轨迹，大部分都是几何扭转之类的低能过程，气相的质子转移（其实应该是氢自由基转移）还不是什么特别优先的过程，所以还得等一堆低能的其他反应。不现实。

以及你确定你要的不是IRC而必须是动力学？

作者
Author: laomalao 时间: 2018-8-8 18:03
IRC分析就可以了

作者
Author: Accelerator 时间: 2018-8-9 06:05

liyuanhe211 发表于 2018-8-8 15:56
你想通过“质子转移的轨迹”说明什么问题？

要用“分子模拟”的软件做，分子模拟即使用周期性盒子也可以 ...

想用MD研究这个过程的话可以从质子转移的过渡态开始吧
虽然不知道这样做能得到什么额外的信息

作者
Author: sobereva 时间: 2018-8-9 06:54

alonewolfyang 发表于 2018-8-8 08:30
个人感觉，既然说到分子受力，为何还要用DFT？分子力学不行吗？

一般的分子力场（非反应力场）没法描述诸如质子转移这样有键的形成和断裂的过程

作者
Author: Daniel_Arndt 时间: 2018-8-10 05:51
你的描述让我想起来14年的一篇Nature Chemistry用TeraChem这个软件来模拟原始大气产生氨基酸的。它的SI里面好像是有一些质子迁移的东西。他们用的方法好像是周期性地压缩体系的体积，得到短时间的高温，再解除对体积的限制，让温度降下来，然后再循环。但按照他们文中的描述的话，数据处理时有点麻烦，他们自己好像是写了一堆python脚本才完成的。
我的感觉是，如果想从过渡态开始研究质子迁移的轨迹，就得保证找到了质子迁移的决速步的过渡态，而且其能垒不能太低，然后才可以在过渡态附近采样，再往反应物、产物两个方向计算轨迹。这方面就得参考Texas Tech University的Hase、Texas A&M University的Singleton、大连化物所张东辉院士等人的工作。
如果你想从反应物开始研究的话，那么就得参考Emory University的Bowman的一些工作。只不过这样的话，最后要看的轨迹的数目有点小多。

作者
Author: 算一算 时间: 2018-8-10 18:32

alonewolfyang 发表于 2018-8-8 08:30
个人感觉，既然说到分子受力，为何还要用DFT？分子力学不行吗？

是这样的，这个过程质子受分子各部位的电子密度的影响很大，所以不能用简单的分子力学来描述。也正如7楼所说，想描述断键过程的话，用DFT免不了的。

作者
Author: 算一算 时间: 2018-8-10 18:37

sobereva 发表于 2018-8-8 09:13
直接在Gaussian中通过ADMP或BOMD关键词做从头算动力学模拟即可

谢谢这个建议。我去Gaussian官网看了一下BOMD，没太明白它这个步长是怎么弄的……它原话是这么写的：
StepSize=n
Sets the dynamic step size to n*0.0001, in the appropriate units. The default step size is 0.25 amu1/2*Bohr except for GradientOnly calculations where it is 0.025 femtoseconds.
请问这个amu开平方再乘以Bohr是描述什么物理量的？如果我选后面那种，就是gradientonly的话，我想把步长设为1飞秒，那我是不是应该写stepsize=10000?

作者
Author: 算一算 时间: 2018-8-10 18:47

liyuanhe211 发表于 2018-8-8 15:56
你想通过“质子转移的轨迹”说明什么问题？

要用“分子模拟”的软件做，分子模拟即使用周期性盒子也可以 ...

是我们头儿丢给我的任务，他就是要看质子运动的路线，让我把每一步的坐标都搞出来。我暂时也不知道他到底想用这个说明什么问题。
你说的从能量极小的几何结构出发受力是0，这个我清楚，我会选择受力不是0的出发点的。
你说的“简单的从头算动力学的问题也不是“电子结构是能量最低的”，” 那它是怎么选择电子结构的呢？我以为软件在每一步都是会在当时的几何结构上优化出能量最低的波函数，如果不是那样的话，软件是怎么做的，它是根据前一步的波函数做初猜？然后尽量让电子结构保持在同一个势能面上是吗？

作者
Author: 算一算 时间: 2018-8-10 18:52

laomalao 发表于 2018-8-8 18:03
IRC分析就可以了

你和楼上都提到IRC了，但是我跟我们头儿商量了一下，他坚持让我做动力学分析，然后把质子运动路线给他看。虽然他也不清楚具体该怎么做。

作者
Author: 算一算 时间: 2018-8-10 18:53

Accelerator 发表于 2018-8-9 06:05
想用MD研究这个过程的话可以从质子转移的过渡态开始吧
虽然不知道这样做能得到什么额外的信息

嗯嗯，我也在考虑从过渡态开始呢

作者
Author: 算一算 时间: 2018-8-10 18:55

Daniel_Arndt 发表于 2018-8-10 05:51
你的描述让我想起来14年的一篇Nature Chemistry用TeraChem这个软件来模拟原始大气产生氨基酸的。它的SI里面 ...

这个信息量好大，好的，我一个一个找找看，我不怕数据多，就怕自己方法用的不好。

作者
Author: 算一算 时间: 2018-8-10 18:59
谢谢大家热情回复，刚发现每天的评分是有配额的，那我就明天继续弄吧。不过这个评分也不足以表达我的感激。你们提供的各种信息都对我很有帮助，真心谢谢

作者
Author: yflchx 时间: 2018-8-10 19:37

算一算发表于 2018-8-10 18:37
谢谢这个建议。我去Gaussian官网看了一下BOMD，没太明白它这个步长是怎么弄的……它原话是这么写的：
St ...

这个跟轨线积分的算法有关系。Gaussian里面的BOMD默认用Hessian-based predictor-corrector（预测-校正）方法，此时的步长单位是sqrt(amu)*bohr，跟构建即时势能面有关。默认步长为0.25 sqrt(amu)*bohr，相当于stepsize=2500。

这种方法需要update Hessian，一般选取update=5左右比较理想（权衡计算耗时和能量守恒）。当然，如果模拟时间过长（比如若干个ps），能量守恒也会明显变差。

我没测试过GradientOnly，不确定GradientOnly是不是用的Velocity-Verlet。如果设定步长为1fs，太大了吧？

作者
Author: sobereva 时间: 2018-8-11 06:49

算一算发表于 2018-8-10 18:37
谢谢这个建议。我去Gaussian官网看了一下BOMD，没太明白它这个步长是怎么弄的……它原话是这么写的：
St ...

论坛里就有过讨论，可以首页google框搜索
amu1/2*Bohr是质权坐标的单位。高斯里有固定距离步长的算法，也有固定时间步长的算法

作者
Author: 算一算 时间: 2018-8-12 04:20

yflchx 发表于 2018-8-10 19:37
这个跟轨线积分的算法有关系。Gaussian里面的BOMD默认用Hessian-based predictor-corrector（预测-校正） ...

谢谢批评我的步长。其实我不知道该选多大的步长合适。我隐约觉得，如果是velocity verlet或者leapfrog的话，也许步长大一点也没事……

作者
Author: yflchx 时间: 2018-8-12 08:32
本帖最后由 yflchx 于 2018-8-12 08:35 编辑

算一算发表于 2018-8-12 04:20
谢谢批评我的步长。其实我不知道该选多大的步长合适。我隐约觉得，如果是velocity verlet或者leapfrog的 ...

可以测试，选择步长要确保：

1.能量守恒、角动量守恒如何？

2.相同的初始条件，轨迹之情形基本不受步长影响。

注：BOMD对于稍微大些的体系，非常耗时。

作者
Author: 算一算 时间: 2018-8-12 08:36

sobereva 发表于 2018-8-11 06:49
论坛里就有过讨论，可以首页google框搜索
amu1/2*Bohr是质权坐标的单位。高斯里有固定距离步长的算法， ...

搜到了，谢谢。另外我还想问一下，4楼和5楼都提到了用IRC解决这个问题，但是我觉得用MD可以加入键的振动，所以更真实一点。但是我有一点担忧就是，如果用MD算的话，每一步都是重新去找最稳定的电子结构，那么有可能出现在短短几步内在不同的PES之间切换的局面，而现实中，短时间内不会在不同的势能面间切换，所以这就导致MD的结果不那么真实，从这个角度讲，或许IRC更真实一点。您对此怎么看？

作者
Author: 算一算 时间: 2018-8-12 08:41

yflchx 发表于 2018-8-12 08:32
可以测试，选择步长要确保：

1.能量守恒、角动量守恒如何？

谢谢这么迅速的回复。我刚才找了一个很简单的小分子试探了一下BOMD，发现那个ETOT一直是在变的，这个是我的问题吗，还是高斯的问题？我原以为这个计算会是constant energy的，看到总能量一直在变，有点懵

作者
Author: yflchx 时间: 2018-8-12 09:35
本帖最后由 yflchx 于 2018-8-12 09:42 编辑

算一算发表于 2018-8-12 08:41
谢谢这么迅速的回复。我刚才找了一个很简单的小分子试探了一下BOMD，发现那个ETOT一直是在变的，这个是我 ...

所谓总能量守恒是指总能量守恒到一个很小的能量范围，并非固定数值不变。

比如：短时，1ps左右的基本应该守恒到 10^（-5）hartree，也就是说总能量浮动在这个范围之内。

既然准备用Gaussian里面的BOMD，H. B. Schlegel组的文章（BOMD计算相关的）应该先读几篇。

如果总能量浮动过大，减小步长试试。减小步长、提高SCF收敛标准、增加DFT积分格点都是解决途径。

作者
Author: k64_cc 时间: 2018-8-13 12:00
质子转移的时间步长常用0.5fs，因为氢原子质量太小，不支持太大的时间步长。
以及MD能跨越的势垒也就是3/2kbT的大小。势垒低于这个数，就平衡好然后祈祷；势垒如果高于这个数，可以在TS做起点跑一堆轨迹分析动力学。
如果选择从TS开始跑轨迹，那么请别用Gaussian直接跑BOMD，好歹也要先预平衡一下……

作者
Author: 算一算 时间: 2018-8-16 21:54

k64_cc 发表于 2018-8-13 12:00
质子转移的时间步长常用0.5fs，因为氢原子质量太小，不支持太大的时间步长。
以及MD能跨越的势垒也就是3/2 ...

预平衡是什么意思?

欢迎光临计算化学公社 (http://bbs.keinsci.com/)