# Program Listing for File model.hpp

↰ Return to documentation for file (`cif++/model.hpp`)

    /*-
     * SPDX-License-Identifier: BSD-2-Clause
     *
     * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
     *
     * Redistribution and use in source and binary forms, with or without
     * modification, are permitted provided that the following conditions are met:
     *
     * 1. Redistributions of source code must retain the above copyright notice, this
     *    list of conditions and the following disclaimer
     * 2. Redistributions in binary form must reproduce the above copyright notice,
     *    this list of conditions and the following disclaimer in the documentation
     *    and/or other materials provided with the distribution.
     *
     * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
     * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
     * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
     * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
     * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
     * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
     * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
     * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
     * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
     * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
     */

    #pragma once

    #include "cif++/atom_type.hpp"
    #include "cif++/datablock.hpp"
    #include "cif++/point.hpp"
    #include "cif++/row.hpp"

    #include <memory>
    #include <numeric>

    #if __cpp_lib_format
    # include <format>
    #endif

    namespace cif::mm
    {

    class atom;
    class residue;
    class monomer;
    class polymer;
    class structure;

    // --------------------------------------------------------------------

    class atom
    {
      private:
        struct atom_impl : public std::enable_shared_from_this<atom_impl>
        {
            atom_impl(const datablock &db, std::string_view id)
                : m_db(db)
                , m_cat(db["atom_site"])
                , m_id(id)
            {
                auto r = row();
                if (r)
                    tie(m_location.m_x, m_location.m_y, m_location.m_z) = r.get("Cartn_x", "Cartn_y", "Cartn_z");
            }

            // constructor for a symmetry copy of an atom
            atom_impl(const atom_impl &impl, const point &loc, const std::string &sym_op)
                : atom_impl(impl)
            {
                m_location = loc;
                m_symop = sym_op;
            }

            atom_impl(const atom_impl &i) = default;

            int compare(const atom_impl &b) const;

            // bool getAnisoU(float anisou[6]) const;

            int get_charge() const;

            void moveTo(const point &p);

            // const compound *compound() const;

            std::string get_property(std::string_view name) const;
            int get_property_int(std::string_view name) const;
            float get_property_float(std::string_view name) const;

            void set_property(const std::string_view name, const std::string &value);

            row_handle row()
            {
                return m_cat[{ { "id", m_id } }];
            }

            const row_handle row() const
            {
                return m_cat[{ { "id", m_id } }];
            }

            row_handle row_aniso()
            {
                row_handle result{};
                auto cat = m_db.get("atom_site_anisotrop");
                if (cat)
                    result = cat->operator[]({ { "id", m_id } });
                return result;
            }

            const row_handle row_aniso() const
            {
                row_handle result{};
                auto cat = m_db.get("atom_site_anisotrop");
                if (cat)
                    result = cat->operator[]({ { "id", m_id } });
                return result;
            }

            const datablock &m_db;
            const category &m_cat;
            std::string m_id;
            point m_location;
            std::string m_symop = "1_555";
        };
      public:
        atom() {}

        atom(std::shared_ptr<atom_impl> impl)
            : m_impl(impl)
        {
        }

        atom(const atom &rhs)
            : m_impl(rhs.m_impl)
        {
        }

        atom(const datablock &db, const row_handle &row)
            : atom(std::make_shared<atom_impl>(db, row["id"].as<std::string>()))
        {
        }

        atom(const atom &rhs, const point &symmmetry_location, const std::string &symmetry_operation)
            : atom(std::make_shared<atom_impl>(*rhs.m_impl, symmmetry_location, symmetry_operation))
        {
        }

        explicit operator bool() const { return (bool)m_impl; }

        atom &operator=(const atom &rhs) = default;

        std::string get_property(std::string_view name) const
        {
            if (not m_impl)
                throw std::logic_error("Error trying to fetch a property from an uninitialized atom");
            return m_impl->get_property(name);
        }

        int get_property_int(std::string_view name) const
        {
            if (not m_impl)
                throw std::logic_error("Error trying to fetch a property from an uninitialized atom");
            return m_impl->get_property_int(name);
        }

        float get_property_float(std::string_view name) const
        {
            if (not m_impl)
                throw std::logic_error("Error trying to fetch a property from an uninitialized atom");
            return m_impl->get_property_float(name);
        }

        void set_property(const std::string_view name, const std::string &value)
        {
            if (not m_impl)
                throw std::logic_error("Error trying to modify an uninitialized atom");
            m_impl->set_property(name, value);
        }

        template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
        void set_property(const std::string_view name, const T &value)
        {
            set_property(name, std::to_string(value));
        }

        const std::string &id() const { return impl().m_id; }

        cif::atom_type get_type() const { return atom_type_traits(get_property("type_symbol")).type(); }

        point get_location() const { return impl().m_location; }

        void set_location(point p)
        {
            if (not m_impl)
                throw std::logic_error("Error trying to modify an uninitialized atom");
            m_impl->moveTo(p);
        }

        void translate(point t)
        {
            set_location(get_location() + t);
        }

        void rotate(quaternion q)
        {
            auto loc = get_location();
            loc.rotate(q);
            set_location(loc);
        }

        void rotate(quaternion q, point p)
        {
            auto loc = get_location();
            loc.rotate(q, p);
            set_location(loc);
        }

        void translate_and_rotate(point t, quaternion q)
        {
            auto loc = get_location();
            loc += t;
            loc.rotate(q);
            set_location(loc);
        }

        void translate_rotate_and_translate(point t1, quaternion q, point t2)
        {
            auto loc = get_location();
            loc += t1;
            loc.rotate(q);
            loc += t2;
            set_location(loc);
        }

        const row_handle get_row() const { return impl().row(); }

        const row_handle get_row_aniso() const { return impl().row_aniso(); }

        bool is_symmetry_copy() const { return impl().m_symop != "1_555"; }

        std::string symmetry() const { return impl().m_symop; }

        bool is_water() const
        {
            auto comp_id = get_label_comp_id();
            return comp_id == "HOH" or comp_id == "H2O" or comp_id == "WAT";
        }

        int get_charge() const { return impl().get_charge(); }

        float get_occupancy() const { return get_property_float("occupancy"); }

        // specifications

        std::string get_label_asym_id() const { return get_property("label_asym_id"); }
        int get_label_seq_id() const { return get_property_int("label_seq_id"); }
        std::string get_label_atom_id() const { return get_property("label_atom_id"); }
        std::string get_label_alt_id() const { return get_property("label_alt_id"); }
        std::string get_label_comp_id() const { return get_property("label_comp_id"); }
        std::string get_label_entity_id() const { return get_property("label_entity_id"); }

        std::string get_auth_asym_id() const { return get_property("auth_asym_id"); }
        std::string get_auth_seq_id() const { return get_property("auth_seq_id"); }
        std::string get_auth_atom_id() const { return get_property("auth_atom_id"); }
        std::string get_auth_alt_id() const { return get_property("pdbx_auth_alt_id"); }
        std::string get_auth_comp_id() const { return get_property("auth_comp_id"); }
        std::string get_pdb_ins_code() const { return get_property("pdbx_PDB_ins_code"); }

        bool is_alternate() const
        {
            if (auto alt_id = get_label_alt_id(); alt_id.empty() or alt_id == ".")
                return false;
            return true;
        }

        std::string pdb_id() const
        {
            return get_label_comp_id() + '_' + get_auth_asym_id() + '_' + get_auth_seq_id() + get_pdb_ins_code();
        }

        bool operator==(const atom &rhs) const
        {
            if (m_impl == rhs.m_impl)
                return true;

            if (not(m_impl and rhs.m_impl))
                return false;

            return &m_impl->m_db == &rhs.m_impl->m_db and m_impl->m_id == rhs.m_impl->m_id;
        }

        bool operator!=(const atom &rhs) const
        {
            return not operator==(rhs);
        }

        bool is_back_bone() const
        {
            auto atomID = get_label_atom_id();
            return atomID == "N" or atomID == "O" or atomID == "C" or atomID == "CA";
        }

        void swap(atom &b)
        {
            std::swap(m_impl, b.m_impl);
        }

        int compare(const atom &b) const { return impl().compare(*b.m_impl); }

        bool operator<(const atom &rhs) const
        {
            return compare(rhs) < 0;
        }

        friend std::ostream &operator<<(std::ostream &os, const atom &atom);

      private:
        friend class structure;

        const atom_impl &impl() const
        {
            if (not m_impl)
                throw std::runtime_error("Uninitialized atom, not found?");
            return *m_impl;
        }

        std::shared_ptr<atom_impl> m_impl;
    };

    inline void swap(atom &a, atom &b)
    {
        a.swap(b);
    }

    inline float distance(const atom &a, const atom &b)
    {
        return distance(a.get_location(), b.get_location());
    }

    inline float distance_squared(const atom &a, const atom &b)
    {
        return distance_squared(a.get_location(), b.get_location());
    }

    // --------------------------------------------------------------------

    enum class EntityType
    {
        Polymer,
        NonPolymer,
        Macrolide,
        Water,
        Branched
    };

    // --------------------------------------------------------------------

    class residue
    {
      public:
        friend class structure;

        residue(structure &structure, const std::string &compoundID,
            const std::string &asymID, int seqID,
            const std::string &authAsymID, const std::string &authSeqID,
            const std::string &pdbInsCode)
            : m_structure(&structure)
            , m_compound_id(compoundID)
            , m_asym_id(asymID)
            , m_seq_id(seqID)
            , m_pdb_strand_id(authAsymID)
            , m_pdb_seq_num(authSeqID)
            , m_pdb_ins_code(pdbInsCode)
        {
        }

        residue(structure &structure, const std::vector<atom> &atoms);

        residue(const residue &rhs) = delete;
        residue &operator=(const residue &rhs) = delete;

        residue(residue &&rhs) = default;
        residue &operator=(residue &&rhs) = default;

        virtual ~residue() = default;
        std::string get_entity_id() const;

        EntityType entity_type() const;

        const std::string &get_asym_id() const { return m_asym_id; }
        int get_seq_id() const { return m_seq_id; }

        const std::string get_pdb_strand_id() const { return m_pdb_strand_id; }
        const std::string get_pdb_seq_num() const { return m_pdb_seq_num; }
        std::string get_pdb_ins_code() const { return m_pdb_ins_code; }

        const std::string &get_compound_id() const { return m_compound_id; }
        void set_compound_id(const std::string &id) { m_compound_id = id; }

        structure *get_structure() const { return m_structure; }

        std::vector<atom> &atoms()
        {
            return m_atoms;
        }

        const std::vector<atom> &atoms() const
        {
            return m_atoms;
        }

        void add_atom(atom &atom);

        std::vector<atom> unique_atoms() const;

        atom get_atom_by_atom_id(const std::string &atomID) const;

        atom get_atom_by_atom_id(const std::string &atomID, const std::string &altID) const;

        std::vector<atom> get_atoms_by_id(const std::string &atomID) const;

        bool is_entity() const;

        bool is_water() const { return m_compound_id == "HOH"; }

        bool has_alternate_atoms() const;

        bool has_alternate_atoms_for(const std::string &atomID) const;

        std::set<std::string> get_alternate_ids() const;

        std::set<std::string> get_atom_ids() const;

        std::tuple<point, float> center_and_radius() const;

        friend std::ostream &operator<<(std::ostream &os, const residue &res);

        bool operator==(const residue &rhs) const
        {
            return this == &rhs or (m_structure == rhs.m_structure and
                                       m_seq_id == rhs.m_seq_id and
                                       m_asym_id == rhs.m_asym_id and
                                       m_compound_id == rhs.m_compound_id and
                                       m_pdb_seq_num == rhs.m_pdb_seq_num);
        }

        virtual atom create_new_atom(atom_type inType, const std::string &inAtomID, point inLocation);

      protected:
        residue() {}

        structure *m_structure = nullptr;
        std::string m_compound_id, m_asym_id;
        int m_seq_id = 0;
        std::string m_pdb_strand_id, m_pdb_seq_num, m_pdb_ins_code;
        std::vector<atom> m_atoms;
    };

    // --------------------------------------------------------------------

    class monomer : public residue
    {
      public:
        monomer(const monomer &rhs) = delete;
        monomer &operator=(const monomer &rhs) = delete;

        monomer(monomer &&rhs);

        monomer &operator=(monomer &&rhs);

        monomer(const polymer &polymer, std::size_t index, int seqID, const std::string &authSeqID,
            const std::string &pdbInsCode, const std::string &compoundID);

        bool is_first_in_chain() const;
        bool is_last_in_chain() const;

        const monomer &prev() const; // Return previous monomer in polymer
        const monomer &next() const; // Return next monomer in polymer

        // convenience
        bool has_alpha() const;
        bool has_kappa() const;

        // Assuming this is really an amino acid...

        float phi() const;
        float psi() const;
        float alpha() const;
        float kappa() const;
        float tco() const;
        float omega() const;

        // torsion angles
        std::size_t nr_of_chis() const;
        float chi(std::size_t i) const;

        bool is_cis() const;

        bool is_complete() const;

        bool has_alternate_backbone_atoms() const;

        atom CAlpha() const { return get_atom_by_atom_id("CA"); }
        atom C() const { return get_atom_by_atom_id("C"); }
        atom N() const { return get_atom_by_atom_id("N"); }
        atom O() const { return get_atom_by_atom_id("O"); }
        atom H() const { return get_atom_by_atom_id("H"); }

        bool is_bonded_to(const monomer &rhs) const
        {
            return this != &rhs and are_bonded(*this, rhs);
        }

        static bool are_bonded(const monomer &a, const monomer &b, float errorMargin = 0.5f);

        static bool is_cis(const monomer &a, const monomer &b);

        static float omega(const monomer &a, const monomer &b);

        float chiral_volume() const;

        bool operator==(const monomer &rhs) const
        {
            return m_polymer == rhs.m_polymer and m_index == rhs.m_index;
        }

        atom create_new_atom(atom_type inType, const std::string &inAtomID, point inLocation) override;

      private:
        const polymer *m_polymer;
        std::size_t m_index;
    };

    // --------------------------------------------------------------------

    class polymer : public std::vector<monomer>
    {
      public:
        polymer(structure &s, const std::string &entityID, const std::string &asymID, const std::string &auth_asym_id);

        polymer(const polymer &) = delete;
        polymer &operator=(const polymer &) = delete;

        structure *get_structure() const { return m_structure; }

        std::string get_asym_id() const { return m_asym_id; }
        std::string get_pdb_strand_id() const { return m_pdb_strand_id; }
        std::string get_entity_id() const { return m_entity_id; }

      private:
        structure *m_structure;
        std::string m_entity_id;
        std::string m_asym_id;
        std::string m_pdb_strand_id;
    };

    // --------------------------------------------------------------------
    // sugar and branch, to describe glycosylation sites

    class branch;

    class sugar : public residue
    {
      public:
        sugar(branch &branch, const std::string &compoundID,
            const std::string &asymID, int authSeqID);

        sugar(sugar &&rhs);
        sugar &operator=(sugar &&rhs);
        int num() const
        {
            int result;
            auto r = std::from_chars(m_pdb_seq_num.data(), m_pdb_seq_num.data() + m_pdb_seq_num.length(), result);
            if ((bool)r.ec)
                throw std::runtime_error("The auth_seq_id should be a number for a sugar");
            return result;
        }

        std::string name() const;

        atom get_link() const { return m_link; }

        void set_link(atom link) { m_link = link; }

        std::size_t get_link_nr() const
        {
            std::size_t result = 0;
            if (m_link)
                result = m_link.get_property_int("auth_seq_id");
            return result;
        }

        atom add_atom(row_initializer atom_info);

      private:
        branch *m_branch;
        atom m_link;
    };

    class branch : public std::vector<sugar>
    {
      public:
        branch(structure &structure, const std::string &asym_id, const std::string &entity_id);

        branch(const branch &) = delete;
        branch &operator=(const branch &) = delete;

        branch(branch &&) = default;
        branch &operator=(branch &&) = default;
        void link_atoms();

        std::string name() const;

        float weight() const;

        std::string get_asym_id() const { return m_asym_id; }
        std::string get_entity_id() const { return m_entity_id; }

        structure &get_structure() { return *m_structure; }
        structure &get_structure() const { return *m_structure; }

        sugar &get_sugar_by_num(int nr);

        const sugar &get_sugar_by_num(int nr) const
        {
            return const_cast<branch *>(this)->get_sugar_by_num(nr);
        }

        sugar &construct_sugar(const std::string &compound_id);

        sugar &construct_sugar(const std::string &compound_id, const std::string &atom_id,
            int linked_sugar_nr, const std::string &linked_atom_id);

      private:
        friend sugar;

        std::string name(const sugar &s) const;

        structure *m_structure;
        std::string m_asym_id, m_entity_id;
    };

    enum class occupancy_policy
    {
        ALL = 0,

        MAX = 1,

        MIN = 2,

        UNOCCUPIED = 3
    };

    struct structure_open_options
    {
        bool skip_hydrogen = false;
        bool skip_hetatom = false;
        bool skip_water = false;
        occupancy_policy occupancy_mode = occupancy_policy::ALL;
        std::vector<std::string> asyms;
        std::optional<float> min_b_factor;
        std::optional<float> max_b_factor;
    };

    // --------------------------------------------------------------------

    class structure
    {
      public:
        structure(file &p, std::size_t modelNr = 1, structure_open_options options = {});

        structure(datablock &db, std::size_t modelNr = 1, structure_open_options options = {});

        structure(structure &&s) = default;
        // structures cannot be copied.

        structure(const structure &) = delete;
        structure &operator=(const structure &) = delete;
        ~structure() = default;

        std::size_t get_model_nr() const { return m_model_nr; }

        const std::vector<atom> &atoms() const { return m_atoms; }

        EntityType get_entity_type_for_entity_id(const std::string entityID) const;
        EntityType get_entity_type_for_asym_id(const std::string asymID) const;

        const std::list<polymer> &polymers() const { return m_polymers; }
        std::list<polymer> &polymers() { return m_polymers; }

        polymer &get_polymer_by_asym_id(const std::string &asymID);
        const polymer &get_polymer_by_asym_id(const std::string &asymID) const
        {
            return const_cast<structure *>(this)->get_polymer_by_asym_id(asymID);
        }

        const std::list<branch> &branches() const { return m_branches; }
        std::list<branch> &branches() { return m_branches; }

        branch &get_branch_by_asym_id(const std::string &asymID);
        const branch &get_branch_by_asym_id(const std::string &asymID) const;

        const std::vector<residue> &non_polymers() const { return m_non_polymers; }

        bool has_atom_id(const std::string &id) const;
        atom get_atom_by_id(const std::string &id) const;

        atom get_atom_by_label(const std::string &atomID, const std::string &asymID,
            const std::string &compID, int seqID, const std::string &altID = "");

        atom get_atom_by_position(point p) const;

        atom get_atom_by_position_and_type(point p, std::string_view type, std::string_view res_type) const;

        residue &create_residue(const std::vector<atom> &atoms);

        residue &get_residue(const std::string &asymID)
        {
            return get_residue(asymID, 0, "");
        }

        const residue &get_residue(const std::string &asymID) const
        {
            return get_residue(asymID, 0, "");
        }

        residue &get_residue(const std::string &asymID, int seqID, const std::string &authSeqID);

        const residue &get_residue(const std::string &asymID, int seqID, const std::string &authSeqID) const
        {
            return const_cast<structure *>(this)->get_residue(asymID, seqID, authSeqID);
        }

        residue &get_residue(const std::string &asymID, const std::string &compID, int seqID, const std::string &authSeqID);

        const residue &get_residue(const std::string &asymID, const std::string &compID, int seqID, const std::string &authSeqID) const
        {
            return const_cast<structure *>(this)->get_residue(asymID, compID, seqID, authSeqID);
        }

        residue &get_residue(const atom &atom)
        {
            return get_residue(atom.get_label_asym_id(), atom.get_label_comp_id(), atom.get_label_seq_id(), atom.get_auth_seq_id());
        }

        const residue &get_residue(const atom &atom) const
        {
            return get_residue(atom.get_label_asym_id(), atom.get_label_comp_id(), atom.get_label_seq_id(), atom.get_auth_seq_id());
        }

        // Actions. Originally a lot more actions were expected here

        void remove_atom(atom &a)
        {
            remove_atom(a, true);
        }

        void swap_atoms(atom a1, atom a2);
        void move_atom(atom a, point p);

        void change_residue(residue &res, const std::string &newcompound,
            const std::vector<std::tuple<std::string, std::string>> &remappedAtoms);

        void remove_residue(const std::string &asym_id, int seq_id, const std::string &auth_seq_id);

        std::string create_non_poly_entity(const std::string &mon_id);

        std::string create_non_poly(const std::string &entity_id, const std::vector<atom> &atoms);

        std::string create_non_poly(const std::string &entity_id, std::vector<row_initializer> atoms);

        std::string create_non_poly(const std::string &compound_id, bool skip_hydrogen);

        void create_water(row_initializer atom);

        std::string create_link(atom a1, atom a2, const std::string &link_type, const std::string &role);

        branch &create_branch();

        // /// \brief Create a new (sugar) branch with one first NAG containing atoms constructed from \a atoms
        // branch &create_branch(std::vector<row_initializer> atoms);

        // /// \brief Extend an existing (sugar) branch identified by \a asymID with one sugar containing atoms constructed from \a atom_info
        // ///
        // /// \param asym_id      The asym id of the branch to extend
        // /// \param atom_info    Array containing the info for the atoms to construct for the new sugar
        // /// \param link_sugar   The sugar to link to, note: this is the sugar number (1 based)
        // /// \param link_atom    The atom id of the atom linked in the sugar
        // branch &extend_branch(const std::string &asym_id, std::vector<row_initializer> atom_info,
        //  int link_sugar, const std::string &link_atom);

        void remove_branch(branch &branch);

        void remove_residue(residue &res);

        void translate(point t);

        void rotate(quaternion t);

        void translate_and_rotate(point t, quaternion q);

        void translate_rotate_and_translate(point t1, quaternion q, point t2);

        void cleanup_empty_categories();

        category &get_category(std::string_view name) const
        {
            return m_db[name];
        }

        datablock &get_datablock() const
        {
            return m_db;
        }

        void validate_atoms() const;

        template <typename... Args>
        atom &emplace_atom(Args &...args)
        {
            return emplace_atom(atom{ std::forward<Args>(args)... });
        }

        atom &emplace_atom(atom &&atom);

        void reorder_atoms();

      private:
        friend polymer;
        friend residue;

        void load_atoms_for_model(structure_open_options options);

        std::string insert_compound(const std::string &compoundID, bool is_entity);

        std::string create_entity_for_branch(branch &branch);

        void load_data();

        void remove_atom(atom &a, bool removeFromResidue);
        void remove_sugar(sugar &sugar);

        datablock &m_db;
        std::size_t m_model_nr;
        std::vector<atom> m_atoms;
        std::vector<std::size_t> m_atom_index;
        std::list<polymer> m_polymers;
        std::list<branch> m_branches;
        std::vector<residue> m_non_polymers;
    };

    } // namespace cif::mm