计算化学公社
标题: 求助金颗粒与多肽动力学模拟添加金属原子力场问题 [打印本页]
作者Author: 不去重蹈 时间: 2021-3-5 20:43
标题: 求助金颗粒与多肽动力学模拟添加金属原子力场问题
老师您好，请教一个关于力场方面的问题。
背景：实验得纳米金颗粒能和多肽吸附（例如，CCY等肽），我们想在gromacs 2019.3中模拟计算纳米金颗粒和多肽CCY的结合能，比较不同的肽与金颗粒结合能力的强弱。
金原子的力场用 “在Gromacs中模拟金纳米线拉伸过程 http://sobereva.com/153 中 ” 的参数。在gromacs\top\gromos53a6.ff 路径下，将参数添加进ions.itp 文件中 f'f'nonbonded.itp
先是生成了CCY多肽的gro和itp文件。然后生成金颗粒的力场时遇到了报错，
gmx pdb2gmx -f .\Au13fixed.pdb -o .\Au13fixed.gro
报错
Fatal error:
Residue 'AU' not found in residue topology database
问题：适合多肽的力场很好选择，但是在构建Au13的力场时，请问应该怎么操作才能成功得到金的 top 和 itp 文件
附件：Au13的pdb文件。
作者Author: sobereva 时间: 2021-3-5 22:11
看我博文里的例子top里是什么样的。当前体系的金团簇部分明显不是直接靠pdb2gmx来产生拓扑信息的。
先拿pdb2gmx产生多肽的拓扑文件，然后自行把金的信息加入其中。
作者Author: 不去重蹈 时间: 2021-3-11 15:50
谢谢sobereva老师的指导，力场问题已经解决，由于gromos53a6力场不能识别部分氢原子（如HA等），改用amber99sb-lidn力场，做动力学后有两个问题。
1: 分子在穿过边界时，键被拉长了，请问周边应该怎么设置？
2: 金颗粒跑散了。请问金颗粒该怎么限制或固定，金颗粒才不散开？
附图，
第一帧，第n帧，最后一帧。动力学使用mdp文件
作者Author: sobereva 时间: 2021-3-12 02:17
1 用dynamic bonds方式显示，或者让每帧重新判断连接关系，看下文
VMD初始化文件(vmd.rc)我的推荐设置
http://sobereva.com/545（http://bbs.keinsci.com/thread-16834-1-1.html）
2 并不需要任何固定。跑散了肯定是拓扑文件有问题
作者Author: 不去重蹈 时间: 2021-3-12 11:16
谢谢sobereva老师，
请问拓扑文件改怎么修改？
作者Author: lyj714 时间: 2021-3-12 12:22
本帖最后由 lyj714 于 2021-3-12 12:27 编辑
我怀疑你金的参数不对，你这个top中也看不出你具体非键参数写的什么(ffnonbonded.itp)。如果你是照搬sob的那篇博文的参数那么你就错的离谱，因为sob写的那个参数只能用于gromos力场(用的是C6和C12)，你这里用的amber系力场，理应看那篇文献中的参数，进行合适的转换得到ε和σ
给你个本人根据那篇原文献转的各种参数，供参考:
https://github.com/liuyujie714/temp/blob/main/param.txt
作者Author: 不去重蹈 时间: 2021-3-12 14:46
谢谢老师！因为换用了力场，所以没有用sob老师的参数。我用的是朋友给的参数，查看ffnonbonded.itp，AU 79 11967 0.0000 A 2.93400e-01 1.63300e-01，应该不对。改用您给的再试一试。
请问金颗粒（Au13.pdb）文件中需要给Au原子之间添加键或者做什么约束吗？建模时13个金原子之间没有任何关系，是不是动力学时13个金原不是一个整体，是13个独立的单体？
作者Author: sobereva 时间: 2021-3-13 06:03
不需要任何约束
合理的参数下自然而然就能靠LJ势维持住簇的状态
作者Author: 牧生 时间: 2021-6-17 09:57
本帖最后由 牧生 于 2021-6-17 10:00 编辑
请问一下大佬，我用excel去推算你给这个txt的参数，从文本中抄过来的前三个参数，以Ag为例：
>> DEBUG:
12-6 LJ parameters from [JPCC 2008,112,17281–17290]:
metal r0 epsilon A B
Ag 2.955 4.56 2021000 6072,
Al 2.925 4.02 1577000 5035,
Au 2.951 5.29 2307000 6987,
For amber/charmm/cvff/oplsaa forcefield of gromacs:
>> DEBUG:
This article 12-6 LJ form:
E = eps*[(r0/r)^12 - 2*(r0/r)^6]
Convert it to below form for gromacs:
E = 4*epsilon*[(sigma/r)^12 - (sigma/r)^6]
That is, sigma = r0 * 2^(-1/6)
#metal sigma(nm) epsilon(kJ/mol)
Ag 0.263260571 19.079040000
Al 0.260587875 16.819680000
Au 0.262904212 22.133360000
我根据公式计算，大了十倍。
(, 下载次数 Times of downloads: 52)
作者Author: lyj714 时间: 2021-6-17 10:55
单位，gromacs用的是nm
作者Author: 牧生 时间: 2021-6-17 11:13
好的，谢谢。这下子清晰明白了
作者Author: 牧生 时间: 2021-6-19 21:57
大佬，请教一下，合理的参数可以使得金属只靠LJ势就维持住簇的状态，且在MD的过程中不跑散。
我目前做了个铁，CTAC，水合氢离子的水溶液接触，尽管铁没有跑散，但是其余的分子进入了铁晶胞内，这个能怎么解决啊。。
具体如下：
以ATB得到CTAB阳离子，水合氢离子的itp文件，
atomtypes.atp里面已经有铁的原子量 FE 55.845
γ铁的参数（该参数来自10.1038/s41524-020-00478-1）
(, 下载次数 Times of downloads: 69)
按照培训班的教程，计算得到c6 c12的参数，并写入添加到ffnonbonded.itp里面
FE 26 0 0 A 0.01514608 2.2873928e-06
延展铁板
gmxeditconf -f FE.pdb -o FE_plate_box.gro -box 3.656 3.656 6
往真空部分加入2个CTAC和10个H3O+
gmx insert-molecules -f FE_plate_box.gro -ci CTAB.pdb -oout_box.gro -nmol 2
gmx insert-molecules -f out_box.gro -ci H3O.pdb -box 4 46 -o out_box2.gro -nmol 10
加水填满盒子
gmxsolvate -cp out_box2.gro -p topol.top -o FE_wat.gro
得到tpr
gmx grompp -f em.mdp -c FE_wat.gro -p topol.top -o CTABa.tpr -maxwarn 80
加入反离子，
gmx genion -s CTABa.tpr -p topol.top -o CTABb4em.gro -nname CL -nn 12
铁晶在盒子中间，隐藏水以后得到的图像如下，
(, 下载次数 Times of downloads: 93)
能量最小化
gmxgrompp -f em.mdp -c CTABb4em.gro -p topol.top -o CTABb.tpr -maxwarn 80
gmxmdrun -deffnm CTABb -v -ntomp10 -ntmpi 1 -pinoffset 10 -gpu_id 0
铁晶体不再保持在盒子中央了，
(, 下载次数 Times of downloads: 55)
隐藏水分子以后的图像如下
(, 下载次数 Times of downloads: 68)
铁晶胞不在盒子中央，且非水的分子进入了铁晶体的内部。请教一下怎么解决啊。
使用的文件如下：
(, 下载次数 Times of downloads: 3)
(, 下载次数 Times of downloads: 3)
(, 下载次数 Times of downloads: 4)
(, 下载次数 Times of downloads: 6)
(, 下载次数 Times of downloads: 3)
(, 下载次数 Times of downloads: 3)
(, 下载次数 Times of downloads: 5)
(, 下载次数 Times of downloads: 9)
作者Author: lyj714 时间: 2021-6-20 11:24
应该是系统某些分子间非键参数有问题，不要忽略警告，看具体的关键性警告
作者Author: 牧生 时间: 2021-6-20 11:30
相同的方法，相同的命令，做金板，铝板和酸性水溶液就完全没有一点问题。唯独做铁板，水分子那些就进入铁内部。
作者Author: lyj714 时间: 2021-6-20 11:43
我说了让你看警告信息的，你怎么就不看。这么对比没意义
作者Author: 牧生 时间: 2021-6-20 13:32
对比了一下铝板和铁板与水溶液接触的过程，唯一的区别在于，使用铁板的做em的时候，多了一个note，
NOTE 1 [file em.mdp]:
You have set rlist larger than the interaction cut-off, but you also have
verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.
其余的警告都是因为使用gromos54a7_atb.ff力场而出现的警告，这个忽略即可。
那么，可能问题出在这个note上面，这个该怎么解决呢？verlet-buffer-tolerance和rlist是自动设定的吗？
作者Author: lyj714 时间: 2021-6-20 14:35
本帖最后由 lyj714 于 2021-6-20 14:39 编辑
就是这个力场的原因啊，非要让我点出问题，我引导你看警告不是没有原因的，如果你做铁，还自己添加了铁的原子信息到ffnbonded.itp中中，肯定会有一个atomtype defined previously警告，这就已经体现了错误地方。明显这个力场中已经添加了有一个铁原子类型了，就在ffnobonded.itp中，自己睁大眼睛去看，注释掉所有关于原先这个铁的所有行，包括预先定义的非键参数那些行。
所以我一直都不提倡那些一遇到警告都不仔细看完，就添加什么-maxwarn去忽略的做法，简直害人不浅。有些警告是必须注意的
作者Author: 牧生 时间: 2021-6-20 17:46
本帖最后由 牧生 于 2021-6-20 17:48 编辑
给大佬递冰可乐。
问题解决了。我输入参数是没有问题的，gromos54a7_atb.ff里面，ffnonbonded.itp只有我输入的参数，先前有的铁参数，我已经删掉的。
真正的原因在NOTE里面提到了，
NOTE 1 [file em.mdp]:
You have set rlist larger than the interaction cut-off, but you also have
verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.
于是，将em.mdp里面，将金板适用的这两个0.9，改成1，结果就很正常了。
rcoulomb = 0.9 改成1
rvdw = 0.9 改成1
作者Author: angervlf 时间: 2023-6-27 06:11
我不太理解，gromacs软件里自带的amber力场不也是已经经过转化的嘛，用的也是C6和C12吧，为什么还要进行转化？
作者Author: 牧生 时间: 2023-6-27 09:37
因为文献给出的参数不是标准的参数，所以要转化
作者Author: angervlf 时间: 2023-7-11 04:54
本帖最后由 angervlf 于 2023-7-11 04:57 编辑
是的是的，感谢，sob大佬的那篇博文后面没有看，现在理解了
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |