计算化学公社
标题: Auto-FEP for Gromacs(Shell 脚本) [打印本页]
作者Author: Huschein 时间: 2023-7-28 12:57
标题: Auto-FEP for Gromacs(Shell 脚本)
本帖最后由 Huschein 于 2023-8-1 19:45 编辑
Auto-FEP for Gromacs
引言：在分子动力学领域的学者都知道，自由能计算是一项非常有用的工具，它在药物设计、生物大分子相互作用分析等领域具有非常重要的地位。而自由能微扰（FEP）则是在经典物理领域（当然我们应该抛开量子化学去谈这件事）目前公认的最准确的计算自由能的方法，而Gromacs很早就实现了FEP的计算。虽然Gromacs为我们降低了很大的工作量，但是FEP所需要的批处理流程依旧十分繁杂，并且我最近的工作又恰好需要基于GMX实现的FEP，所以我在这里写了一些shell脚本来方便大家使用GMX做批处理的FEP计算。其中包括基于服务器的运算脚本（基于slurm调度系统）和基于个人电脑的脚本，但个人电脑脚本我尚未测试，如果出现问题请与我联系（邮箱：Hsuchein0126@outlook.com）。
Github链接：
https://github.com/Hsuchein/AutoFEP_For_GMX
准备工作：
1. Gromacs安装（版本新于5.0.0即可）
2. Alchemical-analysis安装，按照它官网教程即可（需要conda配置python2.7的环境），最后安装完后用我所提供的alchemical-analysis替换原文件的，因为源代码有问题，我已经修改过了。（https://github.com/alchemistry/alchemlyb），在做所有工作之前，务必看清楚这里！！！GMX在FEP任务中最后一步跑MD的时候会对每个模拟都输出一个dhdl.xvg这个文件,而alchemical-analysis是一个用于分析dhdl.xvg这个文件的一个python脚本，所以如果你想要处理FEP的输出，你务必要自己去它的官网，去配置好它所需的环境。而我目前的脚本还不支持直接一键自动分析，我目前的脚本是帮助你快速地跑10段FEP模拟，不用你自己手动调mdp文件，一个一个递交任务，所以我这个脚本本质上就是一个遥控器，同时开启10个模拟，没有太多的技术含量。
操作流程：
注意：全部操作都在主文件夹中完成，不需要去子目录中执行脚本命令！！！
1. 主文件夹配置（真的重要，看清楚！！！）
将你配置好的文件（一般建议是Gromacs中完成Genion指令后得到的文件，在此之前你需要自己完成设置盒子大小，添加溶剂，平衡电荷等操作）。脚本的原理就是：你先自己建模，因为模型都是一样的，FEP改变的只是mdp文件，所以你先建模（盒子大小，溶剂，电荷），全部建完之后把下一步（能量最小化）所需的文件全部丢进一个新文件夹。也就是说，需要你把此前最后一步所得的.top .itp .ndx .mdp .gro(注意自己手动修改mdp ndx等文件以适配自己的任务)等文件全部拷贝到新的主目录中，这里切记.gro文件只能有一个，不然脚本无法识别使用哪个.gro文件，然后如果你的力场是自己提供的，请一并复制.ff文件（example中有例子）到主目录中，然后将全部.sh脚本复制到你的主目录中（如果你是个人PC，那就用PC脚本；如果是slurm调度的服务器，就用servers脚本）。然后先把全部的东西都丢到新文件夹（下文也称之为“主目录”）后，再运行脚本，因为每段模拟都需要自己的指令，所以prepare.sh这个脚本，也会对其他脚本做一些修改，所以如果你不把其他脚本都加到主目录，你可能会遇到：
1.提示你找不到相应的.sh脚本
2.你后续任务递交出问题
2. 运行脚本
提示：！！！对于windows用户，bash脚本不能直接在终端运行，必须安装git-bash后用git-bash终端打开命令行，然后再用命令行运行脚本！！！
i)PC的话是一键运行，但是本人没有测试过，不知道是否会报错。并且PC是跑完一段模拟才能跑下一段，效率会低一些。
ii)servers脚本分为多段，我将逐个介绍
Prepare.sh：在你当前主目录下配置10个窗口的子文件夹目录（我设定的窗口是10，当然可以自己改,注意10个连接窗口会有11个模拟要跑），然后依次修改mdp文件的参数，并且复制文件到子文件夹目录。
Minimization.sh: 依次运行10个窗口的能量最小化任务,这个模拟时间极短,所以就不提交到服务器后台运算了,直接在当前终端运算即可。
Npt.sh/npt_sub.sh：npt_sub.sh会先检查每个窗口中min的模拟是否成功(是否输出min.gro)，如果成功则会执行gmx grompp的任务，随后调用npt.sh来递交任务到服务器后台运算，这里我给每个任务都写了JOB NAME，可以通过查看后台进程来看是否都运行上了（实测的时候有几段模拟会递交不上），如果出现递交不上的情况，那么直接进入相应的子目录运行npt.sh脚本即可。注意：npt_sub.sh中递交任务的部分（第25行）可能需要手动修改以适配自己的任务，同时也注意自己修改npt.sh中mdrun的命令以发挥最大的效能。
Md.sh/md_sub.sh：与npt同理
Collect.sh：在主目录下运行，它会先检查当前目录是否存在DHDL FILES这个文件夹，如果存在则清空里面的内容（所以有需要注意备份），如果不存在则创建，然后从每个窗口子目录中复制dhdl文件到该目录下，最后按照alchemical-analysis的脚本分析即可。
规划：目前打算规划把analysis也做成一键分析，看看有无时间吧。
备注：
1.其实这个shell脚本很简单，大家都可以写，但是也可以方便很多不会shell脚本又想要批量化处理模拟的同行们。
2.如果模拟有出现错误，我的建议是先把文件（你的.gro .ndx .top）名字改成和我example一样，尤其是ndx，因为我也没打算做一个非常普适的成熟脚本，因为这将需要我考虑各种情况，消耗巨大精力
版本修改纪要：
2023/07/24：AutoFEP4GMX 1.0
作者Author: 冰释之川 时间: 2023-7-28 14:03
本帖最后由 冰释之川 于 2023-7-28 14:17 编辑
太好了，又有搞FEP的工具出来了。但好像脚本不涉及到微扰图的产生以及dual 结构的对齐合并？
另外，请问楼主的脚本相比于 https://github.com/luancarvalhomartins/PyAutoFEP 中有什么优势嘛
(, 下载次数 Times of downloads: 67)
作者Author: Huschein 时间: 2023-7-28 14:53
pyautofep是一个大型项目，可以做各个配体之间的拓扑微扰图，并且可以计算ddg。但是我的autoFEP本身没有什么算法，只是构建了一个机器，可以串联运行Gromacs多个窗口任务，最后可以计算得到单个配体的绝对结合自由能，因为我不是做配体改造的，我只关心绝对结合自由能，本质上还是一个gromacs的运行脚本，最后用alchemlyb来分析输出。
作者Author: DwyaneWan 时间: 2023-7-28 15:50
大佬请问如果要算溶解自由能。 是不是只需要算相应分子在水溶液中的FEP就行了呢？
作者Author: Huschein 时间: 2023-7-28 21:43
理论上是的
作者Author: xptracy 时间: 2023-7-29 11:14
本帖最后由 xptracy 于 2023-7-29 11:54 编辑
请问“用我所提供的alchemical-analysis替换原文件的” 这个是指什么文件的替换，是这个alchemlyb-master\src\alchemlyb 文件夹吗
是把mdp、itp、gro都放到examples和mdpfile文件夹下，然后AutoFEP4PC.sh这个文件是在scripts_for_PC文件夹下运行还是？
作者Author: Huschein 时间: 2023-7-29 12:57
我的github里面有一个alchemical-analysis.py，你把他替换源代码alchemical-analysis.py，因为官方那个我不知道是什么问题，反正我用会报错，改了一下就不会了。然后你运行的话就是你自己新建一个文件夹作为主文件夹，然后各个脚本都在主文件夹下运行，包括所有的配置文件也先复制到主文件夹下
作者Author: xptracy 时间: 2023-7-29 14:59
github有这个alchemical-analysis.py文件，没找到要替换的源代码是哪个？
在建的文件夹下需要放什么文件？找不到文件和位置
(, 下载次数 Times of downloads: 66)
(, 下载次数 Times of downloads: 60)
作者Author: Huschein 时间: 2023-7-29 19:03
本帖最后由 Huschein 于 2023-7-29 19:14 编辑
已更新文章，如果还有不懂再问我，你这个存在几个显著性问题：
1.一看就知道配置主目录的时候没有把.sh文件加进去，那后面模拟全白搭
2.你当前目录下没有.gro结尾的文件，因为FIND指令找不到.gro结尾的文件
作者Author: xptracy 时间: 2023-7-29 21:22
本帖最后由 xptracy 于 2023-7-29 21:24 编辑
(, 下载次数 Times of downloads: 63)
(, 下载次数 Times of downloads: 64)
(, 下载次数 Times of downloads: 68)
(, 下载次数 Times of downloads: 66)
请教下您是这些文件吗，还有缺什么吗，有两个错误？
alchemical-analysis替换的原文件是在哪个位置？
作者Author: Huschein 时间: 2023-7-29 22:11
你去对应的子目录查有没有min.tpr这个文件，你这压根就没有执行minimization，2021版本我觉得gmx mdrun -deffnm是可以的，难道说非要-s指定tpr？ 你先看看有没有min.tpr吧，有的话可能我要改一下脚本适配更低版本的gmx，如果没有的话那根本就是你在prepare的时候出了问题，我的建议是，你自己先拿我的example去跑，确定是版本或者代码的问题再找我
作者Author: xptracy 时间: 2023-7-29 22:49
(, 下载次数 Times of downloads: 65)
嗯嗯，我把example的都放进来了，运行也是这两个错误
AutoFEP_For_GMX-main里没有min.tpr，21版之前数据是有生成min.tpr的
作者Author: Huschein 时间: 2023-7-30 09:10
你要去各个Windows窗口里看有没有min.tpr，不是在主目录里看
作者Author: xptracy 时间: 2023-7-30 09:28
(, 下载次数 Times of downloads: 71)
我是在主目录下运行bash AutoFEP4PC.sh，子目录下没有min.tpr
作者Author: Huschein 时间: 2023-7-30 18:30
留一个联系方式吧，qq或微信，我具体看看
作者Author: yyf123 时间: 2024-9-9 19:33
博主您好，最近想用fep计算配体和蛋白的结合自由能，请问gromacs的fep可以使用gpu加速吗，速度会损失多少, 谢谢您
作者Author: yyf123 时间: 2024-9-9 20:45
我还想请教一下，integrator = md和sd的区别在哪里，对精度的影响大吗
作者Author: Huschein 时间: 2024-9-14 02:25
md是正常动力学 sd是随机动力学，sd采样更广，这样可以覆盖窗口之间的空隙，你做完FEP是要算窗口之间重叠的，如果重叠的不够那就说明会有跳跃，最好是每个窗口能重叠上
作者Author: yyf123 时间: 2024-9-14 11:02
好的，谢谢您！
作者Author: Jeffrey0521 时间: 2024-10-21 16:47
博主您好，想问下这个fep的始末状态分别是什么呢？就是一个化合物凭空出现吗？另外，选择10个窗口的依据是什么呢，梯度又是参考什么来设置的呢？求解，十分感谢！！
作者Author: Huschein 时间: 2024-10-21 20:03
我现在就是算绝对结合能 始末就是没有到出现
窗口你可以参考文献，现在一般主流是9个或者11个
梯度是什么？
作者Author: Jeffrey0521 时间: 2024-10-22 12:59
我意思是各个window的范德华、库仑、bonded的变化快慢是怎么确定的？
作者Author: Huschein 时间: 2024-10-22 20:47
文献
作者Author: Jeffrey0521 时间: 2024-10-24 12:46
本帖最后由 Jeffrey0521 于 2024-10-24 12:48 编辑
但我查到的文件大多都是fep+做的，这些能作为参考吗？想问你具体参考的是哪些文献呢，无论怎么操作他的首尾窗口的邻接矩阵有的重合度也不太好
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |