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
%ls /tmp/testmor/pdb


#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 =, topo=topos, param=params, outdir= path+'prot',ionize=False)
#Add sodium in the receptor
sod = Molecule('pdb/sod.pdb')

#Use a POPC membrane created with vmd and C36
memb = Molecule('pdb/membrane80by80C36.pdb')
mcenter = np.mean(memb.get('coords'),axis=0)
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')
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)

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

mol =, topo=topos, param=params, outdir=os.path.join(path,'build'), saltconc=0.15)


from htmd.protocols.equilibration_v3 import Equilibration
from htmd.mdengine.acemd.acemd import GroupRestraint

# Use a 10A flat bottom potential to prevent the ligand from diffusing from original position during equilibration
width = np.array([10, 10, 10])
flatbot = GroupRestraint('segname L and noh', width, [(5, '0ns')])

md = Equilibration()
md.runtime = 40
md.timeunits = 'ns'
md.temperature = 300
md.restraints = [flatbot] + md.defaultEquilRestraints('20ns')
md.useconstantratio = True
md.write(os.path.join(path,'build'), os.path.join(path,'equil'))
# Visualize the flat bottom potential box
mol.view('not water')
fbcentre = mol.get('coords', sel='segid L').mean(axis=0).squeeze()
b = VMDBox(np.vstack((fbcentre - width/2, fbcentre + width/2)).T.flatten())
mdx = LocalGPUQueue()
mdx.submit(os.path.join(path, 'equil'))


# Read in the last frame of the equilibration
mol = Molecule(os.path.join(path,'equil','structure.psf')),'equil','output.xtc'))
from htmd.protocols.production_v6 import Production
from htmd.mdengine.acemd.acemd import GroupRestraint

# Apply a flat bottom potential to prevent the ligand from entering from periodic image of the protein
width = np.array([70, 70, 60])

# Center the box at residue 218 which is on the upper side of the protein
fbcentre = mol.get('coords', sel='protein and resid 218').mean(axis=0).squeeze()
flatbot = GroupRestraint('segname L and noh', width, [(5, '0ns')], fbcentre=fbcentre)

md = Production()
md.runtime = 50
md.timeunits = 'ns'
md.temperature = 300
md.restraints = flatbot
md.write(os.path.join(path,'equil'), os.path.join(path,'prod'))
mol.view('not water')
b = VMDBox(np.vstack((fbcentre - width/2, fbcentre + width/2)).T.flatten())