pdbtbx 0.12.0

#![allow(dead_code)]

use std::collections::HashMap;
use std::fmt;

use doc_cfg::doc_cfg;
#[cfg(feature = "rayon")]
use rayon::prelude::*;

use crate::structs::hierarchy::*;
use crate::transformation::TransformationMatrix;
use crate::{reference_tables, PDBError};
use crate::{structs::*, Context};

#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[derive(Debug, Clone, PartialEq)]
/// A PDB struct is generated by opening a PDB or mmCIF file. It contains
/// all information present in this file, like its atoms, bonds, hierarchy
/// , and metadata. The struct can be used to access, interact with, and
/// edit this data.
///
/// ```rust
/// use pdbtbx::*;
/// let (mut pdb, _errors) = pdbtbx::open("example-pdbs/1ubq.pdb").unwrap();
///
/// pdb.remove_atoms_by(|atom| atom.element() == Some(&Element::H)); // Remove all H atoms
///
/// let mut avg_b_factor = 0.0;
/// for atom in pdb.atoms() { // Iterate over all atoms in the structure
///     avg_b_factor += atom.b_factor();
/// }
/// avg_b_factor /= pdb.atom_count() as f64;
///
/// println!("The average B factor of the protein is: {}", avg_b_factor);
/// pdbtbx::save(&pdb, "dump/1ubq_no_hydrogens.pdb", pdbtbx::StrictnessLevel::Loose);
/// ```
pub struct PDB {
    /// The identifier as posed in the PDB Header or mmCIF entry.id, normally a 4 char string like '1UBQ'.
    pub identifier: Option<String>,
    /// The remarks above the PDB file, containing the remark-type-number and a line of free text.
    remarks: Vec<(usize, String)>,
    /// The Scale needed to transform orthogonal coordinates to fractional coordinates. This is inversely related to the unit cell.
    pub scale: Option<TransformationMatrix>,
    /// The OrigX needed to transform orthogonal coordinates to submitted coordinates. In normal cases this is equal to the identity transformation.
    pub origx: Option<TransformationMatrix>,
    /// The MtriXs needed to transform the Models to the full asymmetric subunit, if needed to contain the non-crystallographic symmetry.
    mtrix: Vec<MtriX>,
    /// The unit cell of the crystal, containing its size and shape. This is the size and shape of the repeating element in the crystal.
    pub unit_cell: Option<UnitCell>,
    /// The Symmetry or space group of the crystal. This is the way in which the protein is placed inside the unit cell.
    pub symmetry: Option<Symmetry>,
    /// The Models making up this PDB, containing all chain, residues, conformers, and atoms.
    models: Vec<Model>,
    /// Bonds in this PDB.
    bonds: Vec<(usize, usize, Bond)>,
}

/// # Creators
/// Creator functions for a PDB file
impl PDB {
    /// Create an empty PDB struct.
    pub const fn new() -> PDB {
        PDB {
            identifier: None,
            remarks: Vec::new(),
            scale: None,
            origx: None,
            mtrix: Vec::new(),
            unit_cell: None,
            symmetry: None,
            models: Vec::new(),
            bonds: Vec::new(),
        }
    }
}

/// # Remarks
/// Functionality for working with remarks.
impl PDB {
    /// Get the number of remark records in the PDB file.
    pub fn remark_count(&self) -> usize {
        self.remarks.len()
    }

    /// Get an iterator of references to the remarks, containing the remark-type-number and a line of free text.
    pub fn remarks(&self) -> impl DoubleEndedIterator<Item = &(usize, String)> + '_ {
        self.remarks.iter()
    }

    /// Get a parallel iterator of references to the remarks, containing the remark-type-number and a line of free text.
    #[doc_cfg(feature = "rayon")]
    pub fn par_remarks(&self) -> impl ParallelIterator<Item = &(usize, String)> + '_ {
        self.remarks.par_iter()
    }

    /// Get an iterator of references to the remarks, containing the remark-type-number and a line of free text.
    pub fn remarks_mut(&mut self) -> impl DoubleEndedIterator<Item = &mut (usize, String)> + '_ {
        self.remarks.iter_mut()
    }

    /// Get a parallel iterator of references to the remarks, containing the remark-type-number and a line of free text.
    #[doc_cfg(feature = "rayon")]
    pub fn par_remarks_mut(&mut self) -> impl ParallelIterator<Item = &mut (usize, String)> + '_ {
        self.remarks.par_iter_mut()
    }

    /// Add a remark
    ///
    /// ## Arguments
    /// * `remark_type` - the remark-type-number
    /// * `remark_text` - the free line of text, containing the actual remark
    ///
    /// ## Fails
    /// It fails if the text if too long, the text contains invalid characters or the remark-type-number is not valid (wwPDB v3.30).
    pub fn add_remark(&mut self, remark_type: usize, remark_text: String) -> Result<(), PDBError> {
        let context = Context::show(format!("REMARK {remark_type:3} {remark_text}"));
        if !reference_tables::valid_remark_type_number(remark_type) {
            return Err(PDBError::new(
                crate::ErrorLevel::InvalidatingError,
                "Remark-type-number invalid",
                "The given remark-type-number is not valid, see wwPDB v3.30 for valid remark-type-numbers",
                context));
        }
        if !valid_text(&remark_text) {
            return Err(PDBError::new(
                crate::ErrorLevel::InvalidatingError,
                "Remark text invalid",
                "The given remark text contains invalid characters.",
                context,
            ));
        }

        // As the text can only contain ASCII len() on strings is fine (it returns the length in bytes)
        let res = if remark_text.len() > 70 {
            Err(PDBError::new(
                crate::ErrorLevel::LooseWarning,
                "Remark text too long",
                format!("The given remark text is too long, the maximal length is 68 characters, the given string is {} characters.", remark_text.len()),
                context))
        } else {
            Ok(())
        };

        self.remarks.push((remark_type, remark_text));
        res
    }

    /// Delete the remarks matching the given predicate.
    pub fn delete_remarks_by<F>(&mut self, predicate: F)
    where
        F: Fn(&(usize, String)) -> bool,
    {
        self.remarks.retain(|r| !predicate(r));
    }
}

/// # MtriX
/// Functionality for working with the MtriX records form the PDB. The MtriX are needed
/// to transform the Models to the full asymmetric subunit, if needed to contain the
/// non-crystallographic symmetry.
impl PDB {
    /// Get an iterator of references to the MtriX records for this PDB.
    pub fn mtrix(&self) -> impl DoubleEndedIterator<Item = &MtriX> + '_ {
        self.mtrix.iter()
    }

    /// Get a parallel iterator of references to the MtriX records for this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_mtrix(&self) -> impl ParallelIterator<Item = &MtriX> + '_ {
        self.mtrix.par_iter()
    }

    /// Get an iterator of mutable references to the MtriX records for this PDB.
    pub fn mtrix_mut(&mut self) -> impl DoubleEndedIterator<Item = &mut MtriX> + '_ {
        self.mtrix.iter_mut()
    }

    /// Get a parallel iterator of mutable references to the MtriX records for this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_mtrix_mut(&mut self) -> impl ParallelIterator<Item = &mut MtriX> + '_ {
        self.mtrix.par_iter_mut()
    }

    /// Add a MtriX to this PDB.
    pub fn add_mtrix(&mut self, mtrix: MtriX) {
        self.mtrix.push(mtrix);
    }

    /// Delete the MtriX matching the given predicate.
    pub fn delete_mtrix_by<F>(&mut self, predicate: F)
    where
        F: Fn(&MtriX) -> bool,
    {
        self.mtrix.retain(|m| !predicate(m));
    }
}

impl<'a> PDB {
    /// Adds a Model to this PDB.
    pub fn add_model(&mut self, new_model: Model) {
        self.models.push(new_model);
    }

    /// Get the number of Models making up this PDB.
    pub fn model_count(&self) -> usize {
        self.models.len()
    }

    /// Get the number of Chains making up this PDB.
    pub fn chain_count(&self) -> usize {
        if self.models.is_empty() {
            0
        } else {
            self.models[0].chain_count()
        }
    }

    /// Get the number of Residues making up this PDB.
    pub fn residue_count(&self) -> usize {
        if self.models.is_empty() {
            0
        } else {
            self.models[0].residue_count()
        }
    }

    /// Get the number of Residues making up this PDB in parallel.
    #[doc_cfg(feature = "rayon")]
    pub fn par_residue_count(&self) -> usize {
        if self.models.is_empty() {
            0
        } else {
            self.models[0].par_residue_count()
        }
    }

    /// Get the number of Conformers making up this PDB.
    pub fn conformer_count(&self) -> usize {
        if self.models.is_empty() {
            0
        } else {
            self.models[0].conformer_count()
        }
    }

    /// Get the number of Conformers making up this PDB in parallel.
    #[doc_cfg(feature = "rayon")]
    pub fn par_conformer_count(&self) -> usize {
        if self.models.is_empty() {
            0
        } else {
            self.models[0].par_conformer_count()
        }
    }

    /// Get the number of Atoms making up this PDB.
    pub fn atom_count(&self) -> usize {
        if self.models.is_empty() {
            0
        } else {
            self.models[0].atom_count()
        }
    }

    /// Get the number of Atoms making up this PDB in parallel.
    #[doc_cfg(feature = "rayon")]
    pub fn par_atom_count(&self) -> usize {
        if self.models.is_empty() {
            0
        } else {
            self.models[0].par_atom_count()
        }
    }

    /// Get the number of Chains making up this PDB. Includes all models.
    pub fn total_chain_count(&self) -> usize {
        self.models
            .iter()
            .fold(0, |acc, item| acc + item.chain_count())
    }

    /// Get the number of Chains making up this PDB in parallel. Includes all models.
    #[doc_cfg(feature = "rayon")]
    pub fn par_total_chain_count(&self) -> usize {
        self.models.par_iter().map(Model::chain_count).sum()
    }

    /// Get the number of Residues making up this PDB. Includes all models.
    pub fn total_residue_count(&self) -> usize {
        self.models
            .iter()
            .fold(0, |acc, item| acc + item.residue_count())
    }

    /// Get the number of Residues making up this PDB in parallel. Includes all models.
    #[doc_cfg(feature = "rayon")]
    pub fn par_total_residue_count(&self) -> usize {
        self.models.par_iter().map(Model::par_residue_count).sum()
    }

    /// Get the number of Conformer making up this PDB. Includes all models.
    pub fn total_conformer_count(&self) -> usize {
        self.models
            .iter()
            .fold(0, |acc, item| acc + item.conformer_count())
    }

    /// Get the number of Conformer making up this PDB in parallel. Includes all models.
    #[doc_cfg(feature = "rayon")]
    pub fn par_total_conformer_count(&self) -> usize {
        self.models.par_iter().map(Model::par_conformer_count).sum()
    }

    /// Get the number of Atoms making up this PDB. Includes all models.
    pub fn total_atom_count(&self) -> usize {
        self.models
            .iter()
            .fold(0, |acc, item| acc + item.atom_count())
    }

    /// Get the number of Atoms making up this PDB in parallel. Includes all models.
    #[doc_cfg(feature = "rayon")]
    pub fn par_total_atom_count(&self) -> usize {
        self.models.par_iter().map(Model::par_atom_count).sum()
    }

    /// Get a reference to a specific Model from the list of Models making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Model
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn model(&self, index: usize) -> Option<&Model> {
        self.models.get(index)
    }

    /// Get a mutable reference to a specific Model from the list of Models making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Model
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn model_mut(&mut self, index: usize) -> Option<&mut Model> {
        self.models.get_mut(index)
    }

    /// Get a reference to a specific Chain from the list of Chains making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Chain
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn chain(&self, index: usize) -> Option<&Chain> {
        self.chains().nth(index)
    }

    /// Get a mutable reference to a specific Chain from the list of Chains making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Chain
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn chain_mut(&mut self, index: usize) -> Option<&mut Chain> {
        self.chains_mut().nth(index)
    }

    /// Get a reference to a specific Residue from the Residues making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Residue
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn residue(&self, index: usize) -> Option<&Residue> {
        self.residues().nth(index)
    }

    /// Get a mutable reference to a specific Residue from the Residues making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Residue
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn residue_mut(&mut self, index: usize) -> Option<&mut Residue> {
        self.residues_mut().nth(index)
    }

    /// Get a reference to a specific Conformer from the Conformers making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Conformer
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn conformer(&self, index: usize) -> Option<&Conformer> {
        self.conformers().nth(index)
    }

    /// Get a mutable reference to a specific Conformer from the Conformers making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Conformer
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn conformer_mut(&mut self, index: usize) -> Option<&mut Conformer> {
        self.conformers_mut().nth(index)
    }

    /// Get a reference to a specific Atom from the Atoms making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Atom
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn atom(&self, index: usize) -> Option<&Atom> {
        self.atoms().nth(index)
    }

    /// Get a mutable reference to a specific Atom from the Atoms making up this PDB.
    ///
    /// ## Arguments
    /// * `index` - the index of the Atom
    ///
    /// ## Fails
    /// It fails and returns `None` when the index is outside bounds.
    pub fn atom_mut(&mut self, index: usize) -> Option<&mut Atom> {
        self.atoms_mut().nth(index)
    }

    /// Get a reference to the specified atom. Its uniqueness is guaranteed by including the
    /// `insertion_code`, with its full hierarchy. The algorithm is based
    /// on binary search so it is faster than an exhaustive search, but the
    /// full structure is assumed to be sorted. This assumption can be enforced
    /// by using `pdb.full_sort()`.
    pub fn binary_find_atom(
        &'a self,
        serial_number: usize,
        alternative_location: Option<&str>,
    ) -> Option<AtomConformerResidueChainModel<'a>> {
        self.models().next().and_then(|m| {
            m.binary_find_atom(serial_number, alternative_location)
                .map(|res| res.extend(m))
        })
    }

    /// Get a mutable reference to the specified atom. Its uniqueness is guaranteed by
    /// including the `insertion_code`, with its full hierarchy. The algorithm is based
    /// on binary search so it is faster than an exhaustive search, but the
    /// full structure is assumed to be sorted. This assumption can be enforced
    /// by using `pdb.full_sort()`.
    pub fn binary_find_atom_mut(
        &'a mut self,
        serial_number: usize,
        alternative_location: Option<&str>,
    ) -> Option<AtomConformerResidueChainModelMut<'a>> {
        self.models_mut().next().and_then(|m| {
            let model: *mut Model = m;
            m.binary_find_atom_mut(serial_number, alternative_location)
                .map(|res| res.extend(model))
        })
    }

    /// Find all hierarchies matching the given search. For more details see [Search].
    /// ```
    /// use pdbtbx::*;
    /// let (pdb, errors) = ReadOptions::new().set_level(StrictnessLevel::Loose).read("example-pdbs/1ubq.pdb").unwrap();
    /// let selection = pdb.find(
    ///     Term::ChainId("A".to_owned())
    ///     & Term::ConformerName("GLY".to_owned())
    ///     & Term::AtomSerialNumber(750)
    /// );
    /// ```
    /// Find all hierarchies matching the given information
    pub fn find(
        &'a self,
        search: Search,
    ) -> impl DoubleEndedIterator<Item = AtomConformerResidueChainModel<'a>> + '_ {
        self.models()
            .map(move |m| (m, search.clone().add_model_info(m)))
            .filter(|(_m, search)| !matches!(search, Search::Known(false)))
            .flat_map(move |(m, search)| m.find(search).map(move |h| h.extend(m)))
    }

    /// Find all hierarchies matching the given search. For more details see [Search].
    pub fn find_mut(
        &'a mut self,
        search: Search,
    ) -> impl DoubleEndedIterator<Item = AtomConformerResidueChainModelMut<'a>> + '_ {
        self.models_mut()
            .map(move |m| {
                let search = search.clone().add_model_info(m);
                (m, search)
            })
            .filter(|(_m, search)| !matches!(search, Search::Known(false)))
            .flat_map(move |(m, search)| {
                let m_ptr: *mut Model = m;
                m.find_mut(search).map(move |h| h.extend(m_ptr))
            })
    }

    /// Get an iterator of references to Models making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn models(&self) -> impl DoubleEndedIterator<Item = &Model> + '_ {
        self.models.iter()
    }

    /// Get a parallel iterator of references to Models making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_models(&self) -> impl ParallelIterator<Item = &Model> + '_ {
        self.models.par_iter()
    }

    /// Get an iterator of mutable references to Models making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn models_mut(&mut self) -> impl DoubleEndedIterator<Item = &mut Model> + '_ {
        self.models.iter_mut()
    }

    /// Get a parallel iterator of mutable references to Models making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_models_mut(&mut self) -> impl ParallelIterator<Item = &mut Model> + '_ {
        self.models.par_iter_mut()
    }

    /// Get an iterator of references to Chains making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn chains(&self) -> impl DoubleEndedIterator<Item = &Chain> + '_ {
        self.models().flat_map(Model::chains)
    }

    /// Get a parallel iterator of references to Chains making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_chains(&self) -> impl ParallelIterator<Item = &Chain> + '_ {
        self.par_models().flat_map(Model::par_chains)
    }

    /// Get a iterator of mutable references to Chains making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn chains_mut(&mut self) -> impl DoubleEndedIterator<Item = &mut Chain> + '_ {
        self.models_mut().flat_map(Model::chains_mut)
    }

    /// Get a parallel iterator of mutable references to Chains making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_chains_mut(&mut self) -> impl ParallelIterator<Item = &mut Chain> + '_ {
        self.par_models_mut().flat_map(Model::par_chains_mut)
    }

    /// Get an iterator of references to Residues making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn residues(&self) -> impl DoubleEndedIterator<Item = &Residue> + '_ {
        self.models().flat_map(Model::residues)
    }

    /// Get a parallel iterator of references to Residues making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_residues(&self) -> impl ParallelIterator<Item = &Residue> + '_ {
        self.par_models().flat_map(Model::par_residues)
    }

    /// Get an iterator of mutable references to Residues making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn residues_mut(&mut self) -> impl DoubleEndedIterator<Item = &mut Residue> + '_ {
        self.models_mut().flat_map(Model::residues_mut)
    }

    /// Get a parallel iterator of mutable references to Residues making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_residues_mut(&mut self) -> impl ParallelIterator<Item = &mut Residue> + '_ {
        self.par_models_mut().flat_map(Model::par_residues_mut)
    }

    /// Get an iterator of references to Conformers making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn conformers(&self) -> impl DoubleEndedIterator<Item = &Conformer> + '_ {
        self.models().flat_map(Model::conformers)
    }

    /// Get a parallel iterator of references to Conformers making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_conformers(&self) -> impl ParallelIterator<Item = &Conformer> + '_ {
        self.par_models().flat_map(Model::par_conformers)
    }

    /// Get an iterator of mutable references to Conformers making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn conformers_mut(&mut self) -> impl DoubleEndedIterator<Item = &mut Conformer> + '_ {
        self.models_mut().flat_map(Model::conformers_mut)
    }

    /// Get a parallel iterator of mutable references to Conformers making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_conformers_mut(&mut self) -> impl ParallelIterator<Item = &mut Conformer> + '_ {
        self.par_models_mut().flat_map(Model::par_conformers_mut)
    }

    /// Get an iterator of references to Atom making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn atoms(&self) -> impl DoubleEndedIterator<Item = &Atom> + '_ {
        self.models().flat_map(Model::atoms)
    }

    /// Get a parallel iterator of references to Atom making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_atoms(&self) -> impl ParallelIterator<Item = &Atom> + '_ {
        self.par_models().flat_map(Model::par_atoms)
    }

    /// Get an iterator of mutable references to Atom making up this PDB.
    /// Double ended so iterating from the end is just as fast as from the start.
    pub fn atoms_mut(&mut self) -> impl DoubleEndedIterator<Item = &mut Atom> + '_ {
        self.models_mut().flat_map(Model::atoms_mut)
    }

    /// Get a parallel iterator of mutable references to Atom making up this PDB.
    #[doc_cfg(feature = "rayon")]
    pub fn par_atoms_mut(&mut self) -> impl ParallelIterator<Item = &mut Atom> + '_ {
        self.par_models_mut().flat_map(Model::par_atoms_mut)
    }

    /// Get an iterator of references to a struct containing all atoms with their hierarchy making up this PDB.
    pub fn atoms_with_hierarchy(
        &'a self,
    ) -> impl DoubleEndedIterator<Item = hierarchy::AtomConformerResidueChainModel<'a>> + '_ {
        self.models()
            .flat_map(|m| m.atoms_with_hierarchy().map(move |h| h.extend(m)))
    }

    /// Get an iterator of mutable references to a struct containing all atoms with their hierarchy making up this PDB.
    pub fn atoms_with_hierarchy_mut(
        &'a mut self,
    ) -> impl DoubleEndedIterator<Item = hierarchy::AtomConformerResidueChainModelMut<'a>> + '_
    {
        self.models_mut().flat_map(|m| {
            let model: *mut Model = m;
            m.atoms_with_hierarchy_mut().map(move |h| h.extend(model))
        })
    }

    /// Remove all Atoms matching the given predicate. The predicate will be run on all Atoms.
    /// As this is done in place this is the fastest way to remove Atoms from this PDB.
    pub fn remove_atoms_by<F>(&mut self, predicate: F)
    where
        F: Fn(&Atom) -> bool,
    {
        for residue in self.residues_mut() {
            residue.remove_atoms_by(&predicate);
        }
    }

    /// Remove all Conformers matching the given predicate. The predicate will be run on all Conformers.
    /// As this is done in place this is the fastest way to remove Conformers from this PDB.
    pub fn remove_conformers_by<F>(&mut self, predicate: F)
    where
        F: Fn(&Conformer) -> bool,
    {
        for chain in self.chains_mut() {
            chain.remove_conformers_by(&predicate);
        }
    }

    /// Remove all Residues matching the given predicate. The predicate will be run on all Residues.
    /// As this is done in place this is the fastest way to remove Residues from this PDB.
    pub fn remove_residues_by<F>(&mut self, predicate: F)
    where
        F: Fn(&Residue) -> bool,
    {
        for chain in self.chains_mut() {
            chain.remove_residues_by(&predicate);
        }
    }

    /// Remove all Residues matching the given predicate. The predicate will be run on all Residues.
    /// As this is done in place this is the fastest way to remove Residues from this PDB.
    pub fn remove_chains_by<F>(&mut self, predicate: F)
    where
        F: Fn(&Chain) -> bool,
    {
        for model in self.models_mut() {
            model.remove_chains_by(&predicate);
        }
    }

    /// Remove all Chains matching the given predicate. The predicate will be run on all Chains.
    /// As this is done in place this is the fastest way to remove Chains from this PDB.
    pub fn remove_models_by<F>(&mut self, predicate: F)
    where
        F: Fn(&Model) -> bool,
    {
        self.models.retain(|model| !predicate(model));
    }

    /// Remove the Model specified.
    ///
    /// ## Complexity
    /// * **Time**: amortized O(n) where n is the number of models, best case O(1) if model is the last one.
    /// * **Memory**: O(1)
    ///
    /// ## Arguments
    /// * `index` - the index of the Model to remove
    ///
    /// ## Panics
    /// Panics if the index is out of bounds.
    pub fn remove_model(&mut self, index: usize) {
        self.models.remove(index);
    }

    /// Remove all Models except for models
    /// specified by idxs.
    ///
    /// ## Complexity
    /// * **Time**: O(n) where n is the number of models.
    /// * **Memory**: O(k) where k is the number of models to keep.
    ///
    /// ## Arguments
    /// * `idxs` - the indices of the Models to keep
    ///
    /// ## Returns
    /// `None` if any of the indices are out of bounds.
    /// `Some(usize)` the number of models removed.
    pub fn remove_models_except(&mut self, idxs: &[usize]) -> Option<usize> {
        if self.models.is_empty() {
            #[cfg(debug_assertions)]
            eprint!("remove_models_except: no models to remove");
            return None;
        }

        let start_count = self.model_count();

        // bounds check the idxs
        if idxs.iter().max()? >= &start_count {
            return None;
        }

        let retained: Vec<_> = self
            .models
            .drain(..)
            .enumerate()
            .filter(|(i, _)| idxs.contains(i))
            .map(|(_, m)| m)
            .collect();

        let removed = start_count - retained.len();
        self.models.clear();
        self.models.extend(retained);

        Some(removed)
    }

    /// Remove all Models except for the first one.
    ///
    pub fn remove_all_models_except_first(&mut self) -> Option<usize> {
        self.remove_models_except(&[0])
    }

    /// Remove the Model specified. It returns `true` if it found a matching Model and removed it.
    /// It removes the first matching Model from the list.
    ///
    /// ## Arguments
    /// * `serial_number` - the serial number of the Model to remove
    pub fn remove_model_serial_number(&mut self, serial_number: usize) -> bool {
        let index = self
            .models
            .iter()
            .position(|a| a.serial_number() == serial_number);

        if let Some(i) = index {
            self.remove_model(i);
            true
        } else {
            false
        }
    }

    /// Remove the Model specified. It returns `true` if it found a matching Model and removed it.
    /// It removes the first matching Model from the list.
    /// Done in parallel.
    ///
    /// ## Arguments
    /// * `serial_number` - the serial number of the Model to remove
    #[doc_cfg(feature = "rayon")]
    pub fn par_remove_model_serial_number(&mut self, serial_number: usize) -> bool {
        let index = self
            .models
            .par_iter()
            .position_first(|a| a.serial_number() == serial_number);

        if let Some(i) = index {
            self.remove_model(i);
            true
        } else {
            false
        }
    }

    /// Remove all empty Models from this PDB, and all empty Chains from the Model, and all empty Residues from the Chains.
    pub fn remove_empty(&mut self) {
        self.models.iter_mut().for_each(Model::remove_empty);
        self.models.retain(|m| m.chain_count() > 0);
    }

    /// Remove all empty Models from this PDB, and all empty Chains from the Model, and all empty Residues from the Chains.
    /// Done in parallel.
    #[doc_cfg(feature = "rayon")]
    pub fn par_remove_empty(&mut self) {
        self.models.par_iter_mut().for_each(Model::remove_empty);
        self.models.retain(|m| m.chain_count() > 0);
    }

    /// This renumbers all numbered structs in the PDB.
    /// So it renumbers models, atoms, residues, chains and [`MtriX`]s.
    pub fn renumber(&mut self) {
        let mut model_counter = 1;
        for model in self.models_mut() {
            model.set_serial_number(model_counter);
            model_counter += 1;

            let mut counter = 1;
            for atom in model.atoms_mut() {
                atom.set_serial_number(counter);
                counter += 1;
            }
            let mut counter_i = 1;
            for residue in model.residues_mut() {
                residue.set_serial_number(counter_i);
                residue.remove_insertion_code();
                counter_i += 1;

                if residue.conformer_count() > 1 {
                    counter = 0;
                    for conformer in residue.conformers_mut() {
                        conformer.set_alternative_location(&number_to_base26(counter));
                        counter += 1;
                    }
                } else if let Some(conformer) = residue.conformer_mut(0) {
                    conformer.remove_alternative_location();
                }
            }
            counter = 0;
            for chain in model.chains_mut() {
                chain.set_id(&number_to_base26(counter));
                counter += 1;
            }
        }
    }

    /// Apply a transformation to the position of all atoms making up this PDB, the new position is immediately set.
    pub fn apply_transformation(&mut self, transformation: &TransformationMatrix) {
        for atom in self.atoms_mut() {
            atom.apply_transformation(transformation);
        }
    }

    /// Apply a transformation to the position of all atoms making up this PDB, the new position is immediately set.
    /// Done in parallel.
    #[doc_cfg(feature = "rayon")]
    pub fn par_apply_transformation(&mut self, transformation: &TransformationMatrix) {
        self.par_atoms_mut()
            .for_each(|atom| atom.apply_transformation(transformation));
    }

    /// Joins two PDBs. If one has multiple models it extends the models of this PDB with the models of the other PDB. If this PDB does
    /// not have any models it moves the models of the other PDB to this PDB. If both have one model it moves all chains/residues/atoms
    /// from the model of the other PDB to the model of this PDB. Effectively the same as calling join on those models.
    pub fn join(&mut self, mut other: PDB) {
        #[allow(clippy::unwrap_used)]
        if self.model_count() > 1 || other.model_count() > 1 {
            self.models.extend(other.models);
        } else if self.model_count() == 0 {
            self.models = other.models;
        } else if other.model_count() == 0 {
            // There is nothing to join
        } else if let Some(model) = self.model_mut(0) {
            model.join(other.models.remove(0));
        }
    }

    /// Sort the Models of this PDB.
    pub fn sort(&mut self) {
        self.models.sort();
    }

    /// Sort the Models of this PDB in parallel.
    #[doc_cfg(feature = "rayon")]
    pub fn par_sort(&mut self) {
        self.models.par_sort();
    }

    /// Sort all structs in this PDB.
    pub fn full_sort(&mut self) {
        self.sort();
        for model in self.models_mut() {
            model.sort();
        }
        for chain in self.chains_mut() {
            chain.sort();
        }
        for residue in self.residues_mut() {
            residue.sort();
        }
        for conformer in self.conformers_mut() {
            conformer.sort();
        }
    }

    /// Sort all structs in this PDB in parallel.
    #[doc_cfg(feature = "rayon")]
    pub fn par_full_sort(&mut self) {
        self.par_sort();
        self.par_models_mut().for_each(Model::par_sort);
        self.par_chains_mut().for_each(Chain::par_sort);
        self.par_residues_mut().for_each(Residue::par_sort);
        self.par_conformers_mut().for_each(Conformer::par_sort);
    }

    /// Create an R star tree of Atoms which can be used for fast lookup of
    /// spatially close atoms. See the crate rstar for documentation
    /// on how to use the tree. (<https://crates.io/crates/rstar>)
    ///
    /// Keep in mind that this creates a tree that is separate from
    /// the original PDB, so any changes to one of the data
    /// structures is not seen in the other data structure (until
    /// you generate a new tree of course).
    #[doc_cfg(feature = "rstar")]
    pub fn create_atom_rtree(&self) -> rstar::RTree<&Atom> {
        rstar::RTree::bulk_load(self.atoms().collect())
    }

    /// Create an R star tree of structs containing Atoms and their hierarchies
    /// which can be used for fast lookup of
    /// spatial close atoms. See the crate rstar for documentation
    /// on how to use the tree. (<https://crates.io/crates/rstar>)
    ///
    /// Keep in mind that this creates a tree that is separate from
    /// the original PDB, so any changes to one of the data
    /// structures is not seen in the other data structure (until
    /// you generate a new tree of course).
    #[doc_cfg(feature = "rstar")]
    pub fn create_hierarchy_rtree(
        &'a self,
    ) -> rstar::RTree<hierarchy::AtomConformerResidueChainModel<'a>> {
        rstar::RTree::bulk_load(self.atoms_with_hierarchy().collect())
    }

    /// Finds the square bounding box around the PDB. The first tuple
    /// is the bottom left point, lowest value for all dimensions
    /// for all points. The second tuple is the top right point, the
    /// highest value for all dimensions for all points.
    pub fn bounding_box(&self) -> ((f64, f64, f64), (f64, f64, f64)) {
        let mut min = [f64::MAX, f64::MAX, f64::MAX];
        let mut max = [f64::MIN, f64::MIN, f64::MIN];
        for atom in self.atoms() {
            if atom.x() < min[0] {
                min[0] = atom.x();
            }
            if atom.y() < min[1] {
                min[1] = atom.y();
            }
            if atom.z() < min[2] {
                min[2] = atom.z();
            }
            if atom.x() > max[0] {
                max[0] = atom.x();
            }
            if atom.y() > max[1] {
                max[1] = atom.y();
            }
            if atom.z() > max[2] {
                max[2] = atom.z();
            }
        }
        ((min[0], min[1], min[2]), (max[0], max[1], max[2]))
    }

    /// Get the bonds in this PDB file. Runtime is `O(bonds_count * 2 * atom_count)` because it
    /// has to iterate over all atoms to prevent borrowing problems.
    pub fn bonds(&self) -> impl DoubleEndedIterator<Item = (&Atom, &Atom, Bond)> + '_ {
        self.bonds.iter().map(move |(a, b, bond)| {
            (
                self.atoms()
                    .find(|atom| atom.counter() == *a)
                    .expect("Could not find an atom in the bonds list"),
                self.atoms()
                    .find(|atom| atom.counter() == *b)
                    .expect("Could not find an atom in the bonds list"),
                *bond,
            )
        })
    }

    /// Add a bond of the given type to the list of bonds in this PDB.
    /// The atoms are selected by serial number and alternative location.
    /// It uses `binary_find_atom` in the background so the PDB should be sorted.
    /// If one of the atoms could not be found it returns `None` otherwise it
    /// will return `Some(())`.
    pub fn add_bond(
        &mut self,
        atom1: (usize, Option<&str>),
        atom2: (usize, Option<&str>),
        bond: Bond,
    ) -> Option<()> {
        self.bonds.push((
            self.binary_find_atom(atom1.0, atom1.1)?.atom().counter(),
            self.binary_find_atom(atom2.0, atom2.1)?.atom().counter(),
            bond,
        ));
        Some(())
    }

    /// Add a bond of the given type to the list of bonds in this PDB.
    /// The raw counters of the atoms are given.
    pub(crate) fn add_bond_counters(&mut self, atom1: usize, atom2: usize, bond: Bond) {
        self.bonds.push((atom1, atom2, bond));
    }

    /// Returns a HashMap with the chains in contact within a given distance.
    ///
    /// # Arguments
    ///
    /// * `distance` - A f64 value representing the maximum distance between two atoms for them to be considered in contact.
    ///
    /// # Returns
    ///
    /// A HashMap with the chains in contact. The keys are the chain IDs and the values are vectors with the IDs of the chains in contact with the key chain.
    pub fn chains_in_contact(&self, distance: f64) -> HashMap<String, Vec<String>> {
        let mut chains = HashMap::new();
        for chain1 in self.chains() {
            for chain2 in self.chains() {
                if chain1.id() == chain2.id() {
                    continue;
                }
                for atom1 in chain1.atoms() {
                    for atom2 in chain2.atoms() {
                        if atom1.distance(atom2) < distance {
                            let chain1_id = chain1.id().to_owned();
                            let chain2_id = chain2.id().to_owned();
                            let entry = chains.entry(chain1_id).or_insert_with(Vec::new);
                            if !entry.contains(&chain2_id) {
                                entry.push(chain2_id)
                            }
                            break;
                        }
                    }
                }
            }
        }
        chains
    }

    /// Returns a vector of unique conformer names present in the PDB file.
    ///
    /// # Arguments
    ///
    /// * `self` - A reference to the PDB struct.
    ///
    /// # Returns
    ///
    /// * `Vec<String>` - A vector of unique conformer names.
    pub fn unique_conformer_names(&self) -> Vec<String> {
        let mut names = Vec::new();
        for conformer in self.conformers() {
            let name = conformer.name().to_owned();
            if let Some(index) = names.binary_search(&name).err() {
                names.insert(index, name);
            }
        }
        names
    }
}

impl fmt::Display for PDB {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "PDB Models: {}", self.models.len())
    }
}

impl Default for PDB {
    fn default() -> Self {
        Self::new()
    }
}

impl Extend<Model> for PDB {
    /// Extend the Models on this PDB by the given iterator of Models.
    fn extend<T: IntoIterator<Item = Model>>(&mut self, iter: T) {
        self.models.extend(iter);
    }
}

impl FromIterator<Model> for PDB {
    fn from_iter<T: IntoIterator<Item = Model>>(iter: T) -> Self {
        let mut pdb = Self::default();
        pdb.extend(iter);
        pdb
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod tests {
    use std::path::Path;

    use crate::ReadOptions;

    use super::*;

    #[test]
    fn remove_model() {
        let pdb_dir = Path::new(env!("CARGO_MANIFEST_DIR")).join("example-pdbs");

        for entry in std::fs::read_dir(pdb_dir).unwrap() {
            let entry = entry.unwrap();
            let path = entry.path();
            if path.extension().unwrap() != "pdb" {
                continue;
            }
            let (pdb, _) = ReadOptions::default()
                .set_level(crate::StrictnessLevel::Loose)
                .read(path.to_str().unwrap())
                .unwrap();

            let model_count = pdb.model_count();
            let mut test_pdb = pdb.clone();
            test_pdb.remove_model(0);
            assert_eq!(test_pdb.model_count(), model_count - 1);

            let mut test_pdb = pdb.clone();
            test_pdb.remove_all_models_except_first();

            assert_eq!(test_pdb.model_count(), 1);
            assert_eq!(test_pdb.model(0).unwrap(), pdb.model(0).unwrap());

            let mut test_pdb = pdb.clone();
            test_pdb.remove_models_except(&[0]);
            assert_eq!(test_pdb.model_count(), 1);
            assert_eq!(test_pdb.model(0).unwrap(), pdb.model(0).unwrap());
        }
    }

    #[test]
    fn remove_model_except() {
        let path = Path::new(env!("CARGO_MANIFEST_DIR"))
            .join("example-pdbs")
            .join("models.pdb");

        let (pdb, _) = ReadOptions::default()
            .set_level(crate::StrictnessLevel::Loose)
            .read(path.to_str().unwrap())
            .unwrap();

        // Has 6 models
        let model_count = pdb.model_count();
        assert_eq!(model_count, 6);
        let mut test_pdb = pdb.clone();

        // Removing at an index higher than the number is not allowed
        assert_eq!(None, test_pdb.remove_models_except(&[6]));

        // Make sure nothing is altered
        assert_eq!(pdb, test_pdb);

        // Only keep models 1, 3, 5 (in a not sorted order because it should be able to handle that)
        assert_eq!(Some(3), test_pdb.remove_models_except(&[5, 1, 3]));

        // Debug output shows the order of the models (each model contains a single atom `ILX` where X is the model number)
        eprintln!(
            "{}",
            pdb.models().fold("Original".to_string(), |acc, m| {
                acc + ", " + m.residue(0).unwrap().name().unwrap()
            })
        );
        eprintln!(
            "{}",
            test_pdb.models().fold("Removed".to_string(), |acc, m| {
                acc + ", " + m.residue(0).unwrap().name().unwrap()
            })
        );

        // Make sure the retained models are all retained and are the only retained models
        assert_eq!(3, test_pdb.model_count());
        assert!(test_pdb.models().any(|m| m == pdb.model(1).unwrap()));
        assert!(test_pdb.models().any(|m| m == pdb.model(3).unwrap()));
        assert!(test_pdb.models().any(|m| m == pdb.model(5).unwrap()));
    }

    #[test]
    fn sort_atoms() {
        let a = Atom::new(false, 0, "", 0.0, 0.0, 0.0, 0.0, 0.0, "", 0).unwrap();
        let b = Atom::new(false, 1, "", 0.0, 0.0, 0.0, 0.0, 0.0, "", 0).unwrap();
        let mut model = Model::new(0);
        model.add_atom(b, "A", (0, None), ("LYS", None));
        model.add_atom(a, "A", (0, None), ("LYS", None));
        let mut pdb = PDB::new();
        pdb.add_model(model);
        assert_eq!(pdb.atom(0).unwrap().serial_number(), 1);
        assert_eq!(pdb.atom(1).unwrap().serial_number(), 0);
        pdb.full_sort();
        assert_eq!(pdb.atom(0).unwrap().serial_number(), 0);
        assert_eq!(pdb.atom(1).unwrap().serial_number(), 1);
    }

    #[test]
    fn binary_lookup() {
        let mut model = Model::new(0);
        model.add_atom(
            Atom::new(false, 1, "", 0.0, 0.0, 0.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", Some("A")),
        );
        model.add_atom(
            Atom::new(false, 1, "", 1.0, 0.0, 0.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", Some("B")),
        );
        model.add_atom(
            Atom::new(false, 1, "", 2.0, 0.0, 0.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        let mut pdb = PDB::new();
        pdb.add_model(model);
        pdb.full_sort();

        assert_eq!(pdb.binary_find_atom(1, Some("A")).unwrap().atom().x(), 0.0);
        assert_eq!(pdb.binary_find_atom(1, Some("B")).unwrap().atom().x(), 1.0);
        assert_eq!(pdb.binary_find_atom(1, None).unwrap().atom().x(), 2.0);
    }

    #[test]
    #[cfg(feature = "rstar")]
    fn spatial_lookup() {
        let mut model = Model::new(0);
        model.add_atom(
            Atom::new(false, 0, "", 0.0, 0.0, 0.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        model.add_atom(
            Atom::new(false, 1, "", 1.0, 1.0, 1.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        model.add_atom(
            Atom::new(false, 2, "", 0.0, 1.0, 1.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        let mut pdb = PDB::new();
        pdb.add_model(model);
        let tree = pdb.create_atom_rtree();
        assert_eq!(tree.size(), 3);
        assert_eq!(
            tree.nearest_neighbor(&(1.0, 1.0, 1.0))
                .unwrap()
                .serial_number(),
            1
        );
        assert_eq!(
            tree.locate_within_distance((1.0, 1.0, 1.0), 1.0)
                .fold(0, |acc, _| acc + 1),
            2
        );
        let mut neighbors = tree.nearest_neighbor_iter(&pdb.atom(0).unwrap().pos());
        assert_eq!(neighbors.next().unwrap().serial_number(), 0);
        assert_eq!(neighbors.next().unwrap().serial_number(), 2);
        assert_eq!(neighbors.next().unwrap().serial_number(), 1);
        assert_eq!(neighbors.next(), None);
    }

    #[test]
    #[cfg(feature = "rstar")]
    fn spatial_lookup_with_hierarchy() {
        let mut model = Model::new(0);
        model.add_atom(
            Atom::new(false, 0, "", 0.0, 0.0, 0.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        model.add_atom(
            Atom::new(false, 1, "", 1.0, 1.0, 1.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        model.add_atom(
            Atom::new(false, 2, "", 0.0, 1.0, 1.0, 0.0, 0.0, "", 0).unwrap(),
            "B",
            (0, None),
            ("MET", None),
        );
        let mut pdb = PDB::new();
        pdb.add_model(model);
        let tree = pdb.create_hierarchy_rtree();
        assert_eq!(tree.size(), 3);
        assert_eq!(
            tree.nearest_neighbor(&(1.0, 1.0, 1.0))
                .unwrap()
                .atom()
                .serial_number(),
            1
        );
        assert_eq!(
            tree.locate_within_distance((1.0, 1.0, 1.0), 1.0)
                .fold(0, |acc, _| acc + 1),
            2
        );
        let mut neighbors = tree.nearest_neighbor_iter(&pdb.atom(0).unwrap().pos());
        let a = neighbors.next().unwrap();
        let b = neighbors.next().unwrap();
        let c = neighbors.next().unwrap();
        assert_eq!(neighbors.next(), None);

        assert_eq!(a.atom().serial_number(), 0);
        assert_eq!(b.atom().serial_number(), 2);
        assert_eq!(c.atom().serial_number(), 1);
        assert_eq!(a.chain().id(), "A");
        assert_eq!(b.chain().id(), "B");
        assert_eq!(c.chain().id(), "A");
    }

    #[test]
    #[cfg(feature = "serde")]
    fn serialization() {
        use serde_json;

        let mut model = Model::new(0);
        model.add_atom(
            Atom::new(false, 0, "", 0.0, 0.0, 0.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        model.add_atom(
            Atom::new(false, 1, "", 1.0, 1.0, 1.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        model.add_atom(
            Atom::new(false, 2, "", 0.0, 1.0, 1.0, 0.0, 0.0, "", 0).unwrap(),
            "B",
            (0, None),
            ("MET", None),
        );
        let pdb = PDB::new();

        let json = serde_json::to_string(&pdb).unwrap();
        let parsed = serde_json::from_str(&json).unwrap();
        assert_eq!(pdb, parsed);
    }

    #[test]
    fn bounding_box() {
        let mut model = Model::new(0);
        model.add_atom(
            Atom::new(false, 0, "", -1.0, 0.0, 2.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        model.add_atom(
            Atom::new(false, 1, "", 1.0, 2.0, -1.0, 0.0, 0.0, "", 0).unwrap(),
            "A",
            (0, None),
            ("MET", None),
        );
        model.add_atom(
            Atom::new(false, 2, "", 2.0, -1.0, 0.5, 0.0, 0.0, "", 0).unwrap(),
            "B",
            (0, None),
            ("MET", None),
        );
        let mut pdb = PDB::new();
        pdb.add_model(model);
        assert_eq!(((-1., -1., -1.), (2., 2., 2.)), pdb.bounding_box());
    }

    #[test]
    fn chains_in_contact() {
        let path = Path::new(env!("CARGO_MANIFEST_DIR"))
            .join("example-pdbs")
            .join("1yyf.pdb");
        let (pdb, _) = ReadOptions::default()
            .set_level(crate::StrictnessLevel::Loose)
            .read(path.to_str().unwrap())
            .unwrap();

        let chainmap = pdb.chains_in_contact(5.0);

        let my_map: HashMap<String, Vec<String>> = [
            ("B".to_string(), vec!["A".to_string(), "C".to_string()]),
            ("A".to_string(), vec!["B".to_string(), "D".to_string()]),
            ("D".to_string(), vec!["A".to_string(), "C".to_string()]),
            ("C".to_string(), vec!["B".to_string(), "D".to_string()]),
        ]
        .iter()
        .cloned()
        .collect();

        assert_eq!(chainmap, my_map);
    }
    #[test]
    fn test_unique_conformer_names() {
        let path = Path::new(env!("CARGO_MANIFEST_DIR"))
            .join("example-pdbs")
            .join("1ubq.pdb");
        let (pdb, _) = ReadOptions::default()
            .set_level(crate::StrictnessLevel::Loose)
            .read(path.to_str().unwrap())
            .unwrap();

        let reslist = pdb.unique_conformer_names();
        let expected_reslist: Vec<String> = vec![
            "MET".to_string(),
            "GLN".to_string(),
            "ILE".to_string(),
            "PHE".to_string(),
            "VAL".to_string(),
            "LYS".to_string(),
            "THR".to_string(),
            "LEU".to_string(),
            "GLY".to_string(),
            "GLU".to_string(),
            "PRO".to_string(),
            "SER".to_string(),
            "ASP".to_string(),
            "ASN".to_string(),
            "ALA".to_string(),
            "ARG".to_string(),
            "TYR".to_string(),
            "HIS".to_string(),
            "HOH".to_string(),
        ];

        assert!(reslist.iter().all(|x| expected_reslist.contains(x)));
        assert!(expected_reslist.iter().all(|x| reslist.contains(x)));
        assert_eq!(reslist.len(), 19);
    }
}