本帖最后由 欢乐多 于 2025-4-22 11:33 编辑
chatGPT辅助生成molclus懒人脚本：一键完成对复杂天然产物NMR和ECD计算 前言
复杂天然产物往往需要进行NMR和ECD的计算，由于多数天然产物专业人员量子化学计算机操作能力有限，所以NMR和ECD的计算困难重重。为了便于计算，尽量减少计算操作，笔者在chatGPT的帮助下，设计了一套计算复杂天然产物NMR和ECD的脚本，脚本经历多次修改和测试，终于实现了一键计算操作，此脚本能自动处理计算结果，计算完成后直接获得NMR和ECD结果，因此拿出来与诸位先进朋友分享交流，还请多多指教，批评指正。 本脚本及其配套的相关输入文件主要是在前人基础进行再次修改和创作，在使用本脚本进行科学研究，研究结果发表时，如若能正确引用，笔者将万分感激！
注：脚本在2025年3月进行了更新，简单的介绍请看
Molclus 懒人脚本简介.pdf
(205.65 KB, 下载次数 Times of downloads: 38)
，视频介绍请B站《Molclus懒人脚本：一键完成复杂天然产物 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。（脚本在2025年3月进行了更新，不再需要提供原始的traj.xyz文件，更新的6molclus112中放入由Chemdraw3D 保存的.pdb格式的结构，比如004rs.pdb，此外还需要放入ECD实验得到的.csv文件和实验的碳谱数据.txt，比如004rs_nmr.txt和004rs.csv，便于计算后的数据分析。）
- 修改00autorunall.sh中第27行，37行中xtb程序crest在本机中的实际绝对路径，例子中为：/root/soft/xtb/crest
- 修改05_set.ini中第15行，orca在本机中的实际路径，例子为：orca_path= "/root/soft/orca504/orca"
- （脚本在2025年3月进行了更新，已不需要手动修改，脚本自动识别甲基编号）修改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元环船椅式构象都不能搜索到，根本不能指望出现多元环的环状构象的变化。比如，一个普通的共轭芳香分子，11和20位羰基并不在一个平面上，理论上共轭芳香环存在至少两种构象，结构如下：
分子动力学构象搜索，可以搜索到多种芳香环的构象，下图是用VMD对齐显示出1-1000个构象的叠加图，11和20位的羰基出现多次的上下摆动，如下：
另外通过这种叠加图的显示，我们也可以清楚的知道，构象搜索是否充分，如果我们发现侧链上下不对称，侧链只在环的一面摆动，那么构象搜索就不充分，一般这种情况很少发生。
而非分子动力学只是对侧链可旋转的键进行搜索，并不能对共轭大芳香环进行构象搜索，下图是用VMD对齐显示出1-1000个构象的叠加图，可以看到，侧链搜索比较充分，环的上下都有出现，而共轭芳香环构象并没有发生该有的变化，如下：
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
计算完成之后，进行实验图谱与计算图谱的拟合，参考B站视频，一般实验图谱是.csv格式的数据，笔者以前使用origin进行拟合，目前有了chatgpt之后，改用python中matlab功能进行数据分析，脚本是
sim_ecd.py
(6.15 KB, 下载次数 Times of downloads: 7)
，更加方便。具体如下，将实验数据改名为*.csv, 比如名字是511rr.csv, 511sr.csv, 当前目录如果有511rr_traj.xyz，511sr_traj.xyz，运行脚本
005_ana_ecd.sh
(5.7 KB, 下载次数 Times of downloads: 6)
，就会产生拟合图谱：
可根据需要调节10_ecd01.txt，10_ecd02.txt，10_ecd03.txt，sim_ecd.py中的参数，首次用sim_ecd.py时，需要修改
首行：#!/root/soft/anaconda3/envs/tjenv/bin/python（自己python的路径，笔者用的conda；若改脚本在window系统中使用，该行为空）
11行：def normalize_y(data, column_name='Y', scale=10): 10为纵坐标的高度，
66行 ：data1 = data1.iloc[35:445]（意思是提取实验数据35行-445行），
134行： input_folder = r"/root/6molclus112/"（改为自己目录的实际路径）。
最后，在10step文件中还有每一步的.sh文件，可以分别对每一步进行运算。在molclus112目录中还有一个01.sh脚本，还可将每步命令放入其中，分步测试。 脚本还保留了.pbs文件的部分命令行，如果使用openpbs作业提交管理功能，只需将00autorunall.sh 修改为00autorunall.pbs，并修改其中的节点就可使用。
5. 总结
对于该脚本00autorunall.sh 中每一行的含义，感兴趣者可以输入到chatGPT中去了解命令的具体含义，同时也可以按照自己需要进行适当修改和改进。
6. 更新提示
2025年3月，在使用的过程中持续的对脚本进行更新，文件12M，将文件拆开上传
6molclus112_part_ab
(4.84 MB, 下载次数 Times of downloads: 8)
和
6molclus112_part_aa
(7 MB, 下载次数 Times of downloads: 6)
，下载后在Linux中合并这两个压缩包cat 6molclus112_part_aa 6molclus112_part_ab > 6molclus112.zip，就可以正常解压。 如果网站下载出现问题，可以百度网盘下载
通过网盘分享的文件：6molclus112_202503sample.zip
更新的脚本如下： 6molclus112中放入由Chemdraw3D 保存的.pdb格式的结构，比如004rs.pdb，此外还需要放入ECD实验得到的.csv文件和实验的碳谱数据.txt，比如004rs_nmr.txt和004rs.csv，以便于对计算结果与实验结果进行分析，碳谱数据按照SVM方法的对比原理，对碳谱的实验值和计算值的相似度进行打分。脚本最明显的更新是可以直接得到发表文章和演讲展示所需要的图表、图片和数据。
脚本的整个运行过程包括：动力学构象搜索——xtb优化——gaussian 优化和振动分析——orca计算单点能——gaussian计算NMR和ECD——Multiwfn 和python脚本结合进行实验和计算值的比对——给出高清的3D结构图和ECD对比图——各个构象的能量和权重比例Excel——各个构象XYZ分子坐标word文档。最终结果可以直接放在PPT、文章或者支持材料部分，非常方便。如果有其他需求，理解脚本中的命令后，可以灵活调节。
使用这个更新的脚本之前，需要安装sobtop、acanda、obabel、python、Gromacs、vmd、ImageMagick等程序，DeepSeek非常有助于指导我们安装和使用这些程序。
仔细检查以上程序可用，就可开始计算了。 首选运行./11cal_nmr_ecd.sh这个脚本，脚本大体计算流程如下：obabel对004rs.pdb转换成004rs.mol2——sobtop计算电荷，为Gromacs准备动力学搜索文件——Gromacs动力学模拟搜索构象——VMD 提取搜索轨迹到004rs_traj.xzy——接下来进行优化等计算……..，计算完成，我们还会收到邮件提醒，会得到如下图文件。
上图是计算任务完成的高斯及相关文件
计算任务完成后，进行提取和分析数据，运行./12ana_nmr.sh，主要是提取三组基组泛函的NMR数据，我们可以进行DP4+的碳谱对比，确定分子相对构型。
上图是提取的计算NMR数据，有碳谱化学位移值和屏蔽常数值
上图是提取的是word格式的xyz文件
上图是提取的各个构象能量及权重
运行./13ana_ecd.sh，提取ECD实验和计算图谱的对比结果，各个构象的三维结构图，并保存到PPT中，我们可以打开PPT非常方便的浏览多个分子的三种基组泛函和实验的对比图，当我们有多个化合物时这个脚本用起来就相当的丝滑和优雅。
上图是提取到PPT中的实验和计算ECD对比图及结构图
上图是高清的3D结构图
运行./14ana_svm.sh，将实验的碳谱和计算的碳谱用SVM的对比方法打分，与DP4一样，辅助确定化合物的相对构型。 在使用过程中请注意这些脚本的路径问题，具体每个脚本的命令含义，DeepSeek讲解非常清楚。
欢迎使用2025年3月的更新，谢谢！！！！ 下面的内容就是2024年对脚本做的优化过程，感兴趣这个脚本的发展经过的朋友可以了解一下：
依然保存原来版本molclus112，增加了一个新的版本6molclus112
6molclus112.zip
(2.76 MB, 下载次数 Times of downloads: 42)
新版本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