Preparation of the \(\mu\) opioid receptor with ligand

This is a complex build system as it has several components, the protein, a sodium ion, the ligand and of course the membrane.

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/'
Please cite HTMD: Doerr et al.(2016)JCTC,12,1845. https://dx.doi.org/10.1021/acs.jctc.6b00049

HTMD Documentation at: https://www.htmd.org/docs/latest/

You are on the latest HTMD version (unpackaged : /shared/sdoerr/Work/htmdacellera/htmd).
%ls /tmp/testmor/pdb
4dkl.pdb  ff.prm  ff.rtf  membrane80by80C36.pdb  QM-min.pdb  sod.pdb

Build

#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)
2017-02-20 10:33:51,140 - htmd.molecule.molecule - INFO - Removed 5120 atoms. 2262 atoms remaining in the molecule.
2017-02-20 10:33:51,274 - htmd.builder.builder - INFO - Created segment P0 between resid 65 and 263.
2017-02-20 10:33:51,275 - htmd.builder.builder - INFO - Created segment P1 between resid 270 and 352.
2017-02-20 10:33:51,673 - htmd.builder.charmm - INFO - Writing out segments.
2017-02-20 10:33:52,150 - htmd.builder.builder - INFO - One disulfide bond was added
2017-02-20 10:33:52,258 - htmd.builder.charmm - INFO - Starting the build.
2017-02-20 10:33:52,326 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 1 atoms due to bad angles.
2017-02-20 10:33:52,327 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 78 atoms.
2017-02-20 10:33:52,327 - htmd.builder.charmm - WARNING - Please check /tmp/testmor/01_prepare/prot/log.txt for further information.
2017-02-20 10:33:52,328 - htmd.builder.charmm - INFO - Finished building.
Bond between A: [serial 3005 resid 140 resname CYS chain B segid P0]
             B: [serial 3615 resid 217 resname CYS chain B segid P0]
prot.view()
#Add sodium in the receptor
sod = Molecule('pdb/sod.pdb')
sod.set('segid','S1')
prot.append(sod)

#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

#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)

#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)))

#Build
mol = charmm.build(mol, topo=topos, param=params, outdir=path+'/build', saltconc=0.15)
2017-02-20 10:41:20,248 - htmd.molecule.molecule - INFO - Removed 311 residues from appended Molecule due to collisions.
2017-02-20 10:41:20,896 - htmd.builder.solvate - INFO - Using water pdb file at: /shared/sdoerr/Work/htmdacellera/htmd/builder/wat.pdb
2017-02-20 10:41:21,793 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2
Solvating: 100% (8/8) [############################################] eta 00:00 /
2017-02-20 10:41:43,603 - htmd.builder.solvate - INFO - After removing water molecules colliding with other molecules, 10037 water molecules were added to the system.
2017-02-20 10:41:52,504 - htmd.builder.charmm - INFO - Writing out segments.
2017-02-20 10:43:02,342 - htmd.builder.builder - INFO - One disulfide bond was added
2017-02-20 10:43:02,454 - htmd.builder.charmm - INFO - Starting the build.
Bond between A: [serial 1212 resid 140 resname CYS chain  segid P0]
             B: [serial 2448 resid 217 resname CYS chain  segid P0]
2017-02-20 10:43:02,930 - htmd.builder.charmm - INFO - Finished building.
2017-02-20 10:43:05,586 - htmd.builder.ionize - INFO - Adding 14 anions + 0 cations for neutralizing and 70 ions for the given salt concentration.
2017-02-20 10:43:05,939 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A
2017-02-20 10:43:05,940 - htmd.builder.ionize - INFO - Min distance between ions: 5A
2017-02-20 10:43:05,940 - htmd.builder.ionize - INFO - Placing 84 ions.
2017-02-20 10:43:43,330 - htmd.builder.charmm - INFO - Writing out segments.
2017-02-20 10:44:55,163 - htmd.builder.charmm - INFO - Starting the build.
2017-02-20 10:44:55,631 - htmd.builder.charmm - INFO - Finished building.

Equilibrate

from htmd.protocols.equilibration_v2 import Equilibration
md = Equilibration()
md.runtime = 10000000
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')
# Visualize the flat bottom potential box
mol.view('not water')
b = VMDBox([-39, 10, -29, 21, 40, 50])
mdx = AcemdLocal()
mdx.submit(path+'/equil')
mdx.wait()

Production

from htmd.protocols.production_v5 import Production
md = Production()
md.runtime = 50
md.timeunits = 'ns'
md.temperature = 300
md.fb_reference = 'protein and resid 293'
md.fb_selection = 'segname L and noh'
md.fb_k = 5
md.fb_box = [-25, 25, -25, 25, -10, 45]
md.write(path +'/equil','gen/s1')
mol.view('not water')
b = VMDBox([-25, 25, -25, 25, -10, 45])