Source Code for Module lib.structure.statistics

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2011-2013 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Module docstring. 
 23  """Module for handling all types of structural statistics.""" 
 24   
 25  # Python module imports. 
 26  from numpy import float64, mean, sqrt, std, zeros 
 27  from numpy.linalg import norm 
 28   
 29   
 30 -def atomic_rmsd(coord, verbosity=0): 
 31      """Determine the RMSD for the given atomic coordinates. 
 32   
 33      This is the per atom RMSD to the mean structure. 
 34   
 35   
 36      @keyword coord:     The array of molecular coordinates.  The first dimension corresponds to the model, the second the atom, the third the coordinate. 
 37      @type coord:        rank-3 numpy array 
 38      @return:            The RMSD value. 
 39      @rtype:             float 
 40      """ 
 41   
 42      # Init. 
 43      M = len(coord) 
 44      N = len(coord[0]) 
 45      model_rmsd = zeros(M, float64) 
 46      mean_str = zeros((N, 3), float64) 
 47   
 48      # Calculate the mean structure. 
 49      calc_mean_structure(coord, mean_str) 
 50   
 51      # Loop over the models. 
 52      for i in range(M): 
 53          # Loop over the atoms. 
 54          for j in range(N): 
 55              # The vector connecting the mean to model atom. 
 56              vect = mean_str[j] - coord[i][j] 
 57   
 58              # The atomic RMSD. 
 59              model_rmsd[i] += norm(vect)**2 
 60   
 61          # Normalise, and sqrt. 
 62          model_rmsd[i] = sqrt(model_rmsd[i] / N) 
 63   
 64          # Print out. 
 65          if verbosity: 
 66              print("Model %2s RMSD:  %s" % (i, model_rmsd[i])) 
 67   
 68      # Calculate the mean. 
 69      rmsd_mean = mean(model_rmsd) 
 70   
 71      # Calculate the normal non-biased standard deviation. 
 72      try: 
 73          rmsd_sd = std(model_rmsd, ddof=1) 
 74   
 75      # Handle old numpy versions not having the ddof argument. 
 76      except TypeError: 
 77          rmsd_sd = std(model_rmsd) * sqrt(M / (M - 1.0)) 
 78   
 79      # Printout. 
 80      if verbosity: 
 81          print("\nGlobal RMSD:  %s +/- %s" % (rmsd_mean, rmsd_sd)) 
 82   
 83      # Return the average RMSD. 
 84      return rmsd_mean 
 85   
 86   
 87 -def calc_mean_structure(coord=None, mean=None): 
 88      """Average the coordinates. 
 89   
 90      @keyword coord:     The list of coordinates of all models to superimpose.  The first index is the models, the second is the atomic positions, and the third is the xyz coordinates. 
 91      @type coord:        list of numpy rank-2, Nx3 arrays 
 92      @keyword mean:      The data storage for the mean structure. 
 93      @type mean:         numpy rank-2, Nx3 array 
 94      """ 
 95   
 96      # The number of atoms. 
 97      N = len(coord[0]) 
 98      M = len(coord) 
 99   
100      # Clear the mean data structure. 
101      for i in range(N): 
102          mean[i] = [0.0, 0.0, 0.0] 
103   
104      # Loop over the atoms. 
105      for i in range(N): 
106          # Loop over the models. 
107          for j in range(M): 
108              mean[i] += coord[j][i] 
109   
110          # Average. 
111          mean[i] = mean[i] / M 
112