HTMD在准备蛋白体系方面比较智能有明显的优势。本文在前面的基础上提供一个构建膜蛋白体系的实例。
序
ACEMD是一款针对NVIDIA的GPU进行动力学优化的引擎,HTMD是基于Python的分子可编程环境,用于准备、处理、模拟、可视化和分析分子系统。 HTMD在准备蛋白体系方面比较智能有明显的优势。本文在前面的基础上提供一个构建膜蛋白体系的实例。 Equilibration protocol for GPCR #!/usr/bin/python3 # coding: utf-8 # In[1]:导入模块,加载文件 from htmd.ui import * from htmd.home import home #get the files shutil.copytree(home()+'/data/mor','/tmp/testmor/pdb') os.chdir('/tmp/testmor') path='./01_prepare/' # In[2]:列出目录下文件 %ls /tmp/testmor/pdb # In[3]:准备体系的坐标和拓扑 #Protein 4dkl is taken from opm topos = charmm.defaultTopo() + ['pdb/ff.rtf'] params = charmm.defaultParam() + ['pdb/ff.prm'] prot = Molecule('pdb/4dkl.pdb') prot.filter('protein and noh and chain B or water within 5 of (chain B and protein)') pcenter = np.mean(prot.get('coords','protein'), axis=0) prot = autoSegment(prot, sel='protein') prot = charmm.build(prot, topo=topos, param=params, outdir= path+'prot',ionize=False) # In[4]:开启可视化,可视化体系 from htmd.config import config config(viewer='ngl') prot.view() # In[5]:向体系中添加离子 #Add sodium in the receptor sod = Molecule('pdb/sod.pdb') sod.set('segid','S1') prot.append(sod) # In[6]:蛋白体系嵌入膜 #Use a POPC membrane created with vmd and C36 memb = Molecule('pdb/membrane80by80C36.pdb') mcenter = np.mean(memb.get('coords'),axis=0) memb.moveBy(pcenter-mcenter) mol = prot.copy() mol.append(memb, collisions=True) # Append membrane and remove colliding atoms # In[7]:膜蛋白体系中加入配体 #Add ligand, previously parametrized using gaamp lig = Molecule('pdb/QM-min.pdb') lig.set('segid','L') lcenter = np.mean(lig.get('coords'),axis=0) newlcenter = [np.random.uniform(-10, 10), np.random.uniform(-10, 10), 43] lig.rotateBy(uniformRandomRotation(), lcenter) lig.moveBy(newlcenter - lcenter) mol.append(lig) # In[8]:体系置于水盒子 #Add water coo = mol.get('coords','lipids or protein') m = np.min(coo,axis=0) + [0,0,-5] M = np.max(coo,axis=0) + [0,0,20] mol = solvate(mol, minmax=np.vstack((m,M))) # In[9]:建立体系的charmm力场相关拓扑和参数文件 mol = charmm.build(mol, topo=topos, param=params, outdir=path+'/build', saltconc=0.15) 运行代码块产生的体系相关坐标和拓扑文件。 # In[10]:平衡体系 from htmd.protocols.equilibration_v2 import Equilibration md = Equilibration() md.runtime = 100000 md.temperature = 300 md.fb_reference = 'protein and resid 293' md.fb_selection = 'segname L and noh' md.fb_box = [-39, 10, -29, 21, 47, 50] md.fb_k = 5 md.useconstantratio = True md.write(path+'/build',path+'/equil') # In[11]:可视化体系 # Visualize the flat bottom potential box from htmd.config import config config(viewer='ngl') mol.view('not water') # In[12]:运行平衡 mdx = AcemdLocal() mdx.submit(path+'/equil') mdx.wait() 运行代码,产生的提交平衡任务的脚本以及平衡的输入参数文件和结果文件。 运行代码,产生的提交动力学任务的脚本以及动力学的输入参数文件和结果文件。 # In[14]:可视化动力学模拟结果 mol.view('not water') 参考资料 https://software.acellera.com/docs/latest/htmd/tutorials/mu-opioid-receptor-gpcr-equilibration.html .https://software.acellera.com/docs/latest/htmd/userguide/introduction.html https://software.acellera.com/docs/latest/htmd/userguide/building.html http://gainstrong.net/works/hudong/