Extend getDihedral to work on many coordinate sets

[jamesmkrieger]

We can then extract several centers from an ensemble and get many dihedrals from them:

    atoms = ens.getAtoms()

    ul_selstr = 'resnum 1 to 116 244 to 353'
    ll_selstr = 'resnum 117 to 243 354 to 380'
    chainA_selstr = 'chain A and '
    chainC_selstr = 'chain C and '

    ens.setAtoms(atoms.select(chainA_selstr + ll_selstr))
    center_A_LL = calcCenter(ens.getCoordsets())

    ens.setAtoms(atoms.select(chainA_selstr + ul_selstr))
    center_A_UL = calcCenter(ens.getCoordsets())

    ens.setAtoms(atoms.select(chainC_selstr + ul_selstr))
    center_C_UL = calcCenter(ens.getCoordsets())

    ens.setAtoms(atoms.select(chainC_selstr + ll_selstr))
    center_C_LL = calcCenter(ens.getCoordsets())

    ens.setAtoms(atoms)

    disp_torsions = getDihedral(center_A_LL, center_A_UL, 
                                center_C_UL, center_C_LL)

I'll send the ensemble to you, Karolina

[jamesmkrieger]

With the existing code before this PR, we get this error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 20
     16 center_C_LL = calcCenter(ens.getCoordsets())
     18 ens.setAtoms(atoms)
---> 20 disp_torsions = getDihedral(center_A_LL, center_A_UL, 
     21                             center_C_UL, center_C_LL)

File ~/software/scipion3/software/em/prody-github/ProDy/prody/measure/measure.py:210, in getDihedral(coords1, coords2, coords3, coords4, radian)
    207 a3 = coords4 - coords3
    209 v1 = cross(a1, a2)
--> 210 v1 = v1 / (v1 * v1).sum(-1)**0.5
    211 v2 = cross(a2, a3)
    212 v2 = v2 / (v2 * v2).sum(-1)**0.5

ValueError: operands could not be broadcast together with shapes (8,3) (8,) 

[karolamik13]

Everything is working fine when you add 'getDihedral' to all.
Without it, I was not able to run the code.

Python 3.9.18 | packaged by conda-forge | (main, Aug 30 2023, 03:49:32)
Type 'copyright', 'credits' or 'license' for more information
IPython 8.18.1 -- An enhanced Interactive Python. Type '?' for help.

In [1]: from prody import *

In [2]: ens = loadEnsemble('old_plus_apo_and_a2.ens.npz')
...: atoms = ens.getAtoms()
...:
...: ul_selstr = 'resnum 1 to 116 244 to 353'
...: ll_selstr = 'resnum 117 to 243 354 to 380'
...: chainA_selstr = 'chain A and '
...: chainC_selstr = 'chain C and '
...:
...: ens.setAtoms(atoms.select(chainA_selstr + ll_selstr))
...: center_A_LL = calcCenter(ens.getCoordsets())
...:
...: ens.setAtoms(atoms.select(chainA_selstr + ul_selstr))
...: center_A_UL = calcCenter(ens.getCoordsets())
...:
...: ens.setAtoms(atoms.select(chainC_selstr + ul_selstr))
...: center_C_UL = calcCenter(ens.getCoordsets())
...:
...: ens.setAtoms(atoms.select(chainC_selstr + ll_selstr))
...: center_C_LL = calcCenter(ens.getCoordsets())
...:
...: ens.setAtoms(atoms)
...:
...: disp_torsions = getDihedral(center_A_LL, center_A_UL,
...: center_C_UL, center_C_LL)

In [3]: disp_torsions
Out[3]:
array([-39.84047146, -31.05299694, -28.57622459, -20.08363918,
-19.12753077, -15.325695 , -17.3971086 , -3.19208338])