Add inline and Tutorial documentation for internal_coords; change edr…

…a .aks to .atomkeys; improve residue.internal_coord init; rebuild-test -quick better matches normal; more tests; typos; NEWS entry

Bio/PDB/PICIO.py

@@ -1,4 +1,4 @@
-# Copyright 2019-2021 by Robert T. Miller.  All rights reserved.
+# Copyright 2019-2022 by Robert T. Miller.  All rights reserved.
 # This file is part of the Biopython distribution and governed by your
 # choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
 # Please see the LICENSE file that should have been included as part of this
@@ -217,38 +217,41 @@ def default_hedron(ek: Tuple, ric: IC_Residue) -> None:
         Adds Hedron to current Chain.internal_coord, see ic_data for default
         values and reference database source.
         """
-        aks = []
+        atomkeys = []
         hkey = None
 
         atmNdx = AtomKey.fields.atm
         resNdx = AtomKey.fields.resname
         resPos = AtomKey.fields.respos
-        aks = [ek[i].akl for i in range(3)]
+        atomkeys = [ek[i].akl for i in range(3)]
 
-        atpl = tuple([aks[i][atmNdx] for i in range(3)])
-        res = aks[0][resNdx]
+        atpl = tuple([atomkeys[i][atmNdx] for i in range(3)])
+        res = atomkeys[0][resNdx]
         if (
-            aks[0][resPos] != aks[2][resPos]  # hedra crosses amide bond so not reversed
+            atomkeys[0][resPos]
+            != atomkeys[2][resPos]  # hedra crosses amide bond so not reversed
             or atpl == ("N", "CA", "C")  # or chain start tau
             or atpl in ic_data_backbone  # or found forward hedron in ic_data
             or (res not in ["A", "G"] and atpl in ic_data_sidechains[res])
         ):
             hkey = ek
-            rhcl = [aks[i][resNdx] + aks[i][atmNdx] for i in range(3)]
+            rhcl = [atomkeys[i][resNdx] + atomkeys[i][atmNdx] for i in range(3)]
             try:
                 dflts = hedra_defaults["".join(rhcl)][0]
             except KeyError:
-                if aks[0][resPos] == aks[1][resPos]:
-                    rhcl = [aks[i][resNdx] + aks[i][atmNdx] for i in range(2)]
-                    rhc = "".join(rhcl) + "X" + aks[2][atmNdx]
+                if atomkeys[0][resPos] == atomkeys[1][resPos]:
+                    rhcl = [atomkeys[i][resNdx] + atomkeys[i][atmNdx] for i in range(2)]
+                    rhc = "".join(rhcl) + "X" + atomkeys[2][atmNdx]
                 else:
-                    rhcl = [aks[i][resNdx] + aks[i][atmNdx] for i in range(1, 3)]
-                    rhc = "X" + aks[0][atmNdx] + "".join(rhcl)
+                    rhcl = [
+                        atomkeys[i][resNdx] + atomkeys[i][atmNdx] for i in range(1, 3)
+                    ]
+                    rhc = "X" + atomkeys[0][atmNdx] + "".join(rhcl)
                 dflts = hedra_defaults[rhc][0]
         else:
             # must be reversed or fail
             hkey = ek[::-1]
-            rhcl = [aks[i][resNdx] + aks[i][atmNdx] for i in range(2, -1, -1)]
+            rhcl = [atomkeys[i][resNdx] + atomkeys[i][atmNdx] for i in range(2, -1, -1)]
             dflts = hedra_defaults["".join(rhcl)][0]
 
         process_hedron(
@@ -959,7 +962,7 @@ def write_PIC(
         * "classic",    # classic_b | chi
         * "hedra",      # all hedra including bond lengths
         * "primary",    # all primary dihedra
-        * "secondary",  # all secondary dihedra
+        * "secondary",  # all secondary dihedra (fixed angle from primary dihedra)
         * "all",        # hedra | primary | secondary
         * "initAtoms",  # XYZ coordinates of initial Tau (N-Ca-C)
         * "bFactors"

Bio/PDB/ic_data.py

@@ -1,4 +1,4 @@
-# Copyright 2019-21 by Robert T. Miller.  All rights reserved.
+# Copyright 2019-22 by Robert T. Miller.  All rights reserved.
 # This file is part of the Biopython distribution and governed by your
 # choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
 # Please see the LICENSE file that should have been included as part of this

Bio/PDB/ic_rebuild.py

@@ -1,4 +1,4 @@
-# Copyright 2019-2021 by Robert T. Miller.  All rights reserved.
+# Copyright 2019-2022 by Robert T. Miller.  All rights reserved.
 # This file is part of the Biopython distribution and governed by your
 # choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
 # Please see the LICENSE file that should have been included as part of this
@@ -169,7 +169,7 @@ def IC_duplicate(entity) -> Structure:
     Calls :meth:`.Chain.atom_to_internal_coordinates` if needed.
 
     :param Entity entity: Biopython PDB Entity (will fail for Atom)
-    :returns: Biopython PDBStructure, no Atom objects
+    :returns: Biopython PDBStructure, no Atom objects except initial coords
     """
     sp = StringIO()
     hasInternalCoords = False
@@ -354,10 +354,10 @@ def compare_residues(
         Whether to print mismatch info, default False
     :param bool quick: default False.
         Only check atomArrays are identical, aCoordMatchCount=0 if different
-    :param float rtol, atol: default 1e-03, 1e-05 or round to 3 places.
+    :param float rtol, atol: default 1e-03, 1e-03 or round to 3 places.
         Numpy allclose parameters; default is to round atom coordinates to 3
         places and test equal.  For 'quick' will use defaults above for
-        comparing ataomArrays
+        comparing atomArrays
     :returns dict:
         Result counts for Residues, Full ID match Residues, Atoms,
         Full ID match atoms, and Coordinate match atoms; report string;
@@ -380,30 +380,38 @@ def compare_residues(
 
     if quick:
         if isinstance(e0, Chain):
-            if numpy.allclose(
+            if e0.internal_coord.atomArray is not None and numpy.allclose(
                 e0.internal_coord.atomArray,
                 e1.internal_coord.atomArray,
                 rtol=1e-03 if rtol is None else rtol,
-                atol=1e-05 if atol is None else atol,
+                atol=1e-03 if atol is None else atol,
             ):
-                cmpdict["pass"] = True
                 cmpdict["aCoordMatchCount"] = numpy.size(e0.internal_coord.atomArray, 0)
+                if cmpdict["aCoordMatchCount"] > 0:
+                    cmpdict["pass"] = True
+                else:
+                    cmpdict["pass"] = False
             else:
                 cmpdict["pass"] = False
         else:
             cmpdict["pass"] = True
             for c0, c1 in zip_longest(e0.get_chains(), e1.get_chains()):
-                if numpy.allclose(
-                    c0.internal_coord.atomArray,
-                    c1.internal_coord.atomArray,
-                    rtol=1e-03 if rtol is None else rtol,
-                    atol=1e-05 if atol is None else atol,
-                ):
-                    cmpdict["aCoordMatchCount"] += numpy.size(
-                        c0.internal_coord.atomArray, 0
-                    )
-                else:
-                    cmpdict["pass"] = False
+                if c0.internal_coord.atomArray is not None:
+                    if numpy.allclose(
+                        c0.internal_coord.atomArray,
+                        c1.internal_coord.atomArray,
+                        rtol=1e-03 if rtol is None else rtol,
+                        atol=1e-03 if atol is None else atol,
+                    ):
+                        cmpdict["aCoordMatchCount"] += numpy.size(
+                            c0.internal_coord.atomArray, 0
+                        )
+                    else:
+                        cmpdict["pass"] = False
+                    cmpdict["aCount"] += numpy.size(c0.internal_coord.atomArray, 0)
+            if cmpdict["aCoordMatchCount"] < cmpdict["aCount"]:
+                cmpdict["pass"] = False
+
     else:
         for r0, r1 in zip_longest(e0.get_residues(), e1.get_residues()):
             if 2 == r0.is_disordered() == r1.is_disordered():