科学网

 找回密码
  注册

tag 标签: gromacs

相关帖子

版块 作者 回复/查看 最后发表

没有相关内容

相关日志

上海超算
richor 2018-9-7 09:03
上海超算上gromacs-4.5.4版本的grompp有时候会在生成*.tpr文件的时候说, ThisrunwillgenerateroughlyxxMbofdata xx is a very large number. 这时候,考虑 在本地生成。 上海超算可能有问题。not very sure.
个人分类: 个人|1 次阅读|0 个评论
如何修改gromacs包中的力场
richor 2018-7-7 15:31
【静电修改】 直接去改top文件即可。 【vdw修改】 由于vdw是nonbond相互作用,不能通过直接修改top文件来达到。 需要在include**/forcefield.itp和 之间加入重新定义的原子类型。 比如修改 某原子类型 与 Na 离子的 vdw 相互作用(离子通道蛋白)。若不是所有的该原子类型全部都改,那么需要单独定义一个新的原子类型,进而需要重新定义与该原子相关的 ffnonbonded.itp 和 ffbonded.itp 中的相互作用参数。 一个例子脚本如下: add_bond.sh #needfileffbonded.itpkcsA_vdw.top,obtainadd.top lbond=6337 lpair=12298 langle=27727 ldihedral1=38499 ldihedral2=54241 lmap=55111 echoprocessingthekcsA_vdw.topfile awk-vlbond=$lbond-vlpair=$lpair'{if(NRlbond+1NRlpair-2)print}'kcsA_vdw.topbond.top awk-vlangle=$langle-vldihedral1=$ldihedral1'{if(NRlangle+1NRldihedral1-2)print}'kcsA_vdw.topangle.top awk-vldihedral1=$ldihedral1-vldihedral2=$ldihedral2'{if(NRldihedral1+1NRldihedral2-2)print}'kcsA_vdw.topdihedral1.top awk-vldihedral2=$ldihedral2-vlmap=$lmap'{if(NRldihedral2+1NRlmap-2)print}'kcsA_vdw.topdihedral2.top echoextractthechangelines forcxin899900901902906907908909236823692370237123752376237723783837383838393840384438453846384753065307530853095313531453155316 do grep$cxbond.topbond_change.top done forcxin899900901902906907908909236823692370237123752376237723783837383838393840384438453846384753065307530853095313531453155316 do grep$cxangle.topangle_change.top done forcxin899900901902906907908909236823692370237123752376237723783837383838393840384438453846384753065307530853095313531453155316 do grep$cxdihedral1.topdihedral1_change.top done forcxin899900901902906907908909236823692370237123752376237723783837383838393840384438453846384753065307530853095313531453155316 do grep$cxdihedral2.topdihedral2_change.top done mvbond_change.topbond.top mvangle_change.topangle.top mvdihedral1_change.topdihedral2.top mvdihedral2_change.topdihedral2.top echoobtaintheindexandnameofatom,thensubstituethenumberswithname lstart=71 lend=6336 awk-vlstart=$lstart-vlend=$lend'{if(NR=lstartNR=lend){if($1!=;)print$1,$2}}'kcsA_vdw.toprefer.top awk'BEGIN{printf(sed\\47)}{printf(s/%s/%s/g;,$1,$2)}END{printf(\\47bond.toptemp\\nmvtempbond.top)}'refer.toprefer.sh echodothesubstitutionofbond,mayneedalittletime bashrefer.sh awk'BEGIN{printf(sed\\47)}{printf(s/%s/%s/g;,$1,$2)}END{printf(\\47angle.toptemp\\nmvtempangle.top)}'refer.toprefer.sh echodothesubstitutionofangle,mayneedalittletime bashrefer.sh awk'BEGIN{printf(sed\\47)}{printf(s/%s/%s/g;,$1,$2)}END{printf(\\47dihedral1.toptemp\\nmvtempdihedral1.top)}'refer.toprefer.sh echodothesubstitutionofdihedral1,mayneedalittletime bashrefer.sh awk'BEGIN{printf(sed\\47)}{printf(s/%s/%s/g;,$1,$2)}END{printf(\\47dihedral2.toptemp\\nmvtempdihedral2.top)}'refer.toprefer.sh echodothesubstitutionofdihedral2,mayneedalittletime bashrefer.sh rmrefer.top rmrefer.sh echofinishthesubstituion awk'/CX|HX/{printf(%s%s,$1,$2);gsub(/CX/,CA);gsub(/HX/,HP);printf(%s%s\\n,$1,$2)}'bond.toptemp sort-ntemp|uniqbond.top awk'/CX|HX/{printf(%s%s%s,$1,$2,$3);gsub(/CX/,CA);gsub(/HX/,HP);printf(%s%s%s\\n,$1,$2,$3)}'angle.toptemp sort-ntemp|uniqangle.top awk'/CX|HX/{printf(%s%s%s%s,$1,$2,$3,$4);gsub(/CX/,CA);gsub(/HX/,HP);printf(%s%s%s%s\\n,$1,$2,$3,$4)}'dihedral1.toptemp sort-ntemp|uniqdihedral1.top awk'/CX|HX/{printf(%s%s%s%s,$1,$2,$3,$4);gsub(/CX/,CA);gsub(/HX/,HP);printf(%s%s%s%s\\n,$1,$2,$3,$4)}'dihedral2.toptemp sort-ntemp|uniqdihedral2.top echoextractingthelinesfromffbonded.itp lbondt=1 lconstr=187 langlet=227 ldihed1=701 ldihed2=1303 awk-vlbondt=$lbondt-vlconstr=$lconstr'{if(NRlbondt+1NRlconstr)print$1,$2,$3,$4,$5}'ffbonded.itpbond.itp awk-vlbondt=$lbondt-vlconstr=$lconstr'{if(NRlbondt+1NRlconstr)print$2,$1,$3,$4,$5}'ffbonded.itpbond.itp sort-nbond.itp|uniqtemp mvtempbond.itp awk-vlanglet=$langlet-vldihed1=$ldihed1'{if(NRlangletNRldihed1)print$1,$2,$3,$4,$5,$6,$7,$8}'ffbonded.itpangle.itp awk-vlanglet=$langlet-vldihed1=$ldihed1'{if(NRlangletNRldihed1)print$3,$2,$1,$4,$5,$6,$7,$8}'ffbonded.itpangle.itp sort-nangle.itp|uniqtemp mvtempangle.itp awk-vldihed1=$ldihed1-vldihed2=$ldihed2'{if(NRldihed1NRldihed2)print$1,$2,$3,$4,$5,$6,$7,$8}'ffbonded.itpdihedral1.itp awk-vldihed1=$ldihed1-vldihed2=$ldihed2'{if(NRldihed1NRldihed2)print$4,$3,$2,$1,$5,$6,$7,$8}'ffbonded.itpdihedral1.itp sort-ndihedral1.itp|uniqtemp mvtempdihedral1.itp awk-vldihed2=$ldihed2'{if(NRldihed2+1)print$1,$2,$3,$4,$5,$6,$7}'ffbonded.itpdihedral2.itp awk-vldihed2=$ldihed2'{if(NRldihed2+1)print$4,$3,$2,$1,$5,$6,$7}'ffbonded.itpdihedral2.itp sort-ndihedral2.itp|uniqtemp mvtempdihedral2.itp echoprocessingbondtypes echo add.top wc-lbond.toplt.txt lt=`awk'{print$1}'lt.txt` i=0 until ;do pair=`awk-vi=$i'{if(NR==i+1)printf(%s%s,$3,$4)}'bond.top` grep$pairbond.itptemp awk-vi=$i'{if(NR==i+1)print$1,$2}'bond.toptemp1 pastetemptemp1to_add awk'{print$6,$7,$3,$4,$5}'to_addadd.top i=$(($i+1)) done rmbond.top rmbond.itp echoprocessingangletypes echo add.top wc-langle.toplt.txt lt=`awk'{print$1}'lt.txt` i=0 until ;do pair=`awk-vi=$i'{if(NR==i+1)printf(%s%s%s,$4,$5,$6)}'angle.top` grep$pairangle.itptemp pair1=`cattemp` if ;then awk-vi=$i'{if(NR==i+1)print$1,$2,$3}'angle.toptemp1 pastetemptemp1to_add awk'{print$9,$10,$11,$4,$5,$6,$7,$8}'to_addadd.top fi i=$(($i+1)) done rmangle.top rmangle.itp echoprocessingdihedraltypes echo add.top wc-ldihedral1.toplt.txt lt=`awk'{print$1}'lt.txt` i=0 until ;do pair=`awk-vi=$i'{if(NR==i+1)printf(%s%s%s%s,$5,$6,$7,$8)}'dihedral1.top` grep$pairdihedral1.itptemp pair1=`cattemp` if ;then awk-vi=$i'{if(NR==i+1)print$1,$2,$3,$4}'dihedral1.toptemp1 pastetemptemp1to_add awk'{print$9,$10,$11,$12,$5,$6,$7,$8}'to_addadd.top fi pair=`awk-vi=$i'{if(NR==i+1)printf(X%s%sX,$6,$7)}'dihedral1.top` grep$pairdihedral1.itptemp pair1=`cattemp` if ;then awk-vi=$i'{if(NR==i+1)printf(X%s%sX\\n,$2,$3)}'dihedral1.toptemp1 pastetemptemp1to_add awk'{print$9,$10,$11,$12,$5,$6,$7,$8}'to_addadd.top fi i=$(($i+1)) done rmdihedral1.top rmdihedral1.itp echo add.top wc-ldihedral2.toplt.txt lt=`awk'{print$1}'lt.txt` i=0 until ;do pair=`awk-vi=$i'{if(NR==i+1)printf(%s%s%s%s,$5,$6,$7,$8)}'dihedral2.top` grep$pairdihedral2.itptemp pair1=`cattemp` if ;then awk-vi=$i'{if(NR==i+1)print$1,$2,$3,$4}'dihedral2.toptemp1 pastetemptemp1to_add awk'{print$8,$9,$10,$11,$5,$6,$7}'to_addadd.top fi pair=`awk-vi=$i'{if(NR==i+1)printf(X%s%sX,$6,$7)}'dihedral2.top` grep$pairdihedral2.itptemp pair1=`cattemp` if ;then awk-vi=$i'{if(NR==i+1)printf(X%s%sX\\n,$2,$3)}'dihedral2.toptemp1 pastetemptemp1to_add awk'{print$8,$9,$10,$11,$5,$6,$7}'to_addadd.top fi i=$(($i+1)) done rmdihedral2.top rmdihedral2.itp rmtemp1 rmlt.txt 这里面包含了awk,grep,sed以及bash的if, until, for语句,可以用来参考。 for循环的使用,完全可以替换掉until: http://gaodr.blog.163.com/blog/static/10461500820093793933887/ 建议使用类C的语法: for((i=0;i5;i++)) do ... done
个人分类: 分子模拟|3168 次阅读|0 个评论
leapfrog error
richor 2018-6-13 11:40
leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1 。 为什么? http://gromacs.org_gmx-users.maillist.sys.kth.narkive.com/cXo18Y84/dr-lemkul-s-umbrella-sampling-tutorial-grompp-note-on-leap-frog-nose-hoover 这里说,可以不用管的。
个人分类: 分子模拟|2 次阅读|0 个评论
氨基酸在PDB文件中的原子命名规则
supernovazjx 2018-6-10 15:23
做分子动力学模拟时,经常需要自定义氨基酸残基,或添加非氨基酸的残基,例如Fmoc基团等等。此时如果不清楚PDB中的残基命名规则,自己的残基就只能用C1,C2,H3,N4等等来命名,十分残念又难以记忆和理解,往往带来构建失败的结局。因此,了解PDB中的原子约定命名规则是很重要的。 笔者了解PDB的命名规则纯属偶然,乃观察PDB中下载的蛋白时偶然发现。 下面说明具体的命名规则。 氨基和羧基上的原子都采用本名,C, N, O, H, etc. 其它原子除 H 外,所有原子命名均采用“原子名+后缀 ”形式。 整体命名方法类似于图论中求解最大流问题时所采用的标号法。首先α-C被命名为CA。其后按照成键关系逐级递推,名字后缀依次为 B-G-D-E-Z-H-...(此为希腊字母表顺序:α(A),β(B),γ(G),δ(D),ε(E),ζ(Z),η(H)……)。以上图中的色氨酸为例,和CA相邻的原子是β-C,故被命名为CB。和CB相邻的原子是γ-C,命名为CG(先不考虑 H)。和CG相邻的有两个δ-C,排到了 D,此处将其命名为CD1和CD2。再向后和CD1相邻的是N原子,和CD2相邻的是两个C原子,故将它们命名为NE1、CE2和CE3。NE1再后面没有未标号的重原子了,停止;CE2和CE3后面是两个C原子,故命名为CZ2和CZ3(标号延续上一原子)。最后还剩下一个C,标为CH2。 对于氢原子,命名采用“H+所连原子后缀 ”形式。 例如连接在CD1上的H命名为HD1,连接在CB上的H命名为HB1和HB2, etc.。 对于自己构建的残基也是同样的道理,按照规则命名后,方便自行建立拓扑文件,即可进行后续的动力学模拟。
个人分类: 分子动力学模拟|15621 次阅读|0 个评论
[转载]多糖分子防止冰晶形成保持食品味道鲜美的机理研究
caixin5120 2018-5-21 09:55
多糖分子抗冻性研究 【分迪科技提供药物设计服务与平台】 期待与您合作! 公司网址: http://www.moldesigner.com/ 公司电话: 028-85160035 公司邮箱:sales@moldesinger.com 抗冻性在食品的保存和运输过程中非常重要的研究课题之一,研究发现一些多糖分子具有抗冻特性,能够防止冰晶结构对食品的挤压造成的细胞或组织损伤,进而保持食品或海鲜产品的鲜美味道。 传统实验方法只能从宏观上得到多糖分子是否具有抗冻性以及抗冻性强弱的结果。但是通过分子对接和分子动力学模拟方法可以从理论上解释多糖分子具有抗冻性质的内在原因: (1)通过分子对接技术可以得到多糖分子在冰晶表面的结合模式,分析结合的主要驱动力; (2)分子动力学模拟可以研究多糖使冰晶融化的过程,分析多糖分子在冰晶结构上的运动状态以及对冰晶结构的影响; (3)通过分析模拟过程中多种参数的变化情况,推测不同多糖分子抗冻性强弱的关系。 下图为抗冻分子(多糖或抗冻蛋白)与冰结合的分子动力学模拟结构:左边是完美的冰晶(模拟的初始条件),右边是经过20ns的模拟冰块表面被多糖融化的结构。 更多成功案例欢迎访问: http://www.moldesigner.com/category/company-case 本文版权属于成都 分迪科技 有限公司,转载请注明出处,商业使用需取得分迪科技书面同意! 扫描下方二维码关注 分迪科技 微信公众号 ,了解更多前沿资讯! \0 本文 转载自: http://www.moldesigner.com/1730.html \0
个人分类: 食品环境|1177 次阅读|0 个评论
classical MD or QM parallel CPU 效率
richor 2018-3-23 11:29
也就是多大的体系选择多少个核的问题(no GPU) 【gromacs 5.0】 benchmark : http://www.gromacs.org/@api/deki/files/240/=gromacs-5.0-benchmarks.pdf 可以看出 5 万原子最优值在 500 个核左右,也就是 100 atoms / thread ,差不多是最优的,实际上,我们的集群 communicate between nodes 要慢一些,那么 200 atoms /thread ,应该是可以的。 几个例子,参考一下(体系设置的比较好的情况): 体系:几个氨基酸( 90 原子),一个离子,其余为水, 6500 个原子, with PME ; 4 threads , 60 ns/day; 12 threads, 120 ns/day; 24 threads, 160ns/day(考虑到只有12个物理核心). ib.q的速度要更快一点。 体系:大蛋白( 1500 原子),若干离子,其余为水, 50,000 个原子, with PME ; 12 threads , 20 ns/day 。 体系:12-mer,4000水,15000原子,with PME; 20 ns/day, macbook air i5双核, same for hp63. 北京超算,元服务器,60ns/day, 20个核。 imac: 4 线程,竟然可以达到 30ns/day, for nvt. 【NAMD】 大体系效率太低,个人倾向于gromacs. http://www.wag.caltech.edu/home/duin/FFgroup/ff_meeting_2march2004.pdf 【cp2k】 benchmark suite: https://www.cp2k.org/performance 根据 suite 里的 H2O-256 ,得到的一个比较全的 test report : http://www.hpcadvisorycouncil.com/pdf/CP2K_Analysis_and_Profiling_AMD.pdf 64 个水, 200 个原子,用 HECToR ,最快 40s/5fs ,大概 1 ps/3 个小时,一天可以 10 ps 。最优处于 300 threads ,大概对应 1 atom/ thread 。 广州超算: 体系:若干个氨基酸( ~70 原子)和几个离子,其余为水, 500 个原子; 2 x 24 threads, 2 ps/day , 3x24 threads, 2ps/day, not improved. diis step, 3.6 ps; opt step:cg+ls Lenovo i5(v2.4.0): 200 s shh_super(v2.6.1): 7.5 s diis step, 3s, 2 nodes ; 5s, 1 nodes 【gaussian】 上海超算:单节点,24个核,体系 16个原子,2 min per opt step with b3lyp/6-31+g(d,p) 4个核i5 win下,体系25个原子,10 min per step with b3lyp/6-31+g(d,p) tcb,单个核, 30 min per step 组里centos服务器: 注意log文件末尾给出的是cpu time: xeon E5-2620 2.00 GHz: 单核2 min per cycle, 2核1 min per cycle, 6-31+g(d,p), 31个原子 6-31g(d,p), 31个原子,6核 10s per cycle
个人分类: 分子模拟|419 次阅读|0 个评论
Gromacs -- Extending Simulations
DaisyJing 2017-12-6 11:56
gmx mdrun -deffnm md_1 -cpi md_1.cpt -maxh 24 -append 经测试,如此续跑是在 md_1.cpt 的基础上延长MD,输出文件在原来的输出文件上继续写,比如log文件和xtc文件。 非常适用于模拟非正常结束和集群有运行时间限制的情况! http://www.gromacs.org/Documentation/How-tos/Extending_Simulations 1. Version 4 and Newer 1.1. Changing .mdp file options 2. Version 3.3.3 and Before 3. Exact vs binary identical continuation It is possible to extend or continue a simulation that has completed, and even those that have crashed (see Doing Restarts ). This is actually good technique for the handling of longer simulations to reduce time lost due to crashes of the computer(s) being utilised. It is only possible to continue seamlessly from a checkpoint file, or the point at which the last full precision coordinates and velocities are available for the system. Even though some coordinate file formats can contain velocities (e.g. .gro format ) these are not precise enough for an exact restart. Therefore, you have to ensure that full-precision data is written at sufficiently short lengths of time (several times per day?) to avoid loss of too much time due to crashes. How to achieve this varies with the version of GROMACS you are using. Version 4 and Newer A simulation that has terminated, but not completed (e.g. because of queue limits, power failure or the use of the -maxh option of mdrun ) can be continued without needing to use tpbconv (which is now called gmx convert-tpr as of version 5.0 ). You may or may not wish to use -append in your mdrun command in this case. Otherwise, a simulation that has completed can be extended using tpbconv , mdrun and checkpoint files ( .cpt ). First, the number of steps or time has to be changed in the .tpr file, then the simulation is continued from the last checkpoint with mdrun . This will produce a simulation that will be the same as if a continuous run was made (but see reproducibility for more discussion). tpbconv -s previous.tpr -extend timetoextendby -o next.tpr mdrun -s next.tpr -cpi previous.cpt You might want to use the - append option of mdrun to append the new output to the old files. Note that this will only work when the old output files have not been modified by the user. Appending is the default behavior as of version 4.5. If you would like to change the default filenames while running a lengthy simulation in manageable parts, then cyclically running commands such as the following will work: (with suitable values for name and time) mdrun -deffnm ${name} -cpi ${name} -append tpbconv -s ${name} -o ${name}_new.tpr -extend ${time} mv ${name}_new.tpr ${name}.tpr If you felt the need to archive your checkpoint and run input files, then you could do that, too. If you used -noappend , then mdrun will add numerical suffixes to a series of files based on your name, just as described in mdrun -h . When running in a queuing system, it is useful to set the number of steps you want for the total simulation with grompp or tpbconv and use the -maxh option of mdrun to gracefully stop the simulation before the queue time ends. With this procedure you can simply continue by letting mdrun read checkpoint files and no other tools are required. However, if your queueing system permits job suspension, the -maxh mechanism will be unaware of the time spent suspended, and you may simulate for less wall time than you would expect. The time can also be extended using the -until and -nsteps options with tpbconv . A simulation can be continued without the checkpoint file , which will be non-binary identical and will have small errors that, for most situations, are negligible. The reason for the errors is that the trajectory and energy files do not store all the state variables of the thermostats and barostats . If this is the case, you must make use of the version 3.3.3 procedure below. Changing .mdp file options If you wish/need to change .mdp file options, then either grompp -f new.mdp -c old.tpr -o new.tpr -t old.cpt mdrun -s new.tpr or grompp -f new.mdp -c old.tpr -o new.tpr mdrun -s new.tpr -cpi old.cpt should work. The former is necessary under GROMACS 4.x if the thermodynamic ensemble has changed. (Someone said If your old.cpt is for a run that has finished, then use tpbconv -extend after grompp and before mdrun . but mabraham disagrees. A run finishing is judged by the contents of the .cpt in the context of the .tpr. So, if the latter is changed, then the run isn't finished.)
个人分类: MD|3245 次阅读|0 个评论
模拟小窝(SimulationWorld)网络课堂介绍
WTianSD 2016-11-30 15:12
1. . 模拟小窝在网易云课堂推出系列分子动力学模拟的课程,包含基本理论、模拟软件、数据分析软件的介绍(目前包含分子模拟入门知识,LAMMPS软件,GROMACS软件,VMD软件)等等。课堂是一线科研人员精心制作的计算模拟技术入门和提升课程。 课程系列以引导如何自学为目标,解决学什么、怎么学的问题;课程以授“渔”为宗旨。入门课程旨在缩短新手的入门时间,降低入门的门槛,实现新人尽快掌握模拟所需要的基本知识和基本技能 课题地址: http://study.163.com/u/simuly 另外课程部分内容可以在优酷 模拟小窝 自频道观看: http://i.youku.com/i/UMjg3MDcwODM1Mg== 2. SimulationWorld微信公众号(扫描关注): 共享分子模拟的技巧,分享模拟经验和相关知识技能,打造国内利用分子模拟进行科学研究的交流平台;构建学习型互联网+社区,助力中国科学技术发展。
个人分类: 模拟技术|7363 次阅读|0 个评论
安装GROMACS with PLUMED
vessel 2015-9-11 09:07
PLUMED是一款计算free energy的免费软件。通常PLUMED作为一种plug-in软件,与常用的MD engine结合计算体系自由能。其中包含了多种常用方法如umbrella sampling, metadynamics等等。先介绍如何与gromacs结合使用。 gromacs是一款常用的MD软件,首先下载gromacs并安装,以gromacs-5.0.4为例。 1.首先安装PLUED 下载plumed-2.2b $ tar -xzvf plumed-2.2b-tar.gz $ cd plumed-2.2b $ ./configure --prefix=$HOME/install/plumed $ make $ make install 2.安装GROMACS,并patch with plumed $ tar -xzvf gromacs-5.0.4.tar.gz $ cd gromacs-5.0.4 $ plumed patch -p (or: patch.sh -p) $ mkdir build $ cd build $ cmake .. -DGMX_DOUBLE=ON -DGMX_X11=OFF -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_C_COMPILER=mpicc -DGMX_MPI=ON -DGMX_PREFER_STATIC_LIBS=ON -DCMAKE_INSTALL_PREFIX=$HOME/install/gromacs (cmake要求2.8之后版本;intel mpich 编译器) $ make $ make install test run: mdrun_mpi_d -version
个人分类: 科研|10634 次阅读|0 个评论
Installation instruction of gromacs-5.0.2
chrishengbee 2014-11-16 19:43
recently i installed an open-source molecular dynamics simulation software -gromacs in my desktop and laptop. since i met some trouble in my installation and i found some settings are very effective to the final performances, i decide to write it down. for ubuntu-14.04 64bit users: step1: install g++ sudo apt-get install g++ step2: install cmake here is the instruction of installing cmake in ubuntu: http://www.cmake.org/cmake/help/install.html it should be noted that gromacs requires cmake version 2.8.8 or higher. step3:install fftw noticed that in official installation guide, they mention that you can use cmake -DGMX_BUILD_OWN_FFTW=ON to download and build FFTW from source automatically for you, and i strongly recommend you to install fftw in this way instead of confirguring and making it manually. to be continued
个人分类: gromacs|0 个评论
Gromacs操作说明
niuyingli 2014-10-27 14:31
Gromacs 操作说明 ( TPBD 作为溶质分子, acetone 作为溶剂分子, amber 力场为例) 一、分别得到 2 种分子的含有 resp 电荷的 mol2 文件(以 acetone 分子为例) ① 优化结构(高斯优化 b3lyp ) #p b3lyp/6-31g*geom=connectivity opt ② g09-c01 版本得到 amber12 下的 antechamber 可以读入的 esp 电荷 # hf/6-31g*Guess=Read Geom=Check iop(6/50=1,6/42=17,6/41=10) pop=mk scf=tight acetone-esp ( 随便起个名字 ) 0 1 (电荷和多重度) antechamber.esp (注意:写上这个文件名) ③ 得到含 resp 电荷的 mol2 文件 antechamber -iantechamber.esp -fi gesp -o acetone.mol2 -fo mol2 -c resp -s 2 二、分别得到 2 种分子的 frcmod 文件(以 acetone 分子为例) parmchk -iacetone.mol2 -f mol2-o acetone.frcmod 三、利用 Leap 模块创建 prmtop 文件和 inpcrd 文件(需要每个分子的 frcmod 文件和 mol2 文件,为了区分两种分子,在 tpbd 的 mol2 文件中将所有的 MOL 改为 TOL ,以 acetone 分子为例) tleap -f tleap.sh 【 tleap.sh 的文件如下】 source leaprc.gaff mods=loadamberparamsacetone.frcmod MOL=loadmol2 acetone.mol2 check MOL saveamberparm MOL acetone.prmtopacetone.inpcrd quit 四、利用 dirac2 上的 acpype 工具将 amber 的力场文件( prmtop 文件和 inpcrd 文件)变为 gromacs 的力场文件( top 文件和 gro 文件) acpype -pacetone.prmtop -x acetone.inpcrd 五、建一个 6*6*6 的空盒子,把溶质分子 TPBD 放在盒子中心 editconf -fTOL_GMX.gro -box 6 6 6 -c -o tpbd_c_666.gro 【在 dirac2 上 vmdtpbd_c_666.gro ,就可以用 vmd 打开 gro 文件,选正交的显示方式,如下: 想要看到 tpbd 分子是不是在盒子中心,要输命令 Extensions → TK Console ,在里面输入 pbcbox 回车,如下所示 六、建立溶液盒子: 方法一:直接往里面注入丙酮分子,可用命令 genbox -cp tpbd_c_666.gro -ci MOL_GMX.gro -o test4.gro-nmol 1400 实际上只装进去 1299 个丙酮分子 方法二:先建溶剂盒子,再建溶液盒子 ① 建一个 2*2*2 的丙酮溶剂盒子(一般建模的时候先平衡溶剂,也可不平衡,后来以溶质和溶剂的整体来平衡) genbox -ci MOL_GMX.gro -box 2 -o sol222.gro -nmol 200 实际上只装进去 52 个丙酮分子 ② 将 2*2*2 的溶剂盒子装进 6*6*6 的溶质在中心的空盒子里(目的是把 3*3*3 个溶剂盒子装进去,密度比较平均) genbox -cptpbd_c_666.gro -cs sol222.gro -o tpbd_sol.gro 实际上装进去 1395 个丙酮分子, 1 个溶质分子在中心,如下图所示 七、将 2 个分子的 top 文件合并为 1 个 top 文件 ( 1 )将 2 个分子的 itp 文件的头和尾剪掉,剩下各自的主体 itp 文件并保存,如下所示: nrexcl = 3 stands forexcluding non-bonded interactions between atoms that are no further than 3bonds away. ( 2 )创建 2 个分子合起来的头文件 head.itp ,把重复的部分去掉; ( 3 )创建总的 top 文件 system.top ,将 head.itp , acetone.itp 和 tpbd.itp 加入进来,尾文件直接写在总的 top 文件内,如下所示: 八、能量最小化 $ ls acetone.itp em.mdp head.itp system.top tpbd.itp tpbd_sol.gro 以下是 em.mdp 的内容: ;title = cpp = /usr/bin/cpp define = -DFLEXIBLE integrator = steep nsteps = 1000 constraints = none ; ; Energy minimizing stuff emtol = 1000.0 ( 默认是 10) emstep = 0.01 nstcomm = 1 ns_type = grid rlist = 1.0 rcoulomb = 1.0 rvdw = 1.0 Tcoupl = no Pcoupl = no gen_vel = no 预处理,打包输入文件 grompp -psystem.top -f em.mdp -c tpbd_sol.gro -o em.tpr (产生 em.tpr 文件和 mdout.mdp 文件)如下所示: 同时会出现 1 个 NOTE: NOTE 1 : You are using aplain Coulomb cut-off, which might produce artifacts. You might want toconsider using PME electrostatics. 提交脚本以后,出现 confout.gro 、 ener.edr 、 traj.trr 和 md.log ( 脚本里的 jobname 要与 .tpr 前面的名称一致 ) 九、在 NPT 系综下跑 MD 使得溶液盒子合理化,先跑 1ns ,再跑 10ns 。以下以 1ns 为例: 将 em 文件夹下的能量最小化后生成的 confout.gro 重命名为 em.gro 以及 top 文件相关的文件拷过来, $ ls acetone.itp em.gro head.itp mdout.mdp md.tpr NPT_md.mdp system.top tpbd.itp NPT_md.mdp 预处理,打包输入文件 grompp -fNPT_md.mdp -p system.top -c em.gro -maxwarn 2 -o md.tpr 会报错: Fatal error: Too many warnings(1), grompp terminated. If you are sureall warnings are harmless, use the -maxwarn option. 同时出现 1 个 NOTE 和 1 个 WARNING: NOTE 1 : The Berendsenthermostat does not generate the correct kinetic energy distribution. You mightwant to consider using the V-rescale thermostat. WARNING 1 : You are not usingcenter of mass motion removal (mdp option comm-mode), numerical rounding errorscan lead to build up of kinetic energy of the center of mass. (WARNING 是说数值误差累积会引起质心平移,比如一开始是 0 0 0,10ns 以后变成 0 0 0.1 ,程序会自动挪回去, comm-mode=none 是指不限制质心平移,时间短几个 ns 是没关系的 ) 忽略掉 WARNING (默认是 0 , Maximum number of warning allowed ) : grompp -fNPT_md.mdp -p system.top -c em.gro -maxwarn 2 -o md.tpr (产生 md.tpr 文件和 mdout.mdp 文件)如下所示: 十、查看轨迹文件 dirac2 上 vmd confout.gro traj.trr ,可以查看每一步的轨迹,右三角可以播放和暂停 * 注意 : ① 如果初始的盒子建立的太密(或者初始构型不好),直接跑 NPT 的时候就会报错(查看 gromacs.out 最后), Fatal error: A charge groupmoved too far between two domain decomposition steps This usually meansthat your system is not well equilibrated. 这时有两种处理办法: ( 1 )重新建盒子,建的不要太密; ( 2 )先将盒子跑一个 NVT 平衡(把 NPT 的文件中压强耦合部分注释掉)或者 NVE ,拿出平衡后的构型跑 NPT 。 ② 如果跑 NPT 的时候,溶质分子跑到了一边,可以通过慢慢升温的办法,在温度和速度部分设定低的温度,慢慢升上去; ③ 摞盒子的办法切出一个球 genconf -f confout.gro-nbox 2 2 2 -o 8.gro 以下在 Graphics → Representaions 里用了涂层的办法, Create Rep 好几层,如果想让一层不显示的话,就双击变红色,即不显示的意思,可以在方框内写“ residue 0” 或“ index 0to 49 ),然后 Apply 显示单个残基或者多个残基。 Vmd 中,想看到显示的画面中的分子 MOL 的编号,可以按下数字 1 ,然后点那个分子即可,在旁边显示 MOL15 ,对应的是 Residue14,vmd 的 residue 和原子编号 index 都是从 0 开始。 以下表格表示从 2*2*2 的盒子里选出目标溶质分子,首先找出所有溶质分子,只有下层的 1 个溶质分子能当作目标分子, 1396 个分子为原胞中的分子数目。 2792 ( 8376 ) 4188 ( 9772 ) 0 ( 5584 ) 1396 ( 6980 ) 保存切出的球,先选中构型文件,变绿色;其中 T 表示 Top ,如果想令其他的为 T ,那么就点一下 T ,变红色即不是 Top 的意思, D 表示 Display ,如果不显示按 D ; File → Save Coordinates… ,在 Selected atoms 的方框旁边的▼,下拉,选择 sameresidue as within 30 of residue 8376 → Save (保存为 pdb 文件即可) 也可以用 tcl 命令保存: ④ 创建索引文件,定义组 make_ndx ‐ f confout.gro ‐ o index.ndx 可见系统自动分为了 3 组,可以选择 1 或者 2 到产生的 .ndx 文件中,下面的选项可以根据比如原子编号进行分组,可以直接输入 a 1-50 ,则 1-50 号的原子就会出现在索引文件 .ndx 下面了。 ⑤ 盒子边长有没有变化看 confout.gro 最后的三个数字 ⑥ 用 g_energy 进行看下 NVT 是否平衡,可以选择总能量、温度、压强 则会将上面几项写入 energy.xvg 文件中, 在 dirac2 上用 gnuplot 打开,可以看到 10ps 的过程,时间 - 能量( kJ/mol )、时间 - 温度( K )、时间 - 压强( bar, 1bar=10 5 Pa )关系图分别如下: 一般从几个 ns 以后开始分析,前面波动太大(压强不需要考虑) ⑦ rdf 分析看溶液微观结构 ( 1 ) g_rdf -ftraj.trr -s md.tpr -o rdf.xvg reference group 可以选 TOL , group 可以选 MOL ,表示在 TOL 附近 MOL 出现的几率;如果有尖峰说明溶液分层。 ( 2 ) g_rdf -ftraj.trr -s md.tpr -cn rdf_cn.xvg 看 0~r 内在 TOL 附近出现的 MOL 的数目。 ⑧ rmsd 分析溶液是否平衡 g_rms -f traj.trr-s md.tpr -o rmsd.xvg
个人分类: Gromacs|0 个评论
Gromacs中添加CHARMM力场
热度 1 wuyuexunxian 2013-8-20 16:09
CHARMM 力场提供了很多生物分子的力场参数,如蛋白质、脂质、核酸等, Gromacs 中自带了 CHARMM27 力场,当然也可在 Gromacs 官网上下载到已转换成 Gromacs 支持的 CHARMM36 力场。不过使用 Gromacs 做 MD 模拟的朋友如果想模拟一个 Gromacs 中没有定义力场参数的分子,如糖,就需要手动为 Gromacs 添加相关力场参数了。这些力场参数多数从文献中得到,并且一般都是 CHARMM 格式的,即 .prm 和 .rtf 格式,如何将它们转换成 Gromacs 支持的格式呢?这篇文章会详细讲解转换过程,并且给出 python 脚本以供参考。 首先分别讲解一下 CHARMM 软件和 Gromacs 软件如何完整定义一个分子的力场。在 CHARMM 软件中,主要通过 2 个文件来实现: .rtf 和 .prm 。 .rtf 文件定义了原子类型、各类型原子的质量、原子偏电荷、分子的键连接以及氢键等,它的格式如图 1 。 MASS 是用来定义原子类型的,每一行都定义了一种原子类型。 图 1. TIP3 水分子的 .rtf 文件, CHARMM 中 “ ! ” 后面的内容是注释 MASS 后面跟着原子编号、原子类型名、原子质量和元素符号。 GROUP 语句用来定义一个残基(也可以是单独的分子): ATOM 定义了残基中各原子的名字、所属的原子类型和该原子所带的偏电荷; BOND 定义残基中原子的键连接情况,以原子对形式列出,如图 1 中的水分子定义了 3 个键 OH2-H1 、 OH2-H2 和 H1-H2 (最后一个键是 SHAKE 限制算法用到的); ANGLE 定义键角; DONOR 和 ACCEPTOR 分别用来定义氢键供体和受体。 .prm 文件中定义了与能量相关的参数,如键伸缩能,键角、二面角的力常数以及范德华参数,图 2 给出了 TIP3 水分子的 .prm 文件。 BONDS 数据块是用来定义键伸缩振动能量的,用简谐振动来近似: BONDS 数据块 每行有 4 个字段,前两个是成键原子的类型,后两个分别是公式中的 K b ( kcal/Å 2 /mol )和 b 0 ( Å ) 。 ANGLES 数据块用来定义键角振动,同样是用简谐振动来近似: ANGLES 数据 每行有 5 或 7 个字段,后者是加了 Urey-Bradley 项的。前 3 项是形成键角的原子的类型,后面 4 项分别对应公式中的 K θ ( kcal/ 度 2 /mol )、 θ 0 (度)、 K UB ( kcal/Å 2 /mol ) 和 b 1-3,0 ( Å )。 DIHEDRALS 数据块用来定义二面角振动,用余弦函数来近似: DIHEDRALS 数据块 每行有 7 个字段,前 4 项是形成二面角的原子的类型,后面 3 项分别对应公式中的 K φ ( kcal/mol ) 、 n (无单位)和 δ (度)。 图 2. TIP3 水分子的 .prm 文件 IMPROPER 数据块用来保证分子中原子共平面的,用简谐振动近似: IMPROPER 数据块 每行有 7 个字段,前 4 项是形成二面角的原子的类型,第 5 项对应公式中的 K ω ( kcal/ 度 2 /mol ),第 6 项一般是 0 ,暂时不知道干嘛用的,第 7 项对应公式中的 ω 0 (度)。 NONBONDED 数据块用来定义非键作用,其中静电力可通过原子偏电荷求得,不需要额外参数,范德华力通常用 L-J 形式来描述: 通常力场中只定义单个原子的 ε 和 r min ,成对原子的范德华力参数可以通过组合规则( combination rule )获得,常用的组合规则是 ε i,j = sqrt( ε i × ε j ) , r min i,j = (r min i + r min j )/2 。另外 CHARMM 力场将 1-4 作用单独分出来描述,也采用 L-J 形式,只是 ε 和 r min 值有所不同。 NONBONDED 数据块每行有 7 个字段组成,第 1 个字段是原子类型,第 2 、 5 字段标识为 ignored ,通常是 0 ,暂不知道干嘛用的,第 3 , 4 字段分别为 vdw 力的 ε ( kcal/mol )和 r min /2 ( Å ),第 6 , 7 字段分别为 1-4 作用的 ε ( kcal/mol )和 r min /2 ( Å )。这里可以提前说明一下, Gromacs 中 1-4 作用是以原子对的形式定义的,叫做 pairtypes ,转换时需要将 CHARMM 中的 1-4 作用通过组合规则两两组合生成 pairtypes ,这在后面还会详细讲到。 下面讲一下 Gromacs 中力场文件的格式。 Gromacs 主要通过 5 个文件来定义: .rtp 、 .hdb 、 atomtypes.atp 、 ffbonded.itp 和 ffnonbonded.itp 。 .rtp 文件用来定义残基(或分子)中原子所属的类型、原子偏电荷、键连接以及共平面原子,其余如键角和二面角是在这里不是必须的,程序一般通过原子类型生成,它将出现在 ffbonded.itp 文件中。图 3 是 .rtp 文件的一个示例,里面的注释很详尽,唯一需要说 图 3. 甘氨酸的 .rtp 文件, Gromacs 中 “ ; ” 后的内容是注释 明的是 chargegroup 项,同一个 chargegroup 中的原子偏电荷变化是以相同比例进行的,对于一个新加入的分子可简单的将每个原子划分到不同的 chargegroup 中。 .hdb 文件是用来给残基(或分子)加氢的,我们知道用 X- 衍射得到的蛋白结构中一般没有氢, pdb2gmx 命令会查询 .hdb 文件中定义的加氢规则为蛋白加氢。如果你的小分子中氢原子已在 .rtp 文件中定义了,就不需要用 .hdb 文件来加氢。因为一般力场中都会详细定义蛋白各种参数,我们只需添加一些小分子,而小分子中的氢原子多在 .rtp 文件中定义,很少用到 .hdb 文件,所以这里不再讲 .hdb 文件格式。 atomtypes.atp 文件中定义了力场中用到的所有原子类型,格式很简单,第一列是原子类型名,第二列是原子质量,如果你添加的分子用到了新的原子类型,就需要在这里添加相关信息。 ffbonded.itp 文件定义了力场中的键作用:键伸缩振动、键角、二面角振动( func=9 ) 图 4. ffbonded.itp 文件格式 以及 IMPROPER (在 Gromacs 中称为 dihedraltypes ( func=2 ))。它们定义及算法与 CHARMM 中相同,只是公式的形式以及单位有差异: 公式的差异表现在键长、键角、 IMPROPER 的简谐振动公式在 Gromcas 中多了一个系数 1/2 ;单位的差异表现在 CHARMM 中能量单位为 kcal ,而 Gromcas 中使用 kJ , CHARMM 中长度单位为 Å ,而 Gromacs 中使用 nm 。参数从 CHARMM 转换到 Gromacs 规则如下: (1) 键伸缩振动: K gromacs = 2 × K charmm × cal2j × 100 (即 K gromacs = 836.8 × K charmm ,其中 cal2j =4.184 是 cal 转换到 J 的转换系数); b 0 gromacs = b 0 charmm /10 。 (2) 键角振动: K θ gromacs = 2 × K θ charmm × cal2j (即 K θ gromacs = 8.368 × K θ charmm ); θ 0 不变; K UB gromacs = 2 × K UB charmm × cal2j × 100 (即 K UB gromacs = 836.8 × K UB charmm ); b U gromacs = b U charmm /10 。 (3) 二面角振动: K φ gromacs = K φ charmm × cal2j (即 K φ gromacs = 4.184 × K φ charmm ); n 和 φ 0 不变 。 (4) IMPROPER:K ξ gromacs = 2 × K ξ charmm × cal2j (即 K ξ gromacs = 8.368 × K ξ charmm ); ξ 0 不变 。 值得注意的是这些参数在文件中位值不同,一定要弄清哪一列放哪个参数! ffnonbonded.itp 定义了力场中的非键作用,因静电力可通过原子偏电荷和 coulomb 公式求的,所以这个文件主要定义范德华力(在 atomtypes 数据块中)。与 CHARMM 类似, Gromacs 也对 1-4 作用进行单独处理,处理方式与 CHARMM 完全相同, 1-4 作用定义在 pairtypes 数据块中。 Gromacs 中描述范德华力的 L-J 公式与 CHARMM 中也有差异: 在转换之前我们需将 CHARMM 和 Gromacs 中的 L-J 公式化成相同的形式: 非键作用的参数转换规则如下: ε gromacs = - ε charmm × cal2j (即 ε gromacs = -4.184 ×ε charmm ,不知道问什么 CHARMM 中加了负号) σ gromacs = 2 ×σ charmm /10/2 1/6 (即 σ gromacs = 0.1781797436 ×σ charmm , CHARMM 中记录的是 r min /2 ) 对于 1-4 作用,转换规则同上,只是 Gromacs 中通过组合规则(上文讲 CHARMM 时有提到)以原子对的形式记录在 pairtypes 数据块中,这需要自己计算添加。图 5 给出了一个转换示例,注意 pairtypes 是通过两两组合计算得到的,不区分先后,共 6 对。 图 5. 非键作用转换示例 对于简单分子可以手动转换并添加到 Gromacs 力场中,但是当分子较为复杂时手动转换将会变得非常繁琐并且容易出错,这时候就需要用脚本来进行批量处理。附件提供了用 python 写的转换脚本: cvt_bd.py: 键伸缩项转换,输入文件 bonds ,输出文件 bonds.out cvt_agl.py: 键角转换,输入文件 angles ,输出文件 angles.out cvt_dihedral.py: 二面角转换,输入文件 dihedrals ,输出文件 dihedrals.out cvt_improper.py: improper 项转换,输入文件 impropers ,输出文件 impropers.out cvt_nb.py: 非键作用转换,输入文件 nobonded ,输出文件 vdw.out, pair.out 使用方法: 将 .prm 文件按数据块分成 bonds (记录 BONDS 数据块)、 angles (记录 ANGLES 数据块)、 dihedrals (记录 DIHEDRALS 数据块)、 impropers (记录 IMPROPER 数据块)和 nonbonded (记录 NONBONDED 数据块) 5 个文件,然后按下面步骤处理: ( 1 )删除关键字 BONDS 、 ANGLES 、 DIHEDRALS 、 IMPROPER 和 NONBONDED 所在的行。 ( 2 )删除“!”及其后面的注释内容,这可以在 vim 中使用命令: :%s/!.*// ( 3 )删除行末多余的空格,可用 vim 命令: :%s/ *$// ( 4 )删除空行, vim 命令: :g/^$/d 处理后的文件应该是每行具有相同列数( angles 和 nonbonded 各行列数不全相同),格式一致的文本。然后运行上面的 python 脚本,将输出文件拷贝到 ffbonded.itp 和 ffnonbonded.itp 文件中对应的数据块内: ( 1 ) bonds.out 内容拷贝到 ffbonded.itp 文件 数据域中 ( 2 ) angles.out 内容拷贝到 ffbonded.itp 文件 数据域中 ( 3 ) dihedrals.out 内容拷贝到 ffbonded.itp 文件 ( func=9 )数据域中 ( 4 ) impropers.out 内容拷贝到 ffbonded.itp 文件 ( func=2 )数据域中 ( 5 ) vdw.out 内容拷贝到 ffnonbonded.itp 文件 数据域中 ( 6 ) pair.out 内容拷贝到 ffnonbonded.itp 文件 数据域中 ( 7 ) vdw.out 的 1 、 3 列内容拷贝到 atomtypes.atp 文件中 最后根据 .rtf 文件中记录的残基(或分子)原子类型、偏电荷、键连接和 IMPROPER 信息,按照 Gromacs 的 .rtp 文件格式制做该残基(或分子)的 .rtp 文件,将该文件放到相应的力场文件目录下(即 ffbonded.itp 所在目录)。最后提醒大家,分子的 pdb 文件一定要与 .rtp 文件原子顺序和命名完全一致! 以上所有操作都正确的话就可以用 Gromacs 模拟该分子了, Happy simulating ! convert_script.rar
个人分类: 分子动力学模拟|22535 次阅读|1 个评论
GROMACS 4.6 beta1 released
热度 1 albumns 2012-11-30 15:48
From Gromacs community We've finally got the first beta release of GROMACS 4.6 ready for you to try out! We've put a lot of very hard work into it, and we hope you'll like the good things we've done. Things won't be perfect yet, so we'll be looking forward to your help finding the things we haven't done well enough yet! Remember, if you want the big performance gains that will be available in 4.6, then you'll want to know things will build and work well on your hardware, and the best way of doing that is helping us over the next few weeks. At the same time, we discourage you from doing work with this code whose scientific reliability you need to trust - this is very much a draft version of the software! You can find the manual here ftp://ftp.gromacs.org/pub/manual/gromacs-manual-4.6-beta1.pdf and the source code here ftp://ftp.gromacs.org/pub/gromacs/gromacs-4.6-beta1.tar.gz It would be great for us if some of you want to try out the new code on lots of different hardware and operation systems and report build problems, inconsistencies, strange or lacking documentation and in worst case pure bugs. To tempt you to do so here's a bit of a carrot corresponding to the new features: * A brand-new native GPU implementation layer. Gromacs now does heterogeneous parallalization using both CPUs and modern NVIDIA GPUs at the same time, the GPU port also works in parallel using both multiple cards in a node or multiple nodes, and it's smoking fast. There's lots of heroic work by Szilard Pall and Berk Hess here, and special thanks to NVIDIA and Mark Berger for their assistance in making this happen. * Gromacs can now use OpenMP parallelization for better scaling inside nodes, in particular when doing the FFT part on the CPU while the GPU does the normal nonbonded interactions. * Automatic load balancing between direct-space and PME nodes, and lots of improvements in domain decomposition load balancing and scaling. * We have a brand new set of classical nonbonded interaction kernels, and Gromacs can now use either SSE2, SSE4.1, 128-bit AVX with FMA support (AMD) or 256-bit AVX (Intel), all of them in both single and double precision. The performance difference depends on your system and parallelization, but it is quite large in many cases - we have seen 40% improvement on ion channels running on modern AMD machines! Did we mention that the classical C kernels are faster too since we can now do force-only interactions for most steps? * There are new kernels using analytical switch/shift functions that are quite a bit faster, and a new CPU-implementation of verlet kernels that guarantee buffered interactions (no atoms drifting in/out of the neighbor list range) that conserve energy extremely well. * There is a large new module to do advanced free energy calculations, thanks to Michael Shirts. Trust us, you need the full manual to decipher all the possibilities… * Gromacs has switched completely to CMake for configuration and building. To be honest, we do expect some hiccups from this, but it has enabled us to provide much more automation and advanced features as part of the setup - and Gromacs now works on Windows out-of-the-box. Please test as many parts of the build system as you can! * All raw assembly has been replaced by machine intrinsics in C. This does wonders for readability, but it means the compiler and compiler flags matter. On x86, you will typically get 5-10% better performance from icc than gcc. Happy simulating! The GROMACS development team
个人分类: 好文转载|3463 次阅读|1 个评论
Dedicated stategies for GPCR MD simulation
albumns 2012-7-27 02:43
As the most important membrane protein, GPCR MD simulation is a quite hot topic nowadays. However, how to assemble protein/membrane system and which forcefield to assign for the whole system are two main critical and tough task for GPCR simulations. 1. Assemble protein/membrane system Many tools are now available for protein/membrane building including: VMD (http://www.ks.uiuc.edu/Research/vmd/) CHARMM-GUI (http://www.charmm-gui.org/) Desmond System Builder (http://www.deshawresearch.com/resources_desmond.html) g_membed (http://wwwuser.gwdg.de/~ggroenh/membed.html) InflateGro (http://moose.bio.ucalgary.ca/index.php?page=Translate_lipdis) VMD can support well for POPC and POPE membrane system building with CHARMM27 and CHARMM36 FF. H owever, one have to add additional solvent and i ons into the syst em by tcl script from VMD tutorial. It will also n eed some script to merge membrane and protein. Moreover, since the lipids are not pre-equilirated, it would be necessary to equ ilibrate the whole system at least 20 ns before MD production. CHARMM-GUI aims to provide more convenient way for NAMD or CHARMM sim ulation. It can gene rate a embeded protein/memb rane system with OPM position through web page interface. It even can helps to assi gn CHARMM CGFF for the ligand. H owever, the re are also some o bvious week nees for it: there are so many atom clashed between lipids that we can hardly believe the lipids are pre-eq uilib rated which claimed by the author; the input file provided by CHARMM-GUI is not good enough for membrane protein simulation s ince the whole quilibration step on ly contains no m ore than 3 ns and obvious GPCR helix movement are of ten observed w ithin such short time which shouldn't expected at this time scale level. So, one have to improve the protocol by himself. Desmond System builder tool is incorporated in Schrodinger Maestro GUI and it provides very friendly interface for users. It can build a OPM based position for pr otein/membrane system very easily by c licking some bottoms . It can also assign CHARMM36 FF by VIPARR tool in D esmond. Both g_membed and InflateGro are tools w ithin Groma cs and both of them can embed the prot ein into a pre-equlibrated membrane system which save lot of ti me for equilibration . Although g_membed is a little bi t diff icult than InflateGro, but the output seems to be much better than the later one. 2. Force filed It is said that CHARMM36 FF is the best FF for lipids which i s currently the only FF can reproduce lipid gel phase property. Howe ver, recen t Lipid 11 FF from latest Amber 12 is also claimed to be as good as CHARMM36 FF , al though related p aper is being revi ewed th ese days. Both full atom FF are quite good option for protein system simulation. There are also other FF and me thods including Gromos FF which is a united atoms FF and nowadays coarse gain MD which use dum my sphere to represent groups and acce lerate the simulation dra matically (24 core workstation can even achieve up to several microsecond/day ). The demerit for those methods are also obvious: we gain w hat we paid. 3. Top ology for Ligand T his is always a head ache problem for many people in MD simulation. In Desmond, ligand to pology can be recognized automatically with OPLS_2005 FF. H owe ver, OPLS_2005 is only go od e nough for t ens of ns MD simulation, it is rather poor if su b mit to micro se cond MD. If we would like to use CHARMM36 FF for protein/membrane system bound with ligand, we can generate to pol ogy from S wissparam (http://swissparam.ch/) and convert them into Desmond Viparr format by script from Desmond 2012 ($SCHRO DINGER/desmond-v31023/data/viparr/converters ). This tool sometimes doesn't work well, it is said that it can be full y supported by D esmond in the next version of Desmond which would be released at the beginning of next year . What's n eed to mention is that CHARMM C G FF is also reachable in Des mond 2012 ($SCHRODING ER/desmond-v31023/data/viparr/ff/cgenff_base_v2b7 ), one can build manually with those molecul ar templates if the target on e is not so complicated. CHARMM CG FF also could be obtained from (https://www.paramchem.org/). If the ligand structure is not so complicated, it may work well. H owever, it some time s may not re cognize ligand bond order and so on correctly. In this way, one have to go to CHARMM forum for helps. Amber GAFF is definitely ext re mely at tractive and the primary choice for a ligand bound system. When the latest LIPID 11 FF in Amber 12 come out, Ambe r should be the first choic e for many people especially those work with ligand. 4. Efficiency Alt hough the hardware of compu ter deve lop s so fast nowada ys that CPU update one generation almost each year, the e ffici ency of MD simulations seems don't imp rove so much these days. F or instan ce: no matter how many CPU we use , for a typical mem brane protein simulation (13 2 lipids, 300aa protein, 50,000 at oms in all ) with full atom FF and typical cutoff (9-10) with PME: Gromacs can get up to 2 0 ns/day (double precision), Amber 12 ns/day, NAMD 4ns/day. Desmond is an exception since the parralization is much more superior tha n any other MD tools, it can up to 100 ns/day with 512 CPU. I t can even up to several microseco nd /day in A nton with full atom FF. It woul d be a wise option to use either D esmond or Gromacs , if one would like to run hundre ds of ns with full atom FF. Of course, it is also accepta ble for Amber and NAMD CPU performance if the simulation on ly last for tens of ns. GPU technology is developing quite fast in recent one or two years and it also bring exc iting news for computational work especially NVIDIA CUDA accel erations. F or instance, with two GTX590, A mber 12 can get up to 20ns/day while 24 core i7 3.6 GH z CPU can only get 4ns/day ( with intel com piler , gnu is even m uch slower). NAMD on the other hand, can get 5ns/day with CUDA acc eleration w hile 24 core i7 3.6 GHz CPU can only get 0.5 ns/day. C urrently, GPU calculation is not supported in Desmond, but this feature is e x pected to be available in next version.
个人分类: 科研笔记|7374 次阅读|0 个评论
compile gromacs GPU version from source code
热度 1 albumns 2012-4-24 14:37
1. download source code: git clone git://git.gromacs.org/gromacs.git 2. install openMM and CUDA library into /soft/gromacs-gpu/openmm (download from https://simtk.org/project/xml/downloads.xml?group_id=506 ) /soft/gromacs-gpu/cuda configuring .cshrc as following etenv OPENMM_ROOT_DIR /soft/openmm setenv LD_LIBRARY_PATH /soft/gromacs4.6-gpu/cuda/lib64:/soft/gromacs4.6-gpu/cuda/lib:/soft/gromacs-gpu /openmm/lib 3. compile GROMACS CUDA from source code mkdir cuda cd cuda cmake ../. -DGMX_OPENMM=ON -DGMX_THREADS=OFF -DCMAKE_INSTALL_PREFIX=/soft/gromacs4.6-gpu make -j24 make install 4. configure .cshrc source /soft/gromacs4.6-gpu/bin/GMXRC.csh done!
5639 次阅读|2 个评论
GPU加速的Gromacs 4.5.1 进行分子动力学模拟
yuanhui80 2011-6-4 11:56
计算机技术的快速发展,算法及相应软件的不断更新,使得当前利用我们手上的普通电脑来模拟相对较大的生物大分子体系和多聚体分子成为了可能,而且这种趋势会越来越明显,尤其是多核心CPU的出现及分子动力学模拟大规模并行化计算能力的提高,让从事生物学研究的人们有可能利用手中的个人电脑对感兴趣的蛋白分子进行有目的的模拟,并充分与生物学实验有机地结合在一起,这是一件非常有意思和好玩的事情。 尽管如此,许多生物大分子体系还是非常巨大的,比如我的一个体系:腺病毒六邻体(Hexon)三聚体蛋白,我想对其进行温控的分子动力学模拟,以动态分析其总表位构成、高温变性机制及高温变性在免疫原性上的反应。该体系共含有约940*3 = 2820个氨基酸残基,再加上一个立方体的水盒子,总体系约200,000个原子数,在QX9650四核心CPU,64位Linux系统下,Gromacs每模似10ns的时间要花上约27天,非常耗时,在一个约500ns的总体设计中,这种计算能力是无法忍受的。 GPU加速的Gromacs为我们带来了非常振奋的好消息,官方称利用Nvidia的CUDA技术,可以将MD模拟提高原单CPU的十倍以上,以下是我利用Nvidia GTX460 2G 进行分子动力学模拟的全过程,现拿出来和大家进行分享。 第一天: Nvida GTX460 2G大显存 按gmx网站( www. gromacs .org/gpu )上的说明,可以 模拟 200,000左右的 原子 ,我模拟的 体系则好 是190,000 硬件:Dell T3400工作站 X38主板 CPU: QX9650 4G ECC内存 GTX460 2G显存 软件 :Ubuntu9.10 64位, CUDA3.1, OpenMM2.0, FFTW3.2.2, CMake, Gromacs 4.5.1, 按照官网要求,独立安装CPU的Gromacs4.5.1(CMake 编译 ),再 下载 预编译好的 mdrun-gpu beta2 for Ubuntu9.10 64位 设好环境变量,运行~ 但是运行后,提示: Fatal error : reading tpx file (md.tpr) version 73 with version 71 program For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors 大概的意思是说:预编译好的mdrun-gpu跑不了由4.5.1版 grompp 程序 生成的tpr 文件。 第二天: 采取从头编译的方法解决了上述问题,因为预编译好的mdrun-gpu与4.5.1里的程序版本号不同,所以会出现不兼容现象, 按照提示,顺利编译4.5.1版的mdrun-gpu成功, —————————————————————————— export OPENMM_ROOT_DIR=path_to_custom_openmm_installation cmake -DGMX_OPENMM=ON make mdrun make install-mdrun —————————————————————————— 但是新的问题来了, 运行出现错误提示: mdrun-gpu: error while loading shared libraries: libopenmm_api_wrapper.so: cannot open shared object file: No such file or directory 很奇怪! 环境变量也设好了,没有问题 在openmm目录下找不到libopenmm_api_wrapper.so文件 第三天: 我将操作系统换成RHEL5.5系统 再利用相同的安装方法,顺利解决上述问题, 但不明白其中原因,不过我想还是有办法解决的(先不管它)! 第四天: 总结一下我所遇到的问题,及解决办法: 1,版本问题 —————————————————————————— Fatal error: reading tpx file (md.tpr) version 73 with version 71 program For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors 这里是说版本不兼容 2,Openmm不支持多组的温度耦合 —————————————————————————— Fatal error: OpenMM does not support multiple temperature coupling groups. For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors 3,不能按以前的mpd来设置 —————————————————————————— Fatal error: OpenMM uses a single cutoff for both Coulomb and VdW interactions. Please set rcoulomb equal to rvdw. For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors 4,GPU加速的gmx现在只支持amber力场及charmm力场 —————————————————————————— Fatal error: The combination rules of the used force-field do not match the one supported by OpenMM: sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j). Switch to a force-field that uses these rules in order to simulate this system using OpenMM. 5,GPU加速的gmx不支持G96里的 interaction ,实际上还是力场问题 —————————————————————————— Fatal error: OpenMM does not support (some) of the provided interaction type(s) (G96Angle) For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors 6,在Ubuntu9.10里面用cmake编译gromacs4.5.1会遇到找不到libopenmm_api_wrapper.so文件的问题,换成RHEL5.5可以解决 —————————————————————————— error while loading shared libraries: libopenmm_api_wrapper.so: cannot open shared object file: No such file or directory 第五天: mdrun-gpu终于跑起来了,mdp文件是用的官网提供的 bench里面的,不过还是有一些warning: It is also possible to optimize the transforms for the current problem by performing some calcula- tions at the start of the run. This is not done by default since it takes a couple of minutes, but for large runs it will save time. Turn it on by specifying optimize_fft = yes WARNING: OpenMM does not support leap-frog, will use velocity-verlet integrator. WARNING: OpenMM supports only Andersen thermostat with the md/md-vv/md-vv-avek integrators. Pre-simulation ~15s memtest in progress...done, no errors detected starting mdrun 'Protein in water' 1000000 steps, 2000.0 ps. NODE (s) Real (s) (%) Time: 33.080 99.577 33.2 (Mnbf/s) (MFlops) (ns/day) (hour/ns) Performance: 0.000 0.074 47.609 0.504 gcq#330: "Go back to the rock from under which you came" (Fiona Apple) ———————————————————————————————— 最后这个Performance,有点看不懂,单从(ns/day) 这一点看,性能是l四核心CPU的五倍,但实际运行,性能仅是CPU的2倍, (MFlops) 一项,竟然是 0.074 CPU的 (MFlops) 是12GFlops 总体上看,GPU加速的GMX,性能提升,至少可以达到传统四核心CPU的2倍, imp模型官网上说可达到10倍以上,继续更新中。。。。。。 第六天: 190,000个原子的体系,共设了10ns,performance显示是5 ns/day,理论上两天就算完了, 实际上得到28号才算完(10月18号下午1点开始),这个结果和performance明显不符~~~ 5000000 steps, 10000.0 ps. step 417300, will finish Thu Oct 28 10:39:59 2010 /10月18号开始,显示10月28号结束 Received the TERM signal, stopping at the next step step 417378, will finish Thu Oct 28 10:39:46 2010 Post-simulation ~15s memtest in progress...done, no errors detected NODE (s) Real (s) (%) Time: 13633.960 71173.931 19.2 3h47:13 (Mnbf/s) (MFlops) (ns/day) (hour/ns) Performance: 0.000 0.003 5.290 4.537 gcq#47: "I Am Testing Your Grey Matter" (Red Hot Chili Peppers) 第七天: 相同的体系,相同的设置,以下是用QX9650 四核心CPU 跑的performance: 性能虽然比GTX460弱,但也只是多算了三天时间 —————————————————————————————— Back Off! I just backed up md.trr to ./#md.trr.1# Back Off! I just backed up md.edr to ./#md.edr.1# WARNING: This run will generate roughly 3924 Mb of data starting mdrun 'Good gRace! Old Maple Actually Chews Slate in water' 5000000 steps, 10000.0 ps. step 0 NOTE: Turning on dynamic load balancing step 500, will finish Mon Nov 1 09:57:48 2010vol 0.74 imb F 2% /10月19号早上开始,显示的是11月1号结束 Received the TERM signal, stopping at the next NS step step 550, will finish Mon Nov 1 10:08:34 2010 Average load imbalance: 2.4 % Part of the total run time spent waiting due to load imbalance: 1.1 % Steps where the load balancing was limited by -rdd, -rcon and/or -dds: Y 0 % Parallel run - timing based on wallclock. NODE (s) Real (s) (%) Time: 123.856 123.856 100.0 2:03 (Mnbf/s) (GFlops) (ns/day) (hour/ns) Performance: 156.395 11.835 0.769 31.220 gcq#358: "Now it's filled with hundreds and hundreds of chemicals" (Midlake) 第八天: 跑官网上的bench: GTX460 的成绩是102ms/day,与c2050并没有想象的那么大差距! Pre-simulation ~15s memtest in progress...done, no errors detected starting mdrun 'Protein' -1 steps, infinite ps. step 285000 performance: 102.1 ns/day Received the TERM signal, stopping at the next step step 285028 performance: 102.1 ns/day Post-simulation ~15s memtest in progress...done, no errors detected NODE (s) Real (s) (%) Time: 481.290 482.224 99.8 8:01 (Mnbf/s) (MFlops) (ns/day) (hour/ns) Performance: 0.000 0.002 102.335 0.235 总结: 新一代GPU加速的Gromacs分子动力学模拟为我们展示了GPU将来在分子动力学领域应用美好前景,但目前还不成熟。从以上测试中我们可以看出在隐性溶剂水模型的MD模拟计算中,GPU加速的计算性能是传统四核心CPU的至少10倍以上,但是在显性溶剂水模型中,GPU加速未见得多么明显;另外最为重要的一点是,当前版本Gromacs 4.5.1 对于GPU加速MD计算有很多限制,如支持力场有限,许多特性还不支持,模拟的可重复性差(与CPU模拟相比)等,不过在足够长的模拟时间下,还是会生成重复性较好的具有统计学意义的模拟轨迹。相关资料请参考: www.gromacs.org/GPU
12854 次阅读|0 个评论
Gromacs 4.0.7并行 含QM/MM功能安装
热度 1 albumns 2010-4-11 04:17
Gromacs 4.0.7并行带QM/MM安装 平台 SUSE Linux Enterprise Desktop 10 SP3 gcc4.1.2 mpich2-1.2.1p1 ifrot 10.1 fftw 3.2.2 解压缩mpich2-1.2.1p1.tar.gz进入此目录,运行: ./configure make make install 运行 touch /etc/mpd.conf chmod 700 /etc/mpd.conf 将下面加入mpd.conf: secretword=secretword (比如secretword=ltwd) 解压fftw3.2.2压缩后进入目录,安装到/soft/fftw下 ./configure --enable-float --enable-threads make make install 把libmopac.a复制到/soft/fftw/lib和/lib下 配置环境变量 setenv CPPFLAGS -I/soft/fftw/include setenv LDFLAGS -L/soft/fftw/lib 解压gromacs4.0.7,进入目录 ./configure --prefix=/soft/gromacs --enable-mpi -enable-fortran --with-qmmm-mopac --enable-shared make make install 配置环境变量 setenv LIBS -lmopac setenv LD_LIBRARY_PATH /soft/gromacs/lib source /soft/gromacs/bin/completion.csh set path=(/soft/gromacs/bin $path) configure中的其他选项 CC C compiler command 一般这个环境变量就是gcc CFLAGS C compiler flags 编译时的参数,一般是-O3 LDFLAGS linker flags, e.g. -Llib dir if you have libraries in a nonstandard directory lib dir 库文件目录 LIBS libraries to pass to the linker, e.g. -llibrary 设的时候不用-l xxx,无需引号 CPPFLAGS C/C++/Objective C preprocessor flags, e.g. -Iinclude dir if you have headers in a nonstandard directory include dir F77 Fortran 77 compiler command 一般这个环境变量就是gfortran或ifort FFLAGS Fortran 77 compiler flags 编译时的参数,一般是-O3 CCAS assembler compiler command (defaults to CC) CCASFLAGS assembler compiler flags (defaults to CFLAGS) CPP C preprocessor CXX C++ compiler command 一般是g++ CXXFLAGS C++ compiler flags CXXCPP C++ preprocessor XMKMF Path to xmkmf, Makefile generator for X Window System Optional Features: --disable-FEATURE do not include FEATURE (same as --enable-FEATURE=no) --enable-FEATURE include FEATURE --enable-shared build shared libraries --disable-float use double instead of single precision --enable-double same effect as --disable-float --enable-fortran use fortran (default on sgi,ibm,sun,axp) --enable-mpi compile for parallel runs using MPI --disable-threads don't try to use multithreading --enable-mpi-environment=VAR only start parallel runs when VAR is set --disable-ia32-3dnow don't build 3DNow! assembly loops on ia32 --disable-ia32-sse don't build SSE/SSE2 assembly loops on ia32 --disable-x86-64-sse don't build SSE assembly loops on X86_64 --disable-ppc-altivec don't build Altivec loops on PowerPC --disable-ia64-asm don't build assembly loops on ia64 --disable-cpu-optimization no detection or tuning flags for cpu version --disable-software-sqrt no software 1/sqrt (disabled on sgi,ibm,ia64) --enable-prefetch-forces prefetch forces in innerloops --enable-all-static make completely static binaries --disable-dependency-tracking speeds up one-time build --enable-dependency-tracking do not reject slow dependency extractors --enable-static build static libraries --enable-fast-install optimize for fast installation --disable-libtool-lock avoid locking (might break parallel builds) --disable-largefile omit support for large files Optional Packages: --with-PACKAGE use PACKAGE --without-PACKAGE do not use PACKAGE (same as --with-PACKAGE=no) --with-fft= FFT library to use. fftw3 is default, fftpack built in. --with-external-blas Use system BLAS library (add to LIBS). Automatic on OS X. --with-external-lapack Use system LAPACK library (add to LIBS). Automatic on OS X. --without-qmmm-gaussian Interface to mod. Gaussian0x for QM-MM (see website) --with-qmmm-gamess use modified Gamess-UK for QM-MM (see website) --with-qmmm-mopac use modified Mopac 7 for QM-MM (see website) --with-gnu-ld assume the C compiler uses GNU ld --with-pic try to use only PIC/non-PIC objects --with-tags include additional configurations --with-dmalloc use dmalloc, as in http://www.dmalloc.com/dmalloc.tar.gz --with-x use the X Window System --with-motif-includes=DIR Motif include files are in DIR --with-motif-libraries=DIR Motif libraries are in DIR --without-gsl do not link to the GNU scientific library, prevents certain analysis tools from being built --with-xml Link to the xml2 library, experimental
个人分类: 科研笔记|8665 次阅读|3 个评论

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-6-2 12:22

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部