计算化学公社
标题: 有机大环分子在gmx中的力场选择与文件生成问题 [打印本页]
作者Author: Rosefinch 时间: 2022-3-12 16:08
标题: 有机大环分子在gmx中的力场选择与文件生成问题
我想做一个较大分子在石墨基底的上的自组装模拟来计算吸附能以及分子间作用力，之前的小分子都是采用的OPLS力场模拟的（参考的一篇文献），但是这个大环分子的原子个数有500多个，用sob老师给的ligpargen网站无法生成，TPPmktop显示有的原子是未知的化学环境，也无法生成。想请教一下如何生成拓扑文件，或者是否有其他的力场也适用于我这个体系的模拟，感谢！
(, 下载次数 Times of downloads: 30)
作者Author: 牧生 时间: 2022-3-12 17:46
本帖最后由 牧生 于 2022-3-12 18:18 编辑
考虑一下 http://sobereva.com/soft/Sobtop
What is Sobtop?
Very very very easily generating forcefield parameters and GROMACS topology file. Applicable to almost any system!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
作者Author: Rosefinch 时间: 2022-3-13 13:15
谢谢您的回复，我利用sobtop生成了完全基于GAFF力场的文件，但是在二面角那里没有生成参数，显示missing and treated as free rotor，请问我该如何填补参数，或者是否需要填补
作者Author: Rosefinch 时间: 2022-3-13 13:16
还想请问一下，如果我的有机分子采用了GAFF力场，那么两层石墨烯基底是否也需要用sobtop来生成GAFF力场的参数？
作者Author: 牧生 时间: 2022-3-13 21:13
本帖最后由 牧生 于 2022-3-13 21:24 编辑
1、如果你能看得懂这个帖子，就按照这个来 http://sobereva.com/635 但这个也有一个问题，就是你用的分子中，原子数较多，量化计算过程耗时间比较长。
2、如果想要简单易行，就参照sobtop的第二个例子 ，如果有些参数missing，还得要回到上一个的操作
3、原则上，都建议使用相同的力场，amber和gaff是兼容的。
作者Author: tjuptz 时间: 2022-3-13 21:57
sobtop没二面角的话就用acpype获得拓扑，别让它优化，自己算电荷
作者Author: sobereva 时间: 2022-3-14 03:05
看具体缺哪些，没结构和拓扑文件别人没法说
缺的二面角可以自己弄个模型体系做扫描并拟合得到，论坛上ztop可以拟合（以后sobtop也可能直接内置拟合功能），对于不重要的二面角可以凑合一个。如果二面角是刚性的就更简单了，用优化后的结构当平衡角度，力常数随便猜一个够大的即可。
作者Author: Rosefinch 时间: 2022-3-14 22:56
好的，谢谢sob老师，我那无法识别的二面角都集中在下图这6个红圈内，应该是可以视为刚性的。
作者Author: sobereva 时间: 2022-3-15 02:08
这肯定都是刚性的
在Sobtop的bonded_param.dat的$HARM_DIH里手动加入缺失的力场项，力常数随便设一个足够大的能维持刚性的值即可。
作者Author: Rosefinch 时间: 2022-3-25 12:28
还想请问一下sob老师：
1. 对于这种大环分子该如何去获取它的原子电荷？总共500+的原子，只有C、S、O、H四种元素，该如何选择单点能的基组？
2. 对于$HARM_DIH中手动设置的力常数，当我设置过大(比如30、50、100)时，做完能量最小化后，与这些二面角相关的环发生了较大的弯折，并且emtol最后也无法收敛得很好，导致后期平衡相出现过多LINCS WARNING导致模拟崩溃。之后我将力常数修改为了2.5kcal/mol/rad^2，进行完EM后，上图红圈所示的环出现了较小的弯折，emtol可以收敛到200一下，但是平衡相跑了一些步数后依旧出现了上述问题，请问这是因为我在$HARM_DIH中力常数设置有问题？还是因为分子没有设定原子电荷以及.mdp中没有规定coul设置造成的无法平衡？（因为此分子自组装主要是靠范德华作用实现，我就想直接在无电荷情况下跑模拟）分子中除H原子外xyz方向限制力大小均为500。
下图是我的两个.mdp文件，GRA组是石墨基底，XFZ是所有有机分子，m3.itp是这个大环分子的力场文件，可以麻烦老师帮我看一下嘛，万分感谢！！
(, 下载次数 Times of downloads: 20)
(, 下载次数 Times of downloads: 25)
(, 下载次数 Times of downloads: 40)
(, 下载次数 Times of downloads: 0)
作者Author: sobereva 时间: 2022-3-25 16:00
1 “选择单点能”语义不明
取出一个大环分子，把烃链短，构造较小的模型体系，然后算RESP电荷。用懒人脚本里的级别就可以
计算RESP原子电荷的超级懒人脚本（一行命令就算出结果）
http://sobereva.com/476（http://bbs.keinsci.com/thread-12858-1-1.html）
然后再构造个模型体系（不考虑所有的环），带着完整的烃链，再计算RESP电荷，级别同上。
以上两次得到的原子电荷恰当拼接构成单个完整环状+烃链的原子电荷。算RESP电荷期间注意恰当设置片段电荷约束，以便于完美拼接，做法参考下文。
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441（http://bbs.keinsci.com/thread-10880-1-1.html）
如果追求傻瓜化，用Sobtop自动调用Openbabel算MMFF94电荷，但质量比RESP有明显差距。
2 跟原子电荷没关系。参数问题
作者Author: Rosefinch 时间: 2022-3-25 22:04
谢谢sob老师，那请问这个力常数一般设置为多少能够刚好保持刚性又不会使环受力太大导致体系崩溃呢？
模拟时会突然弹出大量的报错说是下图中这种类型的C-H键旋转超过30度，请问该如何设置才能避免这类报错呢？（报错在第二张图）
(, 下载次数 Times of downloads: 31)
(, 下载次数 Times of downloads: 52)
作者Author: sobereva 时间: 2022-3-26 08:15
没具体文件没法说
检查拓扑文件里相关参数
用1fs步长，不设约束
作者Author: 牧生 时间: 2022-3-26 09:16
本帖最后由 牧生 于 2022-3-26 09:25 编辑
你这个结构式，确保是全原子模型么？
我认为也有可能是初始结构不合理导致，两个原子如果距离比较近，或被认为成键。画结构时，把两个原子距离拉远一点再优化，不然可能导致程序认错。优化后的位置和距离就合理了。
我有过类似的错误，看第五楼，和你现在的情况有点相似。
http://bbs.keinsci.com/thread-28435-1-1.html
作者Author: Rosefinch 时间: 2022-3-27 00:04
谢谢您的回复，我这个体系因为组装结构的原因导致有的烷烃链上的氢会挨得比较近，我试试修改一下模型。
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |