Source Code for Module lib.structure.statistics

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2011,2013,2015-2016 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 array, float64, mean, ones, 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 quadratic mean. 
 69      rmsd_mean = sqrt(mean(model_rmsd**2)) 
 70   
 71      # Calculate the normal non-biased quadratic standard deviation. 
 72      try: 
 73          rmsd_sd = sqrt(std(model_rmsd**2, ddof=1)) 
 74   
 75      # Handle old numpy versions not having the ddof argument. 
 76      except TypeError: 
 77          rmsd_sd = sqrt(std(model_rmsd**2) * 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, weights=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 weights:   The weights for each structure. 
 93      @type weights:      list of float 
 94      @keyword mean:      The data storage for the mean structure. 
 95      @type mean:         numpy rank-2, Nx3 array 
 96      """ 
 97   
 98      # The number of atoms. 
 99      N = len(coord[0]) 
100      M = len(coord) 
101      if weights is None: 
102          weights = ones(M, float64) 
103      else: 
104          weights = array(weights, float64) 
105   
106      # Clear the mean data structure. 
107      for i in range(N): 
108          mean[i] = [0.0, 0.0, 0.0] 
109   
110      # Loop over the atoms. 
111      for i in range(N): 
112          # Loop over the models. 
113          for j in range(M): 
114              mean[i] += weights[j] * coord[j][i] 
115   
116          # Average. 
117          mean[i] = mean[i] / weights.sum() 
118   
119   
120 -def per_atom_rmsd(coord, verbosity=0): 
121      """Determine the per-atom RMSDs for the given atomic coordinates. 
122   
123      This is the per atom RMSD to the mean structure. 
124   
125   
126      @keyword coord:     The array of molecular coordinates.  The first dimension corresponds to the model, the second the atom, the third the coordinate. 
127      @type coord:        rank-3 numpy array 
128      @return:            The list of RMSD values for each atom. 
129      @rtype:             rank-1 numpy float64 array 
130      """ 
131   
132      # Init. 
133      M = len(coord) 
134      N = len(coord[0]) 
135      model_rmsd = zeros(M, float64) 
136      mean_str = zeros((N, 3), float64) 
137      rmsd = zeros(N, float64) 
138   
139      # Calculate the mean structure. 
140      calc_mean_structure(coord, mean_str) 
141   
142      # Loop over the atoms. 
143      for j in range(N): 
144          # Loop over the models. 
145          for i in range(M): 
146              # The vector connecting the mean to model atom. 
147              vect = mean_str[j] - coord[i][j] 
148   
149              # The atomic RMSD. 
150              rmsd[j] += norm(vect)**2 
151   
152          # Normalise, and sqrt. 
153          rmsd[j] = sqrt(rmsd[j] / M) 
154   
155      # Return the RMSDs. 
156      return rmsd 
157