diff --git a/prody/dynamics/adaptive.py b/prody/dynamics/adaptive.py index 1fb44f7bc..784a706a9 100644 --- a/prody/dynamics/adaptive.py +++ b/prody/dynamics/adaptive.py @@ -14,13 +14,14 @@ import time from numbers import Integral, Number +from decimal import Decimal, ROUND_HALF_UP import numpy as np from prody import LOGGER 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 +718,36 @@ 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() + n_confsA = int(Decimal(n_confs/2).to_integral(rounding=ROUND_HALF_UP)) + + confsA = self._confs[:n_confsA] + if n_confs % 2: + confsB = self._confs[n_confsA:] + else: + confsB = self._confs[n_confsA:] + + RMSDs = np.zeros((n_confs-9)) + n = 0 + for i in range(n_confsA): + for j in range(2): + RMSDs[n] = getRMSD(confsA[i+j], confsB[n_confsA-(i+1)]) + n += 1 + if i == n_confsA - 1: + break + + 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..38b3d1ace 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 @@ -14,6 +14,7 @@ from random import random import os.path import sys +from decimal import Decimal, ROUND_HALF_UP from .adaptive import ONEWAY, ALTERNATING, SERIAL, DEFAULT from .hybrid import Hybrid @@ -537,6 +538,36 @@ 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() + n_confsA = int(Decimal(n_confs/2).to_integral(rounding=ROUND_HALF_UP)) + + confsA = self._confs[:n_confsA] + if n_confs % 2: + confsB = self._confs[n_confsA:] + else: + confsB = self._confs[n_confsA:] + + RMSDs = np.zeros((n_confs-9)) + n = 0 + for i in range(n_confsA): + for j in range(2): + RMSDs[n] = getRMSD(confsA[i+j], confsB[n_confsA-(i+1)]) + n += 1 + if i == n_confsA - 1: + break + + 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,