bowman-lab / enspara

Modeling molecular ensembles with scalable data structures and parallel computing
https://enspara.readthedocs.io
GNU General Public License v3.0
33 stars 16 forks source link

SASA-Exposon calculation problem #221

Closed H-EKE closed 1 year ago

H-EKE commented 1 year ago

Hi,

I am trying to follow the tutorial to calculate exposons. I have tried so far this:

trj = md.load('cent.xtc', top='structure.psf')
sasas = md.shrake_rupley(trj, probe_radius=0.28, mode='residue')

My problem begins here, when I run this 2 commands

sasas = condense_sidechain_sasas(sasas, trj.top)
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-5-26ca3578329b> in <module>
----> 1 sasas = condense_sidechain_sasas(sasas, trj.top)

~/.conda/envs/enspara/lib/python3.6/site-packages/enspara/citation/citation.py in wrapper(*args, **kwargs)
     43         def wrapper(*args, **kwargs):
     44             add_citation(citekey)
---> 45             return func(*args, **kwargs)
     46         return wrapper
     47     return cite_decorator

~/.conda/envs/enspara/lib/python3.6/site-packages/enspara/info_theory/exposons.py in condense_sidechain_sasas(sasas, top)
    144 
    145     assert top.n_residues > 1
--> 146     assert top.n_atoms == sasas.shape[1]
    147 
    148     SELECTION = ('not (name N or name C or name CA or name O or '

AssertionError: 

Did I foget any step? My MD simulation was done with OpenMM and charmm forcefield.

Thanks in advance!

justinrporter commented 1 year ago

Thank you for the very clear report!

The issue is that you've computed residue-level SASAs, which includes the backbone. What you want to do is compute SASAs for each atom, then use condense_sidechain_sasas to combine all the atoms' SASAs into one total SASA for each residue's sidechain.

So instead of

sasas = md.shrake_rupley(trj, probe_radius=0.28, mode='residue')

you want

sasas = md.shrake_rupley(trj, probe_radius=0.28, mode='atom')

Check out the exposons page in the documentation!

Good luck!

justinrporter commented 1 year ago

Reopen this issue if you still have problems!