Fixed issue 93 (#94)

* Fixed rmsd naming convention consistency
* Update Makefile for faster environment setup
* Fixed reorder inconsistency (Issue 93), added test
* Used github instead of gitlab for pre-commit

.pre-commit-config.yaml

@@ -54,7 +54,7 @@ repos:
         args: ["--line-length=99"]
 
 
-  - repo: https://gitlab.com/pycqa/flake8
+  - repo: https://github.com/pycqa/flake8
     rev: 4.0.1
     hooks:
       - id: flake8

Makefile

@@ -1,6 +1,10 @@
 
 env:
 	conda env create -f ./environment.yml -p ./env --quiet
+	./env/bin/python -m pip install -e .
+
+env-qml:
+	./env/bin/pip install git+https://github.com/qmlcode/qml@develop
 
 test:
 	python -m pytest -vrs tests

rmsd/__init__.py

@@ -1,34 +1,45 @@
 # flake8: noqa
-from . import calculate_rmsd
 from .calculate_rmsd import *
 from .calculate_rmsd import __doc__, __version__
 
 __all__ = [
-    "brute_permutation",
-    "centroid",
-    "check_reflections",
-    "generate_permutations",
-    "get_coordinates",
-    "get_coordinates_pdb",
-    "get_coordinates_xyz",
-    "hungarian",
+    "str_atom",
     "int_atom",
-    "kabsch",
+    "rmsd",
     "kabsch_rmsd",
     "kabsch_rotate",
     "kabsch_fit",
+    "kabsch",
+    "kabsch_weighted",
+    "kabsch_weighted_fit",
+    "kabsch_weighted_rmsd",
     "quaternion_rmsd",
-    "quaternion_rotate",
     "quaternion_transform",
-    "reorder_brute",
+    "makeW",
+    "makeQ",
+    "quaternion_rotate",
+    "centroid",
+    "hungarian_vectors",
+    "reorder_similarity",
     "reorder_distance",
+    "hungarian",
     "reorder_hungarian",
     "reorder_inertia_hungarian",
-    "reorder_similarity",
-    "rmsd",
+    "generate_permutations",
+    "brute_permutation",
+    "reorder_brute",
+    "check_reflections",
+    "rotation_matrix_vectors",
+    "get_cm",
+    "get_inertia_tensor",
+    "get_principal_axis",
     "set_coordinates",
-    "str_atom",
+    "get_coordinates",
+    "get_coordinates_pdb",
+    "get_coordinates_xyz_lines",
+    "get_coordinates_xyz",
+    "main",
 ]
 
 if __name__ == "__main__":
-    calculate_rmsd.main()  # pragma: no cover
+    main()  # pragma: no cover

rmsd/calculate_rmsd.py

@@ -2017,10 +2017,14 @@ def main(args: Optional[List[str]] = None) -> None:
     elif settings.reorder:
 
         assert reorder_method is not None, "Cannot reorder without selecting --reorder method"
-
         q_review = reorder_method(p_atoms, q_atoms, p_coord, q_coord)
-        q_coord = q_coord[q_review]
+
+    # If there is a reorder, then apply before print
+    if q_review is not None:
+
+        q_all_atoms = q_all_atoms[q_review]
         q_atoms = q_atoms[q_review]
+        q_coord = q_coord[q_review]
 
         assert all(
             p_atoms == q_atoms
@@ -2029,11 +2033,6 @@ def main(args: Optional[List[str]] = None) -> None:
     # print result
     if settings.output:
 
-        # print(q_coord)
-        # print(q_swap)
-        # print(q_reflection)
-        # print(q_review)
-
         if q_swap is not None:
             q_coord = q_coord[:, q_swap]
 
@@ -2042,10 +2041,6 @@ def main(args: Optional[List[str]] = None) -> None:
 
         q_coord -= centroid(q_coord)
 
-        if q_review is not None:
-            q_coord = q_coord[q_review]
-            q_all_atoms = q_all_atoms[q_review]
-
         # Rotate q coordinates
         # TODO Should actually follow rotation method
         q_coord = kabsch_rotate(q_coord, p_coord)
@@ -2057,8 +2052,6 @@ def main(args: Optional[List[str]] = None) -> None:
         xyz = set_coordinates(q_all_atoms, q_coord, title=f"{settings.structure_b} - modified")
         print(xyz)
 
-        # print(kabsch_rmsd(q_coord, p_coord, translate=True))
-        # assert False
     else:
 
         if not result_rmsd:

tests/context.py

@@ -1,3 +1,21 @@
+import subprocess
 from pathlib import Path
+from typing import List
 
 RESOURCE_PATH = Path("tests/resources")
+
+
+def call_main(args: List[str]) -> List[str]:
+
+    root_path = Path("./")
+    filename = root_path / "rmsd/calculate_rmsd.py"
+
+    cmd = ["python", f"{filename}", *args]
+
+    proc = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
+    stdout, stderr = proc.communicate()
+
+    if stderr is not None:
+        print(stderr.decode())
+
+    return stdout.decode().strip().split("\n")

tests/resources/issue93/a.xyz

@@ -0,0 +1,22 @@
+20
+
+C       -0.4378260222           -0.3027444364           -1.2372902933
+C        0.1107709466            1.1031415940           -1.5729191901
+C        0.5005103668           -1.1044030188           -0.2962381736
+C        0.3998133345           -0.7811385983            1.1631226837
+C        0.0885767200           -1.5765750897            2.2490552567
+C        0.4849378650            0.3701285431            3.0522986790
+H        0.0195747479           -1.5423304575           -2.8383838780
+H       -1.5595867022           -1.6043371753           -2.4047214628
+H       -1.3882210499           -0.1233250051           -0.7119415800
+H        1.5378933274           -0.9934345286           -0.6538216533
+H        0.2377875820           -2.1655077189           -0.4019180847
+H       -0.1633586912           -2.6334567466            2.2467335441
+H        0.6168935065            1.2177913407            3.7166960247
+H        0.8625703727            1.2817813256            1.1335155509
+H       -0.6039973796            0.7377874297           -3.2455391704
+N       -0.7633488058           -0.9817332422           -2.5064360775
+N        0.6557050530            0.4639055355            1.7015091026
+N        0.1451038795           -0.8519411439            3.4200535140
+O        0.7148263821            1.8041376084           -0.7858085661
+O       -0.1586606797            1.5088838316           -2.8169269008 

tests/resources/issue93/b.xyz

@@ -0,0 +1,22 @@
+20
+
+C       -0.2769166422           -0.7099045701           -1.2865831838
+C       -0.0318894674            0.6967020886           -1.8794695568
+C        0.6647089780           -1.0315870523           -0.0954617436
+C        0.2611763277           -0.4666702466            1.2320987930
+C       -0.0191751268           -1.0867015822            2.4345274135
+C       -0.2291714337            1.0136818561            2.8036493128
+H       -0.8594037813           -2.4661239083           -2.2309783928
+H        0.7080595367           -2.0732382850           -2.5032982174
+H       -1.3132498520           -0.6901271514           -0.9163939649
+H        1.6889557061           -0.7206820514           -0.3611372211
+H        0.6818980659           -2.1236321622            0.0197593401
+H        0.2183924342            1.5973044617            0.7751919144
+H       -0.0088301381           -2.1515499995            2.6500889072
+H       -0.4092099932            1.9792682002            3.2647218979
+H       -0.3655783827           -0.1675785054           -3.4867813105
+N       -0.2277128834           -1.6869155940           -2.3914365186
+N        0.1251882275            0.8821440432            1.4922681597
+N       -0.3220367191           -0.1577068300            3.4064874889
+O       -0.1893763020            0.7629926388           -3.2044168223
+O        0.2560010284            1.6738952911           -1.2175802549

tests/resources/issue93/c.xyz

@@ -0,0 +1,22 @@
+20
+b.xyz - modified
+6      -0.40759914     -0.21976042     -1.24619800
+6      -0.00056312      1.26286681     -1.40816417
+6       0.58710277     -1.02255927     -0.36586236
+6       0.41909751     -0.87621203      1.11560424
+6       0.16739220     -1.81528006      2.09726918
+6       0.33911177      0.05972314      3.12041893
+1      -1.36062251     -1.48928504     -2.58617323
+1       0.21329021     -1.22102654     -2.95535420
+1      -1.38339903     -0.19684427     -0.73756537
+1       1.61529960     -0.76790502     -0.67292786
+1       0.43914697     -2.08630272     -0.59541156
+1       0.36609516      0.83692342      3.87711598
+1       0.66401664      1.21334374      1.32621590
+1       0.02744709     -2.88534987      1.97140463
+1      -0.63332683      1.01562723     -3.13450799
+7      -0.63109905     -0.78238518     -2.59199650
+7       0.53063943      0.32047774      1.79455405
+7       0.12049476     -1.22387595      3.34125656
+8      -0.28077247      1.77523394     -2.60968173
+8       0.50821277      1.92922041     -0.52895719

tests/resources/issue93/d.xyz

@@ -0,0 +1,22 @@
+20
+b.xyz - modified
+6      -0.43789836     -0.30253334     -1.23717746
+6       0.11094785      1.10318590     -1.57305854
+6       0.50063992     -1.10429331     -0.29630787
+6       0.39992447     -0.78124662      1.16310428
+6       0.08894976     -1.57679109      2.24904815
+6       0.48469770      0.37004052      3.05228551
+1       0.01882276     -1.54269123     -2.83801055
+1      -1.56033536     -1.60373840     -2.40443087
+1      -1.38810094     -0.12293143     -0.71155788
+1       1.53799278     -0.99306362     -0.65387426
+1       0.23808480     -2.16541809     -0.40219882
+1      -0.16267214     -2.63374746      2.24668396
+1       0.61636707      1.21773342      3.71669406
+1       0.86185680      1.28183964      1.13344692
+1      -0.60323400      0.73728341     -3.24573304
+7      -0.76377924     -0.98157261     -2.50618485
+7       0.65534023      0.46390036      1.70148391
+7       0.14534613     -0.85216749      3.42003187
+8       0.71471651      1.80437288     -0.78585771
+8      -0.15770200      1.50846860     -2.81734749 

tests/test_centroid.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-import rmsd
+import rmsd as rmsdlib
 
 
 def test_centroid() -> None:
@@ -9,7 +9,7 @@ def test_centroid() -> None:
     a3 = np.array([-22.018, 17.551, 26.0], dtype=float)
 
     atms = np.asarray([a1, a2, a3])
-    centroid = rmsd.centroid(atms)
+    centroid = rmsdlib.centroid(atms)
 
     assert 3 == len(centroid)
 

tests/test_file_reading.py

@@ -5,13 +5,13 @@
 import pytest
 from context import RESOURCE_PATH
 
-import rmsd
+import rmsd as rmsdlib
 
 
 def test_get_coordinates_pdb_hetatm() -> None:
 
     filename = Path(RESOURCE_PATH) / "issue88" / "native.pdb"
-    atoms, _ = rmsd.get_coordinates_pdb(filename)
+    atoms, _ = rmsdlib.get_coordinates_pdb(filename)
 
     assert len(atoms)
     assert atoms[0] == "C"
@@ -21,35 +21,35 @@ def test_get_coordinates_pdb_hetatm() -> None:
 def test_get_coordinates_pdb() -> None:
 
     filename = RESOURCE_PATH / "ci2_1.pdb"
-    atoms, coords = rmsd.get_coordinates_pdb(filename)
+    atoms, coords = rmsdlib.get_coordinates_pdb(filename)
     assert "N" == atoms[0]
     assert [-7.173, -13.891, -6.266] == coords[0].tolist()
 
 
 def test_get_coordinates_wrong() -> None:
     filename = RESOURCE_PATH / "ci2_1.pdb"
     with pytest.raises(ValueError):
-        _, _ = rmsd.get_coordinates(filename, "qqq")
+        _, _ = rmsdlib.get_coordinates(filename, "qqq")
 
 
 def test_get_coordinates_xyz() -> None:
 
     filename = RESOURCE_PATH / "ethane.xyz"
-    atoms, coords = rmsd.get_coordinates_xyz(filename)
+    atoms, coords = rmsdlib.get_coordinates_xyz(filename)
 
     assert "C" == atoms[0]
     assert [-0.98353, 1.81095, -0.0314] == coords[0].tolist()
 
 
 def test_get_coordinates_xyz_bad() -> None:
     filename = RESOURCE_PATH / "bad_format.xyz"
-    _, _ = rmsd.get_coordinates(filename, "xyz")
+    _, _ = rmsdlib.get_coordinates(filename, "xyz")
 
 
 def test_get_coordinates() -> None:
 
     filename = RESOURCE_PATH / "ethane.xyz"
-    atoms, coords = rmsd.get_coordinates(filename, "xyz")
+    atoms, coords = rmsdlib.get_coordinates(filename, "xyz")
 
     assert "C" == atoms[0]
     assert [-0.98353, 1.81095, -0.0314] == coords[0].tolist()
@@ -68,7 +68,7 @@ def test_get_coordinates_gzip(tmp_path: Path) -> None:
     with gzip.open(filename_gzip, "wb") as f:  # type: ignore
         f.write(content_byte)  # type: ignore
 
-    atoms, coords = rmsd.get_coordinates(
+    atoms, coords = rmsdlib.get_coordinates(
         filename_gzip, "pdb", is_gzip=True, return_atoms_as_int=True
     )
     assert 7 == atoms[0]
@@ -88,7 +88,7 @@ def test_get_coordinates_gzip_pdb(tmp_path: Path) -> None:
     with gzip.open(filename_gzip, "wb") as f:  # type: ignore
         f.write(content_byte)  # type: ignore
 
-    atoms, coords = rmsd.get_coordinates(filename_gzip, "xyz", is_gzip=True)
+    atoms, coords = rmsdlib.get_coordinates(filename_gzip, "xyz", is_gzip=True)
 
     assert "C" == atoms[0]
     assert [-0.98353, 1.81095, -0.0314] == coords[0].tolist()
@@ -99,10 +99,10 @@ def test_rmsd_pdb() -> None:
     filename_1 = RESOURCE_PATH / "ci2_1.pdb"
     filename_2 = RESOURCE_PATH / "ci2_2.pdb"
 
-    _, p_coord = rmsd.get_coordinates_pdb(filename_1)
-    _, q_coord = rmsd.get_coordinates_pdb(filename_2)
+    _, p_coord = rmsdlib.get_coordinates_pdb(filename_1)
+    _, q_coord = rmsdlib.get_coordinates_pdb(filename_2)
 
-    pure_rmsd = rmsd.rmsd(p_coord, q_coord)
+    pure_rmsd = rmsdlib.rmsd(p_coord, q_coord)
 
     np.testing.assert_almost_equal(26.9750, pure_rmsd, decimal=3)
 
@@ -112,9 +112,9 @@ def test_rmsd_xyz() -> None:
     filename_1 = RESOURCE_PATH / "ethane.xyz"
     filename_2 = RESOURCE_PATH / "ethane_mini.xyz"
 
-    _, p_coord = rmsd.get_coordinates_xyz(filename_1)
-    _, q_coord = rmsd.get_coordinates_xyz(filename_2)
+    _, p_coord = rmsdlib.get_coordinates_xyz(filename_1)
+    _, q_coord = rmsdlib.get_coordinates_xyz(filename_2)
 
-    pure_rmsd = rmsd.rmsd(p_coord, q_coord)
+    pure_rmsd = rmsdlib.rmsd(p_coord, q_coord)
 
     np.testing.assert_almost_equal(0.33512, pure_rmsd, decimal=3)