本帖最后由 欢乐多 于 2024-7-17 20:56 编辑
chatGPT辅助生成molclus懒人脚本：一键完成对复杂天然产物NMR和ECD计算 前言
复杂天然产物往往需要进行NMR和ECD的计算，由于多数天然产物专业人员量子化学计算机操作能力有限，所以NMR和ECD的计算困难重重。为了便于计算，尽量减少计算操作，笔者在chatGPT的帮助下，设计了一套计算复杂天然产物NMR和ECD的脚本，脚本经历多次修改和测试，终于实现了一键计算操作，此脚本能自动处理计算结果，计算完成后直接获得NMR和ECD结果，因此拿出来与诸位先进朋友分享交流，还请多多指教，批评指正。 本脚本及其配套的相关输入文件主要是在前人基础进行再次修改和创作，在使用本脚本进行科学研究，研究结果发表时，如若能正确引用，笔者将万分感激！
1. molclus懒人脚本简介
2. 脚本的依赖程序
脚本依赖Molclus 1.12，xtb -191025及crest功能，Gaussian 16 A03，ORCA 5.04 ，Shermo 2.4，Multiwfn 3.8 (dev) 等程序，在运行脚本前确保服务器中能使用这些程序。这些程序的安装及使用可参考社长大人Sobereva老师的相关博文： 以上软件可正常使用，最好将软件的路径加入系统环境变量，以达到输入软件名就可直接启动软件，准备完毕，就可启动脚本了，该脚本即是压缩包里00autorunall.sh，即是一个简单的Shell命令的组合，在CentOS 7.6和CentOS Stream 9中均可正常运行。
3. 首次使用需要的设定
脚本依赖的程序可正常使用后。 - 将本压缩包解压，进入该目录，赋予执行权限
- unzip molclus112.zip
- cd molclus112
- chmod +x *
- 给解压之后的目录文件中，提供一个输入文件traj.xyz；如果未提供traj.xyz文件将会出现错误提示coord file does not exit
- 修改00autorunall.sh中第27行，37行中xtb程序crest在本机中的实际绝对路径，例子中为：/root/soft/xtb/crest
- 修改05_set.ini中第15行，orca在本机中的实际路径，例子为：orca_path= "/root/soft/orca504/orca"
- 修改09_H01_nmr.txt，09_H02_nmr.txt，09_H03_nmr.txt中第9-18行的相关甲基或次甲基的编号，比如例子结构甲基氢编号为45-47，36-38，39-41，24-26，30-32，编号上一个命令10为Multiwfn屏蔽效应平均化的功能选项，所以第9-18行输为以下：
- 10
- 45-47
- 10
- 36-38
- 10
- 39-41
- 10
- 24-26
- 10
- 30-32
- 如果系统有发邮件功能，将脚本最后一行“echo " Normal termination of " | mail -s "MASTER complete " tttt@163.com”中邮箱部分，换成自己的邮箱，计算完成会发邮件提示。如何给系统安装发邮件功能，感兴趣者可自行Google。
- 准备就已经完成了。接着就是一行命令./00autorunall.sh即可获得NMR和ECD结果。如需要放进后台运行，输入命令nohup ./00autorunall.sh &> output.txt &，计算结束还会有邮件提示。
4. 脚本的计算流程
为了便于管理，笔者将脚本的计算流程分为10步:
4.1 构象搜索 使用Gentor, Spartn, XTB, Gromacs等均可，主要为脚本开始计算提供traj.xyz文件。压缩包中含一个例子的traj.xyz，请替换为自己的构象搜索结果文件。笔者《Windows系统Gromacs进行分子动力学构象搜索》的B站视频及计算化学论坛的博文中提到如何在Windows系统中用Gromacs进行分子动力学搜索，并产生traj.xyz。 有三个原因推荐用Gromacs的分子动力学的方法进行构象搜索。1，在NMR和ECD的计算中，能量最低的构象，或者接近能量最低的构象非常关键，如果没有找到或者找到的构象不接近实际实验中分子的构象或在自然界中可能存在的构象，那么之后步骤再高明的计算，再高精度的基组和泛函的计算，都会存在不容忽视的误差，因此，强烈推荐用分子动力学方法进行构象搜索，特别是含义环状体系的分子，如果连普通6元环船椅式构象都不能搜索到，怎么指望多元环更复杂的环状构象？2，分子动力学构象搜索可以用XTB，Gromacs，Amber等；3，笔者推荐Gromacs，易安装，易用，关键一个字“快”。
4.2 XTB 快速优化 产生初始能量由低到高排列稳定的构象，为Gaussian进一步优化提供可靠的输入文件。注意，首次使用请修改脚本中crest 的路径为本机的实际路径。crest为XTB中的一个辅助插件，当安装好XTB程序后，还需要下载相应版本的crest.tgz插件，解压后，赋予权限，可以放到XTB文件夹中。在使用中需要输入crest的实际绝对路径，比如脚本的路径在/root/soft/xtb/crest。
4.3 对上步中前15个能量最低的构象调佣Gaussian进行精确耗时的优化 如果不想做任何修改，本脚本默认设置数量为15，15基本不会漏掉搜索到的能量最低的构像，熟悉Molclus程序者可以自由调节03_set.ini中参数。计算结束后，并会输出能量由低到高排列的文件0003_output.xyz。压缩包中含有以下几个本步计算所需的文件，如果需要可自行调节相关参数。 03_set.ini 03_tem1.gjf 03_tem2.gjf 03_tem3.gjf
4.4 对0003_output.xyz进行振动分析 排除虚频的结构，得到0004_input.xyz 和能量由低到高排列的0004_output.xyz。压缩包中含有以下几个本步计算所需的文件： 04_set.ini 04_tem.gjf
4.5 调用ORCA 用较高级别的基组和泛函计算单点能。压缩包中含有以下几个本步计算所需的文件： 05_set.ini 05_tem.inp 注意，首次使用，需要修改05_set.ini中orca的路径
4.6 计算0004_input.xyz文件中所有构象的NMR 由于NMR相对不耗时，而复杂天然产物的NMR往往也不是一种泛函和基组就能确定的，所以选用三组基组和泛函对所有构象NMR进行计算，以便选取与实验结果比较接近的结果。以下文件是本步计算所需： 06_set_nmr.ini 06_tem01_nmr.gjf 06_tem02_nmr.gjf 06_tem03_nmr.gjf
4.7 计算0004_input.xyz文件中前5个构象的ECD 由于ECD较为耗时，在以往计算经验中发现第6之后的构象占比在1%以下，对整体ECD图谱影响甚微，所以选择前五个。确保我们选择了主要的构象之后，我们又选用三套基组和泛函对前5个构象的ECD进行计算，以便选取与实验结果比较接近的结果。以下文件是本步计算所需： 07_set_ecd.ini 07_tem01_ecd.gjf 07_tem02_ecd.gjf 07_tem03_ecd.gjf
4.8 Shermo计算Bloztmann比例 首先脚本会自动创建Shermo的输入文件0008_sher01_input.txt，文件第一列是寻找当前目录振动分析结果文件004_freq_gau*.out，第二例是orca计算结果的单点能。0008_sher02_input.txt是前五个构象的信息。结果输出为，前15个构象的Bloztmann比例0008_re_sher01rate.txt和前5个构象的比例0008_sher02_input.txt
4.9 Multiwfn处理得到权重之后的NMR结果 首先脚本会自动生成0009_nmr01_rate.txt文件，第一列是NMR的Gaussian计算结果文件，第二列是由Shermo计算的Bloztmann比例。 此步计算涉及以下文件， 09_C01_nmr.txt 09_C02_nmr.txt 09_C03_nmr.txt 09_H01_nmr.txt 09_H02_nmr.txt 09_H03_nmr.txt 熟悉Multiwfn的人会知道，以上文件就是Multiwfn的输入命令，由于链上C-C单键可以自由旋转，一般实验所得的甲基或有的次甲基的氢化学位移是一样的，而计算中对于一个特定构象上的三个氢化学位移不一样，为了与实验接近，需要对相应的氢进行平均化，因此需要设定09_H01_nmr.txt，09_H02_nmr.txt，09_H03_nmr.txt文件中相应的氢的编号。 由三套基组和泛函因此生成三个Multiwfn输入文件0009_nmr01_rate.txt，0009_nmr02_rate.txt，0009_nmr03_rate.txt。 此步计算结束，会有相应的三组权重之后H NMR和C NMR结果。 0009_HNMRdata01.txt 0009_CNMRdata01.txt 0009_HNMRdata02.txt 0009_CNMRdata02.txt 0009_HNMRdata03.txt 0009_CNMRdata03.txt
4.10 Multiwfn处理得到权重之后的ECD结果 首先脚本会自动生成00010_ecd01_rate.txt文件，第一列是Gaussian计算ECD结果文件，第二列是由Shermo计算的前五个构象的Bloztmann比例。如有需要可适当调节如下文件的相关参数： 10_ecd01.txt 10_ecd02.txt 10_ecd03.txt 由三套基组和泛函因此生成三个Multiwfn输入文件00010_ecd01_rate.txt，00010_ecd01_rate.txt，00010_ecd01_rate.txt。 此步计算结束，会有相应的三组权重之后ECD结果。 00010_dislin01.png 00010_spectrum01_curve.txt 00010_curveall01.txt 00010_dislin02.png 00010_spectrum02_curve.txt 00010_curveall02.txt
00010_dislin03.png 00010_spectrum03_curve.txt 00010_curveall03.txt
最后，在10step文件中还有每一步的.sh文件，可以分别对每一步进行运算。在molclus112目录中还有一个01.sh脚本，还可将每步命令放入其中，分步测试。 脚本还保留了.pbs文件的部分命令行，如果使用openpbs作业提交管理功能，只需将00autorunall.sh 修改为00autorunall.pbs，并修改其中的节点就可使用。
5. 总结
对于该脚本00autorunall.sh 中每一行的含义，感兴趣者可以输入到chatGPT中去了解命令的具体含义，同时也可以按照自己需要进行适当修改和改进。
6. 更新提示 依然保存原来版本molclus112，增加了一个新的版本6molclus112
6molclus112.zip
(2.76 MB, 下载次数 Times of downloads: 28)
新版本6molclus112中：
6.1 增加了一个python脚本，能够自动提取提取分子中甲基氢编号，因此，不需要提供09_H01_nmr.txt，09_H02_nmr.txt，09_H03_nmr.txt，需要有一个只含有一个分子构象的名字为*_traj2.xyz的文件，比如2R4S_traj2.xyz，后缀名“_traj2.xyz”是必须的。并且更改09_H_nmrstat.py脚本中第一行为自己系统的python路径，比如笔者python3.7路径如下#!/usr/local/python37/bin/python3.7，如何安装并行系统自带python版本和python3.7，请自行谷歌搜索安装。
6.2 增加了一个for循环，天然产物中往往好几个相对构型需要计算，比如2S3R4R, 2S3S4S, 2R3R4S等（温馨提示，几个构型的原子编号最好一致，以供与实验对比，若编号不一致，提取化学位移与实验对比将会非常麻烦），加入for循环后，可将这些分子统统放进一个文件，脚本将进行逐个计算。因此需要有多个构象的 “*_traj.xyz”和只有一个构象的“*_traj2.xyz”，比如
511rr_traj.xyz
511rr_traj2.xyz 511sr_traj.xyz 511sr_traj2.xyz
在当前目录下，即6molclus112的目录下，如果只想获得NMR结果，运行./01_mol_nmr.sh，或者nohup ./01_mol_nmr.sh &
6.4 如果想获得NMR和ECD结果，运行./00_mol_nmrecd.sh 或者nohup ./00_mol_nmrecd.sh &
6.5 脚本增加了记录脚本运行时长，运行结束会将运行时间发送到邮箱进行提醒，这个功能是笔者的最爱，非常的方便。
6.6 注意脚本./01_mol_nmr.sh 中第222行： - head -n 3 "${base_name}_08_sher00_input.txt" > "${base_name}_08_sher01_input.txt"
数字3需要与文件06_set_nmr.ini中第二行计算NMR构象ngeom= 3的数字一致，请根据计算机的繁忙程度自行调节，如果需要快速获得结果的话，可以选3，如果想保险一些，计算资源非常充足，可用7, 10等。此外，如果NMR和ECD计算完成后，若发现Gaussian正常完成，而分布比例有问题，或未能成功产生权重的结果，检查错误出现的位置，比如./01_mol_nmr.sh 中第222行数字与实际的06_set_nmr.ini数字不一样，相应修改后，可单独运行./002_shermo.sh。
6.7 增加了一个python脚本dp4.py，提取*_09_*NMRdata*.txt中权重之后的最终结果到excel, 这样可以一次获得每个构型的最终的化学位移及屏蔽常数。使用python脚本还需要安装一些数据库，如下： - python -m pip install numpy pandas
- python -m pip install xlsxwriter
- python -m pip install Workbook
- python -m pip install numpy
- python -m pip install pandas
- python -m pip install glob
- python -m pip install openpyxl
6.8 增加了一个python脚本nmrstat.py和shell脚本004_motosvm.sh，提供一个实验获得的C谱数据exp.txt，运行./004_motosvm.sh将会给出每个构型与实验的差异，判断是哪个构型。该nmrstat.py来自于SVM-M计算NMR方法。
提取码：qcl0
6.10 每次运行molclus程序时，该程序都会检查当前文件的.out，并删除，所以最新的6molcus112脚本多了一步将产生的out转变log的命令，以避开molclus的清理.out作业，相应的Multiwfn需要识别当前目录的.log文件，初始版本molclus112没有相关命令，因此建议使用较新6molclus112相关脚本。
最后，感谢中科院第一海洋研究所崔博士好友设计完成的09_H_nmrstat.py脚本，此外该脚本还得到其他很多人的帮助和支持，在此表示真诚的感谢！ 本文中相关脚本及其配套的输入文件主要是在前人基础进行再次修改和创作，在使用它进行科学研究，研究结果发表时，如若能正确引用，笔者将万分感激！ 由于本人非计算机出身，才疏学浅，脚本中难免有错误，请各位多多指教，谢谢大家！
| 
[复制链接 Copy URL]
发表于 Post on 2023-11-11 11:03:47
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author