# Copyright 2016 Autodesk Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Interface to quaternion-based atom alignment routines.
This code has been translated from Fortran to Python
"""
import numpy as np
[docs]def qtrfit(mol, pos, ref, massweight=True):
r""" Transform the atomic positions in ``pos`` to minimize the overall distance to those in ``ref``
This function minimizes the distance between pos and ref. If `massweight`=True, we minimize the mass-weighted distance:
..math::
D = \sqrt{ \sum_{i=1}^{\text{N}} m_i \left| \mathbf{pos}_i - \mathbf{ref}_i \right|^2 }
Otherwise, the root-squared distance is used:
..math::
D = \sqrt{ \sum_{i=1}^{\text{N}} \left| \mathbf{pos}_i - \mathbf{ref}_i \right|^2 }
Args:
mol (moldesign.Molecule): molecule that these arrays describe
pos (u.Array): The Nx3 array to be rotated
ref (u.Array): The Nx3 array to match
massweight (bool): whether to use mass-weighted coordinates (the default)
Returns:
u.Array: translated and rotated version of pos with minimum RMSD
References:
TODO: reference to original qtrfit code (and license, if it exists)
"""
N = 1.0 * len(pos)
if massweight:
weights = mol.masses
else:
weights = np.ones(mol.num_atoms)
totweight = weights.sum()
pcenter = (weights[:, None] * pos).sum(axis=0) / totweight
rcenter = (weights[:, None] * ref).sum(axis=0) / totweight
p = pos - pcenter
r = ref - rcenter
# yadda yadda yadda
return p + rcenter