计算分子内质子转移的轨迹，该用什么软件呢？大家来给我的想法提提意见吧

算一算

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

算一算

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

预平衡是什么意思?

k64_cc

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

yflchx

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

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

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

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


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

算一算

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

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

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

算一算

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

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

yflchx

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

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

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

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



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

算一算

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

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

sobereva

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

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

yflchx

算一算 发表于 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，太大了吧？

算一算

谢谢大家热情回复，刚发现每天的评分是有配额的，那我就明天继续弄吧。不过这个评分也不足以表达我的感激。你们提供的各种信息都对我很有帮助，真心谢谢

算一算

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

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

算一算

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

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

算一算

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

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

算一算

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

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

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