Improved inertia alignment method and improved documentation (#122)

* Added better readme and --help information.

* Updated the Inertia Hungarian alignment method to align with all three Inertia axes, not just the principal.

* Set it as the default reordering method to inertia-hungarian.

* Change some CLI arguments to be more consistent with standard practice (remove double letter, single dash arguments).


Co-authored-by: hagenmuenkler <[email protected]>

.gitignore

@@ -2,6 +2,7 @@
 .vim
 env
 *.sqlite3
+*.ipynb_checkpoints
 
 *.py[cod]
 

.pre-commit-config.yaml

@@ -70,7 +70,7 @@ repos:
     hooks:
       - id: mypy
         name: checking types
-        entry: mypy
+        entry: mypy --explicit-package-bases
         language: system
         types: [ python ]
 

Makefile

@@ -47,6 +47,9 @@ build:
 upload:
 	${python} -m twine upload ./dist/*
 
+start-jupyter:
+	${python} -m jupyterlab
+
 ## Version
 
 version:

README.rst

@@ -1,90 +1,83 @@
 Calculate Root-mean-square deviation (RMSD) of Two Molecules Using Rotation
 ===========================================================================
 
-The root-mean-square deviation (RMSD) is calculated, using Kabsch algorithm
-(1976) or Quaternion algorithm (1991) for rotation, between two Cartesian
-coordinates in either ``.xyz`` or ``.pdb`` format, resulting in the minimal
-RMSD.
+The root mean Square Deviation (RMSD) is the most common metric for measuring structural similarity between two structures. It is typically used in molecular biology, chemistry, and computational chemistry.
 
-For more information please read RMSD_ and `Kabsch algorithm`_.
+However, the result can become misleadingly large unless the input data has pre-optimized translation and rotation of the molecules in question.
+This solution will perform this optimization before calculating optimal (minimal) RMSD values.
 
-.. _RMSD: http://en.wikipedia.org/wiki/Root-mean-square_deviation
-.. _Kabsch algorithm: http://en.wikipedia.org/wiki/Kabsch_algorithm
+Additionally, if the atoms in the molecules are not correctly ordered, optimal rotation is impossible to achieve.
+This tool utilizes several ways to solve this problem.
 
-Motivation
-----------
+For more details, see below and read RMSD_ and `Kabsch algorithm`_.
 
-You have molecule A and B and want to calculate the structural difference
-between those two. If you just calculate the RMSD_ straight-forward you might
-get a too big of a value as seen below. You would need to first recenter the
-two molecules and then rotate them unto each other to get the true minimal
-RMSD. This is what this script does.
+.. contents:: Overview
+    :depth: 1
 
-==========  ===========  ==========
-No Changes  Re-centered  Rotated
-----------  -----------  ----------
-|begin|     |translate|  |rotate|
-==========  ===========  ==========
-RMSD 2.50   RMSD 1.07    RMSD 0.25
-==========  ===========  ==========
+Installation
+============
 
-.. |begin| image:: https://raw.githubusercontent.com/charnley/rmsd/master/img/plot_beginning.png
-.. |translate| image:: https://raw.githubusercontent.com/charnley/rmsd/master/img/plot_translated.png
-.. |rotate| image:: https://raw.githubusercontent.com/charnley/rmsd/master/img/plot_rotated.png
+The easiest is to get the program via ``pip``.
 
+.. code-block:: bash
 
-Citation
---------
+    pip install rmsd
 
-- **Implementation**:
-    Calculate Root-mean-square deviation (RMSD) of Two Molecules Using Rotation, GitHub,
-    http://github.com/charnley/rmsd, <git commit hash or version number>
+There is only one Python file, so you can also download `calculate_rmsd.py` and put it in your bin folder.
 
-- **Kabsch algorithm**:
-    Kabsch W., 1976,
-    A solution for the best rotation to relate two sets of vectors,
-    Acta Crystallographica, A32:922-923,
-    doi: http://dx.doi.org/10.1107/S0567739476001873
+.. code-block:: bash
 
-- **Quaternion algorithm**:
-    Michael W. Walker and Lejun Shao and Richard A. Volz, 1991,
-    Estimating 3-D location parameters using dual number quaternions, CVGIP: Image Understanding, 54:358-367,
-    doi: http://dx.doi.org/10.1016/1049-9660(91)90036-o
+    wget -O calculate_rmsd https://raw.githubusercontent.com/charnley/rmsd/master/rmsd/calculate_rmsd.py
+    chmod +x calculate_rmsd
 
-Please cite this project when using it for scientific publications.
+Details
+=======
 
+To calculate the structural difference between two molecules, you might initially compute the RMSD directly (**Figure 1.a**).
+However, this straightforward approach could give you a misleadingly large value.
+To get the true minimal RMSD, you must adjust for translation (**Figure 1.b**) and rotation (**Figure 1.c**). This process aligns the two molecules best, ensuring the RMSD accurately reflects their structural similarity after optimal alignment.
 
-Installation
-------------
+.. list-table::
+   :header-rows: 1
 
-Easiest is to get the program vis PyPi under the package name ``rmsd``,
+   * - 1.a
+     - 1.b
+     - 1.c
 
-.. code-block:: bash
+   * - |fig1.a|
+     - |fig1.b|
+     - |fig1.c|
 
-    pip install rmsd
+   * - RMSD = 2.8
+     - RMSD = 0.8
+     - RMSD = 0.2
 
+Atom reordering methods are used when the atoms in two molecules are not in the same order (**Figure 2.a**).
+While brute-force through all possible atom combinations and calculating the optimal rotation for each is possible, this approach is computationally infeasible for large structures, as it scales $O(N!)$.
+Instead, the implemented algorithms efficiently find the optimal mapping of atoms between the two structures using smarter techniques.
 
-or download the project from GitHub via
+Each method has limitations because finding the best atom mapping depends on properly aligning structures.
+This is usually done by comparing atom-pair distances. If the molecules are aligned, using the Hungarian_ cost minimization for atom distance works well.
+If not, you can align the Inertia_ eigenvectors (**Figure 2.b**) as an approximation to align the molecules.
+Or, use atomic descriptors (**Figure 2.c**), independent of the coordinate system, to reorder the atoms. Note that all reordering methods have limitations and drawbacks, and the actual order might not be found.
 
-.. code-block:: bash
+.. list-table::
+   :header-rows: 1
 
-    git clone https://github.com/charnley/rmsd
+   * - 2.a
+     - 2.b
+     - 2.c
 
-
-There is only one Python file, so you can also download `calculate_rmsd.py` and
-put it in your bin folder.
-
-.. code-block:: bash
-
-    wget -O calculate_rmsd https://raw.githubusercontent.com/charnley/rmsd/master/rmsd/calculate_rmsd.py
-    chmod +x calculate_rmsd
+   * - |fig2.a|
+     - |fig2.b|
+     - |fig2.c|
 
 Usage examples
---------------
+==============
 
 Use ``calculate_rmsd --help`` to see all the features. Usage is pretty straight
 forward, call ``calculate_rmsd`` with two structures in either ``.xyz`` or
-``.pdb``. In this example Ethane has the exact same structure, but is
+``.pdb``. In this example, Ethane has the same structure but is
 translated in space, so the RMSD should be zero.
 
 .. code-block:: bash
@@ -99,33 +92,104 @@ visual comparison. The output will be in XYZ format.
 
     calculate_rmsd --no-hydrogen --print tests/ethane.xyz tests/ethane_mini.xyz
 
-If the atoms are scrambled and not aligned you can use the ``--reorder``
-argument which will align the atoms from structure B unto A. Use
-``--reorder-method`` to select what method for reordering. Choose between
-Hungarian_ (default), distance (very approximate) and brute force (slow).
+If the atoms are scrambled and not aligned, you can use the ``--reorder``
+argument, which will align the atoms from structure B onto A.
 
-.. _Hungarian: https://en.wikipedia.org/wiki/Hungarian_algorithm
+Use ``--reorder-method`` to select the reordering method.
+Choose between
+Inertia_ aligned Hungarian_ distance ``inertia-hungarian`` (default),
+Hungarian_ distance ``hungarian`` (if the structure is already aligned),
+sorted distance ``distance``,
+atomic representation ``qml``,
+and brute force ``brute`` (for reference, don't use this).
+More details on which to use in ``--help``.
 
 .. code-block:: bash
 
     calculate_rmsd --reorder tests/water_16.xyz tests/water_16_idx.xyz
 
+If you want to run multiple calculations simultaneously, it's best not to rely solely on the script.
+Instead, you can use GNU Parallel to handle this efficiently. For example, compare all ``ethane_*`` molecules using two cores and print one file and the RMSD per line.
+Bash is good for stuff like that.
 
-It is also possible to use RMSD as a library in other scripts, see
-``example.py`` and ``tests/*`` for example usage.
-
+.. code-block:: bash
 
-Problems?
----------
+    find tests/resources -name "ethane_*xyz" | parallel -j2 "echo -n '{} ' && calculate_rmsd --reorder --no-hydrogen tests/resources/ethane.xyz {}"
 
-Submit issues or pull requests on GitHub.
+It is also possible to use RMSD as a library in other scripts; see ``tests/*`` for example usage.
 
+Known problems
+==============
 
-A note on PDB
--------------
+Found a bug? Submit issues or pull requests on GitHub.
 
-Protein Data Bank format (PDB) is column-based; however, countless examples of non-standard ``.pdb`` files exist.
+**Note on PDB format.** Protein Data Bank format (PDB) is column-based; however, countless examples of non-standard ``.pdb`` files exist.
 We try to read them, but if you have trouble reading the file, check if the file format is compliant with PDB.
 For example, some hydrogens are noted as ``HG11``, which we assume is not mercury.
 
 - https://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#ATOM
+
+Citation
+========
+
+Please cite this project when using it for scientific publications. And cite the relevant methods implemented.
+
+**Implementation**:
+Calculate Root-mean-square deviation (RMSD) of Two Molecules Using Rotation, GitHub,
+http://github.com/charnley/rmsd, <git commit hash or version number>
+
+
+.. list-table::
+   :header-rows: 1
+
+   * - Method
+     - Argument
+     - Citation
+
+   * - **Kabsch**
+     - ``--rotation-method kabsch`` (Default)
+     - Wolfgang Kabsch (1976),
+       Acta Crystallographica, A32:922-923
+
+       http://dx.doi.org/10.1107/S0567739476001873
+
+   * - **Quaternion**
+     - ``--rotation-method quaternion``
+     - Walker, Shao & Volz (1991),
+       CVGIP: Image Understanding, 54:358-367,
+
+       http://dx.doi.org/10.1016/1049-9660(91)90036-o
+
+   * - **Distance Hungarian Assignment**
+     - ``--reorder-method inertia-hungarian`` (Default)
+     - Crouse (2016). Vol. 52, Issue 4, pp. 1679–1696, IEEE.
+
+       http://dx.doi.org/10.1109/TAES.2016.140952
+
+   * - **FCHL19**
+     - ``--reorder-method qml``
+     - Christensen et al (2020), J. Chem. Phys. 152, 044107
+
+       https://doi.org/10.1063/1.5126701
+
+References
+==========
+
+- https://en.wikipedia.org/wiki/Root-mean-square_deviation
+- https://en.wikipedia.org/wiki/Kabsch_algorithm
+- https://en.wikipedia.org/wiki/Hungarian_algorithm
+- https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linear_sum_assignment.html
+
+.. _RMSD: https://en.wikipedia.org/wiki/Root-mean-square_deviation
+.. _Kabsch algorithm: https://en.wikipedia.org/wiki/Kabsch_algorithm
+.. _Hungarian: https://en.wikipedia.org/wiki/Hungarian_algorithm
+.. _Inertia: https://en.wikipedia.org/wiki/Moment_of_inertia
+
+
+.. |fig1.a| image:: https://raw.githubusercontent.com/charnley/rmsd/master/docs/figures/fig_rmsd_nothing.png
+.. |fig1.b| image:: https://raw.githubusercontent.com/charnley/rmsd/master/docs/figures/fig_rmsd_recentered.png
+.. |fig1.c| image:: https://raw.githubusercontent.com/charnley/rmsd/master/docs/figures/fig_rmsd_rotated.png
+
+.. |fig2.a| image:: https://raw.githubusercontent.com/charnley/rmsd/master/docs/figures/fig_reorder_problem.png
+.. |fig2.b| image:: https://raw.githubusercontent.com/charnley/rmsd/master/docs/figures/fig_reorder_inertia.png
+.. |fig2.c| image:: https://raw.githubusercontent.com/charnley/rmsd/master/docs/figures/fig_reorder_qml.png

docs/notebooks/coord_funcs.py

@@ -0,0 +1,29 @@
+import numpy as np
+from numpy import ndarray
+
+
+def get_rotation_matrix(degrees: float) -> ndarray:
+    """https://en.wikipedia.org/wiki/Rotation_matrix"""
+
+    radians = degrees * np.pi / 180.0
+
+    r11 = np.cos(radians)
+    r12 = -np.sin(radians)
+    r21 = np.sin(radians)
+    r22 = np.cos(radians)
+
+    R = np.array([[r11, r12], [r21, r22]])
+
+    return R
+
+
+def degree2radiant(degrees):
+    return degrees * np.pi / 180.0
+
+
+def rotate_coord(angle, coord, axis=[0, 1]):
+    U = get_rotation_matrix(angle)
+    _xy = np.dot(coord[:, axis], U)
+    _coord = np.array(coord, copy=True)
+    _coord[:, axis] = _xy
+    return _coord