search

dimchris / mdanalysis

Automatically exported from code.google.com/p/mdanalysis
0 stars 0 forks source link

cannot use ligands in HBondAnalysis #75

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
What steps will reproduce the problem?

from MDAnalysis import *
from MDAnalysis.analysis import hbonds
from MDAnalysis.analysis.hbonds import HydrogenBondAnalysis

u = Universe("../in/GM3NKT2.prmtop","CD1dNKT-GM3.mdcrd")

#where lig is the resname of the ligand:
h = HydrogenBondAnalysis(u, 'protein', 'resname lig', distance=3.0, angle=120.0)

What is the expected output? What do you see instead?

Ideally, it should be possible to initialise the HydrogenBondAnalysis module 
with custom atom names to allow for easy analysis of hbonds made by a ligand in 
a protein:ligand system.

In actuality, initialisation of HydrogenBondAnalysis fails because -of course- 
ligand atoms names aren't listen in HydrogenBondAnalysis.acceptors/donors, and 
there is no way of specifying them before initialisation.

What version of the product are you using? On what operating system?
trunk version as of 27-07-2011. OS:Linux RedHat

Please provide any additional information below.

Currently, the following hack appears to fix the problem:

from MDAnalysis import *
from MDAnalysis.analysis import hbonds
from MDAnalysis.analysis.hbonds import HydrogenBondAnalysis

u = Universe("../in/GM3NKT2.prmtop","CD1dNKT-GM3.mdcrd")

#cannot specify non-protein residues, so just put protein twice:
h = HydrogenBondAnalysis(u, 'protein', 'protein', distance=3.0, angle=120.0)

#redefine acceptors and donors to include our custom atoms
h.acceptors = ('CO', 'OH2', 'OW', 'OD1', 'OE1', 'SD', 'OD1', 'OE1', 'OG', 
'OD2', 'OE2', 'OG1', 'SG', 'ND1', 'OH', 'O3C', 'OAA', 'O1H', 'O5', 'O6', 'O3', 
'O4', 'O2', 'O7', 'O5N', 'O8', 'O9', )
h.donors = ('NH', 'OH2', 'OW', 'NE', 'ND2', 'NE2', 'OG', 'OH', 'NH1', 'SG', 
'ND1', 'OG1', 'NH2', 'NE2', 'NZ', 'NE1', 'O3C', 'N2', 'O2', 'O6', 'O3', 'O4', 
'O7', 'O8', 'O9', 'N5', )
#update selection2 to select our ligand
h.selection2 = 'resname LIG'

#now we have the correct atom types listed, and the correct residue is 
#selected, we can finally update the selection
h._update_selection_2()
#and then run
results = h.run()

Original issue reported on code.google.com by Brit.of....@googlemail.com on 28 Jul 2011 at 1:47

GoogleCodeExporter commented 9 years ago 
Thanks for reporting the bug, together with the workaround.

Maybe Dave can have a look (I tentatively set owner to dcaplan)?

If nothing is happening for a week then I'll take a look.

Cheers,
Oliver

Original comment by orbeckst on 28 Jul 2011 at 11:12

GoogleCodeExporter commented 9 years ago 
Btw, you could probably derive a custom HydrogenBondAnalysis class along the 
lines of (untested):

  from MDAnalysis.analysis.hbonds import HydrogenBondAnalysis
  class LigandHydrogenBondAnalysis(HydrogenBondAnalysis):
     donors = HydrogenBondAnalysis.donors + ('O3C', 'N2', 'O2', 'O6', 'O3', 'O4', 'O7', 'O8', 'O9', 'N5')
     acceptors = HydrogenBondAnalysis.acceptors + ('O3C', 'OAA', 'O1H', 'O5', 'O6', 'O3', 'O4', 'O2', 'O7', 'O5N', 'O8', 'O9')

and the use the new class. Still not ideal but a little bit nicer than 
monkey-patching.

Original comment by orbeckst on 28 Jul 2011 at 11:21

GoogleCodeExporter commented 9 years ago 
Here's a patch that adds 2 arguments (donors, acceptors). Maybe Oliver can 
merge it.

Original comment by dcaplan on 29 Jul 2011 at 12:52

Attachments:

GoogleCodeExporter commented 9 years ago 
Added the patch and some smaller cosmetic changes in r826.

See also the updated Online Docs 
http://mdanalysis.googlecode.com/svn/trunk/doc/html/documentation_pages/analysis
/hbonds.html

Please re-open the Issue if things are not working.

Thanks, Oliver

Original comment by orbeckst on 29 Jul 2011 at 4:19

GoogleCodeExporter commented 9 years ago 
Er, the fix was in r829 ...

Original comment by orbeckst on 29 Jul 2011 at 1:22