From 315e7d52c4bf9ef3ed22c796f2f1511d59d7abfb Mon Sep 17 00:00:00 2001 From: James Krieger Date: Fri, 26 Feb 2021 17:35:11 +0000 Subject: [PATCH] add calcConvergenceRMSDs --- prody/dynamics/adaptive.py | 28 +++++++++++++++++++++++++++- prody/dynamics/comd.py | 30 ++++++++++++++++++++++++++++-- 2 files changed, 55 insertions(+), 3 deletions(-) diff --git a/prody/dynamics/adaptive.py b/prody/dynamics/adaptive.py index 1fb44f7bc..6a054a232 100644 --- a/prody/dynamics/adaptive.py +++ b/prody/dynamics/adaptive.py @@ -20,7 +20,7 @@ from prody.atomic import Atomic, AtomMap from prody.utilities import getCoords, createStringIO, importLA, checkCoords, copy -from prody.measure.transform import calcTransformation, superpose, applyTransformation, calcRMSD +from prody.measure.transform import calcTransformation, superpose, applyTransformation, calcRMSD, getRMSD from prody.measure.measure import calcDeformVector, calcDistance from prody.ensemble.ensemble import Ensemble @@ -717,6 +717,32 @@ def getRMSDsB(self): return calcRMSD(self._atomsB, self._confs[:, indices], weights) + def getConvergenceRMSDs(self): + if self._confs is None or self._coords is None: + return None + + indices = self._indices + if indices is None: + indices = np.arange(self._confs.shape[1]) + + weights = self._weights[indices] if self._weights is not None else None + + n_confs = self.numConfs() + confsA = self._confs[::2] + if n_confs % 2: + confsB = self._confs[-2::-2] + else: + confsB = self._confs[-1::-2] + + n_confsA = len(confsA) + n_confsB = len(confsB) + RMSDs = zeros((n_confsA, n_confsB)) + for i in range(n_confsA): + for j in range(n_confsB): + RMSDs[i, j] = getRMSD(confsA[i, indices], confsB[j, indices], weights) + + return RMSDs + def _generate(self, confs, **kwargs): LOGGER.info('Sampling conformers in generation %d ...' % self._cycle) diff --git a/prody/dynamics/comd.py b/prody/dynamics/comd.py index 733e8f258..278b4034e 100644 --- a/prody/dynamics/comd.py +++ b/prody/dynamics/comd.py @@ -1,6 +1,6 @@ from prody.proteins.pdbfile import parsePDB, writePDBStream, parsePDBStream -from prody.measure.transform import calcRMSD, calcTransformation, superpose, applyTransformation -from prody.measure.measure import buildDistMatrix, calcDeformVector, calcDistance +from prody.measure.transform import calcRMSD, calcTransformation, getRMSD, applyTransformation +from prody.measure.measure import buildDistMatrix, calcDeformVector from prody.ensemble.ensemble import Ensemble from prody.utilities import createStringIO, importLA, checkCoords, copy @@ -537,6 +537,32 @@ def getRMSDsB(self): return calcRMSD(self._atomsB, self._confs[:, indices], weights) + def getConvergenceRMSDs(self): + if self._confs is None or self._coords is None: + return None + + indices = self._indices + if indices is None: + indices = np.arange(self._confs.shape[1]) + + weights = self._weights[indices] if self._weights is not None else None + + n_confs = self.numConfs() + confsA = self._confs[::2] + if n_confs % 2: + confsB = self._confs[-2::-2] + else: + confsB = self._confs[-1::-2] + + n_confsA = len(confsA) + n_confsB = len(confsB) + RMSDs = zeros((n_confsA, n_confsB)) + for i in range(n_confsA): + for j in range(n_confsB): + RMSDs[i, j] = getRMSD(confsA[i, indices], confsB[j, indices], weights) + + return RMSDs + def run(self, cutoff=15., n_modes=20, gamma=1., n_confs=50, rmsd=1.0, n_gens=5, solvent='imp', sim=False, force_field=None, temp=303.15, t_steps_i=1000, t_steps_g=7500,