Author: bugman Date: Thu Jun 2 19:09:39 2016 New Revision: 28222 URL: http://svn.gna.org/viewcvs/relax?rev=28222&view=rev Log: Implemented the per-atom RMSD calculation for the structure.rmsd user function.
Modified: trunk/lib/structure/statistics.py trunk/pipe_control/structure/main.py Modified: trunk/lib/structure/statistics.py URL: http://svn.gna.org/viewcvs/relax/trunk/lib/structure/statistics.py?rev=28222&r1=28221&r2=28222&view=diff ============================================================================== --- trunk/lib/structure/statistics.py (original) +++ trunk/lib/structure/statistics.py Thu Jun 2 19:09:39 2016 @@ -1,6 +1,6 @@ ############################################################################### # # -# Copyright (C) 2011-2015 Edward d'Auvergne # +# Copyright (C) 2011-2016 Edward d'Auvergne # # # # This file is part of the program relax (http://www.nmr-relax.com). # # # @@ -115,3 +115,42 @@ # Average. mean[i] = mean[i] / weights.sum() + + +def per_atom_rmsd(coord, verbosity=0): + """Determine the per-atom RMSDs for the given atomic coordinates. + + This is the per atom RMSD to the mean structure. + + + @keyword coord: The array of molecular coordinates. The first dimension corresponds to the model, the second the atom, the third the coordinate. + @type coord: rank-3 numpy array + @return: The list of RMSD values for each atom. + @rtype: rank-1 numpy float64 array + """ + + # Init. + M = len(coord) + N = len(coord[0]) + model_rmsd = zeros(M, float64) + mean_str = zeros((N, 3), float64) + rmsd = zeros(N, float64) + + # Calculate the mean structure. + calc_mean_structure(coord, mean_str) + + # Loop over the atoms. + for j in range(N): + # Loop over the models. + for i in range(M): + # The vector connecting the mean to model atom. + vect = mean_str[j] - coord[i][j] + + # The atomic RMSD. + rmsd[j] += norm(vect)**2 + + # Normalise, and sqrt. + rmsd[j] = sqrt(rmsd[j] / M) + + # Return the RMSDs. + return rmsd Modified: trunk/pipe_control/structure/main.py URL: http://svn.gna.org/viewcvs/relax/trunk/pipe_control/structure/main.py?rev=28222&r1=28221&r2=28222&view=diff ============================================================================== --- trunk/pipe_control/structure/main.py (original) +++ trunk/pipe_control/structure/main.py Thu Jun 2 19:09:39 2016 @@ -44,7 +44,7 @@ from lib.structure.internal.object import Internal from lib.structure.pca import pca_analysis from lib.structure.represent.diffusion_tensor import diffusion_tensor -from lib.structure.statistics import atomic_rmsd +from lib.structure.statistics import atomic_rmsd, per_atom_rmsd from lib.structure.superimpose import fit_to_first, fit_to_mean from lib.warnings import RelaxWarning, RelaxNoPDBFileWarning, RelaxZeroVectorWarning from pipe_control import molmol, pipes @@ -1308,7 +1308,7 @@ cdp.structure.load_xyz(file_path, read_mol=read_mol, set_mol_name=set_mol_name, read_model=read_model, set_model_num=set_model_num, verbosity=verbosity) -def rmsd(pipes=None, models=None, molecules=None, atom_id=None): +def rmsd(pipes=None, models=None, molecules=None, atom_id=None, atomic=False): """Calculate the RMSD between the loaded models. The RMSD value will be placed into the current data pipe cdp.structure.rmsd data structure. @@ -1322,6 +1322,8 @@ @type molecules: None or list of lists of str @keyword atom_id: The atom identification string of the coordinates of interest. This matches the spin ID string format. @type atom_id: str or None + @keyword atomic: A flag which if True will allow for per-atom RMSDs to be additionally calculated. + @type atomic: bool @return: The RMSD value. @rtype: float """ @@ -1332,10 +1334,34 @@ # Assemble the structural coordinates. coord, ids, mol_names, res_names, res_nums, atom_names, elements = assemble_structural_coordinates(pipes=pipes, models=models, molecules=molecules, atom_id=atom_id) + # Per-atom RMSDs. + if atomic: + # Printout. + print("\nCalculating the atomic-level RMSDs.") + + # Calculate the per-atom RMSDs. + rmsd = per_atom_rmsd(coord, verbosity=0) + + # Loop over the atoms. + for i in range(len(res_nums)): + # The spin identification string. + id = generate_spin_id_unique(mol_name=mol_names[i], res_num=res_nums[i], res_name=res_names[i], spin_num=i, spin_name=atom_names[i]) + + # Get the spin container. + spin_cont = return_spin(id) + + # Skip the spin if it doesn't exist. + if spin_cont == None: + continue + + # Store the value. + spin_cont.pos_rmsd = rmsd[i] + # Calculate the RMSD. + print("\nCalculating the global RMSD.") cdp.structure.rmsd = atomic_rmsd(coord, verbosity=1) - # Return the RMSD. + # Return the global RMSD. return cdp.structure.rmsd _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-commits mailing list relax-commits@gna.org To unsubscribe from this list, get a password reminder, or change your subscription options, visit the list information page at https://mail.gna.org/listinfo/relax-commits