在NAMD下利用Amber力场计算MD模拟的方法.doc

在NAMD下利用Amber力场计算MD模拟的方法曾经也是费了好的精力, 研究出了在NAMD中用Amber力场的方法. 呵呵，想起来也过去有些时候了，应winner_dk兄弟的要求， 我把这一过程在这里原原本本的呈献给大家，希望做过的朋友给挑挑毛病，正在做的朋友参考一下，准备做朋友借鉴一下，不会做的朋友学习一下，我也就不算白费功夫了！本贴不要钱，不设权限！好，废话少说，教程开始：我的实验条件：前端机 （自己的烂PC.巨烂，512的内存）XP系统， Cluster（46-cpu）Linux系统，提交作业用小企鹅，NAMD2.6, Amber 8实验对象：咱们在这里用HIV-1整合酶来说明操作过程1. 用Amber的tleap/xleap去生成top和 crd文件， 这两个文件的地位相当于我们在NAMD中熟悉的psf和pdb，所以呢，在计算之前一定要确保它们的正确性！否则后面的一切都是瞎扯！关于leap的使用请大家去自行研究AMber的手册。在这里只有一点提醒大家注意，Amber力场和Charmm对于原子的命名原则不同，所以最好不要用Charmm处理过的pdb去生成top/crd，否则是自讨苦吃，老手还好，新手可能搞好几天也过不去，决不是危言耸听。2. 结构优化，本来有两步的， 我只贴出第一步带约束的，第二步大家自己改吧。 咱们看配置文件吧， 以下的都是我的配置文件，久经考验，肯定能用，大家直接拿去贴到文件中就能用了，有一点说一下，我这文件中很多注释文不对题，我懒得改了，大家以鉴赏的眼光去看吧。# 结构优化 # this simulation do a 3 steps minimization#step1 fix all non-water,minimize water# ADJUSTABLE PARAMETERS #amber onparmfile ./IN_2B4J_loop_adjusted_MG.topambercoor ./IN_2B4J_loop_adjusted_MG.crd 看这里，这里要改成Amber的参数啊outputName ./outputs/IN_2B4J_loop_adjusted_MG_min_1set temperature 300firsttimestep 0# SIMULATION PARAMETERS # Input#paraTypeCharmm on#parameters ./pars/par_all27_prot_na.inp # NOTE: Do not set the initial velocity temperature if you 把Charmm的部分注释掉，或者干脆不写# have also specified a .vel restart file!temperature $temperature# Periodic Boundary conditionsprecise values 82.4379997253 57.8550033569 69.1809976101,however well set the box a little bit bigger than it.# NOTE: Do not set the periodic cell basis if you have also # specified an .xsc restart file! 121.847000122 95.2819976806 88.0310020447if 1 cellBasisVector1 82.68 0. 0. cellBasisVector2 0. 74.47 0. cellBasisVector3 0. 0. 69.33 cellOrigin 42.4445991516 38.8152618408 36.2214126587wrapWater onwrapAll on# Force-Field Parametersexclude scaled1-41-4scaling 1.0#Attention#cutoff 12.cutoff 10.# Attention! #switching on 还有这里，注意要把switchingoff了， Amber的要求switching off#switchdist 10.pairlistdist 13.5# Integrator Parameterstimestep 2.0;# 2fs/steprigidBonds all;# needed for 2fs stepsrigidTolerance 0.00001 ;# for Amber, 0.00000001 for Charmm 还有这里nonbondedFreq 1fullElectFrequency2stepspercycle 20#PME (for full-system periodic electrostatics)if 1 PME yesPMEGridSizeX 80PMEGridSizeY 80PMEGridSizeZ 72#Constant T/P contrl, useless here # Constant Temperature Control #langevin on ;# do langevin dynamics #langevinDamping 5 ;# damping coefficient (gamma) of 5/ps #langevinDamping 1langevinTemp $temperature #langevinHydrogen no ;# dont couple langevin bath to hydrogens # # # Constant Pressure Control (variable volume) # #useGroupPressure yes ;# needed for 2fs steps #useFlexibleCell no;# no for water box, yes for membrane #useConstantArea no;# no for water box, yes for membrane #if 1 #langevinPiston on #langevinPistonTarget1.01325 ;#in bar - 1 atm #langevinPistonPeriod100. #langevinPistonDecay 50. langevinPistonPeriod200. #langevinPistonDecay 100. #langevinPistonTemp $temperature # #restartfreq 5000 ;# 500steps = every 1psdcdfreq 5000xstFreq 5000outputEnergies 100outputPressure 100# Fixed Atoms Constraint (set PDB beta-column to 1)#(lipid and name P) or backbone or (segname B12 and noh)if 1 constraints onconsRef con_IN_2B4J_loop_adjusted_MG.pdb ;# 10 kcal/mol, backbone or name MGconsKFile con_IN_2B4J_loop_adjusted_MG.pdbconsKCol B #margin 1.0# EXTRA PARAMETERS # Put here any custom parameters that are specific to # this job (e.g., SMD, TclForces, etc.)# EXECUTION SCRIPT # Minimizationminimize 20000注意以上粗体字的地方，可能还有些Amber/Charmm直间不同的地方我没发现的欢迎大家指正！3. 动力学：可能很多人习惯什么平衡run，和production run的， 其实没什么必要， 两步搞在一起也没关系的， 我这里还是只贴一个MD的conf文件，剩下的大家自己改一下就好了。# Step 1 约束backbone、Mg离子,优化水 # this simulation do a 3 steps minimization#step1 fix all non-water,minimize water# ADJUSTABLE PARAMETERS #amber onparmfile ./IN_2B4J_loop_adjusted_MG.topambercoor ./IN_2B4J_loop_adjusted_MG.crdoutputName ./outputs/IN_2B4J_loop_adjusted_MG_md_1bincoordinates./outputs/IN_2B4J_loop_adjusted_MG_heating.restart.coorbinvelocities ./outputs/IN_2B4J_loop_adjusted_MG_heating.restart.velextendedSystem./outputs/IN_2B4J_loop_adjusted_MG_heating.restart.xscset temperature 300firsttimestep 0# SIMULATION PARAMETERS # Input#paraTypeCharmm on#parameters ./pars/par_all27_prot_na.inp # NOTE: Do not set the initial velocity temperature if you # have also specified a .vel restart file!#temperature $temperature# Periodic Boundary conditionsprecise values 82.4379997253 57.8550033569 69.1809976101,however well set the box a little bit bigger than it.# NOTE: Do not set the periodic cell basis if you have also # specified an .xsc restart file! 121.847000122 95.2819976806 88.0310020447if 0 cellBasisVector1 81.34 0. 0. cellBasisVector2 0. 79.25 0. cellBasisVector3 0. 0. 74.56 cellOrigin 42.1599960327 40.6278266907 38.5585861206wrapWater onwrapAll on# Force-Field Parametersexclude scaled1-41-4scaling 0.833333 ;#scee=1.2 is default in amber, 1-4scaling=1/scee, and the default 1-4scaling in name is 1.scnb 2 ;#2default in amber, 1-4 vdw interaction are diveded by this value.#Attention#cutoff 12.cutoff 10.# Attention! #switching onswitching off#switchdist 10.pairlistdist 13.5# Integrator Parameterstimestep 2.0;# 2fs/steprigidBonds all;# needed for 2fs stepsrigidTolerance 0.00001 ;# for Amber0.00000001 for Charmm nonbondedFreq 1fullElectFrequency2stepspercycle 20#PME (for full-system periodic electrostatics)if 1 PME yesPMEGridSizeX 80PMEGridSizeY 80PMEGridSizeZ 72#Constant T/P contrl, useless here # Constant Temperature Control #langevin on ;# do langevin dynamics #langevinDamping 5 ;# damping coefficient (gamma) of 5/ps #langevinDamping 1langevinTemp $temperature #langevinHydrogen no ;# dont couple langevin bath to hydrogens # # # Constant Pressure Control (variable volume) # #useGroupPressure yes ;# needed for 2fs steps #useFlexibleCell no;# no for water box, yes for membrane #useConstantArea no;# no for water box, yes for membrane #if 1 #langevinPiston on #langevinPistonTarget1.01325 ;#in bar - 1 atm #langevinPistonPeriod100. #langevinPistonDecay 50. langevinPistonPeriod200. #langevinPistonDecay 100. #langevinPistonTemp $temperature # #restartfreq 5000 ;# 500steps = every 1psdcdfreq 5000xstFreq 5000outputEnergies 100outputPressure 100# Fixed Atoms Constraint (set PDB beta-column to 1)#(lipid and name P) or backbone or (segname B12 and noh)if 0 constraints onconsRef con_IN_short_DNA_2MG_md_2.pdb ;# 2 kcal/mol, backbone and nucleic and (resid 316 to 320 329 to 333)consKFile con_IN_short_DNA_2MG_md_2.pdbconsKCol B margin 2.0# EXTRA PARAMETERS # Put here any custom parameters that are specific to # this job (e.g., SMD, TclForces, etc.)# EXECUTION SCRIPT # Minimizationrun 3000000看，其实和之前的优化部分没什么大区别。4. MD 之后的处理：为什么用Amber？主要是因为它能做小分子的参数和拓扑。 那么Amber还有什么优势吗？ 当然了！MM-PBSA/GBSA！ 这里我简要说一下MM-PBSA/GBSA的问题。首先你要做好复合物, 受体，配体的top/crd, 一共6个文件啊！然后把你的DCD文件转化为crd文件（用VMD），然后就可以用MM-PBSA/GBSA的脚本去算了，当然，你要有Amber，还要有耐心。那个玩艺虽然算得挺快，但是总是挂掉，大家要小心，不要提交了作业就跑了，否则你回来了发现它挂了，多郁闷！关于具体的操作，请大家去看AMber手册，我就不越俎代庖了，毕竟咱不是专业的Amber-or。5. 小结：此方法虽然可行，但是在做的过程中大家必然会遇到各种意想不到的情况（必然的），我也无法写出一个让大家一次就能成功地教程，毕竟软件本身也不是完美的。 大家如果遇到了什么搞不定问题，可以跟贴讨论。感谢版主加分!希望大家多参与讨论. 以上的conf文件可以应付一般的用Amb