计算化学公社

 找回密码 Forget password
 注册 Register
Views: 19186|回复 Reply: 13

[综合交流] 利用Gaussian或ORCA程序把分子结构拉直的几种方法

[复制链接 Copy URL]

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

发表于 Post on 2021-5-4 12:14:17 | 显示全部楼层 Show all |阅读模式 Reading model
利用Gaussian或ORCA程序把分子结构拉直的几种方法

文/Sobereva@北京科音
First release: 2021-May-4   Last update: 2021-Jun-15


1 前言

有人问怎么把一个分子从已有的较为卷曲的构象转化成比较平直的构象。这有一些实际意义,比如通过genmixmem程序构造磷脂双层膜、用packmol程序构造表面活性剂的囊泡,都需要用户提供分子比较平直的构象的结构文件。一个非常简单、直观、效果又好的做法是使用SAMSON可视化程序里的twist工具,笔者在《生成混合组分的磷脂双层膜结构文件的工具genmixmem》(http://sobereva.com/245)里还提供了我录的一段操作视频。但可惜如今SAMSON程序免费版当中已经没有twist工具了,需要买收费版才行(虽然收费版也可以申请免费试用)。还有一个办法是在GaussView等程序里一个一个改二面角,但无疑这非常费时费力,特别是分子可旋转的键很多的时候。

在此文,笔者介绍如何利用Gaussian和ORCA量子化学程序以不同方式实现把分子从卷曲的结构拉直,里面用到的做法对于读者研究其它一些问题可能也会有启发。本文使用的例子是《将Confab或Frog2与Molclus联用对有机体系做构象搜索》(http://bbs.keinsci.com/thread-20063-1-1.html)中的例子,即Actos分子。通过Molclus程序(http://www.keinsci.com/research/molclus.html)搜索出的此体系的能量最低构象如下,我们后文将要把此结构拉直:

1.png

本文涉及的所有文件可以在http://sobereva.com/attach/594/file.rar下载,上面的结构是文件包里的Actos.xyz。Gaussian使用的是G16 B.01版,ORCA是4.2.1版。


2 通过ORCA拉直分子:给两个原子间加外力

笔者写过不少与ORCA的安装和使用有关的文章,参看http://sobereva.com/category/ORCA/,相关常识性问题这里就不再说了。在ORCA中优化时可以给两个原子在连线方向上加外力驱使它们远离,很适合用于拉直分子。对于上面图中展示的Actos分子,显然我们应当把最末端的两个重原子,即8和25号之间拉远。为此,写一个ORCA输入文件:

! AM1 opt nopop
%geom
POTENTIALS
{ C 7 24 3.0 }
end
end
* xyz   0   1
[Actos.xyz里的坐标]
*


这里POTENTIALS下面花括号里的C代表constant force。注意原子序号是从0开始算的。3.0是把俩原子间拉远的力的大小,单位是nN。数值大小不是很重要,一般个位数大小就可以。数值太小的话拉不开,而太大的话会最终把一些化学键过度拉长、键角被过度拉大(不过键伸缩、键角弯曲相对于二面角扭转来说是很刚的,只要外力不是特别大这就不是明显问题)。这里用的AM1半经验方法是ORCA里能直接用的几乎最便宜的方法(虽然ORCA也支持力场计算,但还得搞参数,这里不考虑),注意此方法不支持并行计算。如果嫌AM1太low,也可以用稍微好一点的HF-3c(Hartree-Fock结合额外校正项),对不太大的体系也非常便宜。

用ORCA运行上面的输入文件,三分钟就算完了,最终得到的结构如下(即file文件包里ORCA目录下的Actos.xyz),可见完美地实现了我们的目的,分子被拉得很直

2.png

file文件包里的Actos_trj.xyz是优化轨迹,感兴趣者可以将之拖到VMD程序(http://www.ks.uiuc.edu/Research/vmd/)里看看优化过程。顺带一提,考察加外力对化学体系的影响专门有人做过大量研究,参见Chem. Rev., 112, 5412 (2012)、Angew. Chem. Int. Ed., 58, 5232 (2019)。


3 通过Gaussian拉直分子

在Gaussian程序中不支持直接给原子加外力,但有其它的方法也可以实现把分子拉直,将在下面几节介绍。

3.1 方法1:柔性扫描

Gaussian中做柔性扫描在《详谈使用Gaussian做势能面扫描》(http://sobereva.com/474)里笔者详细介绍过。显然,做两端的两个原子距离逐渐拉长的柔性扫描最终就可以得到被拉得平直的结构。写一个这种任务的Gaussian输入文件,如下所示

%nprocs=1
#T UFF opt(loose,nomicro,modredundant)

Generated by Multiwfn

  0  1
[Actos.xyz里的坐标]

8 25 S 50 0.5


这里有几个需要注意的地方:
(1)为了耗时尽可能低,这里使用了Gaussian支持的最便宜的方法,即分子力场。这里用的是较粗糙但支持元素最广、用着最省事(不需要指定原子类型)的UFF力场。由于本例只需要得到一个较为平直的构象,对质量什么没要求,所以UFF就够了。
(2)柔性扫描每一步相当于一次限制性优化,为了节约时间,用了loose收敛限降低每一步扫描的耗时。
(3)opt里必须写nomicro。因为分子力场的优化在Gaussian里有专用的代码,不受modredundant设置的控制。而加了nomicro后,就会用通用的几何优化代码,虽然远没有分子力场专用的优化代码快,但此时可以设置柔性扫描。
(4)分子力场计算时,至少对于这个体系,并行计算会令耗时不降反升,所以这里靠%nprocs=1刻意要求不用并行计算。
(5)扫描的步长太小的话把分子拉直需要的步数太多,显然不行,而步长太大则容易导致中途报错(此时还没充分拉直),根据我的经验步长设0.5埃左右比较合适。由于事先往往不好估计拉直时两个原子相距多远,所以可以把扫描步数设得比较大(比如本例设50步),当分子已经被拉得很直后,继续扫下一步时程序通常会报错,此时用GaussView载入输出文件查看扫描轨迹,若最后一帧结构就是想要的平直的结构的话就取最后一帧结构即可。如果步数一开始设小了,导致最后没有完全拉直也没关系,拿最后一步的结构继续做拉长的扫描即可。
(6)当前体系是之前我用DFT优化过的结构,已经很合理了,所以没有明确靠geom=connectivity指定原子间连接关系,而直接让Gaussian根据当前结构猜连接关系并用于力场计算。如果初始结构没有优化过,最好让GaussView保存出带有恰当连接关系的输入文件以免自动判断的连接关系有误。

使用Gaussian对上面的输入文件进行计算,算了1分钟左右扫描到第31步后报错,报错前最后一步结构如下所示,可见是我们想要的平直的结构。

3.png


3.2 方法2:利用外电场加外力

Gaussian里可以加外电场,外电场会给带电粒子施加额外的力,在《一篇文章深入揭示外电场对18碳环的超强调控作用》(http://sobereva.com/570)里介绍的笔者的文章中有很多介绍和讨论。因此,可以在UFF力场计算时,把25号原子坐标冻结住,而让8号原子带+1.0原子电荷,然后按照下图所示向X轴正方向加足够强度的电场,这样8号原子就会受电场力往X正方向运动,从而实现拉直分子的目的。

4.png

此例输入文件如下

%nprocs=1
# UFF opt(loose,nomicro,modredundant) field=x-300 nosymm

Generated by Multiwfn

  0  1
N      2.96734200   -1.21688000    2.24577900
C      2.32051000   -1.62618600    3.34057500
C      1.63462400   -2.84168400    3.45382100
C      1.64984500   -3.67567100    2.33272400
C      2.32098400   -3.26708600    1.18384000
C      2.96534100   -2.02795600    1.17397000
C      0.96378000   -3.24289700    4.74461900
C--1.0      1.95041300   -3.88295500    5.73536500
C      3.62426800   -1.48867000   -0.07067000
[略]

25 F


此例C--1.0代表对这个8号原子不手动指定原子类型(两个横杠之间没写内容),原子电荷为+1.0。没设原子电荷的原子的原子电荷默认为0。Gaussian里电场方向和习俗相反,因此field=x后面是负号。电场强度用300(0.03 a.u.)一般就可以,如果发现最终结果怪异,如有些键被拉得太长,结构都扭曲了,应当适当减小电场以减小外力;如果分子拉得还不够直,可以稍微增大电场来加大外力。nosymm必须写,要不然结构被自动弄到标准朝向下后外电场相对于分子的方向就和期望不符了,关于nosymm更多信息见《谈谈Gaussian中的对称性与nosymm关键词的使用》(http://sobereva.com/297)。nomicro还是要写,要不然冻结设定不生效,其电场设置也不会对优化起作用。

用Gaussian计算上面的输入文件,68步后收敛,耗时才十几秒。最终结构如下,可见是我们想要的

5.png


3.3 方法3:令两端的原子彼此间受到静电斥力

这个做法超级简单,而且超级快速(直接用分子力场的专用优化代码),比前述两种方法通常更值得推荐。原理很简单,即给两端的原子都设数值较大的相同符号的原子电荷,二者通过强烈的静电斥力自然就会彼此远离,从而把分子拉直。输入文件如下,只有8号和25号原子设了原子电荷,都设为了+10.0。设多大合适应当看实际情况,比如如果分子特别长,那么原子电荷应设大一些,要不然静电斥力不够强。此例输入文件如下

%nprocs=1
# UFF opt

Generated by Multiwfn

  0  1
[略]
C      0.96378000   -3.24289700    4.74461900
C--10.0      1.95041300   -3.88295500    5.73536500
C      3.62426800   -1.48867000   -0.07067000
[略]
S     -2.51218200   -3.05838400    3.26014600
O--10.0     -4.37962400   -3.91727700   -0.05061200
O     -1.13156500   -5.32519500    2.86775900
[略]


用Gaussian计算,仅花了一秒钟就算完了!结果如下所示,非常平直,特别理想!

6.png

网友还给我发过一个含有200多个原子的弯曲的很长的分子,如下所示

7.png

用这一节的做法也很顺利地优化成了平直的构象,如下所示。计算仅仅花了5秒钟。

8.png


4 总结

本文介绍了笔者想到的四种把分子从弯曲结构拉成平直的方法,思路各有不同。实际当中会遇上各种特殊情况,比如可能需要同时拉直两条链(如磷脂的两条尾巴),应当在理解这些方法思想的基础上根据实际情况恰当选择、灵活运用。

本文的做法适合一般中、小分子体系。如果是大分子体系,比如蛋白质,比较适合GROMACS等专门的基于分子力场计算的程序。在GROMACS里面可以用pull相关的选项做拉伸动力学(见比如《在Gromacs中模拟金纳米线拉伸过程(含视频)》http://sobereva.com/153),也可以将体系一端的原子定义为冻结组(freezegrps)而另一端定义为受常外力的组(acc-grps,结合accelerate设置)从而拉开。

值得一提的是,本文介绍的虽然是把分子拉直的方法,但实际上也可以将这些方法应用到研究环状、笼状、簇状等化学体系对外力造成的拉伸和压缩问题的研究上。


补充:使用Avogadro程序拉直分子的做法

通过Avogadro可视化程序也可以拉直分子,此程序可以在http://avogadro.cc免费下载。在此程序中拉直分子比较直观,直接按住鼠标拖拽即可,这和SAMSON可视化程序里的twist工具非常类似。启动Avogadro后,先载入pdb之类的结构文件,点击下图上方箭头所指的按钮,然后再点击下图左侧箭头所指的Start按钮,之后按住鼠标左键拖动特定原子即可,按照下图的绿色箭头分别把两端的原子往两头拉就能最终拉直。

9.png

Avogadro的这个功能实际上是持续通过指定的力场优化体系几何结构,因此当你拖动某个原子时,其它原子的位置会自发地弛豫来降低能量,因此键长、键角、二面角始终能保持比较合理的状态。但我发现Avogadro的这个功能拉直小分子不错,而拉直第3节最后那个特别长的聚合物则很难实现,而且特别卡顿。

评分 Rate

参与人数
Participants 3
eV +15 收起 理由
Reason
lijiayisjtu + 5 好物!
panernie + 5 好物!
zsu007 + 5 赞!

查看全部评分 View all ratings

北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

393

帖子

0

威望

1440

eV
积分
1833

Level 5 (御坂)

娃娃儿鱼

发表于 Post on 2021-5-4 12:25:55 | 显示全部楼层 Show all
本帖最后由 hdhxx123 于 2021-5-4 12:27 编辑

最后这个分子力场加电荷斥力排开的方法莫名萌,联想到聚电解质刷子伸展和塌缩的感觉
真·探

220

帖子

1

威望

1763

eV
积分
2003

Level 5 (御坂)

发表于 Post on 2021-5-4 14:04:18 | 显示全部楼层 Show all
本帖最后由 lyj714 于 2021-5-4 14:08 编辑

补充,我是GROMACS忠实用户,亲测gromacs中的三种方式均可实现本文例子,产生较好的效果。

三种方式分别为拉伸外加电场外加力(加速组)
值得注意的是对于特定体系,拉伸时间,电场强度,加速的大小需要自己尝试,一般设置大点效果明显
相关输入和输出文件见附件,本人使用的是GROMACS 2019.6版本

注意:
不同GROMACS大版本之间的pull,elec的mdp文件部分有显著的差异,一定要看对应版本的GROMACS手册





file.7z

5.33 MB, 下载次数 Times of downloads: 40

pull,acc,elec

评分 Rate

参与人数
Participants 2
eV +10 收起 理由
Reason
sobereva + 5
ggdh + 5 正宗g粉啊

查看全部评分 View all ratings

375

帖子

2

威望

1525

eV
积分
1940

Level 5 (御坂)

发表于 Post on 2021-5-19 19:09:57 | 显示全部楼层 Show all
这是我的感觉 不知道正不正确
我一般画完了一个分子,先用xtb优化,然后再上gaussian和orca。
感觉那种长分子,一版GFN2-xTB优化出来是长的伸展的,但gaussian和orca B3LYP-D3(BJ)/def2-TZVP优化出来就是蜷缩的,不晓得为啥。

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2021-5-20 17:57:23 | 显示全部楼层 Show all
wuzhiyi 发表于 2021-5-19 19:09
这是我的感觉 不知道正不正确
我一般画完了一个分子,先用xtb优化,然后再上gaussian和orca。
感觉那种长 ...

这和色散作用描述有密切关系,色散吸引作用比较倾向于蜷缩的
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

375

帖子

2

威望

1525

eV
积分
1940

Level 5 (御坂)

发表于 Post on 2021-5-20 18:01:08 | 显示全部楼层 Show all
sobereva 发表于 2021-5-20 17:57
这和色散作用描述有密切关系,色散吸引作用比较倾向于蜷缩的

我也觉得是色散作用,但GFN2-xTB里用DFT-D4来描述色散作用,而B3LYP-D3(BJ)里用DFT-D3(BJ),理论上来说应该差不多。。。。

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2021-5-20 20:25:29 | 显示全部楼层 Show all
wuzhiyi 发表于 2021-5-20 18:01
我也觉得是色散作用,但GFN2-xTB里用DFT-D4来描述色散作用,而B3LYP-D3(BJ)里用DFT-D3(BJ),理论上来说应 ...

但毕竟GFN2-xTB底子差得多
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2021-6-15 22:07:42 | 显示全部楼层 Show all
文末新加入了使用Avogadro程序拉直分子的做法
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

15

帖子

0

威望

193

eV
积分
208

Level 3 能力者

发表于 Post on 2021-6-22 09:15:32 | 显示全部楼层 Show all
sob老师,我写bat文件想利用orca批处理施加外力的分子,程序报错orca_gstep不是批处理程序,有什么好方法可以批量处理吗

1188

帖子

5

威望

2758

eV
积分
4046

Level 6 (一方通行)

发表于 Post on 2021-6-22 09:46:45 | 显示全部楼层 Show all
阿社社社 发表于 2021-6-22 09:15
sob老师,我写bat文件想利用orca批处理施加外力的分子,程序报错orca_gstep不是批处理程序,有什么好方法可 ...

orca所在目录加进PATH试试

15

帖子

0

威望

193

eV
积分
208

Level 3 能力者

发表于 Post on 2021-6-22 10:27:22 | 显示全部楼层 Show all
snljty 发表于 2021-6-22 09:46
orca所在目录加进PATH试试

谢谢,问题解决了

114

帖子

0

威望

2431

eV
积分
2545

Level 5 (御坂)

发表于 Post on 2022-5-11 14:59:19 | 显示全部楼层 Show all
请问社长,avogadro1.2如何像gv6里面改变键长键角二面角的时候设定一端fix,另一端整个groups整体变化呢,我看菜单栏里面只有view-properties选择对应属性并修改参数,但是修改后整个视角都变了。


另外必须吐槽还没开发完的avogadro2.0,c++,qt5开发,居然比1代还容易崩溃。。。

4万

帖子

99

威望

4万

eV
积分
89863

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2022-5-11 20:14:18 | 显示全部楼层 Show all
neocc 发表于 2022-5-11 14:59
请问社长,avogadro1.2如何像gv6里面改变键长键角二面角的时候设定一端fix,另一端整个groups整体变化呢, ...

我从来不用avogadro建模,比gview差太多了
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

114

帖子

0

威望

2431

eV
积分
2545

Level 5 (御坂)

发表于 Post on 2022-5-11 20:29:48 | 显示全部楼层 Show all
sobereva 发表于 2022-5-11 20:14
我从来不用avogadro建模,比gview差太多了

avogadro确实不太行
还是gv靠谱

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-2 00:33 , Processed in 0.633032 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list