# Copyright 2020-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
# package.

"""Unit tests for internal_coords module of Bio.PDB."""

import copy
import re
import unittest
import warnings

try:
    import numpy as np  # noqa F401
except ImportError:
    from Bio import MissingPythonDependencyError

    raise MissingPythonDependencyError(
        "Install NumPy if you want to use Bio.PDB."
    ) from None

from io import StringIO

from Bio.File import as_handle
from Bio.PDB.ic_rebuild import compare_residues
from Bio.PDB.ic_rebuild import IC_duplicate
from Bio.PDB.ic_rebuild import structure_rebuild_test
from Bio.PDB.internal_coords import AtomKey
from Bio.PDB.internal_coords import Dihedron
from Bio.PDB.internal_coords import IC_Chain
from Bio.PDB.internal_coords import IC_Residue
from Bio.PDB.MMCIFParser import MMCIFParser
from Bio.PDB.mmtf import MMTFParser
from Bio.PDB.PDBExceptions import PDBConstructionWarning
from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.PICIO import read_PIC_seq
from Bio.PDB.SCADIO import write_SCAD
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord


class Rebuild(unittest.TestCase):
    """Read PDB and mmCIF structures, convert to/from internal coordinates."""

    PDB_parser = PDBParser(PERMISSIVE=True, QUIET=True)
    CIF_parser = MMCIFParser(QUIET=True)
    MMTF_parser = MMTFParser()
    pdb_1LCD = PDB_parser.get_structure("1LCD", "PDB/1LCD.pdb")
    # cif_1A7G = CIF_parser.get_structure("1A7G", "PDB/1A7G.cif")
    # cif_1A7G2 = CIF_parser.get_structure("1A7G", "PDB/1A7G.cif")
    pdb_2XHE = PDB_parser.get_structure("2XHE", "PDB/2XHE.pdb")
    pdb_2XHE2 = PDB_parser.get_structure("2XHE", "PDB/2XHE.pdb")
    pdb_1A8O = PDB_parser.get_structure("1A8O", "PDB/1A8O.pdb")
    cif_3JQH = CIF_parser.get_structure("3JQH", "PDB/3JQH.cif")
    cif_3JQH2 = CIF_parser.get_structure("3JQH", "PDB/3JQH.cif")
    cif_4CUP = CIF_parser.get_structure("4CUP", "PDB/4CUP.cif")
    cif_4CUP2 = CIF_parser.get_structure("4CUP", "PDB/4CUP.cif")
    cif_4ZHL = CIF_parser.get_structure("4ZHL", "PDB/4ZHL.cif")
    cif_4ZHL2 = CIF_parser.get_structure("4ZHL", "PDB/4ZHL.cif")
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always", PDBConstructionWarning)
        mmtf_1A8O = MMTF_parser.get_structure("PDB/1A8O.mmtf")

    def test_mmtf(self):
        chain = next(self.mmtf_1A8O.get_chains())
        ic_chain = IC_Chain(chain)
        self.assertEqual(len(ic_chain.ordered_aa_ic_list), 70)

    def test_rebuild_1a8o(self):
        """Duplicate tutorial doctests which fail under linux."""
        IC_Chain.MaxPeptideBond = 4.0
        r = structure_rebuild_test(self.pdb_1A8O, False)
        self.assertTrue(r["pass"])

        cic = list(self.pdb_1A8O.get_chains())[0].internal_coord
        distances = cic.distance_plot()
        chirality = cic.dihedral_signs()
        c2 = IC_duplicate(cic.chain)[0]["A"]
        cic2 = c2.internal_coord
        cic2.atomArray = np.zeros((cic2.AAsiz, 4), dtype=np.float64)
        cic2.dihedraAngle[:] = 0.0
        cic2.hedraAngle[:] = 0.0
        cic2.hedraL12[:] = 0.0
        cic2.hedraL23[:] = 0.0
        cic2.copy_initNCaCs(cic)
        cic2.distplot_to_dh_arrays(distances, chirality)
        cic2.distance_to_internal_coordinates()
        c2.internal_to_atom_coordinates()
        np.allclose(cic2.atomArray, cic.atomArray)

    def test_rebuild_multichain_missing(self):
        """Convert multichain missing atom struct to, from internal coords."""
        # 2XHE has regions of missing chain, last residue has only N
        r = structure_rebuild_test(self.pdb_2XHE, False)
        self.assertEqual(r["residues"], 787)
        self.assertEqual(r["rCount"], 835)
        self.assertEqual(r["rMatchCount"], 835)
        self.assertEqual(r["aCount"], 6267)
        self.assertEqual(r["disAtmCount"], 0)
        self.assertEqual(r["aCoordMatchCount"], 6267)
        self.assertEqual(len(r["chains"]), 2)
        self.assertTrue(r["pass"])

    def test_rebuild_disordered_atoms_residues(self):
        """Convert disordered protein to internal coordinates and back."""
        # 3jqh has both disordered residues
        # and disordered atoms in ordered residues

        with warnings.catch_warnings(record=True):
            warnings.simplefilter("always", PDBConstructionWarning)
            r = structure_rebuild_test(self.cif_3JQH, False)
        # print(r)
        self.assertEqual(r["residues"], 26)
        self.assertEqual(r["rCount"], 47)
        self.assertEqual(r["rMatchCount"], 47)
        self.assertEqual(r["aCount"], 217)
        self.assertEqual(r["disAtmCount"], 50)
        self.assertEqual(r["aCoordMatchCount"], 217)
        self.assertEqual(len(r["chains"]), 1)
        self.assertTrue(r["pass"])

        IC_Residue.no_altloc = True
        with warnings.catch_warnings(record=True):
            warnings.simplefilter("always", PDBConstructionWarning)
            r = structure_rebuild_test(self.cif_3JQH2, verbose=False, quick=True)
        self.assertEqual(r["aCoordMatchCount"], 167, msg="no_altloc fail")
        self.assertTrue(r["pass"], msg="no_altloc fail")
        IC_Residue.no_altloc = False

    def test_no_crosstalk(self):
        """Deep copy, change few internal coords, test nothing else changes."""
        # IC_Chain.ParallelAssembleResidues = False
        self.cif_4CUP.atom_to_internal_coordinates()
        cpy4cup = copy.deepcopy(self.cif_4CUP)
        cic0 = self.cif_4CUP.child_list[0].child_list[0].internal_coord
        cic1 = cpy4cup.child_list[0].child_list[0].internal_coord
        alist = [
            "omg",
            "phi",
            "psi",
            "chi1",
            "chi2",
            "chi3",
            "chi4",
            "chi5",
            "tau",
        ]
        delta = 33  # degrees to change
        tdelta = delta / 10.0  # more realistic for bond angle
        targPos = 1
        for ang in alist:
            # skip by 2's with alist along original chain changing angle spec
            ricTarg = cic0.chain.child_list[targPos].internal_coord
            # print(targPos + 1, ricTarg.lc, ang)
            targPos += 2
            try:
                edr = ricTarg.pick_angle(ang)
                andx = edr.ndx
                if ang == "tau":
                    cic0.hedraAngle[andx] += tdelta
                    cic0.hAtoms_needs_update[andx] = True
                    cic0.atomArrayValid[cic0.h2aa[andx]] = False
                    cic0.hAtoms_needs_update[:] = True
                    cic0.atomArrayValid[:] = False
                    cic0.dAtoms_needs_update[:] = True
                else:
                    cic0.dihedraAngle[andx] += delta
                    if cic0.dihedraAngle[andx] > 180.0:
                        cic0.dihedraAngle[andx] -= 360.0
                    cic0.dihedraAngleRads[andx] = np.deg2rad(cic0.dihedraAngle[andx])
                    cic0.dAtoms_needs_update[andx] = True
                    cic0.atomArrayValid[cic0.d2aa[andx]] = False
                    # test Dihedron.bits()
                    pfd = IC_Residue.picFlagsDict
                    if ricTarg.rbase[2] == "P" and ang == "omg":
                        self.assertEqual(edr.bits(), (pfd["omg"] | pfd["pomg"]))
                    else:
                        self.assertEqual(edr.bits(), pfd[ang])

            except AttributeError:
                pass  # skip if residue does not have e.g. chi5
        cic0.internal_to_atom_coordinates()  # move atoms
        cic0.atom_to_internal_coordinates()  # get new internal coords
        # generate hdelta and ddelta difference arrays so can look for what
        # changed
        hdelta = cic0.hedraAngle - cic1.hedraAngle
        hdelta[np.abs(hdelta) < 0.00001] = 0.0
        ddelta = cic0.dihedraAngle - cic1.dihedraAngle
        ddelta[np.abs(ddelta) < 0.00001] = 0.0
        ddelta[ddelta < -180.0] += 360.0  # wrap around circle values
        targPos = 1
        for ang in alist:
            # same skip along original chain looking at hdelta and ddelta
            # if change is as specified, set difference to 0 then we can test
            # for any remaining (spurious) changes
            ricTarg = cic0.chain.child_list[targPos].internal_coord
            # print(targPos + 1, ricTarg.lc, ang)
            targPos += 2
            try:
                andx = ricTarg.pick_angle(ang).ndx
                if ang == "tau":
                    self.assertAlmostEqual(hdelta[andx], tdelta, places=4)
                    hdelta[andx] = 0.0
                    # some other angle has to change to accommodate tau change
                    # N-Ca-Cb is artifact of choices in ic_data
                    # expected change so clear relevant hdelta here
                    adjAngNdx = ricTarg.pick_angle("N:CA:CB").ndx
                    self.assertNotAlmostEqual(hdelta[adjAngNdx], 0.0, places=1)
                    hdelta[adjAngNdx] = 0.0
                else:
                    self.assertAlmostEqual(ddelta[andx], delta, places=4)
                    ddelta[andx] = 0.0
            except AttributeError:
                pass  # if residue does not have e.g. chi5

        hsum = hdelta.sum()
        self.assertEqual(hsum, 0.0)
        dsum = ddelta.sum()
        self.assertEqual(dsum, 0.0)

        # test hedron len12, angle, len23 setters and getters
        hed = list(cic0.hedra.values())[10]
        val = hed.len12 + 0.5
        hed.len12 = val
        self.assertEqual(hed.len12, val)
        val = hed.len23 + 0.5
        hed.len23 = val
        self.assertEqual(hed.len23, val)
        val = hed.angle + 1
        hed.angle = val
        self.assertEqual(hed.angle, val)
        dihed = list(cic0.dihedra.values())[10]
        val = dihed.angle + 196
        dihed.angle = val
        if val > 180.0:
            val -= 360.0
        if val < -180.0:
            val += 360.0
        self.assertEqual(dihed.angle, val)

    def test_model_change_internal_coords(self):
        """Get model internal coords, modify psi and chi1 values and check."""
        mdl = self.pdb_1LCD[1]
        mdl.atom_to_internal_coordinates()
        # other tests show can build with arbitrary internal coords
        # build here so changes below trigger more complicated
        # Atoms_needs_update mask arrays
        mdl.internal_to_atom_coordinates()
        nvt = {}
        nvc1 = {}
        nvpsi = {}
        nvlen = {}
        nvlen2 = {}
        tcount = 0
        l2count = 0
        c1count = 0
        psicount = 0
        lcount = 0
        for r in mdl.get_residues():
            ric = r.internal_coord
            if ric:
                # hedra change
                tau = ric.get_angle("tau")
                if ric.rprev != [] and tau is not None:
                    tcount += 1
                    nv = tau + 0.5
                    ric.set_angle("tau", nv)
                    nvt[str(r)] = nv
                    l2count += 1
                    leng2 = ric.get_length("N:CA")
                    nv = leng2 + 0.05
                    ric.set_length("N:CA", nv)
                    nvlen2[str(r)] = nv

                # sidechain dihedron change
                chi1 = ric.get_angle("chi1")
                if chi1 is not None:
                    c1count += 1
                    nv = chi1 + 90
                    if nv > 180.0:
                        nv -= 360.0
                    # ric.set_angle("chi1", nv)
                    ric.bond_set("chi1", nv)
                    nvc1[str(r)] = nv
                # backbone dihedron change
                psi = ric.get_angle("psi")
                if psi is not None:
                    psicount += 1
                    nv = psi - 90
                    if nv < -180.0:
                        nv += 360.0
                    ric.set_angle("psi", nv)
                    nvpsi[str(r)] = nv
                leng = ric.get_length("CA:CB")
                if leng is not None:
                    lcount += 1
                    nv = leng + 0.05
                    ric.set_length("CA:CB", nv)
                    nvlen[str(r)] = nv

        mdl.internal_to_atom_coordinates()

        # prove not using stored results
        for chn in mdl.get_chains():
            if hasattr(chn, "hedraLen"):
                delattr(chn.internal_coord, "hedraLen")
                delattr(chn.internal_coord, "dihedraLen")
                delattr(chn.internal_coord, "hedraAngle")
                delattr(chn.internal_coord, "dihedraAngle")
                for r in chn.get_residues():
                    r.internal_coord.hedra = {}
                    r.internal_coord.dihedra = {}

        mdl.atom_to_internal_coordinates()
        ttcount = 0
        l2tcount = 0
        c1tcount = 0
        psitcount = 0
        ltcount = 0
        for r in mdl.get_residues():
            ric = r.internal_coord
            if ric:
                tau = ric.get_angle("tau")
                if ric.rprev != [] and tau is not None:
                    ttcount += 1
                    # print(str(r), "tau", tau, nvt[str(r)])
                    self.assertAlmostEqual(tau, nvt[str(r)], places=3)
                    l2tcount += 1
                    l2 = ric.get_length("N:CA")
                    self.assertAlmostEqual(l2, nvlen2[str(r)], places=3)
                chi1 = ric.get_angle("chi1")
                if chi1 is not None:
                    c1tcount += 1
                    # print(str(r), "chi1", chi1, nvc1[str(r)])
                    self.assertAlmostEqual(chi1, nvc1[str(r)], places=3)
                psi = ric.get_angle("psi")
                if psi is not None:
                    psitcount += 1
                    # print(str(r), "psi", psi, nvpsi[str(r)])
                    self.assertAlmostEqual(psi, nvpsi[str(r)], places=3)
                leng = ric.get_length("CA:CB")
                if leng is not None:
                    ltcount += 1
                    self.assertAlmostEqual(leng, nvlen[str(r)], places=3)

        self.assertEqual(tcount, ttcount)
        self.assertEqual(l2count, l2tcount)
        self.assertEqual(c1count, c1tcount)
        self.assertEqual(psicount, psitcount)
        self.assertEqual(lcount, ltcount)
        self.assertGreater(ttcount, 0)
        self.assertGreater(c1count, 0)
        self.assertGreater(psicount, 0)
        self.assertGreater(lcount, 0)

    def test_write_SCAD(self):
        """Check SCAD output plus MaxPeptideBond and Gly CB.

        SCAD tests: scaling, transform mtx, extra bond created (allBonds)
        """
        sf = StringIO()
        write_SCAD(
            self.cif_4CUP2,
            sf,
            10.0,
            pdbid="4cup",
            backboneOnly=True,
            includeCode=False,
        )
        sf.seek(0)
        next_one = False
        with as_handle(sf, mode="r") as handle:
            for aline in handle.readlines():
                if "// (1856_S_CB, 1856_S_CA, 1856_S_C)" in aline:
                    m = re.search(r"\[\s+(\d+\.\d+)\,", aline)
                    if m:
                        # test correctly scaled atom bond length
                        self.assertAlmostEqual(float(m.group(1)), 15.30582, places=3)
                    else:
                        self.fail("scaled atom bond length not found")
                elif '[ 1, "1857M",' in aline:
                    next_one = True
                elif next_one:
                    next_one = False
                    # test last residue transform looks roughly correct
                    # some differences due to sorting issues on different
                    # python versions
                    target = [-12.413, -3.303, 35.771, 1.0]
                    ms = re.findall(  # last column of each row
                        r"\s+(-?\d+\.\d+)\s+\]", aline
                    )
                    if ms:
                        for i in range(3):
                            self.assertAlmostEqual(float(ms[i]), target[i], places=0)
                    else:
                        self.fail("transform not found")

        sf.seek(0)
        IC_Residue.gly_Cbeta = True
        IC_Chain.MaxPeptideBond = 100.0
        chn = self.pdb_2XHE2[0]["A"]
        chn.atom_to_internal_coordinates()
        rt0 = chn.internal_coord.ordered_aa_ic_list[12]
        rt1 = chn.internal_coord.ordered_aa_ic_list[16]
        rt0.set_flexible()
        rt1.set_hbond()

        write_SCAD(
            self.pdb_2XHE2[0]["A"],
            sf,
            10.0,
            pdbid="2xhe",
            # maxPeptideBond=100.0,
            includeCode=False,
            start=10,
            fin=570,
        )
        sf.seek(0)
        allBondsPass = False
        maxPeptideBondPass = False
        glyCbetaFound = False
        startPass = True
        finPass = True
        flexPass = False
        hbPass = False
        with as_handle(sf, mode="r") as handle:
            for aline in handle.readlines():
                # test extra bond created in TRP (allBonds is True)
                if '"Cres", 0, 0, 1, 0, StdBond, "W", 24, "CD2CE3CZ3"' in aline:
                    allBondsPass = True
                # test 509_K-561_E long bond created
                if "509_K" in aline and "561_E" in aline:
                    maxPeptideBondPass = True
                if "(21_G_CB, 21_G_CA, 21_G_C)" in aline:
                    glyCbetaFound = True
                    target = [15.33630, 110.17513, 15.13861]
                    ms = re.findall(r"\s+(-?\d+\.\d+)", aline)
                    if ms:
                        for i in range(3):
                            self.assertAlmostEqual(float(ms[i]), target[i], places=0)
                    else:
                        self.fail("Cbeta internal coords not found")
                if "8_K_CA" in aline:
                    startPass = False
                if "572_N_CA" in aline:
                    finPass = False
                if 'FemaleJoinBond, FemaleJoinBond, "N", 13, "NCAC"' in aline:
                    flexPass = True
                if 'HBond, "R", 16, "CACO"' in aline:
                    hbPass = True

        self.assertTrue(allBondsPass, msg="missing extra ring close bonds")
        self.assertTrue(glyCbetaFound, msg="gly CB not created")
        self.assertTrue(maxPeptideBondPass, msg="ignored maxPeptideBond setting")
        self.assertTrue(startPass, msg="writeSCAD wrote residue before start")
        self.assertTrue(finPass, msg="writeSCAD wrote residue past fin")
        self.assertTrue(flexPass, msg="writeSCAD residue 12 not flexible")
        self.assertTrue(hbPass, msg="writeSCAD residue 16 no hbond")

    def test_i2a_start_fin(self):
        """Test assemble start/fin, default NCaC coordinates, IC_duplicate."""
        chn = self.pdb_1LCD[2]["A"]
        cpy = IC_duplicate(chn)[2]["A"]  # generates internal coords as needed
        cpy.internal_to_atom_coordinates(start=31, fin=45)
        cdict = compare_residues(chn, cpy, quick=True)
        self.assertFalse(cdict["pass"])
        # transform source coordinates to put res 31 tau at origin like
        # fragment
        res = chn[31]
        psi = res.internal_coord.pick_angle("psi")
        cst = np.transpose(psi.cst)
        chn.internal_coord.atomArray[:] = chn.internal_coord.atomArray.dot(cst)
        cdict = compare_residues(chn, cpy, rtol=1e-03, atol=1e-05)
        self.assertEqual(cdict["residues"], 51)
        self.assertEqual(cdict["rMatchCount"], 77)
        self.assertEqual(cdict["aCount"], 497)
        self.assertEqual(cdict["disAtmCount"], 0)
        self.assertEqual(cdict["aCoordMatchCount"], 140)
        self.assertEqual(cdict["aFullIdMatchCount"], 140)
        self.assertEqual(len(cdict["chains"]), 1)
        self.assertEqual(cdict["rCount"], 77)
        self.assertFalse(cdict["pass"])

    def test_distplot_rebuild(self):
        """Build identical structure from distplot and chirality data."""
        # load input chain
        for _chn1 in self.cif_4ZHL.get_chains():
            break
        # create atomArray and compute distplot and dihedral signs array
        _chn1.atom_to_internal_coordinates()
        _c1ic = _chn1.internal_coord
        atmNameNdx = AtomKey.fields.atm
        CaSelect = [
            _c1ic.atomArrayIndex.get(k)
            for k in _c1ic.atomArrayIndex.keys()
            if k.akl[atmNameNdx] == "CA"
        ]
        dplot0 = _chn1.internal_coord.distance_plot(filter=CaSelect)
        self.assertAlmostEqual(
            dplot0[3, 9],
            16.296,
            places=3,
            msg="fail generate distance plot with filter",
        )
        dplot1 = _chn1.internal_coord.distance_plot()
        dsigns = _chn1.internal_coord.dihedral_signs()

        # load second copy (same again) input chain
        for _chn2 in self.cif_4ZHL2.get_chains():
            break
        # create internal coord structures but do not compute di/hedra
        cic2 = _chn2.internal_coord = IC_Chain(_chn2)
        cic2.init_edra()
        # load relevant interatomic distances from chn1 distance plot
        cic2.distplot_to_dh_arrays(dplot1, dsigns)
        # compute di/hedra angles from dh_arrays
        cic2.distance_to_internal_coordinates()

        # clear chn2 atom coordinates
        # cic2.atomArrayValid[:] = False  # done in distance_to_internal_coordinates

        # initialize values but this is redundant to Valid=False above
        cic2.atomArray = np.zeros((cic2.AAsiz, 4), dtype=np.float64)
        cic2.atomArray[:, 3] = 1.0

        # 4zhl has chain breaks so copy initial coords of each segment
        cic2.copy_initNCaCs(_chn1.internal_coord)
        # compute chn2 atom coords from di/hedra data
        cic2.internal_to_atom_coordinates()

        # generate distance plot from second chain, confirm minimal distance
        # from original
        dp2 = cic2.distance_plot()
        dpdiff = np.abs(dplot1 - dp2)
        # print(np.amax(dpdiff))
        self.assertTrue(np.amax(dpdiff) < 0.000001)

    def test_seq_as_PIC(self):
        """Read seq as PIC data, extend chain, set each chi angle, check various."""
        pdb_structure = read_PIC_seq(
            SeqRecord(
                Seq("GAVLIMFPSTCNQYWDEHKR"),
                id="1RND",
                description="my 20aa sequence",
            )
        )
        chn = next(pdb_structure.get_chains())
        cic = chn.internal_coord
        cic.make_extended()
        for chi in range(1, 6):
            for ric in chn.internal_coord.ordered_aa_ic_list:
                targ = "chi" + str(chi)
                rchi = ric.get_angle(targ)
                if rchi is not None:
                    nangl = rchi + 60.0
                    # nangl = nangl if (nangl <= 180.0) else nangl - 360.0
                    ric.set_angle(targ, nangl)
                    # ric.bond_set("chi1", nangl)

        pdb_structure.internal_to_atom_coordinates()

        self.assertEqual(
            len(cic.atomArrayValid), 168, msg="wrong number atoms from sequence"
        )

        # test make_extended and each sequential chi rotation places selected
        # atom from each sidechain where expected.
        posnDict = {
            "1_G_CA": [0.000, 0.000, 0.000, 1.000],
            "2_A_CB": [-0.411, 2.505, 4.035, 1.000],
            "3_V_CG2": [-2.246, -2.942, 4.865, 1.000],
            "4_L_CD2": [-4.673, 2.124, 6.060, 1.000],
            "5_I_CD1": [-4.260, -4.000, 10.546, 1.000],
            "6_M_CE": [-10.035, 4.220, 12.322, 1.000],
            "7_F_CE2": [-5.835, -1.076, 15.391, 1.000],
            "8_P_CD": [-15.060, -2.368, 17.751, 1.000],
            "9_S_OG": [-12.733, 0.388, 24.504, 1.000],
            "10_T_CG2": [-19.084, -0.423, 22.650, 1.000],
            "11_C_SG": [-15.455, 2.166, 27.172, 1.000],
            "12_N_ND2": [-20.734, -3.990, 27.233, 1.000],
            "13_Q_NE2": [-19.362, 4.465, 30.847, 1.000],
            "14_Y_CE2": [-20.685, -3.571, 32.854, 1.000],
            "15_W_CH2": [-23.295, 1.706, 32.883, 1.000],
            "16_D_OD2": [-24.285, -3.683, 40.620, 1.000],
            "17_E_OE2": [-31.719, 2.385, 38.887, 1.000],
            "18_H_NE2": [-25.763, -0.071, 46.356, 1.000],
            "19_K_NZ": [-36.073, -0.926, 42.383, 1.000],
            "20_R_NH2": [-28.792, 3.687, 51.738, 1.000],
        }

        for k, v in posnDict.items():
            atm = AtomKey(k)
            ndx = cic.atomArrayIndex[atm]
            coord = np.round(cic.atomArray[ndx], 3)
            self.assertTrue(np.array_equal(coord, v), msg=f"position error on atom {k}")

        cic.update_dCoordSpace()
        rt = cic.ordered_aa_ic_list[5]  # pick a residue
        chi1 = rt.pick_angle("chi1")  # chi1 coord space puts CA at origin
        rt.applyMtx(chi1.cst)
        coord = rt.residue.child_dict["CA"].coord  # Biopython API Atom coords
        self.assertTrue(
            np.allclose(coord, [0.0, 0.0, 0.0]), msg="dCoordSpace transform error"
        )

        # test Dihedron repr and all that leads into it
        psi = rt.pick_angle("psi")

        self.assertEqual(
            psi.__repr__(),
            "4-6_M_N:6_M_CA:6_M_C:7_F_N MNMCAMCFN 123.0 ('1RND', 0, 'A', (' ', 6, ' '))",
            msg="dihedron __repr__ error for M6 psi",
        )

        m = "Edron rich comparison failed"
        self.assertTrue(chi1 != psi, msg=m)
        self.assertFalse(chi1 == psi, msg=m)
        self.assertTrue(psi < chi1, msg=m)
        self.assertTrue(psi <= chi1, msg=m)
        self.assertTrue(chi1 > psi, msg=m)
        self.assertTrue(chi1 >= psi, msg=m)

        self.assertTrue(
            chi1.cre_class == "NsbCsbCsbCsb",
            msg="covalent radii assignment error for chi1",
        )

        # dihedron atomkeys are all in residue atomkeys as expected, including
        # i+1 N for psi.
        self.assertTrue(all(ak in rt for ak in chi1.atomkeys))
        self.assertFalse(all(ak in rt for ak in psi.atomkeys))

        # test Hedron repr and all that leads into it
        tau = rt.pick_angle("tau")
        self.assertEqual(
            tau.__repr__(),
            "3-6_M_N:6_M_CA:6_M_C MNMCAMC 1.46091 110.97184 1.52499",
            msg="hedron __repr__ error for M11 tau",
        )

        # some specific AtomKey comparisons missed in other tests
        a0, a1 = tau.atomkeys[0], tau.atomkeys[1]
        m = "AtomKey rich comparison failed"
        self.assertTrue(a1 > a0, msg=m)
        self.assertTrue(a1 >= a0, msg=m)
        self.assertTrue(a0 <= a1, msg=m)

    def test_angle_fns(self):
        """Test angle_dif and angle_avg across +/-180 boundaries."""
        arr1 = np.array([179.0, 90.0, 88.0, 1.0])
        arr2 = np.array([-179.0, -90.0, -91.0, -1.0])
        assert (
            Dihedron.angle_dif(arr1, arr2) == np.array([2.0, 180.0, -179.0, -2.0])
        ).all()
        assert Dihedron.angle_avg(np.array([179.0, -179.0])) == 180.0
        assert Dihedron.angle_avg(np.array([1.0, -1.0])) == 0.0
        assert Dihedron.angle_avg(np.array([90.0, -90.0])) == 0.0
        assert Dihedron.angle_avg(np.array([91.0, -91.0])) == 180.0


if __name__ == "__main__":
    runner = unittest.TextTestRunner(verbosity=2)
    unittest.main(testRunner=runner)