From 33ec521b82a04dc973e0bd011738fc53e2c34b0e Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 23 Dec 2021 16:45:10 +0000
Subject: [PATCH 01/73] added parseCFlexModes and writeCFlexModes

---
 prody/dynamics/functions.py | 216 +++++++++++++++++++++++++++++++++++-
 prody/proteins/starfile.py  |  13 ++-
 2 files changed, 222 insertions(+), 7 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 5b5521ba7..33e3c684b 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -1,6 +1,9 @@
 # -*- coding: utf-8 -*-
 """This module defines input and output functions."""
 
+from collections import OrderedDict
+import datetime
+
 import os
 from os.path import abspath, join, isfile, isdir, split, splitext
 
@@ -8,23 +11,26 @@
 
 from prody import LOGGER, SETTINGS, PY3K
 from prody.atomic import Atomic, AtomSubset
-from prody.utilities import openFile, isExecutable, which, PLATFORM, addext
+from prody.utilities import openFile, openSQLite, isExecutable, which, PLATFORM, addext, wrapModes
+from prody.proteins.starfile import parseSTAR, writeSTAR
 
 from .nma import NMA, MaskedNMA
 from .anm import ANM, ANMBase, MaskedANM
+from .analysis import calcCollectivity
 from .gnm import GNM, GNMBase, ZERO, MaskedGNM
 from .exanm import exANM, MaskedExANM
 from .rtb import RTB
 from .pca import PCA, EDA
 from .imanm import imANM
 from .exanm import exANM
-from .mode import Vector, Mode
+from .mode import Vector, Mode, VectorBase
 from .modeset import ModeSet
 from .editing import sliceModel, reduceModel, trimModel
 from .editing import sliceModelByMask, reduceModelByMask, trimModelByMask
 
-__all__ = ['parseArray', 'parseModes', 'parseSparseMatrix',
-           'writeArray', 'writeModes',
+__all__ = ['parseArray', 'parseModes', 'parseCFlexModes',
+           'parseSparseMatrix',
+           'writeArray', 'writeModes', 'writeCFlexModes',
            'saveModel', 'loadModel', 'saveVector', 'loadVector',
            'calcENM']
 
@@ -306,6 +312,208 @@ def parseModes(normalmodes, eigenvalues=None, nm_delimiter=None,
     return nma
 
 
+def parseCFlexModes(run_path, title=None):
+    """Returns :class:`.NMA` containing eigenvectors and eigenvalues 
+    parsed from a ContinuousFlex FlexProtNMA Run directory.
+
+    :arg run_path: path to the Run directory
+    :type run_path: str
+    
+    :arg title: title for :class:`.NMA` object
+    :type title: str
+    """
+    run_name = os.path.split(run_path)[-1]
+
+    star_data = parseSTAR(run_path + '/modes.xmd')
+    star_loop = star_data[0][0]
+    
+    n_modes = star_loop.numRows()
+    
+    row1 = star_loop[0]
+    mode1 = parseArray(row1['_nmaModefile']).reshape(-1)
+    dof = mode1.shape[0]
+
+    vectors = np.zeros((dof, n_modes))
+    vectors[:, 0] = mode1
+
+    for i, row in enumerate(star_loop[1:]):
+        vectors[:, i+1] = parseArray(row['_nmaModefile']).reshape(-1)
+        
+    eigvals = np.zeros(n_modes)
+    
+    log_fname = run_path + '/logs/run.stdout'
+    fi = open(log_fname, 'r')
+    lines = fi.readlines()
+    fi.close()
+    
+    for line in lines:
+        if line.find('Eigenvector number') != -1:
+            j = int(line.strip().split()[-1]) - 1
+        if line.find('Corresponding eigenvalue') != -1:
+            eigvals[j] = float(line.strip().split()[-1])
+        
+    if title is None:
+        title = run_name
+
+    nma = NMA(title)
+    nma.setEigens(vectors, eigvals)
+    return nma
+
+def writeCFlexModes(output_path, modes, write_star=False, scores=None, only_sqlite=False, collectivityThreshold=0.):
+    """Writes *modes* to a set of files that can be recognised by Scipion.
+    A directory called **"modes"** will be created if it doesn't already exist. 
+    Filenames inside will start with **"vec"** and have the mode number as the extension.
+    
+    :arg output_path: path to the directory where the modes directory will be
+    :type output_path: str
+
+    :arg modes: modes to be written to files
+    :type modes: :class:`.Mode`, :class:`.ModeSet`, :class:`.NMA`
+
+    :arg write_star: whether to write modes.xmd STAR file.
+        Default is **False** as qualifyModesStep writes it with scores.
+    :type write_star: bool
+
+    :arg scores: scores from qualifyModesStep for re-writing sqlite
+        Default is **None** and then it writes 0 for each.
+    :type scores: list
+
+    :arg only_sqlite: whether to write only the sqlite file instead of everything.
+        Default is **False** but it can be useful to set it to **True** for updating the sqlite file.
+    :type only_sqlite: bool    
+
+    :arg collectivityThreshold: collectivity threshold below which modes are not enabled
+        Default is 0.
+    :type collectivityThreshold: float
+    """
+    if not isinstance(output_path, str):
+        raise TypeError('output_path should be a string, not {0}'
+                        .format(type(output_path)))
+
+    if not isdir(output_path):
+        raise ValueError('output_path should be a working path')
+
+    if not isinstance(modes, (NMA, ModeSet, VectorBase)):
+        raise TypeError('rows must be NMA, ModeSet, or Mode, not {0}'
+                        .format(type(modes)))
+    if not modes.is3d():
+        raise ValueError('modes must be 3-dimensional')
+
+    if not isinstance(write_star, bool):
+        raise TypeError('write_star should be boolean, not {0}'
+                        .format(type(write_star)))
+
+    if scores is not None:
+        if not isinstance(scores, list):
+            raise TypeError('scores should be a list or None, not {0}'
+                            .format(type(scores)))
+
+    if not isinstance(only_sqlite, bool):
+        raise TypeError('only_sqlite should be boolean, not {0}'
+                        .format(type(only_sqlite)))
+
+    if not isinstance(collectivityThreshold, float):
+        raise TypeError('collectivityThreshold should be float, not {0}'
+                        .format(type(collectivityThreshold)))
+
+    if modes.numModes() == 1:
+        modes = wrapModes(modes)
+
+    modes_dir = output_path + '/modes/'
+    if not isdir(modes_dir):
+        os.mkdir(modes_dir)
+
+    modefiles = []
+    for mode in modes:
+        mode_num = mode.getIndex() + 1
+        modefiles.append(writeArray(modes_dir + 'vec.{0}'.format(mode_num),
+                                     mode.getArrayNx3(), '%12.4e', ''))
+
+    if hasattr(modes, 'getIndices'):
+        order = modes.getIndices()
+        collectivities = list(calcCollectivity(modes))
+        enabled = [1 if eigval > ZERO and collectivities[i] > collectivityThreshold else -1
+                   for i, eigval in enumerate(modes.getEigvals())]
+        if scores is None:
+            scores = [0. for eigval in modes.getEigvals()]
+    else:
+        mode = modes[0]
+        collectivities = [calcCollectivity(mode)]
+        order = [mode.getIndex()]
+        enabled = [1 if mode.getEigval() > ZERO and collectivities[0] > collectivityThreshold else -1]
+        if scores is None:
+            scores = [0.]
+
+    modes_sqlite_fn = output_path + '/modes.sqlite'
+    sql_con = openSQLite(modes_sqlite_fn, 'n')
+    cursor = sql_con.cursor()
+    
+    cursor.execute('''CREATE TABLE Properties(key,value)''')
+    properties = [('self', 'SetOfNormalModes'),
+                  ('_size', str(modes.numModes())),
+                  ('_streamState', '2'),
+                  ('_mapperPath', '{0}, '.format(modes_sqlite_fn))]
+    cursor.executemany('''INSERT INTO Properties VALUES(?,?);''', properties);
+    
+    cursor.execute('''CREATE TABLE Classes(id primary key, label_property, column_name, class_name)''')
+    classes = [(1, 'self', 'c00', 'NormalMode'),
+               (2, '_modeFile', 'c01', 'String'),
+               (3, '_collectivity', 'c02', 'Float'),
+               (4, '_score', 'c03', 'Float')]
+    cursor.executemany('''INSERT INTO Classes VALUES(?,?,?,?);''', classes);
+    
+    cursor.execute('''CREATE TABLE Objects(id primary key, enabled, label, comment, creation, c01, c02, c03)''')
+    
+    star_dict = OrderedDict()
+
+    star_dict['noname'] = OrderedDict() # Data Block with title noname
+    loop_dict = star_dict['noname'][0] = OrderedDict() # Loop 0
+
+    loop_dict['fields'] = OrderedDict()
+    fields = ['_enabled', '_nmaCollectivity', '_nmaModefile', #'_nmaScore',
+              '_order_']
+    for j, field in enumerate(fields):
+        loop_dict['fields'][j] = field
+
+    loop_dict['data'] = OrderedDict()
+
+    now = datetime.datetime.now()    
+    creation = now.strftime("%Y-%m-%d %H:%M:%S")
+    for i, mode in enumerate(modes):
+        loop_dict['data'][i] = OrderedDict()
+        
+        id = mode.getIndex() + 1
+        loop_dict['data'][i]['_order_'] = str(id)
+        
+        enab = enabled[i]
+        loop_dict['data'][i]['_enabled'] = '%2i' % enab
+        if enab != 1:
+            enab = 0
+            
+        label = ''
+        comment = ''
+        
+        c01 = loop_dict['data'][i]['_nmaModefile'] = modefiles[i]
+        
+        collec = collectivities[i]
+        loop_dict['data'][i]['_nmaCollectivity'] = '%8.6f' % collec
+        c02 = float('%6.4f' % collec)
+        
+        c03 = scores[i]
+        loop_dict['data'][i]['_nmaScore'] = '%8.6f' % c03
+        
+        cursor.execute('''INSERT INTO Objects VALUES(?,?,?,?,?,?,?,?)''',
+                       (id, enab, label, comment, creation, c01, c02, c03))
+
+    if write_star:
+        writeSTAR(output_path + '/modes.xmd', star_dict)
+
+    sql_con.commit()
+    sql_con.close()
+
+    return modes_dir
+
+
 def writeArray(filename, array, format='%3.2f', delimiter=' '):
     """Write 1-d or 2-d array data into a delimited text file.
 
diff --git a/prody/proteins/starfile.py b/prody/proteins/starfile.py
index 598e52748..63ca9795b 100644
--- a/prody/proteins/starfile.py
+++ b/prody/proteins/starfile.py
@@ -667,7 +667,7 @@ def parseSTARLines(lines, **kwargs):
     return finalDictionary, prog
 
 
-def writeSTAR(filename, starDict):
+def writeSTAR(filename, starDict, **kwargs):
     """Writes a STAR file from a dictionary containing data
     such as that parsed from a Relion STAR file.
 
@@ -679,7 +679,10 @@ def writeSTAR(filename, starDict):
         This should have nested entries starting with data blocks then loops/tables then
         field names and finally data.
     :type starDict: dict
+
+    kwargs can be given to initiate a :class:`.StarDict` from a dict
     """
+    prog=kwargs.get('prog', 'XMIPP')
 
     star = open(filename, 'w')
 
@@ -688,11 +691,15 @@ def writeSTAR(filename, starDict):
         for loopNumber in starDict[dataBlockKey]:
             star.write('\nloop_\n')
             for fieldNumber in starDict[dataBlockKey][loopNumber]['fields']:
-                star.write('_' + starDict[dataBlockKey][loopNumber]['fields'][fieldNumber] + '\n')
+                if prog == 'XMIPP':
+                    star.write(' ')
+                star.write(starDict[dataBlockKey][loopNumber]['fields'][fieldNumber] + '\n')
             for dataItemNumber in starDict[dataBlockKey][loopNumber]['data']:
+                if prog == 'XMIPP':
+                    star.write('\t')
                 for fieldNumber in starDict[dataBlockKey][loopNumber]['fields']:
                     currentField = starDict[dataBlockKey][loopNumber]['fields'][fieldNumber]
-                    star.write(starDict[dataBlockKey][loopNumber]['data'][dataItemNumber][currentField] + ' ')
+                    star.write(starDict[dataBlockKey][loopNumber]['data'][dataItemNumber][currentField] + '\t')
                 star.write('\n')
 
     star.close()

From 389cb545a5780f679d7b9028b0320c5e6080d24f Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 17 Jan 2022 16:28:18 +0000
Subject: [PATCH 02/73] added outcflex and zeros to apps

---
 prody/apps/prody_apps/nmaoptions.py |  4 ++++
 prody/apps/prody_apps/prody_anm.py  | 11 ++++++++++-
 prody/apps/prody_apps/prody_gnm.py  | 10 +++++++++-
 prody/apps/prody_apps/prody_pca.py  |  3 +++
 prody/dynamics/functions.py         |  9 +++++++--
 5 files changed, 33 insertions(+), 4 deletions(-)

diff --git a/prody/apps/prody_apps/nmaoptions.py b/prody/apps/prody_apps/nmaoptions.py
index 7214f3a53..ce31adf92 100644
--- a/prody/apps/prody_apps/nmaoptions.py
+++ b/prody/apps/prody_apps/nmaoptions.py
@@ -17,6 +17,7 @@
     ('outeig', 'write eigenvalues/vectors', False),
     ('outcov', 'write covariance matrix', False),
     ('outnpz', 'write compressed ProDy data file', False),
+    ('outcflex', 'write continuousflex modes directory and sqlite', False),
     ('outcc', 'write cross-correlations', False),
     ('outhm', 'write cross-correlations heatmap file', False),
     ('outsf', 'write square-fluctuations', False),
@@ -94,6 +95,9 @@ def addNMAOutput(parser):
     parser.add_argument('-z', '--npz', dest='outnpz', action='store_true',
         default=DEFAULTS['outnpz'], help=HELPTEXT['outnpz'])
 
+    parser.add_argument('-y', '--export-continuousflex', dest='outcflex', action='store_true',
+        default=DEFAULTS['outcflex'], help=HELPTEXT['outcflex'])
+
     parser.add_argument('-t', '--extend', dest='extend', type=str,
         metavar='STR', choices=set(['bb', 'all', 'backbone']),
         help=HELPTEXT['extend'])
diff --git a/prody/apps/prody_apps/prody_anm.py b/prody/apps/prody_apps/prody_anm.py
index 5e6d34d06..a5a62e1d8 100644
--- a/prody/apps/prody_apps/prody_anm.py
+++ b/prody/apps/prody_apps/prody_anm.py
@@ -14,6 +14,7 @@
     ('model', 'index of model that will be used in the calculations', 1),
     ('cutoff', 'cutoff distance (A)', 15.),
     ('gamma', 'spring constant', 1.),
+    ('zeros', 'calculate zero modes', False),
 
     ('outbeta', 'write beta-factors calculated from GNM modes', False),
     ('hessian', 'write Hessian matrix', False),
@@ -57,6 +58,7 @@ def prody_anm(pdb, **kwargs):
     nmodes = kwargs.get('nmodes')
     selstr = kwargs.get('select')
     model = kwargs.get('model')
+    zeros = kwargs.get('zeros')
 
     pdb = prody.parsePDB(pdb, model=model)
     if prefix == '_anm':
@@ -72,10 +74,15 @@ def prody_anm(pdb, **kwargs):
 
     anm = prody.ANM(pdb.getTitle())
     anm.buildHessian(select, cutoff, gamma)
-    anm.calcModes(nmodes)
+    anm.calcModes(nmodes, zeros=zeros)
     LOGGER.info('Writing numerical output.')
+
     if kwargs.get('outnpz'):
         prody.saveModel(anm, join(outdir, prefix))
+
+    if kwargs.get('outcflex'):
+        prody.writeCFlexModes(outdir, anm)
+
     prody.writeNMD(join(outdir, prefix + '.nmd'), anm, select)
 
     extend = kwargs.get('extend')
@@ -247,6 +254,8 @@ def addCommand(commands):
     group.add_argument('-m', '--model', dest='model', type=int,
         metavar='INT', default=DEFAULTS['model'], help=HELPTEXT['model'])
 
+    group.add_argument('-w', '--zero-modes', dest='zeros', action='store_true',
+        default=DEFAULTS['zeros'], help=HELPTEXT['zeros'])
 
     group = addNMAOutput(subparser)
 
diff --git a/prody/apps/prody_apps/prody_gnm.py b/prody/apps/prody_apps/prody_gnm.py
index cbc7cbb59..e345db1b0 100644
--- a/prody/apps/prody_apps/prody_gnm.py
+++ b/prody/apps/prody_apps/prody_gnm.py
@@ -13,6 +13,7 @@
     ('model', 'index of model that will be used in the calculations', 1),
     ('cutoff', 'cutoff distance (A)', 10.),
     ('gamma', 'spring constant', 1.),
+    ('zeros', 'calculate zero modes', False),
 
     ('outbeta', 'write beta-factors calculated from GNM modes', False),
     ('kirchhoff', 'write Kirchhoff matrix', False),
@@ -54,6 +55,7 @@ def prody_gnm(pdb, **kwargs):
     nmodes = kwargs.get('nmodes')
     selstr = kwargs.get('select')
     model = kwargs.get('model')
+    zeros = kwargs.get('zeros')
 
     pdb = prody.parsePDB(pdb, model=model)
     if prefix == '_gnm':
@@ -68,13 +70,16 @@ def prody_gnm(pdb, **kwargs):
 
     gnm = prody.GNM(pdb.getTitle())
     gnm.buildKirchhoff(select, cutoff, gamma)
-    gnm.calcModes(nmodes)
+    gnm.calcModes(nmodes, zeros=zeros)
 
     LOGGER.info('Writing numerical output.')
 
     if kwargs.get('outnpz'):
         prody.saveModel(gnm, join(outdir, prefix))
 
+    if kwargs.get('outcflex'):
+        prody.writeCFlexModes(outdir, gnm)
+
     prody.writeNMD(join(outdir, prefix + '.nmd'), gnm, select)
 
     extend = kwargs.get('extend')
@@ -268,6 +273,9 @@ def addCommand(commands):
     group.add_argument('-m', '--model', dest='model', type=int,
         metavar='INT', default=DEFAULTS['model'], help=HELPTEXT['model'])
 
+    group.add_argument('-w', '--zero-modes', dest='zeros', action='store_true',
+        default=DEFAULTS['zeros'], help=HELPTEXT['zeros'])
+
     group = addNMAOutput(subparser)
 
     group.add_argument('-b', '--beta-factors', dest='outbeta',
diff --git a/prody/apps/prody_apps/prody_pca.py b/prody/apps/prody_apps/prody_pca.py
index 7afe478bd..c26f00658 100644
--- a/prody/apps/prody_apps/prody_pca.py
+++ b/prody/apps/prody_apps/prody_pca.py
@@ -122,6 +122,9 @@ def prody_pca(coords, **kwargs):
     if kwargs.get('outnpz'):
         prody.saveModel(pca, join(outdir, prefix))
 
+    if kwargs.get('outcflex'):
+        prody.writeCFlexModes(outdir, pca)
+
     prody.writeNMD(join(outdir, prefix + '.nmd'), pca[:nmodes], select)
 
     extend = kwargs.get('extend')
diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 33e3c684b..6ab637154 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -426,8 +426,13 @@ def writeCFlexModes(output_path, modes, write_star=False, scores=None, only_sqli
     modefiles = []
     for mode in modes:
         mode_num = mode.getIndex() + 1
-        modefiles.append(writeArray(modes_dir + 'vec.{0}'.format(mode_num),
-                                     mode.getArrayNx3(), '%12.4e', ''))
+        if mode.is3d():
+            modefiles.append(writeArray(modes_dir + 'vec.{0}'.format(mode_num),
+                                        mode.getArrayNx3(), '%12.4e', ''))
+        else:
+            modefiles.append(writeArray(modes_dir + 'vec.{0}'.format(mode_num),
+                                        mode.getArray(), '%12.4e', ''))            
+
 
     if hasattr(modes, 'getIndices'):
         order = modes.getIndices()

From aacfea95563301b91fee6c40652b9ac14be8f5b7 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 19 Jan 2022 18:17:44 +0000
Subject: [PATCH 03/73] renaming and apps altloc

---
 prody/apps/prody_apps/nmaoptions.py   |  6 +++---
 prody/apps/prody_apps/prody_anm.py    | 11 ++++++++---
 prody/apps/prody_apps/prody_gnm.py    | 11 ++++++++---
 prody/apps/prody_apps/prody_pca.py    |  4 ++--
 prody/apps/prody_apps/prody_select.py |  6 +++++-
 prody/dynamics/functions.py           |  8 ++++----
 6 files changed, 30 insertions(+), 16 deletions(-)

diff --git a/prody/apps/prody_apps/nmaoptions.py b/prody/apps/prody_apps/nmaoptions.py
index ce31adf92..75c4c7773 100644
--- a/prody/apps/prody_apps/nmaoptions.py
+++ b/prody/apps/prody_apps/nmaoptions.py
@@ -17,7 +17,7 @@
     ('outeig', 'write eigenvalues/vectors', False),
     ('outcov', 'write covariance matrix', False),
     ('outnpz', 'write compressed ProDy data file', False),
-    ('outcflex', 'write continuousflex modes directory and sqlite', False),
+    ('outscipion', 'write continuousflex modes directory and sqlite', False),
     ('outcc', 'write cross-correlations', False),
     ('outhm', 'write cross-correlations heatmap file', False),
     ('outsf', 'write square-fluctuations', False),
@@ -95,8 +95,8 @@ def addNMAOutput(parser):
     parser.add_argument('-z', '--npz', dest='outnpz', action='store_true',
         default=DEFAULTS['outnpz'], help=HELPTEXT['outnpz'])
 
-    parser.add_argument('-y', '--export-continuousflex', dest='outcflex', action='store_true',
-        default=DEFAULTS['outcflex'], help=HELPTEXT['outcflex'])
+    parser.add_argument('-S', '--export-scipion', dest='outscipion', action='store_true',
+        default=DEFAULTS['outscipion'], help=HELPTEXT['outscipion'])
 
     parser.add_argument('-t', '--extend', dest='extend', type=str,
         metavar='STR', choices=set(['bb', 'all', 'backbone']),
diff --git a/prody/apps/prody_apps/prody_anm.py b/prody/apps/prody_apps/prody_anm.py
index a5a62e1d8..8cfec1397 100644
--- a/prody/apps/prody_apps/prody_anm.py
+++ b/prody/apps/prody_apps/prody_anm.py
@@ -12,6 +12,7 @@
 HELPTEXT = {}
 for key, txt, val in [
     ('model', 'index of model that will be used in the calculations', 1),
+    ('altloc', 'alternative location identifiers for residues used in the calculations', "A"),
     ('cutoff', 'cutoff distance (A)', 15.),
     ('gamma', 'spring constant', 1.),
     ('zeros', 'calculate zero modes', False),
@@ -58,9 +59,10 @@ def prody_anm(pdb, **kwargs):
     nmodes = kwargs.get('nmodes')
     selstr = kwargs.get('select')
     model = kwargs.get('model')
+    altloc = kwargs.get('altloc')
     zeros = kwargs.get('zeros')
 
-    pdb = prody.parsePDB(pdb, model=model)
+    pdb = prody.parsePDB(pdb, model=model, altloc=altloc)
     if prefix == '_anm':
         prefix = pdb.getTitle() + '_anm'
 
@@ -80,8 +82,8 @@ def prody_anm(pdb, **kwargs):
     if kwargs.get('outnpz'):
         prody.saveModel(anm, join(outdir, prefix))
 
-    if kwargs.get('outcflex'):
-        prody.writeCFlexModes(outdir, anm)
+    if kwargs.get('outscipion'):
+        prody.writeScipionModes(outdir, anm)
 
     prody.writeNMD(join(outdir, prefix + '.nmd'), anm, select)
 
@@ -254,6 +256,9 @@ def addCommand(commands):
     group.add_argument('-m', '--model', dest='model', type=int,
         metavar='INT', default=DEFAULTS['model'], help=HELPTEXT['model'])
 
+    group.add_argument('-L', '--altloc', dest='altloc', type=str,
+        metavar='STR', default=DEFAULTS['altloc'], help=HELPTEXT['altloc'])
+
     group.add_argument('-w', '--zero-modes', dest='zeros', action='store_true',
         default=DEFAULTS['zeros'], help=HELPTEXT['zeros'])
 
diff --git a/prody/apps/prody_apps/prody_gnm.py b/prody/apps/prody_apps/prody_gnm.py
index e345db1b0..81599d1aa 100644
--- a/prody/apps/prody_apps/prody_gnm.py
+++ b/prody/apps/prody_apps/prody_gnm.py
@@ -11,6 +11,7 @@
 HELPTEXT = {}
 for key, txt, val in [
     ('model', 'index of model that will be used in the calculations', 1),
+    ('altloc', 'alternative location identifiers for residues used in the calculations', "A"),
     ('cutoff', 'cutoff distance (A)', 10.),
     ('gamma', 'spring constant', 1.),
     ('zeros', 'calculate zero modes', False),
@@ -55,9 +56,10 @@ def prody_gnm(pdb, **kwargs):
     nmodes = kwargs.get('nmodes')
     selstr = kwargs.get('select')
     model = kwargs.get('model')
+    altloc = kwargs.get('altloc')
     zeros = kwargs.get('zeros')
 
-    pdb = prody.parsePDB(pdb, model=model)
+    pdb = prody.parsePDB(pdb, model=model, altloc=altloc)
     if prefix == '_gnm':
         prefix = pdb.getTitle() + '_gnm'
 
@@ -77,8 +79,8 @@ def prody_gnm(pdb, **kwargs):
     if kwargs.get('outnpz'):
         prody.saveModel(gnm, join(outdir, prefix))
 
-    if kwargs.get('outcflex'):
-        prody.writeCFlexModes(outdir, gnm)
+    if kwargs.get('outscipion'):
+        prody.writeScipionModes(outdir, gnm)
 
     prody.writeNMD(join(outdir, prefix + '.nmd'), gnm, select)
 
@@ -273,6 +275,9 @@ def addCommand(commands):
     group.add_argument('-m', '--model', dest='model', type=int,
         metavar='INT', default=DEFAULTS['model'], help=HELPTEXT['model'])
 
+    group.add_argument('-L', '--altloc', dest='altloc', type=int,
+        metavar='INT', default=DEFAULTS['altloc'], help=HELPTEXT['altloc'])
+
     group.add_argument('-w', '--zero-modes', dest='zeros', action='store_true',
         default=DEFAULTS['zeros'], help=HELPTEXT['zeros'])
 
diff --git a/prody/apps/prody_apps/prody_pca.py b/prody/apps/prody_apps/prody_pca.py
index 3def0b505..40b11474d 100644
--- a/prody/apps/prody_apps/prody_pca.py
+++ b/prody/apps/prody_apps/prody_pca.py
@@ -123,8 +123,8 @@ def prody_pca(coords, **kwargs):
     if kwargs.get('outnpz'):
         prody.saveModel(pca, join(outdir, prefix))
 
-    if kwargs.get('outcflex'):
-        prody.writeCFlexModes(outdir, pca)
+    if kwargs.get('outscipion'):
+        prody.writeScipionModes(outdir, pca)
 
     prody.writeNMD(join(outdir, prefix + '.nmd'), pca[:nmodes], select)
 
diff --git a/prody/apps/prody_apps/prody_select.py b/prody/apps/prody_apps/prody_select.py
index 1b9aeeed2..757c7194b 100644
--- a/prody/apps/prody_apps/prody_select.py
+++ b/prody/apps/prody_apps/prody_select.py
@@ -35,9 +35,10 @@ def prody_select(selstr, *pdbs, **kwargs):
     prefix = kwargs.get('prefix', None)
     suffix = kwargs.get('suffix', '_selected')
     output = kwargs.get('output', None)
+    altloc = kwargs.get('altloc', None)
 
     for pdb in pdbs:
-        pdb = parsePDB(pdb)
+        pdb = parsePDB(pdb, altloc=altloc)
 
         pdbselect = pdb.select(selstr)
         if pdbselect is None:
@@ -81,6 +82,9 @@ def addCommand(commands):
     group.add_argument('-p', '--prefix', dest='prefix', metavar='STR',
         type=str, help=('output filename prefix (default: PDB filename)'))
 
+    group.add_argument('-L', '--altloc', dest='altloc', metavar='STR',
+        type=str, help=('altloc (default: "A")'))
+
     group.add_argument('-x', '--suffix', dest='suffix', metavar='STR',
         type=str, default='_selected',
         help=('output filename suffix (default: %(default)s)'))
diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index cce5bb395..08821750b 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -28,9 +28,9 @@
 from .editing import sliceModel, reduceModel, trimModel
 from .editing import sliceModelByMask, reduceModelByMask, trimModelByMask
 
-__all__ = ['parseArray', 'parseModes', 'parseCFlexModes',
+__all__ = ['parseArray', 'parseModes', 'parseScipionModes',
            'parseSparseMatrix',
-           'writeArray', 'writeModes', 'writeCFlexModes',
+           'writeArray', 'writeModes', 'writeScipionModes',
            'saveModel', 'loadModel', 'saveVector', 'loadVector',
            'calcENM']
 
@@ -312,7 +312,7 @@ def parseModes(normalmodes, eigenvalues=None, nm_delimiter=None,
     return nma
 
 
-def parseCFlexModes(run_path, title=None):
+def parseScipionModes(run_path, title=None):
     """Returns :class:`.NMA` containing eigenvectors and eigenvalues 
     parsed from a ContinuousFlex FlexProtNMA Run directory.
 
@@ -359,7 +359,7 @@ def parseCFlexModes(run_path, title=None):
     nma.setEigens(vectors, eigvals)
     return nma
 
-def writeCFlexModes(output_path, modes, write_star=False, scores=None, only_sqlite=False, collectivityThreshold=0.):
+def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sqlite=False, collectivityThreshold=0.):
     """Writes *modes* to a set of files that can be recognised by Scipion.
     A directory called **"modes"** will be created if it doesn't already exist. 
     Filenames inside will start with **"vec"** and have the mode number as the extension.

From 283a309c04412b6b186d789b322e00123a3f04dc Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 19 Jan 2022 21:16:43 +0000
Subject: [PATCH 04/73] parseScipionModes allow no eigvals

---
 prody/dynamics/functions.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 08821750b..0867ee07c 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -346,15 +346,22 @@ def parseScipionModes(run_path, title=None):
     lines = fi.readlines()
     fi.close()
     
+    found_eigvals = False
     for line in lines:
         if line.find('Eigenvector number') != -1:
             j = int(line.strip().split()[-1]) - 1
         if line.find('Corresponding eigenvalue') != -1:
             eigvals[j] = float(line.strip().split()[-1])
+            if not found_eigvals:
+                found_eigvals = True
         
     if title is None:
         title = run_name
 
+    if not found_eigvals:
+        LOGGER.warn('No eigenvalues found in {0}/logs/run.stdout'.format())
+        eigvals=None
+
     nma = NMA(title)
     nma.setEigens(vectors, eigvals)
     return nma

From 7d941d8cac5078b10a90e50266c770c1a3778108 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 3 Feb 2022 18:43:43 +0000
Subject: [PATCH 05/73] eigvals in xmd

---
 prody/dynamics/functions.py | 45 ++++++++++++++++++++++++-------------
 1 file changed, 30 insertions(+), 15 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 5210ae258..d12973a3a 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -341,30 +341,38 @@ def parseScipionModes(run_path, title=None):
     vectors = np.zeros((dof, n_modes))
     vectors[:, 0] = mode1
 
+    eigvals = np.zeros(n_modes)
+
+    try:
+        eigvals[0] = float(row1['_nmaEigenval'])
+        found_eigvals = True
+    except:
+        found_eigvals = False
+
     for i, row in enumerate(star_loop[1:]):
         vectors[:, i+1] = parseArray(row['_nmaModefile']).reshape(-1)
-        
-    eigvals = np.zeros(n_modes)
+        if found_eigvals:
+            eigvals[i+1] = float(row['_nmaEigenval'])
     
     log_fname = run_path + '/logs/run.stdout'
     fi = open(log_fname, 'r')
     lines = fi.readlines()
     fi.close()
     
-    found_eigvals = False
-    for line in lines:
-        if line.find('Eigenvector number') != -1:
-            j = int(line.strip().split()[-1]) - 1
-        if line.find('Corresponding eigenvalue') != -1:
-            eigvals[j] = float(line.strip().split()[-1])
-            if not found_eigvals:
-                found_eigvals = True
+    if not found_eigvals:
+        for line in lines:
+            if line.find('Eigenvector number') != -1:
+                j = int(line.strip().split()[-1]) - 1
+            if line.find('Corresponding eigenvalue') != -1:
+                eigvals[j] = float(line.strip().split()[-1])
+                if not found_eigvals:
+                    found_eigvals = True
         
     if title is None:
         title = run_name
 
     if not found_eigvals:
-        LOGGER.warn('No eigenvalues found in {0}/logs/run.stdout'.format())
+        LOGGER.warn('No eigenvalues found')
         eigvals=None
 
     nma = NMA(title)
@@ -445,16 +453,17 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
             modefiles.append(writeArray(modes_dir + 'vec.{0}'.format(mode_num),
                                         mode.getArray(), '%12.4e', ''))            
 
-
     if hasattr(modes, 'getIndices'):
         order = modes.getIndices()
         collectivities = list(calcCollectivity(modes))
+        eigvals = modes.getEigvals()
         enabled = [1 if eigval > ZERO and collectivities[i] > collectivityThreshold else -1
-                   for i, eigval in enumerate(modes.getEigvals())]
+                   for i, eigval in enumerate(eigvals)]
         if scores is None:
             scores = [0. for eigval in modes.getEigvals()]
     else:
         mode = modes[0]
+        eigvals = np.array([mode.getEigval()])
         collectivities = [calcCollectivity(mode)]
         order = [mode.getIndex()]
         enabled = [1 if mode.getEigval() > ZERO and collectivities[0] > collectivityThreshold else -1]
@@ -487,8 +496,8 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
     loop_dict = star_dict['noname'][0] = OrderedDict() # Loop 0
 
     loop_dict['fields'] = OrderedDict()
-    fields = ['_enabled', '_nmaCollectivity', '_nmaModefile', #'_nmaScore',
-              '_order_']
+    fields = ['_enabled', '_nmaCollectivity', '_nmaModefile', '_nmaScore',
+              '_nmaEigenval', '_order_']
     for j, field in enumerate(fields):
         loop_dict['fields'][j] = field
 
@@ -518,6 +527,12 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
         
         c03 = scores[i]
         loop_dict['data'][i]['_nmaScore'] = '%8.6f' % c03
+
+        eigval = eigvals[i]
+        if float('%9.6f' % eigval) > 0:
+            loop_dict['data'][i]['_nmaEigenval'] = '%9.6f' % eigval
+        else:
+            loop_dict['data'][i]['_nmaEigenval'] = '%9.6e' % eigval
         
         cursor.execute('''INSERT INTO Objects VALUES(?,?,?,?,?,?,?,?)''',
                        (id, enab, label, comment, creation, c01, c02, c03))

From 5cca1687421acd862c69568698a9554cb2115784 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 10 Feb 2022 15:06:12 +0000
Subject: [PATCH 06/73] _getEnsembleENMs error typo fix

---
 prody/dynamics/signature.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/dynamics/signature.py b/prody/dynamics/signature.py
index cc71a4808..4dd20ec87 100644
--- a/prody/dynamics/signature.py
+++ b/prody/dynamics/signature.py
@@ -1015,7 +1015,7 @@ def _getEnsembleENMs(ensemble, **kwargs):
             enms = ModeEnsemble()
             enms.addModeSet(ensemble)
         except TypeError:
-            raise TypeError('ensemble must be an Ensemble or a ModeEnsemble instance,'
+            raise TypeError('ensemble must be an Ensemble or a ModeEnsemble instance, '
                             'or a list of NMA, Mode, or ModeSet instances.')
     return enms
 

From 7dc7b689fc7604479f177ad5508931646928c2e3 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 10 Feb 2022 15:07:08 +0000
Subject: [PATCH 07/73] docs fix

---
 INSTALL.rst  | 2 +-
 docs/conf.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/INSTALL.rst b/INSTALL.rst
index 511de3e17..499b5f1f5 100644
--- a/INSTALL.rst
+++ b/INSTALL.rst
@@ -87,7 +87,7 @@ environment variable.
 Recommended Software
 --------------------
 
-* `Scipy`_, when installed, replaces linear algebra module of Numpy.
+* `Scipy`_ , when installed, replaces linear algebra module of Numpy.
   Scipy linear algebra module is more flexible and can be faster.
 * `IPython`_ is a must have for interactive ProDy sessions.
 * `PyReadline`_ for colorful IPython sessions on Windows.
diff --git a/docs/conf.py b/docs/conf.py
index c02bd673d..6ba8006d2 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -170,7 +170,7 @@
 .. _PDB: http://www.pdb.org
 .. _pip: https://pypi.python.org/pypi/pip
 
-.. _MDAnalysis: http://code.google.com/p/mdanalysis
+.. _MDAnalysis: https://www.mdanalysis.org/
 .. _pyparsing: http://pyparsing.wikispaces.com
 .. _Matplotlib: http://matplotlib.org
 .. _Biopython: http://biopython.org

From 344f81a42464d7900e21ba48f9194290a1015c53 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 10 Feb 2022 15:49:40 +0000
Subject: [PATCH 08/73] calcScipionScore

---
 prody/dynamics/analysis.py | 33 +++++++++++++++++++++++++++++++--
 1 file changed, 31 insertions(+), 2 deletions(-)

diff --git a/prody/dynamics/analysis.py b/prody/dynamics/analysis.py
index a1ad2c85c..4d8713f9a 100644
--- a/prody/dynamics/analysis.py
+++ b/prody/dynamics/analysis.py
@@ -23,7 +23,7 @@
            'calcProjection', 'calcCrossProjection',
            'calcSpecDimension', 'calcPairDeformationDist',
            'calcDistFlucts', 'calcHinges', 'calcHitTime', 'calcHitTime',
-           'calcAnisousFromModel']
+           'calcAnisousFromModel', 'calcScipionScore']
            #'calcEntropyTransfer', 'calcOverallNetEntropyTransfer']
 
 def calcCollectivity(mode, masses=None, is3d=None):
@@ -36,7 +36,7 @@ def calcCollectivity(mode, masses=None, is3d=None):
        spin relaxation. *J Chem Phys* **1995** 102:3396-3403.
 
     :arg mode: mode(s) or vector(s)
-    :type mode: :class:`.Mode`, :class:`.Vector`, :class:`.ModeSet`
+    :type mode: :class:`.Mode`, :class:`.Vector`, :class:`.ModeSet`, :class:`.NMA`
 
     :arg masses: atomic masses
     :type masses: :class:`numpy.ndarray`
@@ -707,3 +707,32 @@ def calcAnisousFromModel(model, ):
         anisou[index, 4] = submatrix[0, 2]
         anisou[index, 5] = submatrix[1, 2]
     return anisou
+
+
+def calcScipionScore(modes):
+    """Calculate Scipion continuousflex score, 
+    which is a function of mode number and collectivity order.
+
+    :arg modes: mode(s) or vector(s)
+    :type modes: :class:`.Mode`, :class:`.Vector`, :class:`.ModeSet`, :class:`.NMA`
+    """
+    n_modes = modes.numModes()
+    
+    if n_modes > 1:
+        collectivityList = list(calcCollectivity(modes))
+    else:
+        collectivityList = [calcCollectivity(modes)]
+
+    idxSorted = [i[0] for i in sorted(enumerate(collectivityList),
+                                      key=lambda x: x[1],
+                                      reverse=True)]
+
+    score = np.zeros(n_modes)
+    modeNum = list(range(n_modes))
+
+    for i in range(n_modes):
+        score[idxSorted[i]] = idxSorted[i] + modeNum[i] + 2  
+
+    score = score / (2.0 * n_modes) 
+
+    return score

From e5ae155c8a228b5e76070fd77373aa7b1f329868 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 10 Feb 2022 15:50:31 +0000
Subject: [PATCH 09/73] sliceModel funcs norm option

---
 prody/dynamics/editing.py | 24 ++++++++++++++++++------
 1 file changed, 18 insertions(+), 6 deletions(-)

diff --git a/prody/dynamics/editing.py b/prody/dynamics/editing.py
index 1dfdd801e..9ce1d47fe 100644
--- a/prody/dynamics/editing.py
+++ b/prody/dynamics/editing.py
@@ -266,10 +266,10 @@ def sliceMode(mode, atoms, select):
     return (vec, sel)
 
 
-def sliceModel(model, atoms, select):
+def sliceModel(model, atoms, select, norm=False):
     """Returns a part of the *model* (modes calculated) for *atoms* matching *select*. 
     Note that normal modes are sliced instead the connectivity matrix. Sliced normal 
-    modes (eigenvectors) are not normalized.
+    modes (eigenvectors) are not normalized unless *norm* is **True**.
 
     :arg mode: NMA model instance to be sliced
     :type mode: :class:`.NMA`
@@ -280,6 +280,9 @@ def sliceModel(model, atoms, select):
     :arg select: an atom selection or a selection string
     :type select: :class:`.Selection`, str
 
+    :arg norm: whether to normalize eigenvectors, default **False**
+    :type norm: bool
+
     :returns: (:class:`.NMA`, :class:`.Selection`)"""
 
     if not isinstance(model, NMA):
@@ -292,13 +295,13 @@ def sliceModel(model, atoms, select):
         raise ValueError('number of atoms in model and atoms must be equal')
 
     which, sel = sliceAtoms(atoms, select)
-    nma = sliceModelByMask(model, which)
+    nma = sliceModelByMask(model, which, norm=norm)
 
     return (nma, sel)
 
-def sliceModelByMask(model, mask):
+def sliceModelByMask(model, mask, norm=False):
     """Returns a part of the *model* indicated by *mask*.  Note that
-    normal modes (eigenvectors) are not normalized.
+    normal modes (eigenvectors) are not normalized unless *norm* is **True**.
 
     :arg mode: NMA model instance to be sliced
     :type mode: :class:`.NMA`
@@ -307,6 +310,9 @@ def sliceModelByMask(model, mask):
         the parts being selected 
     :type mask: list, :class:`~numpy.ndarray`
 
+    :arg norm: whether to normalize eigenvectors, default **False**
+    :type norm: bool
+
     :returns: :class:`.NMA`"""
 
     if not isListLike(mask):
@@ -331,7 +337,13 @@ def sliceModelByMask(model, mask):
                     .format(model.getTitle()))
     if model.is3d():
         which = np.repeat(which, 3)
-    nma.setEigens(array[which, :], model.getEigvals())
+
+    evecs = array[which, :]
+    if norm:
+        evecs /= np.array([((evecs[:, i]) ** 2).sum() ** 0.5
+                           for i in range(evecs.shape[1])])
+
+    nma.setEigens(evecs, model.getEigvals())
     return nma
 
 def reduceModel(model, atoms, select):

From 9832d39d56561a97a0909e6bb395921465f11578 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 10 Feb 2022 15:56:55 +0000
Subject: [PATCH 10/73] writeScipionModes use calcScipionScore

---
 prody/dynamics/functions.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index d12973a3a..31be96fd7 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -20,7 +20,7 @@
 
 from .nma import NMA, MaskedNMA
 from .anm import ANM, ANMBase, MaskedANM
-from .analysis import calcCollectivity
+from .analysis import calcCollectivity, calcScipionScore
 from .analysis import calcProjection
 from .gnm import GNM, GNMBase, ZERO, MaskedGNM
 from .exanm import exANM, MaskedExANM
@@ -379,7 +379,9 @@ def parseScipionModes(run_path, title=None):
     nma.setEigens(vectors, eigvals)
     return nma
 
-def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sqlite=False, collectivityThreshold=0.):
+
+def writeScipionModes(output_path, modes, write_star=False, scores=None,
+                      only_sqlite=False, collectivityThreshold=0.):
     """Writes *modes* to a set of files that can be recognised by Scipion.
     A directory called **"modes"** will be created if it doesn't already exist. 
     Filenames inside will start with **"vec"** and have the mode number as the extension.
@@ -460,7 +462,7 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
         enabled = [1 if eigval > ZERO and collectivities[i] > collectivityThreshold else -1
                    for i, eigval in enumerate(eigvals)]
         if scores is None:
-            scores = [0. for eigval in modes.getEigvals()]
+            scores = list(calcScipionScore(modes))
     else:
         mode = modes[0]
         eigvals = np.array([mode.getEigval()])
@@ -468,7 +470,7 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
         order = [mode.getIndex()]
         enabled = [1 if mode.getEigval() > ZERO and collectivities[0] > collectivityThreshold else -1]
         if scores is None:
-            scores = [0.]
+            scores = [calcScipionScore(mode)[0]]
 
     modes_sqlite_fn = output_path + '/modes.sqlite'
     sql_con = openSQLite(modes_sqlite_fn, 'n')

From 394698e1124505a1b73664c25d73e6d8f2b51e82 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 10 Feb 2022 20:28:32 +0000
Subject: [PATCH 11/73] modehunter interpolateModel

---
 prody/dynamics/editing.py | 126 ++++++++++++++++++++++++++++++++++++--
 1 file changed, 122 insertions(+), 4 deletions(-)

diff --git a/prody/dynamics/editing.py b/prody/dynamics/editing.py
index 9ce1d47fe..fc00d7d0a 100644
--- a/prody/dynamics/editing.py
+++ b/prody/dynamics/editing.py
@@ -5,8 +5,9 @@
 
 from prody.atomic import Atomic, AtomGroup, AtomMap, AtomSubset
 from prody.atomic import Selection, SELECT, sliceAtoms, extendAtoms
+from prody.measure import calcDistance
 from prody.utilities import importLA, isListLike
-from prody import _PY3K
+from prody import _PY3K, LOGGER
 
 from .nma import NMA
 from .mode import VectorBase, Mode, Vector
@@ -19,7 +20,8 @@
 
 __all__ = ['extendModel', 'extendMode', 'extendVector',
            'sliceMode', 'sliceModel', 'sliceModelByMask', 'sliceVector',
-           'reduceModel', 'reduceModelByMask', 'trimModel', 'trimModelByMask']
+           'reduceModel', 'reduceModelByMask', 'trimModel', 'trimModelByMask',
+           'interpolateModel']
 
 
 def extendModel(model, nodes, atoms, norm=False):
@@ -37,7 +39,7 @@ def extendModel(model, nodes, atoms, norm=False):
         raise ValueError('model must be an NMA instance')
 
     if model.numAtoms() != nodes.numAtoms():
-        raise ValueError('atom numbers must be the same')
+        raise ValueError('mode and nodes atom numbers must be the same')
 
     indices, atommap = extendAtoms(nodes, atoms, model.is3d())
 
@@ -90,7 +92,7 @@ def extendVector(vector, nodes, atoms):
         raise ValueError('vector must be a Vector instance')
 
     if vector.numAtoms() != nodes.numAtoms():
-        raise ValueError('atom numbers must be the same')
+        raise ValueError('vector and nodes atom numbers must be the same')
 
     indices, atommap = extendAtoms(nodes, atoms, vector.is3d())
     extended = Vector(vec[indices], 'Extended ' + str(vector), vector.is3d())
@@ -485,3 +487,119 @@ def _reduceModel(matrix, system):
         matrix = ss
 
     return matrix
+
+
+def interpolateModel(model, nodes, atoms, norm=False, **kwargs):
+    """Interpolate a coarse grained *model* built for *nodes* to *atoms*.  
+    
+    *model* may be :class:`.ANM`, :class:`.PCA`, or :class:`.NMA`
+    instance. *model* should have all modes calculated.
+    
+    This function will take the part of the normal modes for each node
+    (i.e. Cα atoms) and extend it to nearby atoms. 
+    
+    If *norm* is **True**, extended modes are normalized.
+    
+    Adapted from ModeHunter as described in [JS09]_.
+
+    .. [JS09] Stember JN, Wriggers W. Bend-twist-stretch model for coarse 
+    elastic network simulation of biomolecular motion. *J Chem Phys* **2009** 131:074112.
+    """
+
+    quiet = kwargs.pop('quiet', False)
+
+    try:
+        eigvecs = model._getArray()
+        eigvals = model.getEigvals()
+    except AttributeError:
+        raise TypeError('model must be an NMA instance')
+
+    if model.numAtoms() != nodes.numAtoms():
+        raise ValueError('model and nodes atom numbers must be the same')
+
+    if not model.is3d():
+        raise ValueError('model must be 3D')
+
+    cg_coords = nodes.getCoords()
+    len_cg = nodes.numAtoms()
+
+    fg_coords = atoms.getCoords()
+    len_fg = atoms.numAtoms()
+
+    n_modes = model.numModes()
+
+    eigvecs_fg = np.zeros((3*len_fg, n_modes))
+
+    if not quiet:
+        LOGGER.progress('Interpolating {0} coarse modes to fine modes...'.format(n_modes),
+                        n_modes, '_prody_iterp')
+
+    for k in np.arange(n_modes): # loop through cg eigenvecs
+        coarse_grained_nms = eigvecs[:, k]
+        coarse_grained_nms.shape = (-1,3)
+
+        lmat = np.zeros((len_cg+4,len_cg+4))
+        for i in np.arange(len_cg):
+            for j in np.arange(len_cg):
+                lmat[i,j] = np.linalg.norm(cg_coords[i]-cg_coords[j])
+
+        for i in np.arange(len_cg):
+            lmat[i,len_cg] = 1.0
+            lmat[i,len_cg+1] = cg_coords[i,0]
+            lmat[i,len_cg+2] = cg_coords[i,1]
+            lmat[i,len_cg+3] = cg_coords[i,2]
+            lmat[len_cg,i] = 1.0
+            lmat[len_cg+1,i] = cg_coords[i,0]
+            lmat[len_cg+2,i] = cg_coords[i,1]
+            lmat[len_cg+3,i] = cg_coords[i,2]
+
+        lmat_inv = np.linalg.inv(lmat)
+
+        vxprime = coarse_grained_nms[:,0]
+        for i in np.arange(4):
+            vxprime = np.append(vxprime,0.0)
+
+        vyprime = coarse_grained_nms[:,1]
+        for i in np.arange(4):
+            vyprime = np.append(vyprime,0.0)
+            
+        vzprime = coarse_grained_nms[:,2]
+        for i in np.arange(4):
+            vzprime = np.append(vzprime,0.0)
+
+        vx = np.dot(lmat_inv,vxprime)
+        vy = np.dot(lmat_inv,vyprime)
+        vz = np.dot(lmat_inv,vzprime)
+
+        nmx=np.zeros(len_fg)
+        nmy=np.zeros(len_fg)
+        nmz=np.zeros(len_fg)
+
+        for j in np.arange(len_fg):
+            nmx[j] += vx[len_cg] + fg_coords[j,0]*vx[len_cg+1] + fg_coords[j,1]*vx[len_cg+2] + fg_coords[j,2]*vx[len_cg+3]  
+            nmy[j] += vy[len_cg] + fg_coords[j,0]*vy[len_cg+1] + fg_coords[j,1]*vy[len_cg+2] + fg_coords[j,2]*vy[len_cg+3]  
+            nmz[j] += vz[len_cg] + fg_coords[j,0]*vz[len_cg+1] + fg_coords[j,1]*vz[len_cg+2] + fg_coords[j,2]*vz[len_cg+3]  
+
+            dist_j = calcDistance(cg_coords, fg_coords[j])
+            for i in np.arange(len_cg):
+                nmx[j] += vx[i]*dist_j[i]
+                nmy[j] += vy[i]*dist_j[i]
+                nmz[j] += vz[i]*dist_j[i]
+
+        eigvecs_fg[::3, k] = nmx
+        eigvecs_fg[1::3, k] = nmy
+        eigvecs_fg[2::3, k] = nmz  
+
+        if not quiet:
+            LOGGER.update(k+1, label='_prody_iterp')
+
+    if not quiet:
+        LOGGER.finish()
+
+    if norm:
+        for k in np.arange(np.size(cg_coords)):
+            eigvecs_fg[k] = eigvecs_fg[k]/np.linalg.norm(eigvecs_fg[k])  
+            
+    interpolated = NMA('Interpolated ' + str(model))
+    interpolated.setEigens(eigvecs_fg, eigvals)
+    return interpolated, atoms

From 2b5c8fdfb66538ff1192bd57433989ee69f3b98f Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Fri, 18 Feb 2022 11:29:59 +0000
Subject: [PATCH 12/73] writeScipionModes fix for single modes

---
 prody/dynamics/functions.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 31be96fd7..819675894 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -438,8 +438,10 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
         raise TypeError('collectivityThreshold should be float, not {0}'
                         .format(type(collectivityThreshold)))
 
-    if modes.numModes() == 1:
-        modes = wrapModes(modes)
+    if modes.numModes() == 1 and not isinstance(modes, NMA):
+        old_modes = modes
+        modes = NMA(old_modes)
+        modes.setEigens(old_modes.getArray().reshape(-1, 1))
 
     modes_dir = output_path + '/modes/'
     if not isdir(modes_dir):
@@ -455,7 +457,7 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
             modefiles.append(writeArray(modes_dir + 'vec.{0}'.format(mode_num),
                                         mode.getArray(), '%12.4e', ''))            
 
-    if hasattr(modes, 'getIndices'):
+    if modes.numModes() > 1:
         order = modes.getIndices()
         collectivities = list(calcCollectivity(modes))
         eigvals = modes.getEigvals()

From 13b997a8d7a0e8f5387327a6525f74d838c230ac Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 21 Feb 2022 15:18:25 +0000
Subject: [PATCH 13/73] modehunter license notice

---
 prody/dynamics/editing.py | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/prody/dynamics/editing.py b/prody/dynamics/editing.py
index fc00d7d0a..f033c3fe9 100644
--- a/prody/dynamics/editing.py
+++ b/prody/dynamics/editing.py
@@ -504,6 +504,30 @@ def interpolateModel(model, nodes, atoms, norm=False, **kwargs):
 
     .. [JS09] Stember JN, Wriggers W. Bend-twist-stretch model for coarse 
     elastic network simulation of biomolecular motion. *J Chem Phys* **2009** 131:074112.
+
+    # Legal notice:
+    # 
+    # This software is copyrighted, (c) 2009-10, by Joseph N. Stember and Willy Wriggers
+    # under the following terms:
+    #
+    # The authors hereby grant permission to use, copy, modify, and re-distribute this
+    # software and its documentation for any purpose, provided that existing copyright
+    # notices are retained in all copies and that this notice is included verbatim in
+    # any distributions. No written agreement, license, or royalty fee is required for
+    # any of the authorized uses.
+    #
+    # IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR DIRECT,
+    # INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE
+    # OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF, EVEN IF THE
+    # AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+    #
+    # THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING,
+    # BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
+    # PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE IS PROVIDED ON AN "AS
+    # IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE
+    # MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+    #
+    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     """
 
     quiet = kwargs.pop('quiet', False)

From 4c38bf301b693e1cb6f20b769ebf665387f63959 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 24 Feb 2022 21:37:42 +0000
Subject: [PATCH 14/73] switched parseGromacsModes to mdtraj

---
 prody/dynamics/functions.py | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 819675894..059c785f9 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -791,11 +791,15 @@ def parseGromacsModes(run_path, title="", model='nma', **kwargs):
     :arg eigvec_fname: filename for trr file containing eigenvectors
         Default is ``"eigenvec.trr"`` as this is the default from Gromacs
     :type eigvec_fname: str
+
+    :arg average_pdb: filename for pdb file containing average structure
+        Default is ``"average.pdb"``
+    :type average_pdb: str
     """ 
     try:
-        from MDAnalysis.coordinates import TRR
+        from mdtraj import load_trr
     except ImportError:
-        raise ImportError('Please install MDAnalysis in order to use parseGromacsModes.')
+        raise ImportError('Please install mdtraj in order to use parseGromacsModes.')
 
     if not isinstance(run_path, str):
         raise TypeError('run_path should be a string')
@@ -817,6 +821,10 @@ def parseGromacsModes(run_path, title="", model='nma', **kwargs):
     eigvec_fname = kwargs.get('eigvec_fname', 'eigenvec.trr')
     if not isinstance(eigvec_fname, str):
         raise TypeError('eigvec_fname should be a string')
+
+    average_pdb = kwargs.get('average_pdb', 'average.pdb')
+    if not isinstance(average_pdb, str):
+        raise TypeError('average_pdb should be a string')
     
     vals_fname = run_path + eigval_fname
     fi = open(vals_fname, 'r')
@@ -832,10 +840,10 @@ def parseGromacsModes(run_path, title="", model='nma', **kwargs):
 
     # Parse eigenvectors trr with MDAnalysis, which assumes trajectory and multiplies by 10
     # to get A even though actually they are unit vectors
-    vecs_traj = TRR.TRRReader(run_path + eigvec_fname)
+    vecs_traj = load_trr(run_path+eigvec_fname, top=run_path+average_pdb)
 
     # format vectors appropriately, reversing *10 and skipping initial and average structures
-    vectors = np.array([frame.positions.flatten()/10 for frame in vecs_traj[2:]]).T
+    vectors = np.array([frame.xyz.flatten()/10 for frame in vecs_traj[2:]]).T
 
     result.setEigens(vectors, eigvals)
     return result

From bd2dd80e394d34f53e340add1baea618dc29b772 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Fri, 25 Feb 2022 10:08:48 +0000
Subject: [PATCH 15/73] parseGromacsModes support abs paths

---
 prody/dynamics/functions.py | 57 ++++++++++++++++++++++++++++---------
 1 file changed, 43 insertions(+), 14 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 059c785f9..5d7bbd2b7 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -784,17 +784,17 @@ def parseGromacsModes(run_path, title="", model='nma', **kwargs):
         or ``"pca"``. If it is not changed to ``"pca"`` then ``"nma"`` will be assumed.
     :type model: str
 
-    :arg eigval_fname: filename for xvg file containing eigenvalues
+    :arg eigval_fname: filename or path for xvg file containing eigenvalues
         Default is ``"eigenval.xvg"`` as this is the default from Gromacs
     :type eigval_fname: str
 
-    :arg eigvec_fname: filename for trr file containing eigenvectors
+    :arg eigvec_fname: filename or path for trr file containing eigenvectors
         Default is ``"eigenvec.trr"`` as this is the default from Gromacs
     :type eigvec_fname: str
 
-    :arg average_pdb: filename for pdb file containing average structure
-        Default is ``"average.pdb"``
-    :type average_pdb: str
+    :arg pdb_fname: filename or path for pdb file containing the reference structure
+        Default is ``"average.pdb"`` although this is probably suboptimal
+    :type pdb_fname: str
     """ 
     try:
         from mdtraj import load_trr
@@ -804,6 +804,9 @@ def parseGromacsModes(run_path, title="", model='nma', **kwargs):
     if not isinstance(run_path, str):
         raise TypeError('run_path should be a string')
 
+    if not run_path.endswith('/'):
+        run_path += '/'
+
     if not isinstance(title, str):
         raise TypeError('title should be a string')
 
@@ -814,19 +817,46 @@ def parseGromacsModes(run_path, title="", model='nma', **kwargs):
             LOGGER.warn('model not recognised so using NMA')
         result = NMA(title)
 
+
     eigval_fname = kwargs.get('eigval_fname', 'eigenval.xvg')
     if not isinstance(eigval_fname, str):
         raise TypeError('eigval_fname should be a string')
 
+    if isfile(eigval_fname):
+        vals_fname = eigval_fname
+    elif isfile(run_path + eigval_fname):
+        vals_fname = run_path + eigval_fname
+    else:
+        raise ValueError('eigval_fname should point be a path to a file '
+                         'either relative to run_path or an absolute one')
+
+
     eigvec_fname = kwargs.get('eigvec_fname', 'eigenvec.trr')
     if not isinstance(eigvec_fname, str):
         raise TypeError('eigvec_fname should be a string')
 
-    average_pdb = kwargs.get('average_pdb', 'average.pdb')
-    if not isinstance(average_pdb, str):
-        raise TypeError('average_pdb should be a string')
+    if isfile(eigvec_fname):
+        vecs_fname = eigval_fname
+    elif isfile(run_path + eigvec_fname):
+        vecs_fname = run_path + eigvec_fname
+    else:
+        raise ValueError('eigvec_fname should point be a path to a file '
+                         'either relative to run_path or an absolute one')
+
+
+    pdb_fname = kwargs.get('pdb_fname', 'average.pdb')
+    if not isinstance(pdb_fname, str):
+        raise TypeError('pdb_fname should be a string')
+
+    if isfile(pdb_fname):
+        pdb = eigval_fname
+    elif isfile(run_path + pdb_fname):
+        pdb = run_path + pdb_fname
+    else:
+        raise ValueError('pdb_fname should point be a path to a file '
+                         'either relative to run_path or an absolute one')
+    
     
-    vals_fname = run_path + eigval_fname
     fi = open(vals_fname, 'r')
     lines = fi.readlines()
     fi.close()
@@ -838,12 +868,11 @@ def parseGromacsModes(run_path, title="", model='nma', **kwargs):
 
     eigvals = np.array(eigvals)
 
-    # Parse eigenvectors trr with MDAnalysis, which assumes trajectory and multiplies by 10
-    # to get A even though actually they are unit vectors
-    vecs_traj = load_trr(run_path+eigvec_fname, top=run_path+average_pdb)
+    # Parse eigenvectors trr with mdtraj, which uses nm so doesn't rescale
+    vecs_traj = load_trr(vecs_fname, top=pdb)
 
-    # format vectors appropriately, reversing *10 and skipping initial and average structures
-    vectors = np.array([frame.xyz.flatten()/10 for frame in vecs_traj[2:]]).T
+    # format vectors appropriately, skipping initial and average structures
+    vectors = np.array([frame.xyz.flatten() for frame in vecs_traj[2:]]).T
 
     result.setEigens(vectors, eigvals)
     return result

From ee864f3692036739adf59f73375ba152e67604af Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 1 Mar 2022 15:46:28 +0000
Subject: [PATCH 16/73] cealign non-ca fix

---
 prody/proteins/compare.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/prody/proteins/compare.py b/prody/proteins/compare.py
index 89a096511..79845f5c1 100644
--- a/prody/proteins/compare.py
+++ b/prody/proteins/compare.py
@@ -1394,8 +1394,8 @@ def getCEAlignMapping(target, chain):
                     'It may not be installed properly.')
         return None
 
-    tar_coords = target.getCoords().tolist()
-    mob_coords = chain.getCoords().tolist()
+    tar_coords = target.getCoords(calpha=True).tolist()
+    mob_coords = chain.getCoords(calpha=True).tolist()
 
     if len(tar_coords) < 8:
         LOGGER.warn('target ({1}) is too small to be aligned '

From 4cb7c603c27209348a9e8f275ae3c3ad14a7ea43 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 2 Mar 2022 17:30:05 +0000
Subject: [PATCH 17/73] docs fix

---
 prody/atomic/functions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/atomic/functions.py b/prody/atomic/functions.py
index 13a55cd9c..1f15f8ff4 100644
--- a/prody/atomic/functions.py
+++ b/prody/atomic/functions.py
@@ -614,7 +614,7 @@ def assignBlocks(atoms, res_per_block=None, secstr=False, shortest_block=2):
         secstrs = sel_ca.getSecstrs()
         if secstrs is None:
             raise OSError("Please parse secstr information "
-                          "from PDB header or make sure DSSP or STRIDE "
+                          "from PDB header or use DSSP or STRIDE "
                           "to use secstr for assigning blocks")
         
         blocks.append(0)

From c0308c805f935e80333bcfe49227cb28dd42da34 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 6 Mar 2022 10:34:30 +0000
Subject: [PATCH 18/73] starfile error fix

---
 prody/proteins/cifheader.py | 0
 prody/proteins/starfile.py  | 2 +-
 2 files changed, 1 insertion(+), 1 deletion(-)
 create mode 100644 prody/proteins/cifheader.py

diff --git a/prody/proteins/cifheader.py b/prody/proteins/cifheader.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/prody/proteins/starfile.py b/prody/proteins/starfile.py
index d392b49c9..96fabb182 100644
--- a/prody/proteins/starfile.py
+++ b/prody/proteins/starfile.py
@@ -674,7 +674,7 @@ def parseSTARLines(lines, **kwargs):
                 prog = 'XMIPP'
 
         else:
-            raise TypeError('This file does not conform to the STAR file format.'
+            raise TypeError('This file does not conform to the STAR file format. '
                             'There is a problem with line {0}:\n {1}'.format(lineNumber, line))
 
         lineNumber += 1

From 5dce94811b998b6a4dfa6f362ae79ac98d1f1182 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 7 Mar 2022 16:35:03 +0000
Subject: [PATCH 19/73] writePDB secstr

---
 prody/proteins/header.py  |  4 +--
 prody/proteins/pdbfile.py | 71 ++++++++++++++++++++++++++++++++++++++-
 2 files changed, 72 insertions(+), 3 deletions(-)

diff --git a/prody/proteins/header.py b/prody/proteins/header.py
index c55a7aafd..086725b77 100644
--- a/prody/proteins/header.py
+++ b/prody/proteins/header.py
@@ -470,8 +470,8 @@ def _getSheet(lines):
     for i, line in lines['SHEET ']:
         try:
             chid = line[21]
-            value = (int(line[38:40]), int(line[7:10]),
-                     line[11:14].strip())
+                     # sense           # strand num     # sheet id
+            value = (int(line[38:40]), int(line[7:10]), line[11:14].strip())
         except:
             continue
 
diff --git a/prody/proteins/pdbfile.py b/prody/proteins/pdbfile.py
index a73041f52..5df1f87b2 100644
--- a/prody/proteins/pdbfile.py
+++ b/prody/proteins/pdbfile.py
@@ -1542,6 +1542,75 @@ def writePQRStream(stream, atoms, **kwargs):
 
     s_or_u = np.array(['a']).dtype.char
 
+
+    write = stream.write
+
+    calphas = atoms.ca
+    ssa = calphas.getSecstrs()
+    helix = []
+    sheet = []
+    if ssa is not None:
+        ss_prev = ssa[0]
+        ss_start = 0
+        ss_end = 1
+        for i, ss in enumerate(ssa):
+            if ss != ss_prev:
+                # store prev secstr and prepare for next
+                ss_end = i-1
+                init = calphas[ss_start]
+                end = calphas[ss_end]
+                length = ss_end - ss_start + 1
+
+                entries = [init.getSecindex(), init.getSecid(),
+                           init.getResname(), init.getChid(), 
+                           init.getResnum(), init.getIcode(),
+                           end.getResname(), end.getChid(),
+                           end.getResnum(), end.getIcode(),
+                           init.getSecclass()]
+                
+                if ssa[ss_end] == 'H':
+                    helix.append(["HELIX "] + entries +
+                                 ['', length])
+
+                elif ssa[ss_end] == 'E':
+                    sheet.append(["SHEET "] + entries)
+
+                ss_start = i
+                ss_prev = ss
+
+    format_helix = ('{0:6s} {1:3d} {2:3s} ' +
+                    '{3:3s} {4:1s} {5:4d}{6:1s} ' +
+                    '{7:3s} {8:1s} {9:4d}{10:1s} ' +
+                    '{11:2d} {12:30s} {13:5d}\n').format
+    for line in helix:
+        write(format_helix(*line))
+
+
+    sorted_sheet = sorted(sheet, key=lambda item: (item[2], item[1]))
+    sheet_prev = 'A'
+    num_strands_list = []
+    for i, item1 in enumerate(sorted_sheet):
+        if item1[2] != sheet_prev:
+            num_strands = sorted_sheet[i-1][1]
+            num_strands_list.append(num_strands)
+
+            sheet_prev = item1[2]
+
+            for item2 in sorted_sheet[i-num_strands:i]:
+                item2.append(num_strands)
+
+    num_strands = item1[1]
+    for item2 in sorted_sheet[i-num_strands+1:]:
+        item2.append(num_strands)    
+
+    format_sheet = ('{0:6s} {1:3d} {2:3s}{12:2d} ' +
+                    '{3:3s} {4:1s}{5:4d}{6:1s}' +
+                    '{7:3s} {8:1s}{9:4d}{10:1s}' +
+                    '{11:2d}\n').format
+
+    for i, line in enumerate(sorted_sheet):
+        write(format_sheet(*line))
+
     resnames = atoms._getResnames()
     if resnames is None:
         resnames = ['UNK'] * n_atoms
@@ -1576,7 +1645,7 @@ def writePQRStream(stream, atoms, **kwargs):
               '{8:8.3f} {9:8.3f} {10:8.3f}' +
               '{11:8.4f} {12:7.4f}\n').format
     coords = atoms._getCoords()
-    write = stream.write
+    
     for i, xyz in enumerate(coords):
         write(format(hetero[i], i+1, atomnames[i], altlocs[i],
                      resnames[i], chainids[i], int(resnums[i]),

From 4835ba9dd88556ee219bf15b7cef6252756d7cf6 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 7 Mar 2022 17:53:15 +0000
Subject: [PATCH 20/73] added max block size

---
 prody/atomic/functions.py | 32 +++++++++++++++++++++++++++++++-
 1 file changed, 31 insertions(+), 1 deletion(-)

diff --git a/prody/atomic/functions.py b/prody/atomic/functions.py
index 1f15f8ff4..39c252dd4 100644
--- a/prody/atomic/functions.py
+++ b/prody/atomic/functions.py
@@ -6,6 +6,7 @@
 
 from numpy import load, savez, ones, zeros, array, argmin, where
 from numpy import ndarray, asarray, isscalar, concatenate, arange, ix_
+from numpy import append, unique
 
 from prody.utilities import openFile, rangeString, getDistance, fastin
 from prody import LOGGER
@@ -547,7 +548,7 @@ def extendAtomicData(data, nodes, atoms, axis=None):
 extendData = extendAtomicData
 
 
-def assignBlocks(atoms, res_per_block=None, secstr=False, shortest_block=2):
+def assignBlocks(atoms, res_per_block=None, secstr=False, **kwargs):
     """Assigns blocks to protein from *atoms*
     using a block size of *res_per_block* or 
     secondary structure information if *secstr* is **True**.
@@ -573,6 +574,11 @@ def assignBlocks(atoms, res_per_block=None, secstr=False, shortest_block=2):
         in a block before merging with the previous block
         Default is **2**
     :type shortest_block: int
+
+    :arg longest_block: largest number of residues to be included 
+        in a block before splitting it in half.
+        Default is the length of the protein so it isn't triggered.
+    :type longest_block: int
     """
 
     if not isinstance(atoms, Atomic):
@@ -592,12 +598,19 @@ def assignBlocks(atoms, res_per_block=None, secstr=False, shortest_block=2):
     if not isinstance(secstr, bool):
         raise TypeError('secstr should be a Boolean')
 
+    shortest_block = kwargs.get('shortest_block', 2)
+    shortest_block = kwargs.get('min_size', shortest_block)
     if not isinstance(shortest_block, Integral):
         raise TypeError("shortest_block should be an integer")
 
     sel_ca = atoms.ca
     n_res = sel_ca.numAtoms()
 
+    longest_block = kwargs.get('max_size', n_res)
+    longest_block = kwargs.get('longest_block', longest_block)
+    if not isinstance(longest_block, Integral):
+        raise TypeError("longest_block should be an integer")
+
     blocks = []
 
     if res_per_block:
@@ -624,12 +637,29 @@ def assignBlocks(atoms, res_per_block=None, secstr=False, shortest_block=2):
             if secstr != secstr_prev:
                 secstr_prev = secstr
                 if len(where(array(blocks) == i)[0]) >= shortest_block:
+                    # long enough so next sec str is a new block
                     i += 1
             blocks.append(i)
         
         # include last residue in previous block
         blocks.append(i)
 
+    blocks = array(blocks)
+
+    unique_blocks, lengths = unique(blocks, return_counts=True)
+    max_length = max(lengths)
+
+    while max_length > longest_block:
+        for i in unique_blocks:
+            len_block = len(where(blocks == i)[0])
+
+            if len_block > longest_block:
+                new_block = max(unique_blocks)+1
+                blocks[where(blocks == i)[0][len_block // 2 :]] = new_block
+
+                unique_blocks, lengths = unique(blocks, return_counts=True)
+                max_length = max(lengths)
+
     blocks, amap = extendAtomicData(blocks, sel_ca, atoms)
 
     if amap.getHierView().numResidues() < atoms.getHierView().numResidues():

From 62c8e6dd5af58ab4a9996a0cfd268ddb86fc26d1 Mon Sep 17 00:00:00 2001
From: Ricardo Serrano <rserranogut@hotmail.com>
Date: Wed, 9 Mar 2022 14:02:24 +0100
Subject: [PATCH 21/73] gnm app altloc fix

---
 prody/apps/prody_apps/prody_gnm.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/apps/prody_apps/prody_gnm.py b/prody/apps/prody_apps/prody_gnm.py
index 81599d1aa..e779ad995 100644
--- a/prody/apps/prody_apps/prody_gnm.py
+++ b/prody/apps/prody_apps/prody_gnm.py
@@ -275,7 +275,7 @@ def addCommand(commands):
     group.add_argument('-m', '--model', dest='model', type=int,
         metavar='INT', default=DEFAULTS['model'], help=HELPTEXT['model'])
 
-    group.add_argument('-L', '--altloc', dest='altloc', type=int,
+    group.add_argument('-L', '--altloc', dest='altloc', type=str,
         metavar='INT', default=DEFAULTS['altloc'], help=HELPTEXT['altloc'])
 
     group.add_argument('-w', '--zero-modes', dest='zeros', action='store_true',

From 06e962ae4a0e5ace7e207f29969f7df36622a90b Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 9 Mar 2022 15:11:14 +0000
Subject: [PATCH 22/73] writeScipionModes no is3d error

---
 prody/dynamics/functions.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 5d7bbd2b7..c346489fe 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -418,8 +418,6 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
     if not isinstance(modes, (NMA, ModeSet, VectorBase)):
         raise TypeError('rows must be NMA, ModeSet, or Mode, not {0}'
                         .format(type(modes)))
-    if not modes.is3d():
-        raise ValueError('modes must be 3-dimensional')
 
     if not isinstance(write_star, bool):
         raise TypeError('write_star should be boolean, not {0}'

From 74dd0445e57f45c3477d08f105920456e84c3ba1 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 9 Mar 2022 15:15:47 +0000
Subject: [PATCH 23/73] writeScipionModes eigvals fix

---
 prody/dynamics/functions.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index c346489fe..99ede780f 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -354,12 +354,12 @@ def parseScipionModes(run_path, title=None):
         if found_eigvals:
             eigvals[i+1] = float(row['_nmaEigenval'])
     
-    log_fname = run_path + '/logs/run.stdout'
-    fi = open(log_fname, 'r')
-    lines = fi.readlines()
-    fi.close()
-    
     if not found_eigvals:
+        log_fname = run_path + '/logs/run.stdout'
+        fi = open(log_fname, 'r')
+        lines = fi.readlines()
+        fi.close()
+
         for line in lines:
             if line.find('Eigenvector number') != -1:
                 j = int(line.strip().split()[-1]) - 1

From f7b24caa7c67dd63a838a1dc1799f9422342e227 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 10 Feb 2022 20:28:32 +0000
Subject: [PATCH 24/73] modehunter interpolateModel

---
 prody/dynamics/editing.py | 150 +++++++++++++++++++++++++++++++++++++-
 1 file changed, 146 insertions(+), 4 deletions(-)

diff --git a/prody/dynamics/editing.py b/prody/dynamics/editing.py
index 1dfdd801e..230f083bd 100644
--- a/prody/dynamics/editing.py
+++ b/prody/dynamics/editing.py
@@ -5,8 +5,9 @@
 
 from prody.atomic import Atomic, AtomGroup, AtomMap, AtomSubset
 from prody.atomic import Selection, SELECT, sliceAtoms, extendAtoms
+from prody.measure import calcDistance
 from prody.utilities import importLA, isListLike
-from prody import _PY3K
+from prody import _PY3K, LOGGER
 
 from .nma import NMA
 from .mode import VectorBase, Mode, Vector
@@ -19,7 +20,8 @@
 
 __all__ = ['extendModel', 'extendMode', 'extendVector',
            'sliceMode', 'sliceModel', 'sliceModelByMask', 'sliceVector',
-           'reduceModel', 'reduceModelByMask', 'trimModel', 'trimModelByMask']
+           'reduceModel', 'reduceModelByMask', 'trimModel', 'trimModelByMask',
+           'interpolateModel']
 
 
 def extendModel(model, nodes, atoms, norm=False):
@@ -37,7 +39,7 @@ def extendModel(model, nodes, atoms, norm=False):
         raise ValueError('model must be an NMA instance')
 
     if model.numAtoms() != nodes.numAtoms():
-        raise ValueError('atom numbers must be the same')
+        raise ValueError('mode and nodes atom numbers must be the same')
 
     indices, atommap = extendAtoms(nodes, atoms, model.is3d())
 
@@ -90,7 +92,7 @@ def extendVector(vector, nodes, atoms):
         raise ValueError('vector must be a Vector instance')
 
     if vector.numAtoms() != nodes.numAtoms():
-        raise ValueError('atom numbers must be the same')
+        raise ValueError('vector and nodes atom numbers must be the same')
 
     indices, atommap = extendAtoms(nodes, atoms, vector.is3d())
     extended = Vector(vec[indices], 'Extended ' + str(vector), vector.is3d())
@@ -473,3 +475,143 @@ def _reduceModel(matrix, system):
         matrix = ss
 
     return matrix
+
+
+def interpolateModel(model, nodes, atoms, norm=False, **kwargs):
+    """Interpolate a coarse grained *model* built for *nodes* to *atoms*.  
+    
+    *model* may be :class:`.ANM`, :class:`.PCA`, or :class:`.NMA`
+    instance. *model* should have all modes calculated.
+    
+    This function will take the part of the normal modes for each node
+    (i.e. Cα atoms) and extend it to nearby atoms. 
+    
+    If *norm* is **True**, extended modes are normalized.
+    
+    Adapted from ModeHunter as described in [JS09]_.
+
+    .. [JS09] Stember JN, Wriggers W. Bend-twist-stretch model for coarse 
+    elastic network simulation of biomolecular motion. *J Chem Phys* **2009** 131:074112.
+
+    # Legal notice:
+    # 
+    # This software is copyrighted, (c) 2009-10, by Joseph N. Stember and Willy Wriggers
+    # under the following terms:
+    #
+    # The authors hereby grant permission to use, copy, modify, and re-distribute this
+    # software and its documentation for any purpose, provided that existing copyright
+    # notices are retained in all copies and that this notice is included verbatim in
+    # any distributions. No written agreement, license, or royalty fee is required for
+    # any of the authorized uses.
+    #
+    # IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR DIRECT,
+    # INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE
+    # OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF, EVEN IF THE
+    # AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+    #
+    # THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING,
+    # BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
+    # PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE IS PROVIDED ON AN "AS
+    # IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE
+    # MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+    #
+    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    """
+
+    quiet = kwargs.pop('quiet', False)
+
+    try:
+        eigvecs = model._getArray()
+        eigvals = model.getEigvals()
+    except AttributeError:
+        raise TypeError('model must be an NMA instance')
+
+    if model.numAtoms() != nodes.numAtoms():
+        raise ValueError('model and nodes atom numbers must be the same')
+
+    if not model.is3d():
+        raise ValueError('model must be 3D')
+
+    cg_coords = nodes.getCoords()
+    len_cg = nodes.numAtoms()
+
+    fg_coords = atoms.getCoords()
+    len_fg = atoms.numAtoms()
+
+    n_modes = model.numModes()
+
+    eigvecs_fg = np.zeros((3*len_fg, n_modes))
+
+    if not quiet:
+        LOGGER.progress('Interpolating {0} coarse modes to fine modes...'.format(n_modes),
+                        n_modes, '_prody_iterp')
+
+    for k in np.arange(n_modes): # loop through cg eigenvecs
+        coarse_grained_nms = eigvecs[:, k]
+        coarse_grained_nms.shape = (-1,3)
+
+        lmat = np.zeros((len_cg+4,len_cg+4))
+        for i in np.arange(len_cg):
+            for j in np.arange(len_cg):
+                lmat[i,j] = np.linalg.norm(cg_coords[i]-cg_coords[j])
+
+        for i in np.arange(len_cg):
+            lmat[i,len_cg] = 1.0
+            lmat[i,len_cg+1] = cg_coords[i,0]
+            lmat[i,len_cg+2] = cg_coords[i,1]
+            lmat[i,len_cg+3] = cg_coords[i,2]
+            lmat[len_cg,i] = 1.0
+            lmat[len_cg+1,i] = cg_coords[i,0]
+            lmat[len_cg+2,i] = cg_coords[i,1]
+            lmat[len_cg+3,i] = cg_coords[i,2]
+
+        lmat_inv = np.linalg.inv(lmat)
+
+        vxprime = coarse_grained_nms[:,0]
+        for i in np.arange(4):
+            vxprime = np.append(vxprime,0.0)
+
+        vyprime = coarse_grained_nms[:,1]
+        for i in np.arange(4):
+            vyprime = np.append(vyprime,0.0)
+            
+        vzprime = coarse_grained_nms[:,2]
+        for i in np.arange(4):
+            vzprime = np.append(vzprime,0.0)
+
+        vx = np.dot(lmat_inv,vxprime)
+        vy = np.dot(lmat_inv,vyprime)
+        vz = np.dot(lmat_inv,vzprime)
+
+        nmx=np.zeros(len_fg)
+        nmy=np.zeros(len_fg)
+        nmz=np.zeros(len_fg)
+
+        for j in np.arange(len_fg):
+            nmx[j] += vx[len_cg] + fg_coords[j,0]*vx[len_cg+1] + fg_coords[j,1]*vx[len_cg+2] + fg_coords[j,2]*vx[len_cg+3]  
+            nmy[j] += vy[len_cg] + fg_coords[j,0]*vy[len_cg+1] + fg_coords[j,1]*vy[len_cg+2] + fg_coords[j,2]*vy[len_cg+3]  
+            nmz[j] += vz[len_cg] + fg_coords[j,0]*vz[len_cg+1] + fg_coords[j,1]*vz[len_cg+2] + fg_coords[j,2]*vz[len_cg+3]  
+
+            dist_j = calcDistance(cg_coords, fg_coords[j])
+            for i in np.arange(len_cg):
+                nmx[j] += vx[i]*dist_j[i]
+                nmy[j] += vy[i]*dist_j[i]
+                nmz[j] += vz[i]*dist_j[i]
+
+        eigvecs_fg[::3, k] = nmx
+        eigvecs_fg[1::3, k] = nmy
+        eigvecs_fg[2::3, k] = nmz  
+
+        if not quiet:
+            LOGGER.update(k+1, label='_prody_iterp')
+
+    if not quiet:
+        LOGGER.finish()
+
+    if norm:
+        for k in np.arange(np.size(cg_coords)):
+            eigvecs_fg[k] = eigvecs_fg[k]/np.linalg.norm(eigvecs_fg[k])  
+            
+    interpolated = NMA('Interpolated ' + str(model))
+    interpolated.setEigens(eigvecs_fg, eigvals)
+    return interpolated, atoms

From 9af75888dfae01dfe1a05488bc74508ba2777912 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 13 Mar 2022 11:39:24 +0000
Subject: [PATCH 25/73] revert INSTALL scipy space

---
 INSTALL.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/INSTALL.rst b/INSTALL.rst
index d92f1c330..dcf4f1185 100644
--- a/INSTALL.rst
+++ b/INSTALL.rst
@@ -87,7 +87,7 @@ environment variable.
 Recommended Software
 --------------------
 
-* `Scipy`_ , when installed, replaces linear algebra module of Numpy.
+* `Scipy`_, when installed, replaces linear algebra module of Numpy.
   Scipy linear algebra module is more flexible and can be faster.
 * `IPython`_ is a must have for interactive ProDy sessions.
 * `PyReadline`_ for colorful IPython sessions on Windows.

From 9912aff1d61e3be9b2099ba5e894e3db0c7db302 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 13 Mar 2022 12:07:55 +0000
Subject: [PATCH 26/73] -

---
 prody/dynamics/functions.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 99ede780f..bee1517af 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -33,8 +33,8 @@
 from .editing import sliceModel, reduceModel, trimModel
 from .editing import sliceModelByMask, reduceModelByMask, trimModelByMask
 
-__all__ = ['parseArray', 'parseModes', 'parseScipionModes',
-           'parseSparseMatrix', 'parseGromacsModes',
+__all__ = ['parseArray', 'parseModes', 'parseSparseMatrix', 
+           'parseGromacsModes', 'parseScipionModes',
            'writeArray', 'writeModes', 'writeScipionModes',
            'saveModel', 'loadModel', 'saveVector', 'loadVector',
            'calcENM', 'realignModes']

From b88732cc346af5c0d09f040af023e703703ecb57 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 13 Mar 2022 12:08:37 +0000
Subject: [PATCH 27/73] -

---
 prody/dynamics/functions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index bee1517af..f6bbcae26 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -33,7 +33,7 @@
 from .editing import sliceModel, reduceModel, trimModel
 from .editing import sliceModelByMask, reduceModelByMask, trimModelByMask
 
-__all__ = ['parseArray', 'parseModes', 'parseSparseMatrix', 
+__all__ = ['parseArray', 'parseModes', 'parseSparseMatrix',
            'parseGromacsModes', 'parseScipionModes',
            'writeArray', 'writeModes', 'writeScipionModes',
            'saveModel', 'loadModel', 'saveVector', 'loadVector',

From 2621576afe517f7943eabbbf627e404bd37fec68 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 13 Mar 2022 13:02:59 +0000
Subject: [PATCH 28/73] doc using calcScipionScore

---
 prody/dynamics/functions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index f6bbcae26..fa458adfd 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -397,7 +397,7 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
     :type write_star: bool
 
     :arg scores: scores from qualifyModesStep for re-writing sqlite
-        Default is **None** and then it writes 0 for each.
+        Default is **None** and then it uses :func:`.calcScipionScore`
     :type scores: list
 
     :arg only_sqlite: whether to write only the sqlite file instead of everything.

From 4c7ed9e57f010c00a05054479c7ebc4d130755f7 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 13 Mar 2022 15:22:39 +0000
Subject: [PATCH 29/73] merge clean up

---
 prody/dynamics/analysis.py  | 17 ++++++++++++++---
 prody/dynamics/functions.py | 37 +++++++++++++++++++++++++++++++++++--
 prody/proteins/starfile.py  |  2 +-
 3 files changed, 50 insertions(+), 6 deletions(-)

diff --git a/prody/dynamics/analysis.py b/prody/dynamics/analysis.py
index 234a4e109..91328c2ef 100644
--- a/prody/dynamics/analysis.py
+++ b/prody/dynamics/analysis.py
@@ -23,7 +23,7 @@
            'calcProjection', 'calcCrossProjection',
            'calcSpecDimension', 'calcPairDeformationDist',
            'calcDistFlucts', 'calcHinges', 'calcHitTime', 'calcHitTime',
-           'calcAnisousFromModel', 'calcScipionScore']
+           'calcAnisousFromModel', 'calcScipionScore', 'calcHemnmaScore']
            #'calcEntropyTransfer', 'calcOverallNetEntropyTransfer']
 
 def calcCollectivity(mode, masses=None, is3d=None):
@@ -711,8 +711,17 @@ def calcAnisousFromModel(model, ):
 
 
 def calcScipionScore(modes):
-    """Calculate Scipion continuousflex score, 
-    which is a function of mode number and collectivity order.
+    """Calculate the score from hybrid electron microscopy normal mode analysis (HEMNMA) 
+    [CS14]_ as implemented in the Scipion continuousflex plugin [MH20]_. This score 
+    prioritises modes as a function of mode number and collectivity order.
+
+    .. [CS14] Sorzano COS, de la Rosa-Trevín JM, Tama F, Jonić S.
+       Hybrid Electron Microscopy Normal Mode Analysis graphical interface and protocol.
+       *J Struct Biol* **2014** 188:134-41.
+
+    .. [MH20] Harastani M, Sorzano COS, Jonić S. 
+       Hybrid Electron Microscopy Normal Mode Analysis with Scipion.
+       *Protein Sci* **2020** 29:223-236.
 
     :arg modes: mode(s) or vector(s)
     :type modes: :class:`.Mode`, :class:`.Vector`, :class:`.ModeSet`, :class:`.NMA`
@@ -737,3 +746,5 @@ def calcScipionScore(modes):
     score = score / (2.0 * n_modes) 
 
     return score
+
+calcHemnmaScore = calcScipionScore
\ No newline at end of file
diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index fa458adfd..658d8fc43 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -14,14 +14,13 @@
 from prody.utilities import openFile, openSQLite, isExecutable, which, PLATFORM, addext, wrapModes
 from prody.proteins.starfile import parseSTAR, writeSTAR
 from prody.proteins import alignChains
-from prody.utilities import openFile, isExecutable, which, PLATFORM, addext
-
 from prody.ensemble import PDBEnsemble
 
 from .nma import NMA, MaskedNMA
 from .anm import ANM, ANMBase, MaskedANM
 from .analysis import calcCollectivity, calcScipionScore
 from .analysis import calcProjection
+from .analysis import calcCollectivity
 from .gnm import GNM, GNMBase, ZERO, MaskedGNM
 from .exanm import exANM, MaskedExANM
 from .rtb import RTB
@@ -354,12 +353,21 @@ def parseScipionModes(run_path, title=None):
         if found_eigvals:
             eigvals[i+1] = float(row['_nmaEigenval'])
     
+<<<<<<< HEAD
     if not found_eigvals:
         log_fname = run_path + '/logs/run.stdout'
         fi = open(log_fname, 'r')
         lines = fi.readlines()
         fi.close()
 
+=======
+    log_fname = run_path + '/logs/run.stdout'
+    fi = open(log_fname, 'r')
+    lines = fi.readlines()
+    fi.close()
+    
+    if not found_eigvals:
+>>>>>>> e82d36ae189edc248fb140de4637cd9326114032
         for line in lines:
             if line.find('Eigenvector number') != -1:
                 j = int(line.strip().split()[-1]) - 1
@@ -379,9 +387,13 @@ def parseScipionModes(run_path, title=None):
     nma.setEigens(vectors, eigvals)
     return nma
 
+<<<<<<< HEAD
 
 def writeScipionModes(output_path, modes, write_star=False, scores=None,
                       only_sqlite=False, collectivityThreshold=0.):
+=======
+def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sqlite=False, collectivityThreshold=0.):
+>>>>>>> e82d36ae189edc248fb140de4637cd9326114032
     """Writes *modes* to a set of files that can be recognised by Scipion.
     A directory called **"modes"** will be created if it doesn't already exist. 
     Filenames inside will start with **"vec"** and have the mode number as the extension.
@@ -397,7 +409,11 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
     :type write_star: bool
 
     :arg scores: scores from qualifyModesStep for re-writing sqlite
+<<<<<<< HEAD
         Default is **None** and then it uses :func:`.calcScipionScore`
+=======
+        Default is **None** and then it writes 0 for each.
+>>>>>>> e82d36ae189edc248fb140de4637cd9326114032
     :type scores: list
 
     :arg only_sqlite: whether to write only the sqlite file instead of everything.
@@ -418,6 +434,11 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
     if not isinstance(modes, (NMA, ModeSet, VectorBase)):
         raise TypeError('rows must be NMA, ModeSet, or Mode, not {0}'
                         .format(type(modes)))
+<<<<<<< HEAD
+=======
+    if not modes.is3d():
+        raise ValueError('modes must be 3-dimensional')
+>>>>>>> e82d36ae189edc248fb140de4637cd9326114032
 
     if not isinstance(write_star, bool):
         raise TypeError('write_star should be boolean, not {0}'
@@ -453,7 +474,11 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
                                         mode.getArrayNx3(), '%12.4e', ''))
         else:
             modefiles.append(writeArray(modes_dir + 'vec.{0}'.format(mode_num),
+<<<<<<< HEAD
                                         mode.getArray(), '%12.4e', ''))            
+=======
+                                        mode.getArray(), '%12.4e', ''))
+>>>>>>> e82d36ae189edc248fb140de4637cd9326114032
 
     if modes.numModes() > 1:
         order = modes.getIndices()
@@ -462,7 +487,11 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
         enabled = [1 if eigval > ZERO and collectivities[i] > collectivityThreshold else -1
                    for i, eigval in enumerate(eigvals)]
         if scores is None:
+<<<<<<< HEAD
             scores = list(calcScipionScore(modes))
+=======
+            scores = [0. for eigval in modes.getEigvals()]
+>>>>>>> e82d36ae189edc248fb140de4637cd9326114032
     else:
         mode = modes[0]
         eigvals = np.array([mode.getEigval()])
@@ -470,7 +499,11 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
         order = [mode.getIndex()]
         enabled = [1 if mode.getEigval() > ZERO and collectivities[0] > collectivityThreshold else -1]
         if scores is None:
+<<<<<<< HEAD
             scores = [calcScipionScore(mode)[0]]
+=======
+            scores = [0.]
+>>>>>>> e82d36ae189edc248fb140de4637cd9326114032
 
     modes_sqlite_fn = output_path + '/modes.sqlite'
     sql_con = openSQLite(modes_sqlite_fn, 'n')
diff --git a/prody/proteins/starfile.py b/prody/proteins/starfile.py
index d4d2950a1..ea6011ef4 100644
--- a/prody/proteins/starfile.py
+++ b/prody/proteins/starfile.py
@@ -695,7 +695,7 @@ def writeSTAR(filename, starDict, **kwargs):
         field names and finally data.
     :type starDict: dict
 
-    kwargs can be given to initiate a :class:`.StarDict` from a dict
+    kwargs can be given including the program style to follow (*prog*)
     """
     prog=kwargs.get('prog', 'XMIPP')
 

From 7419183b26bb8a38365c1480a9ca350d8e9a0bb7 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 13 Mar 2022 15:29:53 +0000
Subject: [PATCH 30/73] merge clean up

---
 prody/dynamics/functions.py | 34 ----------------------------------
 1 file changed, 34 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 658d8fc43..34d11dfa7 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -353,21 +353,12 @@ def parseScipionModes(run_path, title=None):
         if found_eigvals:
             eigvals[i+1] = float(row['_nmaEigenval'])
     
-<<<<<<< HEAD
     if not found_eigvals:
         log_fname = run_path + '/logs/run.stdout'
         fi = open(log_fname, 'r')
         lines = fi.readlines()
         fi.close()
 
-=======
-    log_fname = run_path + '/logs/run.stdout'
-    fi = open(log_fname, 'r')
-    lines = fi.readlines()
-    fi.close()
-    
-    if not found_eigvals:
->>>>>>> e82d36ae189edc248fb140de4637cd9326114032
         for line in lines:
             if line.find('Eigenvector number') != -1:
                 j = int(line.strip().split()[-1]) - 1
@@ -387,13 +378,9 @@ def parseScipionModes(run_path, title=None):
     nma.setEigens(vectors, eigvals)
     return nma
 
-<<<<<<< HEAD
 
 def writeScipionModes(output_path, modes, write_star=False, scores=None,
                       only_sqlite=False, collectivityThreshold=0.):
-=======
-def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sqlite=False, collectivityThreshold=0.):
->>>>>>> e82d36ae189edc248fb140de4637cd9326114032
     """Writes *modes* to a set of files that can be recognised by Scipion.
     A directory called **"modes"** will be created if it doesn't already exist. 
     Filenames inside will start with **"vec"** and have the mode number as the extension.
@@ -409,11 +396,7 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
     :type write_star: bool
 
     :arg scores: scores from qualifyModesStep for re-writing sqlite
-<<<<<<< HEAD
         Default is **None** and then it uses :func:`.calcScipionScore`
-=======
-        Default is **None** and then it writes 0 for each.
->>>>>>> e82d36ae189edc248fb140de4637cd9326114032
     :type scores: list
 
     :arg only_sqlite: whether to write only the sqlite file instead of everything.
@@ -434,11 +417,6 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
     if not isinstance(modes, (NMA, ModeSet, VectorBase)):
         raise TypeError('rows must be NMA, ModeSet, or Mode, not {0}'
                         .format(type(modes)))
-<<<<<<< HEAD
-=======
-    if not modes.is3d():
-        raise ValueError('modes must be 3-dimensional')
->>>>>>> e82d36ae189edc248fb140de4637cd9326114032
 
     if not isinstance(write_star, bool):
         raise TypeError('write_star should be boolean, not {0}'
@@ -474,11 +452,7 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
                                         mode.getArrayNx3(), '%12.4e', ''))
         else:
             modefiles.append(writeArray(modes_dir + 'vec.{0}'.format(mode_num),
-<<<<<<< HEAD
-                                        mode.getArray(), '%12.4e', ''))            
-=======
                                         mode.getArray(), '%12.4e', ''))
->>>>>>> e82d36ae189edc248fb140de4637cd9326114032
 
     if modes.numModes() > 1:
         order = modes.getIndices()
@@ -487,11 +461,7 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
         enabled = [1 if eigval > ZERO and collectivities[i] > collectivityThreshold else -1
                    for i, eigval in enumerate(eigvals)]
         if scores is None:
-<<<<<<< HEAD
             scores = list(calcScipionScore(modes))
-=======
-            scores = [0. for eigval in modes.getEigvals()]
->>>>>>> e82d36ae189edc248fb140de4637cd9326114032
     else:
         mode = modes[0]
         eigvals = np.array([mode.getEigval()])
@@ -499,11 +469,7 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None, only_sq
         order = [mode.getIndex()]
         enabled = [1 if mode.getEigval() > ZERO and collectivities[0] > collectivityThreshold else -1]
         if scores is None:
-<<<<<<< HEAD
             scores = [calcScipionScore(mode)[0]]
-=======
-            scores = [0.]
->>>>>>> e82d36ae189edc248fb140de4637cd9326114032
 
     modes_sqlite_fn = output_path + '/modes.sqlite'
     sql_con = openSQLite(modes_sqlite_fn, 'n')

From 50690d73c4a082d41b58cf11644f2786026e3b0f Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 15 Mar 2022 13:33:04 +0000
Subject: [PATCH 31/73] fix missing dir in parseScipionModes

---
 prody/dynamics/functions.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 34d11dfa7..9e9045d0a 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -326,7 +326,10 @@ def parseScipionModes(run_path, title=None):
     :arg title: title for :class:`.NMA` object
     :type title: str
     """
+    if run_path.endswith("/"):
+        run_path = run_path[:-1]
     run_name = os.path.split(run_path)[-1]
+    top_dirs = os.path.split(run_path)[0][:-4] # exclude "Runs"
 
     star_data = parseSTAR(run_path + '/modes.xmd')
     star_loop = star_data[0][0]
@@ -334,7 +337,7 @@ def parseScipionModes(run_path, title=None):
     n_modes = star_loop.numRows()
     
     row1 = star_loop[0]
-    mode1 = parseArray(row1['_nmaModefile']).reshape(-1)
+    mode1 = parseArray(top_dirs + row1['_nmaModefile']).reshape(-1)
     dof = mode1.shape[0]
 
     vectors = np.zeros((dof, n_modes))
@@ -349,7 +352,7 @@ def parseScipionModes(run_path, title=None):
         found_eigvals = False
 
     for i, row in enumerate(star_loop[1:]):
-        vectors[:, i+1] = parseArray(row['_nmaModefile']).reshape(-1)
+        vectors[:, i+1] = parseArray(top_dirs + row['_nmaModefile']).reshape(-1)
         if found_eigvals:
             eigvals[i+1] = float(row['_nmaEigenval'])
     

From 8125f8517820687817ec70f089122b8a2fa6c45c Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 28 Mar 2022 14:50:21 +0200
Subject: [PATCH 32/73] clustenm accept sparse kdtree turbo

---
 prody/dynamics/clustenm.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/prody/dynamics/clustenm.py b/prody/dynamics/clustenm.py
index c2c73592c..2b166edd7 100644
--- a/prody/dynamics/clustenm.py
+++ b/prody/dynamics/clustenm.py
@@ -478,7 +478,8 @@ def _multi_targeted_sim(self, args):
     def _buildANM(self, cg):
 
         anm = ANM()
-        anm.buildHessian(cg, cutoff=self._cutoff, gamma=self._gamma)
+        anm.buildHessian(cg, cutoff=self._cutoff, gamma=self._gamma,
+                         sparse=self._sparse, kdtree=self._kdtree)
 
         return anm
 
@@ -513,7 +514,7 @@ def _sample(self, conf):
         if not self._checkANM(anm_cg):
             return None
 
-        anm_cg.calcModes(self._n_modes)
+        anm_cg.calcModes(self._n_modes, turbo=self._turbo)
 
         anm_ex = self._extendModel(anm_cg, cg, tmp)
         ens_ex = sampleModes(anm_ex, atoms=tmp,
@@ -980,6 +981,12 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
         self._cutoff = cutoff
         self._n_modes = n_modes
         self._gamma = gamma
+        self._sparse = kwargs.get('sparse', False)
+        self._kdtree = kwargs.get('kdtree', False)
+        self._turbo = kwargs.get('turbo', False)
+        if kwargs.get('zeros', False):
+            LOGGER.warn('ClustENM cannot use zero modes so ignoring this kwarg')
+
         self._n_confs = n_confs
         self._rmsd = (0.,) + rmsd if isinstance(rmsd, tuple) else (0.,) + (rmsd,) * n_gens
         self._n_gens = n_gens

From dd65f366e7beeabb0a08ea3f63fcc651e09fd29e Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 28 Mar 2022 14:51:23 +0200
Subject: [PATCH 33/73] anm app accept sparse kdtree turbo

---
 prody/apps/prody_apps/prody_anm.py | 21 +++++++++++++++++++--
 1 file changed, 19 insertions(+), 2 deletions(-)

diff --git a/prody/apps/prody_apps/prody_anm.py b/prody/apps/prody_apps/prody_anm.py
index 8cfec1397..4aadb4b69 100644
--- a/prody/apps/prody_apps/prody_anm.py
+++ b/prody/apps/prody_apps/prody_anm.py
@@ -15,7 +15,10 @@
     ('altloc', 'alternative location identifiers for residues used in the calculations', "A"),
     ('cutoff', 'cutoff distance (A)', 15.),
     ('gamma', 'spring constant', 1.),
+    ('sparse', 'use sparse matrices', False),
+    ('kdtree', 'use kdtree for Hessian', False),    
     ('zeros', 'calculate zero modes', False),
+    ('turbo', 'use memory-intensive turbo option for modes', False),
 
     ('outbeta', 'write beta-factors calculated from GNM modes', False),
     ('hessian', 'write Hessian matrix', False),
@@ -56,6 +59,8 @@ def prody_anm(pdb, **kwargs):
     prefix = kwargs.get('prefix')
     cutoff = kwargs.get('cutoff')
     gamma = kwargs.get('gamma')
+    sparse = kwargs.get('sparse')
+    kdtree = kwargs.get('kdtree')
     nmodes = kwargs.get('nmodes')
     selstr = kwargs.get('select')
     model = kwargs.get('model')
@@ -75,8 +80,8 @@ def prody_anm(pdb, **kwargs):
                 .format(len(select)))
 
     anm = prody.ANM(pdb.getTitle())
-    anm.buildHessian(select, cutoff, gamma)
-    anm.calcModes(nmodes, zeros=zeros)
+    anm.buildHessian(select, cutoff, gamma, sparse=sparse, kdtree=kdtree)
+    anm.calcModes(nmodes, zeros=zeros, turbo=turbo)
     LOGGER.info('Writing numerical output.')
 
     if kwargs.get('outnpz'):
@@ -253,6 +258,18 @@ def addCommand(commands):
         default=DEFAULTS['gamma'], metavar='FLOAT',
         help=HELPTEXT['gamma'] + ' (default: %(default)s)')
 
+    group.add_argument('-C', '--sparse-hessian', dest='sparse', type=bool,
+        default=DEFAULTS['sparse'],
+        help=HELPTEXT['sparse'] + ' (default: %(default)s)')
+
+    group.add_argument('-G', '--use-kdtree', dest='kdtree', type=bool,
+        default=DEFAULTS['kdtree'],
+        help=HELPTEXT['kdtree'] + ' (default: %(default)s)')
+
+    group.add_argument('-y', '--turbo', dest='turbo', type=bool,
+        default=DEFAULTS['turbo'],
+        help=HELPTEXT['turbo'] + ' (default: %(default)s)')
+
     group.add_argument('-m', '--model', dest='model', type=int,
         metavar='INT', default=DEFAULTS['model'], help=HELPTEXT['model'])
 

From 8a965f82ea3e5685fb75b1eee46e6146b44bf3d9 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 28 Mar 2022 14:51:57 +0200
Subject: [PATCH 34/73] assignBlocks shortest default 4

---
 prody/measure/measure.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/prody/measure/measure.py b/prody/measure/measure.py
index 51f4fea0e..047204137 100644
--- a/prody/measure/measure.py
+++ b/prody/measure/measure.py
@@ -891,7 +891,7 @@ def assignBlocks(atoms, res_per_block=None, secstr=False, **kwargs):
 
     :arg shortest_block: smallest number of residues to be included 
         in a block before merging with the previous block
-        Default is **2**
+        Default is **4** as smaller numbers can cause problems for distance matrices.
     :type shortest_block: int
 
     :arg longest_block: largest number of residues to be included 
@@ -922,7 +922,7 @@ def assignBlocks(atoms, res_per_block=None, secstr=False, **kwargs):
     if not isinstance(secstr, bool):
         raise TypeError('secstr should be a Boolean')
 
-    shortest_block = kwargs.get('shortest_block', 2)
+    shortest_block = kwargs.get('shortest_block', 4)
     shortest_block = kwargs.get('min_size', shortest_block)
     if not isinstance(shortest_block, Integral):
         raise TypeError("shortest_block should be an integer")

From 10bcef09f50f6cdfcd4d10b10e881c2eeb556133 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 28 Mar 2022 20:50:25 +0200
Subject: [PATCH 35/73] clustenm threads for sims

---
 prody/dynamics/clustenm.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/prody/dynamics/clustenm.py b/prody/dynamics/clustenm.py
index 2b166edd7..10599d061 100644
--- a/prody/dynamics/clustenm.py
+++ b/prody/dynamics/clustenm.py
@@ -291,6 +291,12 @@ def _prep_sim(self, coords, external_forces=[]):
             properties = {'Precision': 'single'}
         elif self._platform in ['CUDA', 'OpenCL']:
             properties = {'Precision': 'single'}
+        elif self._platform == 'CPU':
+            if self._threads == 0:
+                cpus = cpu_count()
+            else:
+                cpus = self._threads
+            properties = {'Threads': str(cpus)}
 
         simulation = Simulation(modeller.topology, system, integrator,
                                 platform, properties)
@@ -968,10 +974,15 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
         :arg platform: Architecture on which the OpenMM part runs, default is None.
             It can be chosen as 'CUDA', 'OpenCL' or 'CPU'.
             For efficiency, 'CUDA' or 'OpenCL' is recommended.
+            'CPU' is needed for setting threads per simulation.
         :type platform: str
 
         :arg parallel: If it is True (default is False), conformer generation will be parallelized.
         :type parallel: bool
+
+        :arg threads: Number of threads to use for an individual simulation
+            Default of 0 uses all CPUs on the machine.
+        :type threads: int
         '''
 
         if self._isBuilt():
@@ -992,6 +1003,7 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
         self._n_gens = n_gens
         self._platform = kwargs.pop('platform', None)
         self._parallel = kwargs.pop('parallel', False)
+        self._threads = kwargs.pop('threads', 0)
         self._targeted = kwargs.pop('targeted', False)
         self._tmdk = kwargs.pop('tmdk', 15.)
 

From 10928cf1d1e6d0e045aec09f0e342ac889d0302a Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 28 Mar 2022 23:16:25 +0200
Subject: [PATCH 36/73] fix pointer psf iters

---
 prody/atomic/pointer.py | 106 +++++++++++++++++++++-------------------
 1 file changed, 57 insertions(+), 49 deletions(-)

diff --git a/prody/atomic/pointer.py b/prody/atomic/pointer.py
index e556e3938..6740dbdcf 100644
--- a/prody/atomic/pointer.py
+++ b/prody/atomic/pointer.py
@@ -2,7 +2,7 @@
 """This module defines atom pointer base class."""
 
 from numbers import Integral
-from numpy import all, array, concatenate, ones, unique
+from numpy import all, array, concatenate, ones, unique, any
 
 from .atomic import Atomic
 from .bond import Bond
@@ -279,7 +279,7 @@ def _iterBonds(self):
         Use :meth:`setBonds` for setting bonds."""
 
         if self._ag._bonds is None:
-            raise ValueError('bonds are not set, use `setBonds` or `inferBonds`')
+            LOGGER.warning('bonds are not set, use `setBonds` or `inferBonds`')
 
         indices = self._getIndices()
         iset = set(indices)
@@ -288,11 +288,12 @@ def _iterBonds(self):
                 if a in iset and b in iset:
                     yield a, b
         else:
-            for a, bmap in zip(indices, self._ag._bmap[indices]):
-                for b in bmap:
-                    if b > -1 and b in iset:
-                        yield a, b
-                iset.remove(a)
+            if any(self._ag._bmap):
+                for a, bmap in zip(indices, self._ag._bmap[indices]):
+                    for b in bmap:
+                        if b > -1 and b in iset:
+                            yield a, b
+                    iset.remove(a)
 
     def iterBonds(self):
         """Yield bonds formed by the atom.  Use :meth:`setBonds` or 
@@ -308,7 +309,7 @@ def _iterAngles(self):
         Use :meth:`setAngles` for setting angles."""
 
         if self._ag._angles is None:
-            raise ValueError('angles are not set, use `AtomGroup.setAngles`')
+            LOGGER.warning('angles are not set, use `AtomGroup.setAngles`')
 
         indices = self._getIndices()
         iset = set(indices)
@@ -317,11 +318,12 @@ def _iterAngles(self):
                 if a in iset and b in iset and c in iset:
                     yield a, b, c
         else:
-            for a, amap in zip(indices, self._ag._angmap[indices]):
-                for b, c in amap:
-                    if b > -1 and b in iset and c > -1 and c in iset:
-                        yield a, b, c
-                iset.remove(a)
+            if any(self._ag._angmap):
+                for a, amap in zip(indices, self._ag._angmap[indices]):
+                    for b, c in amap:
+                        if b > -1 and b in iset and c > -1 and c in iset:
+                            yield a, b, c
+                    iset.remove(a)
 
     def iterAngles(self):
         """Yield angles formed by the atom.  Use :meth:`setAngles` for setting
@@ -337,7 +339,7 @@ def _iterDihedrals(self):
         Use :meth:`setDihedrals` for setting dihedrals."""
 
         if self._ag._dihedrals is None:
-            raise ValueError('dihedrals are not set, use `AtomGroup.setDihedrals`')
+            LOGGER.warning('dihedrals are not set, use `AtomGroup.setDihedrals`')
 
         indices = self._getIndices()
         iset = set(indices)
@@ -346,12 +348,13 @@ def _iterDihedrals(self):
                 if a in iset and b in iset and c in iset and d in iset:
                     yield a, b, c, d
         else:
-            for a, dmap in zip(indices, self._ag._dmap[indices]):
-                for b, c, d in dmap:
-                    if b > -1 and b in iset and c > -1 and c in iset \
-                    and d > -1 and d in iset:
-                        yield a, b, c, d
-                iset.remove(a)
+            if any(self._ag._dmap):
+                for a, dmap in zip(indices, self._ag._dmap[indices]):
+                    for b, c, d in dmap:
+                        if b > -1 and b in iset and c > -1 and c in iset \
+                        and d > -1 and d in iset:
+                            yield a, b, c, d
+                    iset.remove(a)
 
     def iterDihedrals(self):
         """Yield dihedrals formed by the atom.  Use :meth:`setDihedrals` for setting
@@ -367,7 +370,7 @@ def _iterImpropers(self):
         Use :meth:`setImpropers` for setting impropers."""
 
         if self._ag._impropers is None:
-            raise ValueError('impropers are not set, use `AtomGroup.setImpropers`')
+            LOGGER.warning('impropers are not set, use `AtomGroup.setImpropers`')
 
         indices = self._getIndices()
         iset = set(indices)
@@ -376,12 +379,13 @@ def _iterImpropers(self):
                 if a in iset and b in iset and c in iset and d in iset:
                     yield a, b, c, d
         else:
-            for a, imap in zip(indices, self._ag._imap[indices]):
-                for b, c, d in imap:
-                    if b > -1 and b in iset and c > -1 and c in iset \
-                    and d > -1 and d in iset:
-                        yield a, b, c, d
-                iset.remove(a)
+            if any(self._ag._imap):
+                for a, imap in zip(indices, self._ag._imap[indices]):
+                    for b, c, d in imap:
+                        if b > -1 and b in iset and c > -1 and c in iset \
+                        and d > -1 and d in iset:
+                            yield a, b, c, d
+                    iset.remove(a)
 
     def iterImpropers(self):
         """Yield impropers formed by the atom.  Use :meth:`setImpropers` for setting
@@ -397,7 +401,7 @@ def _iterCrossterms(self):
         Use :meth:`setCrossterms` for setting crossterms."""
 
         if self._ag._crossterms is None:
-            raise ValueError('crossterms are not set, use `AtomGroup.setCrossterms`')
+            LOGGER.warning('crossterms are not set, use `AtomGroup.setCrossterms`')
 
         indices = self._getIndices()
         iset = set(indices)
@@ -406,12 +410,13 @@ def _iterCrossterms(self):
                 if a in iset and b in iset and c in iset and d in iset:
                     yield a, b, c, d
         else:
-            for a, cmap in zip(indices, self._ag._cmap[indices]):
-                for b, c, d in cmap:
-                    if b > -1 and b in iset and c > -1 and c in iset \
-                    and d > -1 and d in iset:
-                        yield a, b, c, d
-                iset.remove(a)
+            if any(self._ag._cmap):
+                for a, cmap in zip(indices, self._ag._cmap[indices]):
+                    for b, c, d in cmap:
+                        if b > -1 and b in iset and c > -1 and c in iset \
+                        and d > -1 and d in iset:
+                            yield a, b, c, d
+                    iset.remove(a)
 
     def iterCrossterms(self):
         """Yield crossterms formed by the atom.  Use :meth:`setCrossterms` for setting
@@ -436,11 +441,12 @@ def _iterDonors(self):
                 if a in iset and b in iset:
                     yield a, b
         else:
-            for a, dmap in zip(indices, self._ag._domap[indices]):
-                for b in dmap:
-                    if b > -1 and b in iset:
-                        yield a, b
-                iset.remove(a)
+            if any(self._ag._domap):
+                for a, dmap in zip(indices, self._ag._domap[indices]):
+                    for b in dmap:
+                        if b > -1 and b in iset:
+                            yield a, b
+                    iset.remove(a)
 
     def iterDonors(self):
         """Yield donors formed by the atom.  Use :meth:`setDonors` for setting
@@ -465,11 +471,12 @@ def _iterAcceptors(self):
                 if a in iset and b in iset:
                     yield a, b
         else:
-            for a, amap in zip(indices, self._ag._acmap[indices]):
-                for b in amap:
-                    if b > -1 and b in iset:
-                        yield a, b
-                iset.remove(a)
+            if any(self._ag._acmap):
+                for a, amap in zip(indices, self._ag._acmap[indices]):
+                    for b in amap:
+                        if b > -1 and b in iset:
+                            yield a, b
+                    iset.remove(a)
 
     def iterAcceptors(self):
         """Yield acceptors formed by the atom.  Use :meth:`setAcceptors` for setting
@@ -494,11 +501,12 @@ def _iterNBExclusions(self):
                 if a in iset and b in iset:
                     yield a, b
         else:
-            for a, nbemap in zip(indices, self._ag._nbemap[indices]):
-                for b in nbemap:
-                    if b > -1 and b in iset:
-                        yield a, b
-                iset.remove(a)
+            if any(self._ag._nbemap):
+                for a, nbemap in zip(indices, self._ag._nbemap[indices]):
+                    for b in nbemap:
+                        if b > -1 and b in iset:
+                            yield a, b
+                    iset.remove(a)
 
     def iterNBExclusions(self):
         """Yield nbexclusions formed by the atom.  Use :meth:`setNBExclusions` for setting

From f47486fa9a68b06bdf354c21baf2fa04bb7f9cc0 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 28 Mar 2022 23:18:22 +0200
Subject: [PATCH 37/73] clustenm num of parallel

---
 prody/dynamics/clustenm.py | 27 +++++++++++++++++++++++++--
 1 file changed, 25 insertions(+), 2 deletions(-)

diff --git a/prody/dynamics/clustenm.py b/prody/dynamics/clustenm.py
index 10599d061..542e0b890 100644
--- a/prody/dynamics/clustenm.py
+++ b/prody/dynamics/clustenm.py
@@ -530,7 +530,12 @@ def _sample(self, conf):
 
         if self._targeted:
             if self._parallel:
-                with Pool(cpu_count()) as p:
+                if isinstance(self._parallel, int):
+                    repeats = self._parallel
+                elif isinstance(self._parallel, bool):
+                    repeats = cpu_count()
+
+                with Pool(repeats) as p:
                     pot_conf = p.map(self._multi_targeted_sim,
                                      [(conf, coords) for coords in coordsets])
             else:
@@ -614,7 +619,12 @@ def _generate(self, confs):
         sample_method = self._sample_v1 if self._v1 else self._sample
 
         if self._parallel:
-            with Pool(cpu_count()) as p:
+            if isinstance(self._parallel, int):
+                repeats = self._parallel
+            elif isinstance(self._parallel, bool):
+                repeats = cpu_count()
+
+            with Pool(repeats) as p:
                 tmp = p.map(sample_method, [conf for conf in confs])
         else:
             tmp = [sample_method(conf) for conf in confs]
@@ -978,6 +988,8 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
         :type platform: str
 
         :arg parallel: If it is True (default is False), conformer generation will be parallelized.
+            This can also be set to a number for how many simulations are run in parallel.
+            Setting 0 or True means run as many as there are CPUs on the machine.
         :type parallel: bool
 
         :arg threads: Number of threads to use for an individual simulation
@@ -1002,7 +1014,18 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
         self._rmsd = (0.,) + rmsd if isinstance(rmsd, tuple) else (0.,) + (rmsd,) * n_gens
         self._n_gens = n_gens
         self._platform = kwargs.pop('platform', None)
+
         self._parallel = kwargs.pop('parallel', False)
+        if not isinstance(self._parallel, (int, bool)):
+            raise TypeError('parallel should be an int or bool')
+
+        if self._parallel == 1:
+            # this is equivalent to not actually being parallel
+            self._parallel = False
+        elif self._parallel == 0:
+            # this is a deliberate choice and is documented
+            self._parallel = True
+
         self._threads = kwargs.pop('threads', 0)
         self._targeted = kwargs.pop('targeted', False)
         self._tmdk = kwargs.pop('tmdk', 15.)

From 17b6b129aa668d0fba1ec515ce7cbdf4bae72d57 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 29 Mar 2022 01:40:43 +0200
Subject: [PATCH 38/73] anm kwargs store_true

---
 prody/apps/prody_apps/prody_anm.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/prody/apps/prody_apps/prody_anm.py b/prody/apps/prody_apps/prody_anm.py
index 29610098b..701d7b297 100644
--- a/prody/apps/prody_apps/prody_anm.py
+++ b/prody/apps/prody_apps/prody_anm.py
@@ -271,15 +271,15 @@ def addCommand(commands):
         default=DEFAULTS['gamma'], metavar='FLOAT',
         help=HELPTEXT['gamma'] + ' (default: %(default)s)')
 
-    group.add_argument('-C', '--sparse-hessian', dest='sparse', type=bool,
+    group.add_argument('-C', '--sparse-hessian', dest='sparse', action='store_true',
         default=DEFAULTS['sparse'],
         help=HELPTEXT['sparse'] + ' (default: %(default)s)')
 
-    group.add_argument('-G', '--use-kdtree', dest='kdtree', type=bool,
+    group.add_argument('-G', '--use-kdtree', dest='kdtree', action='store_true',
         default=DEFAULTS['kdtree'],
         help=HELPTEXT['kdtree'] + ' (default: %(default)s)')
 
-    group.add_argument('-y', '--turbo', dest='turbo', type=bool,
+    group.add_argument('-y', '--turbo', dest='turbo', action='store_true',
         default=DEFAULTS['turbo'],
         help=HELPTEXT['turbo'] + ' (default: %(default)s)')
 

From aeb8ceb74a2b28ac7a8f79d89003e8c0d98104e5 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 29 Mar 2022 17:36:16 +0200
Subject: [PATCH 39/73] missed MRC2015 to 14

---
 prody/proteins/emdfile.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/proteins/emdfile.py b/prody/proteins/emdfile.py
index 434d65687..9cd457878 100644
--- a/prody/proteins/emdfile.py
+++ b/prody/proteins/emdfile.py
@@ -26,7 +26,7 @@ class EMDParseError(Exception):
 
 
 def parseEMD(emd, **kwargs):
-    """Parses an EM density map in EMD/MRC2015 format and 
+    """Parses an EM density map in EMD/MRC2014 format and 
     optionally returns an :class:`.AtomGroup` containing  
     beads built in the density using the TRN algorithm [_TM94]. 
 

From 04ea96c48eb501d88aa5ec081ea7e2879652809f Mon Sep 17 00:00:00 2001
From: Ricardo Serrano <rserranogut@hotmail.com>
Date: Thu, 31 Mar 2022 14:04:23 +0200
Subject: [PATCH 40/73] Fixing problem with parseScipionModes when using GNM

---
 prody/dynamics/functions.py | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 9e9045d0a..e33c18e69 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -12,8 +12,7 @@
 from prody import LOGGER, SETTINGS, PY3K
 from prody.atomic import Atomic, AtomSubset
 from prody.utilities import openFile, openSQLite, isExecutable, which, PLATFORM, addext, wrapModes
-from prody.proteins.starfile import parseSTAR, writeSTAR
-from prody.proteins import alignChains
+from prody.proteins import parseSTAR, writeSTAR, alignChains, parsePDB
 from prody.ensemble import PDBEnsemble
 
 from .nma import NMA, MaskedNMA
@@ -316,7 +315,7 @@ def parseModes(normalmodes, eigenvalues=None, nm_delimiter=None,
     return nma
 
 
-def parseScipionModes(run_path, title=None):
+def parseScipionModes(run_path, title=None, pdb=None):
     """Returns :class:`.NMA` containing eigenvectors and eigenvalues 
     parsed from a ContinuousFlex FlexProtNMA Run directory.
 
@@ -336,6 +335,9 @@ def parseScipionModes(run_path, title=None):
     
     n_modes = star_loop.numRows()
     
+    atoms = parsePDB(pdb)
+    n_atoms = atoms.numAtoms()
+
     row1 = star_loop[0]
     mode1 = parseArray(top_dirs + row1['_nmaModefile']).reshape(-1)
     dof = mode1.shape[0]
@@ -377,7 +379,11 @@ def parseScipionModes(run_path, title=None):
         LOGGER.warn('No eigenvalues found')
         eigvals=None
 
-    nma = NMA(title)
+    if dof == n_atoms * 3:
+        nma = NMA(title)
+    else:
+        nma = GNM(title)
+
     nma.setEigens(vectors, eigvals)
     return nma
 

From 20fbe24f34670e105c9f8bf3d5a8432262cd7f8d Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 14 Apr 2022 17:06:36 +0100
Subject: [PATCH 41/73] parsePDB consistent logging

---
 prody/proteins/pdbfile.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/prody/proteins/pdbfile.py b/prody/proteins/pdbfile.py
index da5eeb553..4cb76a241 100644
--- a/prody/proteins/pdbfile.py
+++ b/prody/proteins/pdbfile.py
@@ -210,11 +210,11 @@ def _parsePDB(pdb, **kwargs):
         filename = fetchPDB(pdb, **kwargs)
         if filename is None:
             try:
-                LOGGER.info("Trying to parse mmCIF file instead")
+                LOGGER.warn("Trying to parse mmCIF file instead")
                 return parseMMCIF(pdb+chain, **kwargs)
             except:
                 try:
-                    LOGGER.info("Trying to parse EMD file instead")
+                    LOGGER.warn("Trying to parse EMD file instead")
                     return parseEMD(pdb+chain, **kwargs)
                 except:                
                     raise IOError('PDB file for {0} could not be downloaded.'

From af6f6a15aec6cc87f6df9ad7d52152b93b1744a0 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 18 Apr 2022 19:54:04 +0100
Subject: [PATCH 42/73] parseScipionModes fix no pdb

---
 prody/dynamics/functions.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index e33c18e69..0859d69f9 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -335,13 +335,17 @@ def parseScipionModes(run_path, title=None, pdb=None):
     
     n_modes = star_loop.numRows()
     
-    atoms = parsePDB(pdb)
-    n_atoms = atoms.numAtoms()
-
     row1 = star_loop[0]
     mode1 = parseArray(top_dirs + row1['_nmaModefile']).reshape(-1)
     dof = mode1.shape[0]
 
+    if pdb is not None:
+        atoms = parsePDB(pdb)
+        n_atoms = atoms.numAtoms()
+    else:
+        # assume standard NMA
+        n_atoms = dof//3
+
     vectors = np.zeros((dof, n_modes))
     vectors[:, 0] = mode1
 

From cc1a42d3f17bb7ec82a0e2c8b63f0f039d841be3 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 10 May 2022 23:57:40 +0200
Subject: [PATCH 43/73] writeScipionModes writes eigenvals

---
 prody/dynamics/functions.py | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 0859d69f9..aa723c4c0 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -499,10 +499,11 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
     classes = [(1, 'self', 'c00', 'NormalMode'),
                (2, '_modeFile', 'c01', 'String'),
                (3, '_collectivity', 'c02', 'Float'),
-               (4, '_score', 'c03', 'Float')]
+               (4, '_score', 'c03', 'Float'),
+               (5, '_eigenval', 'c04', 'Float')]
     cursor.executemany('''INSERT INTO Classes VALUES(?,?,?,?);''', classes);
     
-    cursor.execute('''CREATE TABLE Objects(id primary key, enabled, label, comment, creation, c01, c02, c03)''')
+    cursor.execute('''CREATE TABLE Objects(id primary key, enabled, label, comment, creation, c01, c02, c03, c04)''')
     
     star_dict = OrderedDict()
 
@@ -542,14 +543,14 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
         c03 = scores[i]
         loop_dict['data'][i]['_nmaScore'] = '%8.6f' % c03
 
-        eigval = eigvals[i]
-        if float('%9.6f' % eigval) > 0:
-            loop_dict['data'][i]['_nmaEigenval'] = '%9.6f' % eigval
+        c04 = eigvals[i]
+        if float('%9.6f' % c04) > 0:
+            loop_dict['data'][i]['_nmaEigenval'] = '%9.6f' % c04
         else:
-            loop_dict['data'][i]['_nmaEigenval'] = '%9.6e' % eigval
+            loop_dict['data'][i]['_nmaEigenval'] = '%9.6e' % c04
         
-        cursor.execute('''INSERT INTO Objects VALUES(?,?,?,?,?,?,?,?)''',
-                       (id, enab, label, comment, creation, c01, c02, c03))
+        cursor.execute('''INSERT INTO Objects VALUES(?,?,?,?,?,?,?,?,?)''',
+                       (id, enab, label, comment, creation, c01, c02, c03, c04))
 
     if write_star:
         writeSTAR(output_path + '/modes.xmd', star_dict)

From f3dcdda9203449c7ffaa2772c94536d9c1213435 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 11 May 2022 00:30:48 +0200
Subject: [PATCH 44/73] parseScipionModes sqlite

---
 prody/dynamics/functions.py | 52 +++++++++++++++++++++++++++++++++++--
 prody/proteins/starfile.py  |  2 +-
 2 files changed, 51 insertions(+), 3 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index aa723c4c0..73cf470ac 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -12,7 +12,7 @@
 from prody import LOGGER, SETTINGS, PY3K
 from prody.atomic import Atomic, AtomSubset
 from prody.utilities import openFile, openSQLite, isExecutable, which, PLATFORM, addext, wrapModes
-from prody.proteins import parseSTAR, writeSTAR, alignChains, parsePDB
+from prody.proteins import parseSTAR, writeSTAR, StarDict, alignChains, parsePDB
 from prody.ensemble import PDBEnsemble
 
 from .nma import NMA, MaskedNMA
@@ -330,7 +330,55 @@ def parseScipionModes(run_path, title=None, pdb=None):
     run_name = os.path.split(run_path)[-1]
     top_dirs = os.path.split(run_path)[0][:-4] # exclude "Runs"
 
-    star_data = parseSTAR(run_path + '/modes.xmd')
+    try:
+        star_data = parseSTAR(run_path + '/modes.xmd')
+    except OSError:
+        try:
+            sql_con = openSQLite(run_path + "/modes.sqlite")
+            cursor = sql_con.cursor()
+
+            #names = ['id', 'enabled', 'label', 'comment', 'creation', 'c01', 'c02', 'c03', 'c04']
+
+            star_dict = OrderedDict()
+            star_block_dict = OrderedDict()
+            
+            star_loop_dict = OrderedDict()
+            star_loop_dict["fields"] = OrderedDict([(0, '_enabled'),
+                                               (1, '_nmaCollectivity'),
+                                               (2, '_nmaModefile'),
+                                               (3, '_nmaScore'),
+                                               (4, '_nmaEigenval'),
+                                               (5, '_order_')])
+            star_loop_dict["data"] = OrderedDict()
+
+            for row in cursor.execute("SELECT * FROM Objects;"):
+
+                id_ = row[0]
+                key = id_-1 # sqlite ids start from 1 not 0
+
+                star_loop_dict["data"][key] = OrderedDict()
+
+                star_loop_dict["data"][key]['_order_'] = id_
+
+                star_loop_dict["data"][key]['_enabled'] = row[1]
+
+                star_loop_dict["data"][key]['_nmaCollectivity'] = row[6]
+
+                star_loop_dict["data"][key]['_nmaModefile'] = row[5]
+
+                star_loop_dict["data"][key]['_nmaScore'] = row[7]
+
+                if len(row) > 8:
+                    star_loop_dict["data"][key]['_nmaEigenval'] = row[8]
+
+            star_block_dict[0] = star_loop_dict
+            star_dict[0] = star_block_dict
+
+            star_data = StarDict(star_dict, prog='XMIPP')
+
+        except OSError:
+            raise OSError("neither modes.xmd nor modes.sqlite found")
+
     star_loop = star_data[0][0]
     
     n_modes = star_loop.numRows()
diff --git a/prody/proteins/starfile.py b/prody/proteins/starfile.py
index e535852c2..c24dce22d 100644
--- a/prody/proteins/starfile.py
+++ b/prody/proteins/starfile.py
@@ -379,7 +379,7 @@ def __init__(self, dataBlock, key, indices=None):
         self.data = list(self._dict['data'].values())
         self._n_fields = len(self.fields)
         self._n_rows = len(self.data)
-        self._title = dataBlock._title + ' loop ' + str(key)
+        self._title = str(dataBlock._title) + ' loop ' + str(key)
 
     def numRows(self):
         return self._n_rows

From 9fd48cfaecd08e05fe3e8f7890a7f9f396a997ea Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 11 May 2022 11:04:19 +0200
Subject: [PATCH 45/73] no svd in pca app to allow nmodes control

---
 prody/apps/prody_apps/prody_pca.py | 35 ++++++++++--------------------
 1 file changed, 12 insertions(+), 23 deletions(-)

diff --git a/prody/apps/prody_apps/prody_pca.py b/prody/apps/prody_apps/prody_pca.py
index 1eb89bbf1..dd3675e04 100644
--- a/prody/apps/prody_apps/prody_pca.py
+++ b/prody/apps/prody_apps/prody_pca.py
@@ -96,27 +96,13 @@ def prody_pca(coords, **kwargs):
                 raise ImportError('Please install threadpoolctl to control threads')
 
             with threadpool_limits(limits=nproc, user_api="blas"):
-                if len(dcd) > 1000:
-                    pca.buildCovariance(dcd, aligned=kwargs.get('aligned'), quiet=quiet)
-                    pca.calcModes(nmodes)
-                    ensemble = dcd
-                else:
-                    ensemble = dcd[:]
-                    if not kwargs.get('aligned'):
-                        ensemble.iterpose(quiet=quiet)
-                    pca.performSVD(ensemble)
-                nmodes = pca.numModes()
-        else:
-            if len(dcd) > 1000:
                 pca.buildCovariance(dcd, aligned=kwargs.get('aligned'), quiet=quiet)
                 pca.calcModes(nmodes)
                 ensemble = dcd
-            else:
-                ensemble = dcd[:]
-                if not kwargs.get('aligned'):
-                    ensemble.iterpose(quiet=quiet)
-                pca.performSVD(ensemble)
-            nmodes = pca.numModes()
+        else:
+            pca.buildCovariance(dcd, aligned=kwargs.get('aligned'), quiet=quiet)
+            pca.calcModes(nmodes)
+            ensemble = dcd
 
     else:
         pdb = prody.parsePDB(coords)
@@ -147,10 +133,13 @@ def prody_pca(coords, **kwargs):
                 raise ImportError('Please install threadpoolctl to control threads')
 
             with threadpool_limits(limits=nproc, user_api="blas"):
-                pca.performSVD(ensemble)
+                pca.buildCovariance(ensemble, aligned=kwargs.get('aligned'), quiet=quiet)
+                pca.calcModes(nmodes)
         else:
-            pca.performSVD(ensemble)
+            pca.buildCovariance(ensemble, aligned=kwargs.get('aligned'), quiet=quiet)
+            pca.calcModes(nmodes)
 
+    pca = pca[:nmodes]
 
     LOGGER.info('Writing numerical output.')
     if kwargs.get('outnpz'):
@@ -159,15 +148,15 @@ def prody_pca(coords, **kwargs):
     if kwargs.get('outscipion'):
         prody.writeScipionModes(outdir, pca)
 
-    prody.writeNMD(join(outdir, prefix + '.nmd'), pca[:nmodes], select)
+    prody.writeNMD(join(outdir, prefix + '.nmd'), pca, select)
 
     extend = kwargs.get('extend')
     if extend:
         if pdb:
             if extend == 'all':
-                extended = prody.extendModel(pca[:nmodes], select, pdb)
+                extended = prody.extendModel(pca, select, pdb)
             else:
-                extended = prody.extendModel(pca[:nmodes], select,
+                extended = prody.extendModel(pca, select,
                                              select | pdb.bb)
             prody.writeNMD(join(outdir, prefix + '_extended_' +
                            extend + '.nmd'), *extended)

From f3464503b1645bfce84e9fe54363603e66fde420 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 11 May 2022 11:14:21 +0200
Subject: [PATCH 46/73] -

---
 prody/apps/prody_apps/prody_pca.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/prody/apps/prody_apps/prody_pca.py b/prody/apps/prody_apps/prody_pca.py
index dd3675e04..0e79c91fb 100644
--- a/prody/apps/prody_apps/prody_pca.py
+++ b/prody/apps/prody_apps/prody_pca.py
@@ -139,8 +139,6 @@ def prody_pca(coords, **kwargs):
             pca.buildCovariance(ensemble, aligned=kwargs.get('aligned'), quiet=quiet)
             pca.calcModes(nmodes)
 
-    pca = pca[:nmodes]
-
     LOGGER.info('Writing numerical output.')
     if kwargs.get('outnpz'):
         prody.saveModel(pca, join(outdir, prefix))

From ec52f9b4f8409faba86f93c950968a949bba0509 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Fri, 13 May 2022 17:36:47 +0200
Subject: [PATCH 47/73] cif header fix for 3p3w

---
 prody/proteins/cifheader.py | 152 ++++++++++++------------------------
 prody/proteins/compare.py   |   2 +-
 2 files changed, 50 insertions(+), 104 deletions(-)

diff --git a/prody/proteins/cifheader.py b/prody/proteins/cifheader.py
index 8cc7b9574..ae2097113 100644
--- a/prody/proteins/cifheader.py
+++ b/prody/proteins/cifheader.py
@@ -1005,123 +1005,69 @@ def _getPolymers(lines):
 
     # COMPND double block. 
     # Block 6 has most info. Block 7 has synonyms
-    i = 0
-    fields6 = OrderedDict()
-    fieldCounter6 = -1
-    foundPolyBlock6 = False
-    foundPolyBlockData6 = False
-    donePolyBlock6 = False
-    foundPolyBlock7 = False
-    donePolyBlock7 = False
-    fields7 = OrderedDict()
-    fieldCounter7 = -1
-    start6 = 0
-    stop6 = 0
-    stop7 = 0
-    while not donePolyBlock7 and i < len(lines):
-        line = lines[i]
-        if line.split('.')[0] == '_entity':
-            fieldCounter6 += 1
-            fields6[line.split('.')[1].strip()] = fieldCounter6
-            if not foundPolyBlock6:
-                start6 = i
-                foundPolyBlock6 = True
-
-        if line.split('.')[0] == '_entity_name_com':
-            fieldCounter7 += 1
-            fields7[line.split('.')[1].strip()] = fieldCounter7
-            if not foundPolyBlock7:
-                foundPolyBlock7 = True
-
-        if foundPolyBlock6:
-            if not line.startswith('#'):
-                if not foundPolyBlockData6:
-                    foundPolyBlockData6 = True
-            else:
-                if foundPolyBlock6 and not foundPolyBlock7:
-                    donePolyBlock6 = True
-                    stop6 = i
+    data6 = parseSTARSection(lines, "_entity")
+    data7 = parseSTARSection(lines, "_entity_name_com")
 
-                if foundPolyBlock7:
-                    donePolyBlock7 = True
-                    stop7 = i
+    dict_ = {}
+    for molecule in data6:
+        dict_.clear()
+        for k, value in molecule.items():
+            if k == '_entity.id':
+                dict_['CHAIN'] = ', '.join(entities[value])
 
-        i += 1
+            try:
+                key = _COMPND_KEY_MAPPINGS[k]
+            except:
+                continue
+            val = value.strip()
+            if val == '?':
+                val = ''
+            dict_[key.strip()] = val
 
-    if i < len(lines):
-        star_dict6, _ = parseSTARLines(lines[:2] + lines[start6-1:stop6], shlex=True)
-        loop_dict6 = list(star_dict6.values())[0]
+        chains = dict_.pop('CHAIN', '').strip()
 
-        data6 = list(loop_dict6[0]["data"].values())
+        if not chains:
+            continue
+        for ch in chains.split(','):
+            ch = ch.strip()
+            poly = polymers.get(ch, Polymer(ch))
+            polymers[ch] = poly
+            poly.name = dict_.get('MOLECULE', '').upper()
 
-        star_dict7, _ = parseSTARLines(lines[:2] + lines[stop6:stop7], shlex=True)
-        loop_dict7 = list(star_dict7.values())[0]
+            poly.fragment = dict_.get('FRAGMENT', '').upper()
 
-        if lines[stop6+1].strip() == "loop_":
-            data7 = loop_dict7[0]["data"].values()
-        else:
-            data7 = [loop_dict7["data"]]
+            poly.comments = dict_.get('OTHER_DETAILS', '').upper()
 
-        dict_ = {}
-        for molecule in data6:
-            dict_.clear()
-            for k, value in molecule.items():
-                if k == '_entity.id':
-                    dict_['CHAIN'] = ', '.join(entities[value])
+            val = dict_.get('EC', '')
+            poly.ec = [s.strip() for s in val.split(',')] if val else []
 
-                try:
-                    key = _COMPND_KEY_MAPPINGS[k]
-                except:
-                    continue
-                val = value.strip()
-                if val == '?':
-                    val = ''
-                dict_[key.strip()] = val
+            poly.mutation = dict_.get('MUTATION', '') != ''
+            poly.engineered = dict_.get('ENGINEERED', poly.mutation)
 
-            chains = dict_.pop('CHAIN', '').strip()
+    for molecule in data7:
+        dict_.clear()
+        for k, value in molecule.items():
+            if k.find('entity_id') != -1:
+                dict_['CHAIN'] = ', '.join(entities[value])
 
-            if not chains:
+            try:
+                key = _COMPND_KEY_MAPPINGS[k]
+            except:
                 continue
-            for ch in chains.split(','):
-                ch = ch.strip()
-                poly = polymers.get(ch, Polymer(ch))
-                polymers[ch] = poly
-                poly.name = dict_.get('MOLECULE', '').upper()
-
-                poly.fragment = dict_.get('FRAGMENT', '').upper()
+            dict_[key.strip()] = value.strip()
 
-                poly.comments = dict_.get('OTHER_DETAILS', '').upper()
+        chains = dict_.pop('CHAIN', '').strip()
 
-                val = dict_.get('EC', '')
-                poly.ec = [s.strip() for s in val.split(',')] if val else []
-
-                poly.mutation = dict_.get('MUTATION', '') != ''
-                poly.engineered = dict_.get('ENGINEERED', poly.mutation)
-
-        for molecule in data7:
-            dict_.clear()
-            for k, value in molecule.items():
-                if k.find('entity_id') != -1:
-                    dict_['CHAIN'] = ', '.join(entities[value])
-
-                try:
-                    key = _COMPND_KEY_MAPPINGS[k]
-                except:
-                    continue
-                dict_[key.strip()] = value.strip()
-
-            chains = dict_.pop('CHAIN', '').strip()
+        if not chains:
+            continue
+        for ch in chains.split(','):
+            ch = ch.strip()
+            poly = polymers.get(ch, Polymer(ch))
+            polymers[ch] = poly
 
-            if not chains:
-                continue
-            for ch in chains.split(','):
-                ch = ch.strip()
-                poly = polymers.get(ch, Polymer(ch))
-                polymers[ch] = poly
-
-                val = dict_.get('SYNONYM', '')
-                poly.synonyms = [s.strip().upper() for s in val.split(',')
-                                    ] if val else []
+            val = dict_.get('SYNONYM', '')
+            poly.synonyms = [s.strip().upper() for s in val.split(',')
+                                ] if val else []
 
     return list(polymers.values())
 
diff --git a/prody/proteins/compare.py b/prody/proteins/compare.py
index e0e2c6feb..58ab1660c 100644
--- a/prody/proteins/compare.py
+++ b/prody/proteins/compare.py
@@ -913,7 +913,7 @@ def mapChainOntoChain(mobile, target, **kwargs):
     :arg target: chain to which atoms will be mapped
     :type target: :class:`.Chain`
 
-    :keyword seqid: percent sequence identity, default is **90**. Note that This parameter is 
+    :keyword seqid: percent sequence identity, default is **90**. Note that this parameter is 
         only effective for sequence alignment
     :type seqid: float
 

From 711fb78761936dfb4d35a2b32e13a10dbb41ac37 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 17 May 2022 08:59:53 +0200
Subject: [PATCH 48/73] add missed quiet check

---
 prody/dynamics/pca.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/prody/dynamics/pca.py b/prody/dynamics/pca.py
index 65a936cd9..70ffe0ccc 100644
--- a/prody/dynamics/pca.py
+++ b/prody/dynamics/pca.py
@@ -158,13 +158,16 @@ def buildCovariance(self, coordsets, **kwargs):
                     cov = np.zeros((dof, dof))
                     coordsets = coordsets.reshape((n_confs, dof))
                     mean = coordsets.mean(0)
-                    LOGGER.progress('Building covariance', n_confs,
+                    if not quiet:
+                        LOGGER.progress('Building covariance', n_confs,
                                     '_prody_pca')
                     for i, coords in enumerate(coordsets.reshape(s)):
                         deviations = coords - mean
                         cov += np.outer(deviations, deviations)
-                        LOGGER.update(n_confs, label='_prody_pca')
-                    LOGGER.finish()
+                        if not quiet:
+                            LOGGER.update(n_confs, label='_prody_pca')
+                    if not quiet:
+                        LOGGER.finish()
                     cov /= n_confs
                     self._cov = cov
             else:
@@ -255,10 +258,10 @@ def performSVD(self, coordsets):
                               coordsets._getCoords())
 
         n_confs = deviations.shape[0]
-        if n_confs < 3:
+        if n_confs <= 3:
             raise ValueError('coordsets must have more than 3 coordinate sets')
         n_atoms = deviations.shape[1]
-        if n_atoms < 3:
+        if n_atoms <= 3:
             raise ValueError('coordsets must have more than 3 atoms')
 
         dof = n_atoms * 3

From ba968714112e1273bf9cea8e19fc08c3f4c82ae2 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 18 May 2022 08:34:50 +0200
Subject: [PATCH 49/73] parseScipionModes from metadata file

---
 prody/dynamics/functions.py | 79 ++++++++++++++++++-------------------
 1 file changed, 39 insertions(+), 40 deletions(-)

diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 73cf470ac..92f69266a 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -315,7 +315,7 @@ def parseModes(normalmodes, eigenvalues=None, nm_delimiter=None,
     return nma
 
 
-def parseScipionModes(run_path, title=None, pdb=None):
+def parseScipionModes(metadata_file, title=None, pdb=None):
     """Returns :class:`.NMA` containing eigenvectors and eigenvalues 
     parsed from a ContinuousFlex FlexProtNMA Run directory.
 
@@ -325,59 +325,58 @@ def parseScipionModes(run_path, title=None, pdb=None):
     :arg title: title for :class:`.NMA` object
     :type title: str
     """
-    if run_path.endswith("/"):
-        run_path = run_path[:-1]
+    run_path = os.path.split(metadata_file)[0]
+    top_dirs = os.path.split(run_path)[0][:-4]
     run_name = os.path.split(run_path)[-1]
-    top_dirs = os.path.split(run_path)[0][:-4] # exclude "Runs"
 
-    try:
-        star_data = parseSTAR(run_path + '/modes.xmd')
-    except OSError:
-        try:
-            sql_con = openSQLite(run_path + "/modes.sqlite")
-            cursor = sql_con.cursor()
-
-            #names = ['id', 'enabled', 'label', 'comment', 'creation', 'c01', 'c02', 'c03', 'c04']
-
-            star_dict = OrderedDict()
-            star_block_dict = OrderedDict()
-            
-            star_loop_dict = OrderedDict()
-            star_loop_dict["fields"] = OrderedDict([(0, '_enabled'),
-                                               (1, '_nmaCollectivity'),
-                                               (2, '_nmaModefile'),
-                                               (3, '_nmaScore'),
-                                               (4, '_nmaEigenval'),
-                                               (5, '_order_')])
-            star_loop_dict["data"] = OrderedDict()
+    if metadata_file.endswith('.xmd'):
+        star_data = parseSTAR(metadata_file)
+        
+    elif metadata_file.endswith('.sqlite'):
+        # reconstruct star data from sqlite
 
-            for row in cursor.execute("SELECT * FROM Objects;"):
+        sql_con = openSQLite(metadata_file)
+        cursor = sql_con.cursor()
 
-                id_ = row[0]
-                key = id_-1 # sqlite ids start from 1 not 0
+        star_dict = OrderedDict()
+        star_block_dict = OrderedDict()
+        
+        star_loop_dict = OrderedDict()
+        star_loop_dict["fields"] = OrderedDict([(0, '_enabled'),
+                                            (1, '_nmaCollectivity'),
+                                            (2, '_nmaModefile'),
+                                            (3, '_nmaScore'),
+                                            (4, '_nmaEigenval'),
+                                            (5, '_order_')])
+        star_loop_dict["data"] = OrderedDict()
+
+        for row in cursor.execute("SELECT * FROM Objects;"):
+        
+            id_ = row[0]
+            key = id_ - 1 # sqlite ids start from 1 not 0
 
-                star_loop_dict["data"][key] = OrderedDict()
+            star_loop_dict["data"][key] = OrderedDict()
 
-                star_loop_dict["data"][key]['_order_'] = id_
+            star_loop_dict["data"][key]['_order_'] = id_
 
-                star_loop_dict["data"][key]['_enabled'] = row[1]
+            star_loop_dict["data"][key]['_enabled'] = row[1]
 
-                star_loop_dict["data"][key]['_nmaCollectivity'] = row[6]
+            star_loop_dict["data"][key]['_nmaCollectivity'] = row[6]
 
-                star_loop_dict["data"][key]['_nmaModefile'] = row[5]
+            star_loop_dict["data"][key]['_nmaModefile'] = row[5]
 
-                star_loop_dict["data"][key]['_nmaScore'] = row[7]
+            star_loop_dict["data"][key]['_nmaScore'] = row[7]
 
-                if len(row) > 8:
-                    star_loop_dict["data"][key]['_nmaEigenval'] = row[8]
+            if len(row) > 8:
+                star_loop_dict["data"][key]['_nmaEigenval'] = row[8]
 
-            star_block_dict[0] = star_loop_dict
-            star_dict[0] = star_block_dict
+        star_block_dict[0] = star_loop_dict
+        star_dict[0] = star_block_dict
 
-            star_data = StarDict(star_dict, prog='XMIPP')
+        star_data = StarDict(star_dict, prog='XMIPP')
 
-        except OSError:
-            raise OSError("neither modes.xmd nor modes.sqlite found")
+    else:
+        raise ValueError("Metadata file should be an xmd or sqlite file")
 
     star_loop = star_data[0][0]
     

From 0d07ef5325c809412ac66282443cc7d9dee8273b Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 19 May 2022 08:25:57 +0200
Subject: [PATCH 50/73] add back pca nmd slicing

---
 prody/apps/prody_apps/prody_pca.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/prody/apps/prody_apps/prody_pca.py b/prody/apps/prody_apps/prody_pca.py
index 49af7bad9..0c96b58d5 100644
--- a/prody/apps/prody_apps/prody_pca.py
+++ b/prody/apps/prody_apps/prody_pca.py
@@ -150,15 +150,15 @@ def prody_pca(coords, **kwargs):
     if kwargs.get('outscipion'):
         prody.writeScipionModes(outdir, pca)
 
-    prody.writeNMD(join(outdir, prefix + '.nmd'), pca, select)
+    prody.writeNMD(join(outdir, prefix + '.nmd'), pca[:nmodes], select)
 
     extend = kwargs.get('extend')
     if extend:
         if pdb:
             if extend == 'all':
-                extended = prody.extendModel(pca, select, pdb)
+                extended = prody.extendModel(pca[:nmodes], select, pdb)
             else:
-                extended = prody.extendModel(pca, select,
+                extended = prody.extendModel(pca[:nmodes], select,
                                              select | pdb.bb)
             prody.writeNMD(join(outdir, prefix + '_extended_' +
                            extend + '.nmd'), *extended)

From efd19651339f76625e28e5e73ac0481ff5f8aa2e Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 26 May 2022 11:39:18 +0200
Subject: [PATCH 51/73] interp takes coords

---
 prody/dynamics/editing.py | 46 ++++++++++++++++++++++++++++-----------
 1 file changed, 33 insertions(+), 13 deletions(-)

diff --git a/prody/dynamics/editing.py b/prody/dynamics/editing.py
index 230f083bd..bfc51aaea 100644
--- a/prody/dynamics/editing.py
+++ b/prody/dynamics/editing.py
@@ -6,7 +6,7 @@
 from prody.atomic import Atomic, AtomGroup, AtomMap, AtomSubset
 from prody.atomic import Selection, SELECT, sliceAtoms, extendAtoms
 from prody.measure import calcDistance
-from prody.utilities import importLA, isListLike
+from prody.utilities import importLA, isListLike, getCoords
 from prody import _PY3K, LOGGER
 
 from .nma import NMA
@@ -268,10 +268,10 @@ def sliceMode(mode, atoms, select):
     return (vec, sel)
 
 
-def sliceModel(model, atoms, select):
+def sliceModel(model, atoms, select, norm=False):
     """Returns a part of the *model* (modes calculated) for *atoms* matching *select*. 
     Note that normal modes are sliced instead the connectivity matrix. Sliced normal 
-    modes (eigenvectors) are not normalized.
+    modes (eigenvectors) are not normalized unless *norm* is **True**.
 
     :arg mode: NMA model instance to be sliced
     :type mode: :class:`.NMA`
@@ -282,6 +282,9 @@ def sliceModel(model, atoms, select):
     :arg select: an atom selection or a selection string
     :type select: :class:`.Selection`, str
 
+    :arg norm: whether to normalize eigenvectors, default **False**
+    :type norm: bool
+
     :returns: (:class:`.NMA`, :class:`.Selection`)"""
 
     if not isinstance(model, NMA):
@@ -294,13 +297,13 @@ def sliceModel(model, atoms, select):
         raise ValueError('number of atoms in model and atoms must be equal')
 
     which, sel = sliceAtoms(atoms, select)
-    nma = sliceModelByMask(model, which)
+    nma = sliceModelByMask(model, which, norm=norm)
 
     return (nma, sel)
 
-def sliceModelByMask(model, mask):
+def sliceModelByMask(model, mask, norm=False):
     """Returns a part of the *model* indicated by *mask*.  Note that
-    normal modes (eigenvectors) are not normalized.
+    normal modes (eigenvectors) are not normalized unless *norm* is **True**.
 
     :arg mode: NMA model instance to be sliced
     :type mode: :class:`.NMA`
@@ -309,6 +312,9 @@ def sliceModelByMask(model, mask):
         the parts being selected 
     :type mask: list, :class:`~numpy.ndarray`
 
+    :arg norm: whether to normalize eigenvectors, default **False**
+    :type norm: bool
+
     :returns: :class:`.NMA`"""
 
     if not isListLike(mask):
@@ -333,7 +339,13 @@ def sliceModelByMask(model, mask):
                     .format(model.getTitle()))
     if model.is3d():
         which = np.repeat(which, 3)
-    nma.setEigens(array[which, :], model.getEigvals())
+
+    evecs = array[which, :]
+    if norm:
+        evecs /= np.array([((evecs[:, i]) ** 2).sum() ** 0.5
+                           for i in range(evecs.shape[1])])
+
+    nma.setEigens(evecs, model.getEigvals())
     return nma
 
 def reduceModel(model, atoms, select):
@@ -477,11 +489,19 @@ def _reduceModel(matrix, system):
     return matrix
 
 
-def interpolateModel(model, nodes, atoms, norm=False, **kwargs):
-    """Interpolate a coarse grained *model* built for *nodes* to *atoms*.  
+def interpolateModel(model, nodes, coords, norm=False, **kwargs):
+    """Interpolate a coarse grained *model* built for *nodes* to *coords*.  
     
     *model* may be :class:`.ANM`, :class:`.PCA`, or :class:`.NMA`
-    instance. *model* should have all modes calculated.
+    instance
+
+    :arg nodes: the coordinate set or object with :meth:`getCoords` method
+        that corresponds to the model
+    :type nodes: :class:`.Atomic`, :class:`~numpy.ndarray`
+
+    :arg coords: a coordinate set or an object with :meth:`getCoords` method
+        onto which the model should be interpolated
+    :type coords: :class:`.Atomic`, :class:`~numpy.ndarray`
     
     This function will take the part of the normal modes for each node
     (i.e. Cα atoms) and extend it to nearby atoms. 
@@ -535,8 +555,8 @@ def interpolateModel(model, nodes, atoms, norm=False, **kwargs):
     cg_coords = nodes.getCoords()
     len_cg = nodes.numAtoms()
 
-    fg_coords = atoms.getCoords()
-    len_fg = atoms.numAtoms()
+    fg_coords = getCoords(coords)
+    len_fg = fg_coords.shape[0]//3
 
     n_modes = model.numModes()
 
@@ -614,4 +634,4 @@ def interpolateModel(model, nodes, atoms, norm=False, **kwargs):
             
     interpolated = NMA('Interpolated ' + str(model))
     interpolated.setEigens(eigvecs_fg, eigvals)
-    return interpolated, atoms
+    return interpolated, coords

From 77ffdeb707351f927fea69c5be1d27ef3e37d201 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 7 Aug 2022 16:10:55 +0200
Subject: [PATCH 52/73] support 3-char chids

---
 prody/atomic/fields.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/atomic/fields.py b/prody/atomic/fields.py
index 2451765d4..342010b3e 100644
--- a/prody/atomic/fields.py
+++ b/prody/atomic/fields.py
@@ -125,7 +125,7 @@ def getDocstr(self, meth, plural=True, selex=True):
                        selstr=('altloc A B', 'altloc _'),),
     'anisou':    Field('anisou', float, doc='anisotropic temperature factor',
                        ndim=2),
-    'chain':     Field('chain', DTYPE + '2',  doc='chain identifier',
+    'chain':     Field('chain', DTYPE + '3',  doc='chain identifier',
                        meth='Chid', none=HVNONE, synonym='chid',
                        selstr=('chain A', 'chid A B C', 'chain _')),
     'element':   Field('element', DTYPE + '2', doc='element symbol',

From f7511cc07e4b38b039873a0cdbb7c4d604ce9193 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 7 Aug 2022 16:11:39 +0200
Subject: [PATCH 53/73] cif secstr fix

---
 prody/proteins/cifheader.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/prody/proteins/cifheader.py b/prody/proteins/cifheader.py
index ae2097113..09c5e05e2 100644
--- a/prody/proteins/cifheader.py
+++ b/prody/proteins/cifheader.py
@@ -285,8 +285,8 @@ def _getHelix(lines):
             data = line.split()
 
             try:
-                chid = data[fields["beg_label_asym_id"]]
-                segn = data[fields["beg_auth_asym_id"]]
+                chid = data[fields["beg_auth_asym_id"]]
+                segn = data[fields["beg_label_asym_id"]]
                 value = (int(data[fields["pdbx_PDB_helix_class"]]),
                         int(data[fields["id"]][6:]),
                         data[fields["pdbx_PDB_helix_id"]].strip())
@@ -352,7 +352,7 @@ def _getHelixRange(lines):
             data = line.split()
 
             try:
-                chid = data[fields["beg_label_asym_id"]]
+                chid = data[fields["beg_auth_asym_id"]]
                 Hclass=int(data[fields["pdbx_PDB_helix_class"]])
                 Hnr=int(data[fields["id"]][6:])
             except:
@@ -504,8 +504,8 @@ def _getSheet(lines):
             data = line.split()
 
             try:
-                chid = data[fields3["beg_label_asym_id"]]
-                segn = data[fields3["beg_auth_asym_id"]]
+                chid = data[fields3["beg_auth_asym_id"]]
+                segn = data[fields3["beg_label_asym_id"]]
 
                 strand_id = int(data[fields3["id"]])
                 sheet_id = data[fields3["sheet_id"]]
@@ -675,7 +675,7 @@ def _getSheetRange(lines):
             data = line.split()
 
             try:
-                chid = data[fields3["beg_label_asym_id"]]
+                chid = data[fields3["beg_auth_asym_id"]]
 
                 Snr = int(data[fields3["id"]])
                 sheet_id = data[fields3["sheet_id"]]

From c30b27065f86f3daf82a6d9aee857461f49c7087 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 7 Aug 2022 16:14:34 +0200
Subject: [PATCH 54/73] poly strand id to star section

---
 prody/proteins/cifheader.py | 137 ++++++++++++++----------------------
 1 file changed, 54 insertions(+), 83 deletions(-)

diff --git a/prody/proteins/cifheader.py b/prody/proteins/cifheader.py
index 09c5e05e2..87c72b82f 100644
--- a/prody/proteins/cifheader.py
+++ b/prody/proteins/cifheader.py
@@ -916,92 +916,63 @@ def _getPolymers(lines):
                                   data[fields4["details"]]))
 
     # SEQADV block
-    i = 0
-    fields5 = OrderedDict()
-    fieldCounter5 = -1
-    foundPolyBlock5 = False
-    foundPolyBlockData5 = False
-    donePolyBlock5 = False
-    start5 = 0
-    stop5 = 0
-    while not donePolyBlock5 and i < len(lines):
-        line = lines[i]
-        if line.split('.')[0] == '_struct_ref_seq_dif':
-            fieldCounter5 += 1
-            fields5[line.split('.')[1].strip()] = fieldCounter5
-            if not foundPolyBlock5:
-                foundPolyBlock5 = True
-
-        if foundPolyBlock5:
-            if not line.startswith('#') and not line.startswith('_'):
-                if not foundPolyBlockData5:
-                    start5 = i
-                    foundPolyBlockData5 = True
-            else:
-                if foundPolyBlockData5:
-                    donePolyBlock5 = True
-                    stop5 = i
-
-        i += 1
-
-    if i < len(lines):
-        for i, line in enumerate(lines[start5:stop5]):
-            data = split(line, shlex=True)
-
-            ch = data[fields5["pdbx_pdb_strand_id"]]
-
-            poly = polymers.get(ch, Polymer(ch))
-            polymers[ch] = poly
-            dbabbr = data[fields5["pdbx_seq_db_name"]]
-            resname = data[fields5["mon_id"]]
-            if resname == '?':
-                resname = '' # strip for pdb
+    data5 = parseSTARSection(lines, "_struct_ref_seq_dif")
+
+    for i, data in enumerate(data5):
+        ch = data["_struct_ref_seq_dif.pdbx_pdb_strand_id"]
+
+        poly = polymers.get(ch, Polymer(ch))
+        polymers[ch] = poly
+        dbabbr = data["_struct_ref_seq_dif.pdbx_seq_db_name"]
+        resname = data["_struct_ref_seq_dif.mon_id"]
+        if resname == '?':
+            resname = '' # strip for pdb
+
+        try:
+            resnum = int(data["_struct_ref_seq_dif.pdbx_auth_seq_num"])
+        except:
+            #LOGGER.warn('SEQADV for chain {2}: failed to parse PDB sequence '
+            #            'number ({0}:{1})'.format(pdbid, i, ch))
+            continue
 
-            try:
-                resnum = int(data[fields5["pdbx_auth_seq_num"]])
-            except:
-                #LOGGER.warn('SEQADV for chain {2}: failed to parse PDB sequence '
-                #            'number ({0}:{1})'.format(pdbid, i, ch))
+        icode = data["_struct_ref_seq_dif.pdbx_pdb_ins_code"]
+        if icode == '?':
+            icode = '' # strip for pdb            
+        
+        try:
+            dbnum = int(data["_struct_ref_seq_dif.pdbx_seq_db_seq_num"])
+        except:
+            #LOGGER.warn('SEQADV for chain {2}: failed to parse database '
+            #            'sequence number ({0}:{1})'.format(pdbid, i, ch))
+            continue            
+
+        comment = data["_struct_ref_seq_dif.details"].upper()
+        if comment == '?':
+            comment = '' # strip for pdb 
+
+        match = False
+        for dbref in poly.dbrefs:
+            if not dbref.first[0] <= resnum <= dbref.last[0]:
                 continue
-
-            icode = data[fields5["pdbx_pdb_ins_code"]]
-            if icode == '?':
-                icode = '' # strip for pdb            
-            
-            try:
-                dbnum = int(data[fields5["pdbx_seq_db_seq_num"]])
-            except:
-                #LOGGER.warn('SEQADV for chain {2}: failed to parse database '
-                #            'sequence number ({0}:{1})'.format(pdbid, i, ch))
-                continue            
-
-            comment = data[fields5["details"]].upper()
-            if comment == '?':
-                comment = '' # strip for pdb 
-
-            match = False
-            for dbref in poly.dbrefs:
-                if not dbref.first[0] <= resnum <= dbref.last[0]:
-                    continue
-                match = True
-                if dbref.dbabbr != dbabbr:
-                    LOGGER.warn('SEQADV for chain {2}: reference database '
-                                'mismatch, expected {3} parsed {4} '
-                                '({0}:{1})'.format(pdbid, i+1, ch,
-                                repr(dbref.dbabbr), repr(dbabbr)))
-                    continue
-                dbacc = data[fields5["pdbx_seq_db_accession_code"]]
-                if dbref.accession[:9] != dbacc[:9]:
-                    LOGGER.warn('SEQADV for chain {2}: accession code '
-                                'mismatch, expected {3} parsed {4} '
-                                '({0}:{1})'.format(pdbid, i+1, ch,
-                                repr(dbref.accession), repr(dbacc)))
-                    continue
-                dbref.diff.append((resname, resnum, icode, dbnum, dbnum, comment))
-            if not match:
-                LOGGER.warn('SEQADV for chain {2}: database sequence reference '
-                            'not found ({0}:{1})'.format(pdbid, i+1, ch))
+            match = True
+            if dbref.dbabbr != dbabbr:
+                LOGGER.warn('SEQADV for chain {2}: reference database '
+                            'mismatch, expected {3} parsed {4} '
+                            '({0}:{1})'.format(pdbid, i+1, ch,
+                            repr(dbref.dbabbr), repr(dbabbr)))
+                continue
+            dbacc = data["_struct_ref_seq_dif.pdbx_seq_db_accession_code"]
+            if dbref.accession[:9] != dbacc[:9]:
+                LOGGER.warn('SEQADV for chain {2}: accession code '
+                            'mismatch, expected {3} parsed {4} '
+                            '({0}:{1})'.format(pdbid, i+1, ch,
+                            repr(dbref.accession), repr(dbacc)))
                 continue
+            dbref.diff.append((resname, resnum, icode, dbnum, dbnum, comment))
+        if not match:
+            LOGGER.warn('SEQADV for chain {2}: database sequence reference '
+                        'not found ({0}:{1})'.format(pdbid, i+1, ch))
+            continue
 
     # COMPND double block. 
     # Block 6 has most info. Block 7 has synonyms

From 767704a96bd6fb1efe856c243923586979a1780f Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 7 Aug 2022 16:34:13 +0200
Subject: [PATCH 55/73] modres parseSTARSection too

---
 prody/proteins/cifheader.py | 61 ++++++++++---------------------------
 1 file changed, 16 insertions(+), 45 deletions(-)

diff --git a/prody/proteins/cifheader.py b/prody/proteins/cifheader.py
index 87c72b82f..ce542f49d 100644
--- a/prody/proteins/cifheader.py
+++ b/prody/proteins/cifheader.py
@@ -867,53 +867,24 @@ def _getPolymers(lines):
             last = temp
 
     # MODRES block
-    i = 0
-    fields4 = OrderedDict()
-    fieldCounter4 = -1
-    foundPolyBlock4 = False
-    foundPolyBlockData4 = False
-    donePolyBlock4 = False
-    start4 = 0
-    stop4 = 0
-    while not donePolyBlock4 and i < len(lines):
-        line = lines[i]
-        if line.split('.')[0] == '_pdbx_struct_mod_residue':
-            fieldCounter4 += 1
-            fields4[line.split('.')[1].strip()] = fieldCounter4
-            if not foundPolyBlock4:
-                foundPolyBlock4 = True
-
-        if foundPolyBlock4:
-            if not line.startswith('#') and not line.startswith('_'):
-                if not foundPolyBlockData4:
-                    start4 = i
-                    foundPolyBlockData4 = True
-            else:
-                if foundPolyBlockData4:
-                    donePolyBlock4 = True
-                    stop4 = i
+    data4 = parseSTARSection(lines, "_pdbx_struct_mod_residue")
 
-        i += 1
-
-    if i < len(lines):
-        for line in lines[start4:stop4]:
-            data = split(line, shlex=True)
+    for data in data4:
+        ch = data["_pdbx_struct_mod_residue.label_asym_id"]
 
-            ch = data[fields4["label_asym_id"]]
-
-            poly = polymers.get(ch, Polymer(ch))
-            polymers[ch] = poly
-            if poly.modified is None:
-                poly.modified = []
-
-            iCode = data[fields4["PDB_ins_code"]]
-            if iCode == '?':
-                iCode == '' # PDB one is stripped
-            poly.modified.append((data[fields4["auth_comp_id"]],
-                                  data[fields4["auth_asym_id"]],
-                                  data[fields4["auth_seq_id"]] + iCode,
-                                  data[fields4["parent_comp_id"]],
-                                  data[fields4["details"]]))
+        poly = polymers.get(ch, Polymer(ch))
+        polymers[ch] = poly
+        if poly.modified is None:
+            poly.modified = []
+
+        iCode = data["_pdbx_struct_mod_residue.PDB_ins_code"]
+        if iCode == '?':
+            iCode == '' # PDB one is stripped
+        poly.modified.append((data["_pdbx_struct_mod_residue.auth_comp_id"],
+                                data["_pdbx_struct_mod_residue.auth_asym_id"],
+                                data["_pdbx_struct_mod_residue.auth_seq_id"] + iCode,
+                                data["_pdbx_struct_mod_residue.parent_comp_id"],
+                                data["_pdbx_struct_mod_residue.details"]))
 
     # SEQADV block
     data5 = parseSTARSection(lines, "_struct_ref_seq_dif")

From 6fafab2ddd40362528e7b70aa76e428d99738f8e Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Sun, 7 Aug 2022 17:18:34 +0200
Subject: [PATCH 56/73] fixed updateDefinitions

---
 prody/atomic/flags.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/atomic/flags.py b/prody/atomic/flags.py
index 3cb38a227..110b19706 100644
--- a/prody/atomic/flags.py
+++ b/prody/atomic/flags.py
@@ -690,7 +690,7 @@ def updateDefinitions():
 
     global DEFINITIONS, AMINOACIDS, BACKBONE, TIMESTAMP
     DEFINITIONS = {}
-    user = SETTINGS.get('flag_definitions', {})
+    user = SETTINGS.get('flags_definitions', {})
 
     # nucleics
     nucleic = set()

From 74b378b88977e54500b205eb204f043452ade81f Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Fri, 4 Nov 2022 09:23:18 +0100
Subject: [PATCH 57/73] fix

---
 prody/apps/prody_apps/prody_energy.py | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/prody/apps/prody_apps/prody_energy.py b/prody/apps/prody_apps/prody_energy.py
index 0d9ed4416..1b712ca1f 100644
--- a/prody/apps/prody_apps/prody_energy.py
+++ b/prody/apps/prody_apps/prody_energy.py
@@ -101,6 +101,7 @@ def prody_energy(*pdbs, **kwargs):
             
             if minimise:
                 from openmm.unit import kilojoule_per_mole, angstrom
+                from openmm.app.pdbfile import PDBFile
                 simulation.minimizeEnergy(tolerance=10.0 * kilojoule_per_mole, maxIterations=0)
                 
             state = simulation.context.getState(getEnergy=True, getPositions=True)
@@ -108,9 +109,9 @@ def prody_energy(*pdbs, **kwargs):
             
             if minimise:
                 pos = array(state.getPositions().in_units_of(angstrom)._value)
-                struct = clu.getAtoms().copy()
-                struct.setCoords(pos)
-                writePDB(ag.getTitle() + '_minim.pdb', struct)
+                fo = open(outname + '.pdb', 'w')
+                PDBFile.writeFile(simulation.topology, state.getPositions(), fo)
+                fo.close()
             
             f.write(str(energy) + " kJ/mol\n")
         
@@ -185,4 +186,8 @@ def addCommand(commands):
 
     group_energy.add_argument('-M', '--minimise', dest='minimise', action='store_true',
         default=False, help=('whether to energy minimise (default: %(default)s)'))
-    
\ No newline at end of file
+
+    group_energy.add_argument('-b', '--padding', dest='padding',
+                              type=float, metavar='FLOAT', default=1., 
+                              help=('padding distance in nanometers (default: %(default)s)'))
+

From 82a8837a08162aff1fc10a85e1a93cb60808f59d Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 7 Nov 2022 13:03:33 +0100
Subject: [PATCH 58/73] add boxSize

---
 prody/apps/prody_apps/prody_energy.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/prody/apps/prody_apps/prody_energy.py b/prody/apps/prody_apps/prody_energy.py
index 1b712ca1f..82fc24675 100644
--- a/prody/apps/prody_apps/prody_energy.py
+++ b/prody/apps/prody_apps/prody_energy.py
@@ -19,6 +19,9 @@ def prody_energy(*pdbs, **kwargs):
     
     :arg padding: padding in nanometers
     
+    :arg boxSize: boxSize array with x, y and z in nanometers
+        If not None then this overrides padding
+
     :arg force_field_protein: name of force field for protein in OpenMM
         default depends on solvent type as in ClustENM
     
@@ -51,6 +54,7 @@ def prody_energy(*pdbs, **kwargs):
         return
         
     padding = kwargs.get('padding', 1.0)
+    boxSize = kwargs.get('boxSize', "None")
     
     force_field_protein = kwargs.get('force_field_protein', None)
     force_field_solvent = kwargs.get('force_field_solvent', None)
@@ -96,6 +100,7 @@ def prody_energy(*pdbs, **kwargs):
                 clu._force_field = ('amber14-all.xml', 'amber14/tip3pfb.xml') if force_field is None else force_field
                 
             clu._padding = padding
+            clu._boxSize = eval(boxSize)
             
             simulation = clu._prep_sim(clu._atoms.getCoords())
             
@@ -191,3 +196,7 @@ def addCommand(commands):
                               type=float, metavar='FLOAT', default=1., 
                               help=('padding distance in nanometers (default: %(default)s)'))
 
+    group_energy.add_argument('-B', '--boxSize', dest='boxSize',
+                              type=str, metavar='STR', default="None",
+                              help=('x, y, z of box size in nanometers instead of padding (default: use padding)'))
+

From 99fdbcf8cf11b795929b70e02bfc5dd4141273cd Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 7 Nov 2022 13:24:31 +0100
Subject: [PATCH 59/73] add repeats

---
 prody/apps/prody_apps/prody_energy.py | 42 +++++++++++++++++----------
 1 file changed, 26 insertions(+), 16 deletions(-)

diff --git a/prody/apps/prody_apps/prody_energy.py b/prody/apps/prody_apps/prody_energy.py
index 82fc24675..136ef3f51 100644
--- a/prody/apps/prody_apps/prody_energy.py
+++ b/prody/apps/prody_apps/prody_energy.py
@@ -33,6 +33,8 @@ def prody_energy(*pdbs, **kwargs):
     :arg minimise: whether to energy minimise
     
     :arg select: atom selection string, default is "all"
+
+    :arg repeats: number of repeats
     """
 
     from os.path import isfile
@@ -63,6 +65,8 @@ def prody_energy(*pdbs, **kwargs):
     if model == 0:
         model = None
         
+    repeats = kwargs.get('repeats', 1)
+
     minimise = kwargs.get('minimise', False)
     selstr = kwargs.get('select', 'all')
     
@@ -102,23 +106,25 @@ def prody_energy(*pdbs, **kwargs):
             clu._padding = padding
             clu._boxSize = eval(boxSize)
             
-            simulation = clu._prep_sim(clu._atoms.getCoords())
-            
-            if minimise:
-                from openmm.unit import kilojoule_per_mole, angstrom
-                from openmm.app.pdbfile import PDBFile
-                simulation.minimizeEnergy(tolerance=10.0 * kilojoule_per_mole, maxIterations=0)
+            for j in range(repeats):
                 
-            state = simulation.context.getState(getEnergy=True, getPositions=True)
-            energy = state.getPotentialEnergy()._value
-            
-            if minimise:
-                pos = array(state.getPositions().in_units_of(angstrom)._value)
-                fo = open(outname + '.pdb', 'w')
-                PDBFile.writeFile(simulation.topology, state.getPositions(), fo)
-                fo.close()
-            
-            f.write(str(energy) + " kJ/mol\n")
+                simulation = clu._prep_sim(clu._atoms.getCoords())
+
+                if minimise:
+                    from openmm.unit import kilojoule_per_mole, angstrom
+                    from openmm.app.pdbfile import PDBFile
+                    simulation.minimizeEnergy(tolerance=10.0 * kilojoule_per_mole, maxIterations=0)
+
+                state = simulation.context.getState(getEnergy=True, getPositions=True)
+                energy = state.getPotentialEnergy()._value
+
+                if minimise:
+                    pos = array(state.getPositions().in_units_of(angstrom)._value)
+                    fo = open(outname + '.pdb', 'w')
+                    PDBFile.writeFile(simulation.topology, state.getPositions(), fo)
+                    fo.close()
+
+                f.write(str(energy) + " kJ/mol\n")
         
         f.close()
         LOGGER.info('\nEnergy is written into: ' + outname + '\n')
@@ -140,6 +146,10 @@ def addCommand(commands):
         default=0, metavar='INT',
         help=('index of model that will be used in the calculations (default: all of them)'))
 
+    subparser.add_argument('-r', '--repeats', dest='repeats', type=int,
+        default=1, metavar='INT',
+        help=('number of repeats (default: %(default)s)'))
+
     subparser.set_defaults(usage_example=
     """This command fixes missing atoms, solvates and writes energies in a txt file.
 

From d3a8ff540d7e2368e739691cefbaacf7e39e6995 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 7 Nov 2022 13:25:21 +0100
Subject: [PATCH 60/73] replace imp find prody

---
 docs/apps/docapps.py              | 3 ++-
 prody/apps/evol_apps/__init__.py  | 4 +++-
 prody/apps/prody_apps/__init__.py | 4 +++-
 scripts/evol                      | 4 +++-
 scripts/prody                     | 4 +++-
 5 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/docs/apps/docapps.py b/docs/apps/docapps.py
index 2dda9363d..07401800f 100755
--- a/docs/apps/docapps.py
+++ b/docs/apps/docapps.py
@@ -2,9 +2,10 @@
 
 import os
 import imp
+import importlib
 from subprocess import Popen, PIPE
 
-path = [imp.find_module('prody')[1]]
+path = [importlib.util.find_spec("prody").submodule_search_locations[0]]
 apps = imp.load_module('prody.apps',
                            *imp.find_module('apps', path))
 
diff --git a/prody/apps/evol_apps/__init__.py b/prody/apps/evol_apps/__init__.py
index 0d69ae4a7..32d521467 100644
--- a/prody/apps/evol_apps/__init__.py
+++ b/prody/apps/evol_apps/__init__.py
@@ -1,6 +1,7 @@
 """This module defines some sequence evolution applications."""
 
 import imp
+import importlib
 import sys
 
 try:
@@ -10,7 +11,8 @@
 
 from ..apptools import *
 
-path_prody = imp.find_module('prody')[1]
+#path_prody = imp.find_module('prody')[1]
+path_prody = importlib.util.find_spec("prody").submodule_search_locations[0]
 path_apps = imp.find_module('apps', [path_prody])[1]
 path_apps = imp.find_module('evol_apps', [path_apps])[1]
 
diff --git a/prody/apps/prody_apps/__init__.py b/prody/apps/prody_apps/__init__.py
index 510bc9f7b..fbfb5916c 100644
--- a/prody/apps/prody_apps/__init__.py
+++ b/prody/apps/prody_apps/__init__.py
@@ -1,6 +1,7 @@
 """This module defines structure and dynamics analysis applications."""
 
 import imp
+import importlib
 import sys
 
 try:
@@ -10,7 +11,8 @@
 
 from ..apptools import *
 
-path_prody = imp.find_module('prody')[1]
+#path_prody = imp.find_module('prody')[1]
+path_prody = importlib.util.find_spec("prody").submodule_search_locations[0]
 path_apps = imp.find_module('apps', [path_prody])[1]
 path_apps = imp.find_module('prody_apps', [path_apps])[1]
 
diff --git a/scripts/evol b/scripts/evol
index b723cb561..f4a95d247 100755
--- a/scripts/evol
+++ b/scripts/evol
@@ -1,7 +1,9 @@
 #!/usr/bin/env python
 
 import imp
-path = imp.find_module('prody')[1]
+import importlib
+#path = imp.find_module('prody')[1]
+path = importlib.util.find_spec("prody").submodule_search_locations[0]
 apps = imp.find_module('apps', [path])
 apps = imp.load_module('prody.apps', *apps)
 
diff --git a/scripts/prody b/scripts/prody
index 553fdc9e1..b0129eac4 100755
--- a/scripts/prody
+++ b/scripts/prody
@@ -1,7 +1,9 @@
 #!/usr/bin/env python
 
 import imp
-path = imp.find_module('prody')[1]
+import importlib
+#path = imp.find_module('prody')[1]
+path = importlib.util.find_spec("prody").submodule_search_locations[0]
 apps = imp.find_module('apps', [path])
 apps = imp.load_module('prody.apps', *apps)
 

From bcc955af937d7a5c465adbd36890c637729b8eed Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 16 Nov 2022 12:35:09 +0100
Subject: [PATCH 61/73] no sorting of strands in writing

---
 prody/proteins/pdbfile.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/prody/proteins/pdbfile.py b/prody/proteins/pdbfile.py
index f33a27006..cbf69da61 100644
--- a/prody/proteins/pdbfile.py
+++ b/prody/proteins/pdbfile.py
@@ -1366,7 +1366,9 @@ def writePDBStream(stream, atoms, csets=None, **kwargs):
         torf_all_sheets = isSheet(secstrs)
         sheet_secids = secids[torf_all_sheets]
 
-        for sheet_id in np.unique(sheet_secids):
+        unique_sheet_secids, indices = np.unique(sheet_secids, return_index=True)
+        unique_sheet_secids = unique_sheet_secids[indices.argsort()]
+        for sheet_id in unique_sheet_secids:
             torf_strands_in_sheet = np.logical_and(torf_all_sheets, secids==sheet_id)
             strand_indices = secindices[torf_strands_in_sheet]
             numStrands = len(np.unique(strand_indices))
@@ -1604,21 +1606,20 @@ def writePQRStream(stream, atoms, **kwargs):
         write(format_helix(*line))
 
 
-    sorted_sheet = sorted(sheet, key=lambda item: (item[2], item[1]))
     sheet_prev = 'A'
     num_strands_list = []
-    for i, item1 in enumerate(sorted_sheet):
+    for i, item1 in enumerate(sheet):
         if item1[2] != sheet_prev:
-            num_strands = sorted_sheet[i-1][1]
+            num_strands = sheet[i-1][1]
             num_strands_list.append(num_strands)
 
             sheet_prev = item1[2]
 
-            for item2 in sorted_sheet[i-num_strands:i]:
+            for item2 in sheet[i-num_strands:i]:
                 item2.append(num_strands)
 
     num_strands = len(num_strands_list)
-    for item2 in sorted_sheet[i-num_strands+1:]:
+    for item2 in sheet[i-num_strands+1:]:
         item2.append(num_strands)
 
     format_sheet = ('{0:6s} {1:3d} {2:3s}{12:2d} ' +
@@ -1626,7 +1627,7 @@ def writePQRStream(stream, atoms, **kwargs):
                     '{7:3s} {8:1s}{9:4d}{10:1s}' +
                     '{11:2d}\n').format
 
-    for i, line in enumerate(sorted_sheet):
+    for i, line in enumerate(sheet):
         write(format_sheet(*line))
 
     resnames = atoms._getResnames()

From 94d6cf09d95e719d5c92b4879bd3a7dc83c8bb4f Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 19 Jan 2023 17:20:00 +0100
Subject: [PATCH 62/73] destrict biopython for scipion

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index b295048f6..9a0a07233 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,2 +1,2 @@
 [build-system]
-requires = ["setuptools", "wheel", "numpy>=1.10,<1.24", "pyparsing", "biopython<=1.76", "scipy"]
\ No newline at end of file
+requires = ["setuptools", "wheel", "numpy>=1.10,<1.24", "pyparsing", "biopython", "scipy"]
\ No newline at end of file

From 56343aa30dd83967c7fa6f18d6325b2545f7db2e Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 2 Feb 2023 18:45:47 +0100
Subject: [PATCH 63/73] remove clustenm parallel and threads changes

---
 prody/dynamics/clustenm.py | 39 ++------------------------------------
 1 file changed, 2 insertions(+), 37 deletions(-)

diff --git a/prody/dynamics/clustenm.py b/prody/dynamics/clustenm.py
index f396a6d26..59eabb5a6 100644
--- a/prody/dynamics/clustenm.py
+++ b/prody/dynamics/clustenm.py
@@ -297,12 +297,6 @@ def _prep_sim(self, coords, external_forces=[]):
             properties = {'Precision': 'single'}
         elif self._platform in ['CUDA', 'OpenCL']:
             properties = {'Precision': 'single'}
-        elif self._platform == 'CPU':
-            if self._threads == 0:
-                cpus = cpu_count()
-            else:
-                cpus = self._threads
-            properties = {'Threads': str(cpus)}
 
         simulation = Simulation(modeller.topology, system, integrator,
                                 platform, properties)
@@ -536,12 +530,7 @@ def _sample(self, conf):
 
         if self._targeted:
             if self._parallel:
-                if isinstance(self._parallel, int):
-                    repeats = self._parallel
-                elif isinstance(self._parallel, bool):
-                    repeats = cpu_count()
-
-                with Pool(repeats) as p:
+                with Pool(cpu_count()) as p:
                     pot_conf = p.map(self._multi_targeted_sim,
                                      [(conf, coords) for coords in coordsets])
             else:
@@ -625,12 +614,7 @@ def _generate(self, confs):
         sample_method = self._sample_v1 if self._v1 else self._sample
 
         if self._parallel:
-            if isinstance(self._parallel, int):
-                repeats = self._parallel
-            elif isinstance(self._parallel, bool):
-                repeats = cpu_count()
-
-            with Pool(repeats) as p:
+            with Pool(cpu_count()) as p:
                 tmp = p.map(sample_method, [conf for conf in confs])
         else:
             tmp = [sample_method(conf) for conf in confs]
@@ -990,17 +974,10 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
         :arg platform: Architecture on which the OpenMM part runs, default is None.
             It can be chosen as 'CUDA', 'OpenCL' or 'CPU'.
             For efficiency, 'CUDA' or 'OpenCL' is recommended.
-            'CPU' is needed for setting threads per simulation.
         :type platform: str
 
         :arg parallel: If it is True (default is False), conformer generation will be parallelized.
-            This can also be set to a number for how many simulations are run in parallel.
-            Setting 0 or True means run as many as there are CPUs on the machine.
         :type parallel: bool
-
-        :arg threads: Number of threads to use for an individual simulation
-            Default of 0 uses all CPUs on the machine.
-        :type threads: int
         '''
 
         if self._isBuilt():
@@ -1020,19 +997,7 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
         self._rmsd = (0.,) + rmsd if isinstance(rmsd, tuple) else (0.,) + (rmsd,) * n_gens
         self._n_gens = n_gens
         self._platform = kwargs.pop('platform', None)
-
         self._parallel = kwargs.pop('parallel', False)
-        if not isinstance(self._parallel, (int, bool)):
-            raise TypeError('parallel should be an int or bool')
-
-        if self._parallel == 1:
-            # this is equivalent to not actually being parallel
-            self._parallel = False
-        elif self._parallel == 0:
-            # this is a deliberate choice and is documented
-            self._parallel = True
-
-        self._threads = kwargs.pop('threads', 0)
         self._targeted = kwargs.pop('targeted', False)
         self._tmdk = kwargs.pop('tmdk', 15.)
 

From 8a0aef4a7180a61b554d1824bfc2029b19c58f5e Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 29 May 2023 17:46:37 +0200
Subject: [PATCH 64/73] no adjustText for many structs

---
 prody/dynamics/plotting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/dynamics/plotting.py b/prody/dynamics/plotting.py
index e995594ec..072cfc41d 100644
--- a/prody/dynamics/plotting.py
+++ b/prody/dynamics/plotting.py
@@ -377,7 +377,7 @@ def showProjection(ensemble, modes, *args, **kwargs):
         except ImportError:
             pass
         else:
-            if len(modes) == 2:
+            if len(modes) == 2 and len(texts) < 15:
                 adjust_text(ts)
 
     if len(modes) == 2:

From 3581ab34bf98ed3979d2dc43d52c2be773f7c240 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Wed, 12 Jul 2023 15:52:09 +0200
Subject: [PATCH 65/73] shorter atom numbers error

---
 prody/dynamics/editing.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/prody/dynamics/editing.py b/prody/dynamics/editing.py
index f033c3fe9..2c2663116 100644
--- a/prody/dynamics/editing.py
+++ b/prody/dynamics/editing.py
@@ -39,7 +39,7 @@ def extendModel(model, nodes, atoms, norm=False):
         raise ValueError('model must be an NMA instance')
 
     if model.numAtoms() != nodes.numAtoms():
-        raise ValueError('mode and nodes atom numbers must be the same')
+        raise ValueError('atom numbers must be the same')
 
     indices, atommap = extendAtoms(nodes, atoms, model.is3d())
 
@@ -92,7 +92,7 @@ def extendVector(vector, nodes, atoms):
         raise ValueError('vector must be a Vector instance')
 
     if vector.numAtoms() != nodes.numAtoms():
-        raise ValueError('vector and nodes atom numbers must be the same')
+        raise ValueError('atom numbers must be the same')
 
     indices, atommap = extendAtoms(nodes, atoms, vector.is3d())
     extended = Vector(vec[indices], 'Extended ' + str(vector), vector.is3d())
@@ -539,7 +539,7 @@ def interpolateModel(model, nodes, atoms, norm=False, **kwargs):
         raise TypeError('model must be an NMA instance')
 
     if model.numAtoms() != nodes.numAtoms():
-        raise ValueError('model and nodes atom numbers must be the same')
+        raise ValueError('atom numbers must be the same')
 
     if not model.is3d():
         raise ValueError('model must be 3D')

From a0ca9085b0be3516cc7e8596c771a8f363281f0a Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Thu, 13 Jul 2023 13:10:42 +0200
Subject: [PATCH 66/73] clean up atom num errors

---
 prody/dynamics/editing.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/prody/dynamics/editing.py b/prody/dynamics/editing.py
index bfc51aaea..764249896 100644
--- a/prody/dynamics/editing.py
+++ b/prody/dynamics/editing.py
@@ -39,7 +39,7 @@ def extendModel(model, nodes, atoms, norm=False):
         raise ValueError('model must be an NMA instance')
 
     if model.numAtoms() != nodes.numAtoms():
-        raise ValueError('mode and nodes atom numbers must be the same')
+        raise ValueError('atom numbers must be the same')
 
     indices, atommap = extendAtoms(nodes, atoms, model.is3d())
 
@@ -92,7 +92,7 @@ def extendVector(vector, nodes, atoms):
         raise ValueError('vector must be a Vector instance')
 
     if vector.numAtoms() != nodes.numAtoms():
-        raise ValueError('vector and nodes atom numbers must be the same')
+        raise ValueError('atom numbers must be the same')
 
     indices, atommap = extendAtoms(nodes, atoms, vector.is3d())
     extended = Vector(vec[indices], 'Extended ' + str(vector), vector.is3d())
@@ -547,7 +547,7 @@ def interpolateModel(model, nodes, coords, norm=False, **kwargs):
         raise TypeError('model must be an NMA instance')
 
     if model.numAtoms() != nodes.numAtoms():
-        raise ValueError('model and nodes atom numbers must be the same')
+        raise ValueError('atom numbers must be the same')
 
     if not model.is3d():
         raise ValueError('model must be 3D')

From 2b32e8bf254835777cf21ae95e6a43a48baf63b2 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 7 Aug 2023 17:07:31 +0200
Subject: [PATCH 67/73] fix sliceMode and sliceVector for 1D

---
 prody/dynamics/editing.py | 23 ++++++++++++++++-------
 1 file changed, 16 insertions(+), 7 deletions(-)

diff --git a/prody/dynamics/editing.py b/prody/dynamics/editing.py
index 077a28dbc..bccd5db58 100644
--- a/prody/dynamics/editing.py
+++ b/prody/dynamics/editing.py
@@ -124,10 +124,14 @@ def sliceVector(vector, atoms, select):
 
     which, sel = sliceAtoms(atoms, select)
 
-    vec = Vector(vector.getArrayNx3()[
-                 which, :].flatten(),
-                 '{0} slice {1}'.format(str(vector), select),
-                 vector.is3d())
+    if vector.is3d():
+        vec = Vector(vector.getArrayNx3()[which, :].flatten(),
+                    '{0} slice {1}'.format(str(vector), select),
+                    vector.is3d())
+    else:
+        vec = Vector(vector.getArray()[which].flatten(),
+                    '{0} slice {1}'.format(str(vector), select),
+                    vector.is3d())
     return (vec, sel)
 
 def trimModel(model, atoms, select):
@@ -262,9 +266,14 @@ def sliceMode(mode, atoms, select):
 
     which, sel = sliceAtoms(atoms, select)
 
-    vec = Vector(mode.getArrayNx3()[which, :].flatten() *
-                 mode.getVariance()**0.5,
-                 '{0} slice {1}'.format(str(mode), select), mode.is3d())
+    if mode.is3d():
+        vec = Vector(mode.getArrayNx3()[which, :].flatten() *
+                     mode.getVariance()**0.5,
+                    '{0} slice {1}'.format(str(mode), select), mode.is3d())
+    else:
+        vec = Vector(mode.getArray()[which].flatten() *
+                     mode.getVariance()**0.5,
+                    '{0} slice {1}'.format(str(mode), select), mode.is3d())
     return (vec, sel)
 
 

From cc3876c8cbe19102b9fd8b4663037a06336ca703 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 7 Aug 2023 17:07:46 +0200
Subject: [PATCH 68/73] add slice and interp tests

---
 prody/tests/dynamics/test_editing.py | 64 ++++++++++++++++++++++++++--
 1 file changed, 61 insertions(+), 3 deletions(-)

diff --git a/prody/tests/dynamics/test_editing.py b/prody/tests/dynamics/test_editing.py
index f86185363..4c71b9b49 100644
--- a/prody/tests/dynamics/test_editing.py
+++ b/prody/tests/dynamics/test_editing.py
@@ -1,11 +1,13 @@
 """This module contains unit tests for :mod:`~prody.KDTree` module."""
 
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_array_equal, assert_equal
 
-from numpy import concatenate
+from numpy import concatenate, array
 
 from prody.dynamics import calcGNM, calcANM
-from prody.dynamics import extendModel, extendMode, extendVector
+from prody.dynamics import (extendModel, extendMode, extendVector,
+                            sliceModel, sliceMode, sliceVector,
+                            interpolateModel)
 
 from prody.tests import unittest
 from prody.tests.datafiles import parseDatafile
@@ -28,6 +30,14 @@
 EXT3D = concatenate([list(range(i*3, (i+1)*3)) * len(res)
                      for i, res in enumerate(ATOMS.iterResidues())])
 
+SLC1D = concatenate([[i] for i in ATOMS.ca.getIndices()])
+
+SLC3D = concatenate([list(range(i*3, (i+1)*3))
+                     for i in ATOMS.ca.getIndices()])
+
+ANM_AA, NODES_AA = calcANM(ATOMS, selstr='all', cutoff=8)
+GNM_AA = calcGNM(ATOMS, selstr='all', cutoff=4)[0]
+
 
 class TestExtending(unittest.TestCase):
 
@@ -67,3 +77,51 @@ def testVector1d(self):
         ext = extendVector(vector, NODES, ATOMS)[0]
         assert_array_equal(ext._getArray(), vector._getArray()[EXT1D])
 
+
+class TestSlicing(unittest.TestCase):
+
+    def testModel3d(self):
+
+        slc = sliceModel(ANM_AA, NODES_AA, NODES)[0]
+        assert_array_equal(slc._getArray(), ANM_AA._getArray()[SLC3D, :])
+
+    def testModel1d(self):
+
+        slc = sliceModel(GNM_AA, NODES_AA, NODES)[0]
+        assert_array_equal(slc._getArray(), GNM_AA._getArray()[SLC1D, :])
+
+    def testMode3d(self):
+
+        mode = ANM_AA[0]
+        slc = sliceMode(mode, NODES_AA, NODES)[0]
+        assert_array_equal(slc._getArray(), mode._getArray()[SLC3D] *
+                                            mode.getVariance()**0.5)
+
+    def testMode1d(self):
+
+        mode = GNM_AA[0]
+        slc = sliceMode(mode, NODES_AA, NODES)[0]
+        assert_array_equal(slc._getArray(), mode._getArray()[SLC1D] *
+                                            mode.getVariance()**0.5)
+
+    def testVector3d(self):
+
+        vector = ANM_AA[0] * 1
+        slc = sliceVector(vector, NODES_AA, NODES)[0]
+        assert_array_equal(slc._getArray(), vector._getArray()[SLC3D])
+
+    def testVector1d(self):
+
+        vector = GNM_AA[0] * 1
+        slc = sliceVector(vector, NODES_AA, NODES)[0]
+        assert_array_equal(slc._getArray(), vector._getArray()[SLC1D])
+
+
+class TestInterpolating(unittest.TestCase):
+
+    def testModel3d(self):
+
+        interp, atoms = interpolateModel(ANM, NODES, ATOMS)
+        assert_equal(interp._getArray().shape, ANM._getArray()[EXT3D, :].shape)
+        assert_equal(atoms.numAtoms(), NODES_AA.numAtoms())
+        
\ No newline at end of file

From 922d147b655d673b33843d9ee33359c78aa9042c Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Mon, 7 Aug 2023 17:10:53 +0200
Subject: [PATCH 69/73] -

---
 prody/tests/dynamics/test_editing.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/prody/tests/dynamics/test_editing.py b/prody/tests/dynamics/test_editing.py
index 4c71b9b49..c4f7f75f5 100644
--- a/prody/tests/dynamics/test_editing.py
+++ b/prody/tests/dynamics/test_editing.py
@@ -124,4 +124,3 @@ def testModel3d(self):
         interp, atoms = interpolateModel(ANM, NODES, ATOMS)
         assert_equal(interp._getArray().shape, ANM._getArray()[EXT3D, :].shape)
         assert_equal(atoms.numAtoms(), NODES_AA.numAtoms())
-        
\ No newline at end of file

From f500e7e6e34dbab6392acd51679c7494e771d7fd Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 8 Aug 2023 12:44:11 +0200
Subject: [PATCH 70/73] better GammaED docs

---
 prody/dynamics/gamma.py | 30 +++++++++++++++++++++++++++---
 1 file changed, 27 insertions(+), 3 deletions(-)

diff --git a/prody/dynamics/gamma.py b/prody/dynamics/gamma.py
index 7566aa869..4c5af0345 100644
--- a/prody/dynamics/gamma.py
+++ b/prody/dynamics/gamma.py
@@ -37,7 +37,7 @@ class GammaStructureBased(Gamma):
     """Facilitate setting the spring constant based on the secondary structure
     and connectivity of the residues.
 
-    A recent systematic study [LT10]_ of a large set of NMR-structures analyzed
+    A systematic study [LT10]_ of a large set of NMR-structures analyzed
     using a method based on entropy maximization showed that taking into
     consideration properties such as sequential separation between
     contacting residues and the secondary structure types of the interacting
@@ -301,15 +301,37 @@ def gamma(self, dist2, i, j):
         return gamma
 
 
-class GammaGOdMD(Gamma):
+class GammaED(Gamma):
     """Facilitate setting the spring constant based on 
-    sequence distance and spatial distance as in GOdMD.
+    sequence distance and spatial distance as in ed-ENM [OL10]_.
+    This ENM is refined based on comparison to essential dynamics 
+    from MD simulations and can reproduce flexibility in NMR ensembles.
+
+    It has also been implemented in FlexServ [CJ09]_ and 
+    used in MDdMD [SP12]_ and GOdMD [SP13]_.
 
     The sequence distance-dependent term is Cseq/(S**2)
     for S=abs(i,j) <= Slim
 
     The structure distance-dependent term is (Ccart/dist)**Ex
 
+    .. [OL10] Orellana L, Rueda M, Ferrer-Costa C, Lopez-Blanco JR, Chacón P, Orozco M. 
+       Approaching Elastic Network Models to Molecular Dynamics Flexibility.
+       *J Chem Theory Comput* **2010** 6(9):2910-23.
+
+    .. [CJ09] Camps J, Carrillo O, Emperador A, Orellana L, Hospital A, Rueda M, 
+       Cicin-Sain D, D'Abramo M, Gelpí JL, Orozco M.
+       FlexServ: an integrated tool for the analysis of protein flexibility.
+       *Bioinformatics* **2009** 25(13):1709-10.
+
+    .. [SP12] Sfriso P, Emperador A, Orellana L, Hospital A, Gelpí JL, Orozco M.
+       Finding Conformational Transition Pathways from Discrete Molecular Dynamics Simulations.
+       *J Chem Theory Comput* **2012** 8(11):4707-18.
+
+    .. [SP13] Sfriso P, Hospital A, Emperador A, Orozco M. 
+       Exploration of conformational transition pathways from coarse-grained simulations.
+       *Bioinformatics* **2013** 29(16):1980-6.     
+       
     **Example**:
 
     Let's parse coordinates from a PDB file.
@@ -396,3 +418,5 @@ def gamma(self, dist2, i, j):
         else:
             dist = dist2**0.5
             return (Ccart/dist)**Ex
+
+GammaGOdMD = GammaED

From a42aae1d7466d2598cbe10dd6ba13c4db1f97ff1 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 8 Aug 2023 12:45:07 +0200
Subject: [PATCH 71/73] better mapping docs

---
 prody/proteins/compare.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/prody/proteins/compare.py b/prody/proteins/compare.py
index 6c4171e71..9c624560f 100644
--- a/prody/proteins/compare.py
+++ b/prody/proteins/compare.py
@@ -923,12 +923,12 @@ def mapChainOntoChain(mobile, target, **kwargs):
     :type overlap: float
 
     :keyword mapping: what method will be used if the trivial mapping based on residue numbers 
-        fails. If ``"ce"`` or ``"cealign"``, then the CE algorithm [IS98]_ will be 
+        fails. If ``"ce"`` or ``"cealign"``, then the CE structural alignment [IS98]_ will be
         performed. It can also be a list of prealigned sequences, a :class:`.MSA` instance,
         or a dict of indices such as that derived from a :class:`.DaliRecord`.
         If set to **True** then the sequence alignment from Biopython 
         will be used. If set to **False**, only the trivial mapping will be performed. 
-        Default is **"auto"**
+        Default is **"auto"**, which means try sequence alignment then CE.
     :type mapping: list, str, bool
 
     :keyword pwalign: if **True**, then pairwise sequence alignment will 

From 343f75c3fa0baf5eebbadffc6b1b8aaaafc027a7 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 8 Aug 2023 14:27:35 +0200
Subject: [PATCH 72/73] add custom cutoff to anm app

---
 prody/apps/prody_apps/prody_anm.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/prody/apps/prody_apps/prody_anm.py b/prody/apps/prody_apps/prody_anm.py
index 7574aba7f..f2eb02574 100644
--- a/prody/apps/prody_apps/prody_anm.py
+++ b/prody/apps/prody_apps/prody_anm.py
@@ -13,7 +13,7 @@
 for key, txt, val in [
     ('model', 'index of model that will be used in the calculations', 1),
     ('altloc', 'alternative location identifiers for residues used in the calculations', "A"),
-    ('cutoff', 'cutoff distance (A)', 15.),
+    ('cutoff', 'cutoff distance (A)', '15.'),
     ('gamma', 'spring constant', '1.'),
     ('sparse', 'use sparse matrices', False),
     ('kdtree', 'use kdtree for Hessian', False),    
@@ -57,7 +57,6 @@ def prody_anm(pdb, **kwargs):
 
     selstr = kwargs.get('select')
     prefix = kwargs.get('prefix')
-    cutoff = kwargs.get('cutoff')
     sparse = kwargs.get('sparse')
     kdtree = kwargs.get('kdtree')
     nmodes = kwargs.get('nmodes')
@@ -91,6 +90,19 @@ def prody_anm(pdb, **kwargs):
             raise NameError("Please provide gamma as a float or ProDy Gamma class")
         except TypeError:
             raise TypeError("Please provide gamma as a float or ProDy Gamma class")
+        
+    try:
+        cutoff = float(kwargs.get('cutoff'))
+        LOGGER.info("Using cutoff {0}".format(cutoff))
+    except ValueError:
+        try:
+            import math
+            cutoff = eval(kwargs.get('cutoff'))
+            LOGGER.info("Using cutoff {0}".format(cutoff))
+        except NameError:
+            raise NameError("Please provide cutoff as a float or equation using math")
+        except TypeError:
+            raise TypeError("Please provide cutoff as a float or equation using math")
 
     anm = prody.ANM(pdb.getTitle())
 
@@ -276,7 +288,7 @@ def addCommand(commands):
 
     group = addNMAParameters(subparser)
 
-    group.add_argument('-c', '--cutoff', dest='cutoff', type=float,
+    group.add_argument('-c', '--cutoff', dest='cutoff', type=str,
         default=DEFAULTS['cutoff'], metavar='FLOAT',
         help=HELPTEXT['cutoff'] + ' (default: %(default)s)')
 

From 84affd90a4f5b5f1eb26fac87edd85e3ea1bc5f0 Mon Sep 17 00:00:00 2001
From: James Krieger <jamesmkrieger@gmail.com>
Date: Tue, 8 Aug 2023 14:29:32 +0200
Subject: [PATCH 73/73] add GammaED to __all__

---
 prody/dynamics/gamma.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prody/dynamics/gamma.py b/prody/dynamics/gamma.py
index 4c5af0345..19e64f673 100644
--- a/prody/dynamics/gamma.py
+++ b/prody/dynamics/gamma.py
@@ -6,7 +6,7 @@
 
 from prody.atomic import Atomic
 
-__all__ = ['Gamma', 'GammaStructureBased', 'GammaVariableCutoff', 'GammaGOdMD']
+__all__ = ['Gamma', 'GammaStructureBased', 'GammaVariableCutoff', 'GammaED', 'GammaGOdMD']
 
 
 class Gamma(object):