molecule.h

Go to the documentation of this file.

//
// Project NuriKit - Copyright 2023 SNU Compbio Lab.
// SPDX-License-Identifier: Apache-2.0
//
#ifndef NURI_CORE_MOLECULE_H_
#define NURI_CORE_MOLECULE_H_

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <ostream>
#include <string>
#include <string_view>
#include <type_traits>
#include <utility>
#include <vector>
 
#include <absl/algorithm/container.h>
#include <absl/base/attributes.h>
#include <absl/base/optimization.h>
#include <absl/container/fixed_array.h>
#include <absl/container/flat_hash_map.h>
#include <absl/log/absl_check.h>
#include <boost/container/flat_map.hpp>
#include <Eigen/Dense>
 
#include "nuri/eigen_config.h"
#include "nuri/core/container/index_set.h"
#include "nuri/core/container/property_map.h"
#include "nuri/core/element.h"
#include "nuri/core/graph/graph.h"
#include "nuri/core/graph/subgraph.h"
#include "nuri/meta.h"
#include "nuri/utils.h"
 
namespace nuri {
namespace constants {
  enum Hybridization : int {
    kUnbound = 0,   // Unbound
    kTerminal = 1,  // Terminal
    kSP = 2,
    kSP2 = 3,
    kSP3 = 4,
    kSP3D = 5,
    kSP3D2 = 6,
    kOtherHyb = 7,  // Unknown/other
  };
 
  inline std::ostream &operator<<(std::ostream &os, Hybridization hyb) {
    switch (hyb) {
    case kUnbound:
      return os << "unbound";
    case kTerminal:
      return os << "terminal";
    case kSP:
      return os << "sp";
    case kSP2:
      return os << "sp2";
    case kSP3:
      return os << "sp3";
    case kSP3D:
      return os << "sp3d";
    case kSP3D2:
      return os << "sp3d2";
    case kOtherHyb:
      break;
    }
 
    return os << "other";
  }

  enum BondOrder : int {
    kOtherBond = 0,
    kSingleBond = 1,
    kDoubleBond = 2,
    kTripleBond = 3,
    kQuadrupleBond = 4,
    kAromaticBond = 5,
  };
 
  inline std::ostream &operator<<(std::ostream &os, BondOrder bo) {
    switch (bo) {
    case kOtherBond:
      break;
    case kSingleBond:
      return os << "single";
    case kDoubleBond:
      return os << "double";
    case kTripleBond:
      return os << "triple";
    case kQuadrupleBond:
      return os << "quadruple";
    case kAromaticBond:
      return os << "aromatic";
    }
 
    return os << "other";
  }
 
  constexpr double kBondOrderToDouble[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 1.5 };
}  // namespace constants
 
inline constants::Hybridization clamp_hyb(int hyb) {
  return nuri::clamp(static_cast<constants::Hybridization>(hyb),
                     constants::kTerminal, constants::kSP3D2);
}
 
inline constants::BondOrder clamp_ord(int ord) {
  return nuri::clamp(static_cast<constants::BondOrder>(ord),
                     constants::kSingleBond, constants::kTripleBond);
}
 
enum class AtomFlags : std::uint32_t {
  kAromatic = 0x1,
  kConjugated = 0x2,
  kRing = 0x4,
  kChiral = 0x8,
  kClockWise = 0x10,
};
 
class AtomData {
public:
  AtomData(): AtomData(kPt[0]) { }
 
  AtomData(const Element &element, int implicit_hydrogens = 0,
           int formal_charge = 0,
           constants::Hybridization hyb = constants::kOtherHyb,
           double partial_charge = 0.0, int mass_number = -1,
           bool is_aromatic = false, bool is_in_ring = false,
           bool is_chiral = false, bool is_clockwise = false);

  int atomic_number() const { return element().atomic_number(); }

  double atomic_weight() const { return element().atomic_weight(); }

  std::string_view element_symbol() const { return element().symbol(); }

  std::string_view element_name() const { return element().name(); }

  AtomData &set_element(const Element &element) {
    element_ = &element;
    return *this;
  }

  AtomData &set_element(int atomic_number) {
    set_element(kPt[atomic_number]);
    return *this;
  }

  const Element &element() const noexcept {
    // GCOV_EXCL_START
    ABSL_ASSUME(element_ != nullptr);
    // GCOV_EXCL_STOP
    return *element_;
  }

  const Isotope &isotope() const {
    return ABSL_PREDICT_TRUE(isotope_ == nullptr) ? element().major_isotope()
                                                  : *isotope_;
  }

  AtomData &set_isotope(const Isotope &isotope) {
    isotope_ = &isotope;
    return *this;
  }

  AtomData &set_isotope(int mass_number) {
    isotope_ = element().find_isotope(mass_number);
    return *this;
  }

  internal::Nullable<const Isotope *> explicit_isotope() const {
    return isotope_;
  }
 
  AtomData &set_hybridization(constants::Hybridization hyb) {
    hyb_ = hyb;
    return *this;
  }
 
  constants::Hybridization hybridization() const { return hyb_; }

  AtomData &set_implicit_hydrogens(int implicit_hydrogens) {
    ABSL_DCHECK(implicit_hydrogens >= 0);
    implicit_hydrogens_ = implicit_hydrogens;
    return *this;
  }

  int implicit_hydrogens() const { return implicit_hydrogens_; }
 
  AtomData &set_aromatic(bool is_aromatic) {
    internal::update_flag(flags_, is_aromatic, AtomFlags::kAromatic);
    return *this;
  }
 
  bool is_aromatic() const {
    return internal::check_flag(flags_, AtomFlags::kAromatic);
  }
 
  AtomData &set_conjugated(bool is_conjugated) {
    internal::update_flag(flags_, is_conjugated, AtomFlags::kConjugated);
    return *this;
  }
 
  bool is_conjugated() const {
    return internal::check_flag(flags_, AtomFlags::kConjugated);
  }
 
  AtomData &set_ring_atom(bool is_ring_atom) {
    internal::update_flag(flags_, is_ring_atom, AtomFlags::kRing);
    return *this;
  }
 
  bool is_ring_atom() const {
    return internal::check_flag(flags_, AtomFlags::kRing);
  }
 
  AtomData &set_chiral(bool is_chiral) {
    internal::update_flag(flags_, is_chiral, AtomFlags::kChiral);
    return *this;
  }
 
  bool is_chiral() const {
    return internal::check_flag(flags_, AtomFlags::kChiral);
  }
 
  AtomData &set_clockwise(bool is_clockwise) {
    internal::update_flag(flags_, is_clockwise, AtomFlags::kClockWise);
    return *this;
  }

  bool is_clockwise() const {
    return internal::check_flag(flags_, AtomFlags::kClockWise);
  }
 
  AtomFlags flags() const { return flags_; }
 
  AtomData &add_flags(AtomFlags flags) {
    flags_ |= flags;
    return *this;
  }
 
  AtomData &del_flags(AtomFlags flags) {
    flags_ &= ~flags;
    return *this;
  }
 
  AtomData &reset_flags() {
    flags_ = static_cast<AtomFlags>(0);
    return *this;
  }
 
  AtomData &set_partial_charge(double charge) {
    partial_charge_ = charge;
    return *this;
  }
 
  double partial_charge() const { return partial_charge_; }
 
  AtomData &set_formal_charge(int charge) {
    formal_charge_ = charge;
    return *this;
  }
 
  int formal_charge() const { return formal_charge_; }
 
  std::string_view get_name() const { return internal::get_name(props_); }
 
  template <class ST>
  AtomData &set_name(ST &&name) {
    internal::set_name(props_, std::forward<ST>(name));
    return *this;
  }
 
  template <class KT, class VT>
  AtomData &add_prop(KT &&key, VT &&val) {
    internal::set_key(props_, std::forward<KT>(key), std::forward<VT>(val));
    return *this;
  }
 
  internal::PropertyMap &props() { return props_; }
 
  const internal::PropertyMap &props() const { return props_; }
 
private:
  friend bool operator==(const AtomData &lhs, const AtomData &rhs) noexcept;
 
  const Element *element_;
  int implicit_hydrogens_;
  int formal_charge_;
  constants::Hybridization hyb_;
  AtomFlags flags_;
  double partial_charge_;
  const Isotope *isotope_;
  internal::PropertyMap props_;
};
 
inline bool operator==(const AtomData &lhs, const AtomData &rhs) noexcept {
  return lhs.element() == rhs.element()
         && lhs.hybridization() == rhs.hybridization()
         && lhs.flags_ == rhs.flags_
         && lhs.formal_charge() == rhs.formal_charge();
}
 
enum class BondFlags : std::uint32_t {
  kRing = 0x1,
  kAromatic = 0x2,
  kConjugated = 0x4,
  kConfigSpecified = 0x8,
  kTransConfig = 0x10,
};
 
class BondData {
public:
  BondData(): BondData(constants::kOtherBond) { }
 
  explicit BondData(constants::BondOrder order)
      : order_(order), flags_(static_cast<BondFlags>(0)) { }

  constants::BondOrder order() const { return order_; }

  double approx_order() const { return constants::kBondOrderToDouble[order_]; }

  constants::BondOrder &order() { return order_; }
 
  BondData &set_order(constants::BondOrder order) {
    order_ = order;
    return *this;
  }

  bool is_rotatable() const {
    return order_ <= constants::kSingleBond
           && !internal::check_flag(flags_,
                                    BondFlags::kConjugated | BondFlags::kRing);
  }
 
  bool is_ring_bond() const {
    return internal::check_flag(flags_, BondFlags::kRing);
  }
 
  BondData &set_ring_bond(bool ring) {
    internal::update_flag(flags_, ring, BondFlags::kRing);
    return *this;
  }
 
  bool is_aromatic() const {
    return internal::check_flag(flags_, BondFlags::kAromatic);
  }
 
  BondData &set_aromatic(bool aromatic) {
    internal::update_flag(flags_, aromatic, BondFlags::kAromatic);
    return *this;
  }
 
  bool is_conjugated() const {
    return internal::check_flag(flags_, BondFlags::kConjugated);
  }
 
  BondData &set_conjugated(bool conj) {
    internal::update_flag(flags_, conj, BondFlags::kConjugated);
    return *this;
  }

  bool has_config() const {
    return internal::check_flag(flags_, BondFlags::kConfigSpecified);
  }

  BondData &set_config(bool config) {
    internal::update_flag(flags_, config, BondFlags::kConfigSpecified);
    return *this;
  }

  bool is_trans() const {
    return internal::check_flag(flags_, BondFlags::kTransConfig);
  }

  BondData &set_trans(bool trans) {
    internal::update_flag(flags_, trans, BondFlags::kTransConfig);
    return *this;
  }
 
  BondFlags flags() const { return flags_; }
 
  BondData &add_flags(BondFlags flags) {
    flags_ |= flags;
    return *this;
  }
 
  BondData &del_flags(BondFlags flags) {
    flags_ &= ~flags;
    return *this;
  }
 
  BondData &reset_flags() {
    flags_ = static_cast<BondFlags>(0);
    return *this;
  }
 
  std::string_view get_name() const { return internal::get_name(props_); }
 
  template <class ST>
  BondData &set_name(ST &&name) {
    internal::set_name(props_, std::forward<ST>(name));
    return *this;
  }
 
  template <class KT, class VT>
  BondData &add_prop(KT &&key, VT &&val) {
    internal::set_key(props_, std::forward<KT>(key), std::forward<VT>(val));
    return *this;
  }
 
  internal::PropertyMap &props() { return props_; }
 
  const internal::PropertyMap &props() const { return props_; }
 
private:
  constants::BondOrder order_;
  BondFlags flags_;
  internal::PropertyMap props_;
};
 
class Molecule;
class MoleculeMutator;
 
enum class SubstructCategory {
  kUnknown,
  kResidue,
  kChain,
};
 
namespace internal {
  template <bool is_const = false>
  class Substructure {
  public:
    using GraphType = const_if_t<is_const, Graph<AtomData, BondData>>;
 
    // Should use SubgraphOf<GraphType> here, but Subgraph was used directly for
    // better clangd autocompletion.
    using SubgraphType = Subgraph<AtomData, BondData, is_const>;
 
    using MutableAtom = typename SubgraphType::NodeRef;
    using Atom = typename SubgraphType::ConstNodeRef;
 
    using iterator = typename SubgraphType::iterator;
    using const_iterator = typename SubgraphType::const_iterator;
    using atom_iterator = iterator;
    using const_atom_iterator = const_iterator;
 
    using bond_iterator = typename SubgraphType::edge_iterator;
    using const_bond_iterator = typename SubgraphType::const_edge_iterator;
    using MutableBond = typename SubgraphType::EdgeRef;
    using Bond = typename SubgraphType::ConstEdgeRef;
 
    using BondsWrapper = typename SubgraphType::EdgesWrapper;
    using ConstBondsWrapper = typename SubgraphType::ConstEdgesWrapper;
 
    using MutableNeighbor = typename SubgraphType::AdjRef;
    using Neighbor = typename SubgraphType::ConstAdjRef;
    using neighbor_iterator = typename SubgraphType::adjacency_iterator;
    using const_neighbor_iterator =
        typename SubgraphType::const_adjacency_iterator;
 
    Substructure(const SubgraphType &sub,
                 SubstructCategory cat = SubstructCategory::kUnknown)
        : graph_(sub), cat_(cat) { }
 
    Substructure(SubgraphType &&sub,
                 SubstructCategory cat = SubstructCategory::kUnknown) noexcept
        : graph_(std::move(sub)), cat_(cat) { }
 
    Substructure(const SubgraphType &sub, const std::string &name,
                 SubstructCategory cat = SubstructCategory::kUnknown)
        : graph_(sub), name_(name), cat_(cat) { }
 
    Substructure(SubgraphType &&sub, std::string &&name,
                 SubstructCategory cat = SubstructCategory::kUnknown) noexcept
        : graph_(std::move(sub)), name_(std::move(name)), cat_(cat) { }
 
    template <bool other_const,
              std::enable_if_t<is_const && !other_const, int> = 0>
    Substructure(const Substructure<other_const> &other)
        : graph_(other.graph_), name_(other.name_), id_(other.id_),
          cat_(other.cat_), props_(other.props_) { }
 
    template <bool other_const,
              std::enable_if_t<is_const && !other_const, int> = 0>
    Substructure(Substructure<other_const> &&other) noexcept
        : graph_(std::move(other.graph_)), name_(std::move(other.name_)),
          id_(other.id_), cat_(other.cat_), props_(std::move(other.props_)) { }
 
    template <bool other_const,
              std::enable_if_t<is_const && !other_const, int> = 0>
    Substructure &operator=(const Substructure<other_const> &other) {
      graph_ = other.graph_;
      name_ = other.name_;
      id_ = other.id_;
      cat_ = other.cat_;
      props_ = other.props_;
      return *this;
    }
 
    template <bool other_const,
              std::enable_if_t<is_const && !other_const, int> = 0>
    Substructure &operator=(Substructure<other_const> &&other) noexcept {
      graph_ = std::move(other.graph_);
      name_ = std::move(other.name_);
      id_ = other.id_;
      cat_ = other.cat_;
      props_ = std::move(other.props_);
      return *this;
    }
 
    bool empty() const { return graph_.empty(); }
    int size() const { return graph_.size(); }
    int num_atoms() const { return graph_.num_nodes(); }
    int num_bonds() const { return graph_.num_edges(); }
    int count_heavy_atoms() const {
      return absl::c_count_if(graph_, [](Substructure::Atom atom) {
        return atom.data().atomic_number() != 1;
      });
    }
 
    void clear() noexcept {
      graph_.clear();
      name_.clear();
      id_ = 0;
      cat_ = SubstructCategory::kUnknown;
      props_.clear();
    }
 
    void clear_atoms() noexcept { graph_.clear(); }
 
    void update(IndexSet &&atoms, IndexSet &&bonds) {
      graph_.update(std::move(atoms), std::move(bonds));
    }
 
    void update_atoms(IndexSet &&atoms) noexcept {
      graph_.update_nodes(std::move(atoms));
    }
 
    void reserve_atoms(int n) { graph_.reserve_nodes(n); }
 
    void add_atom(int id) { graph_.add_node(id); }
 
    void add_atoms(const IndexSet &atoms, bool bonds = false) {
      if (bonds) {
        graph_.add_nodes_with_edges(atoms);
      } else {
        graph_.add_nodes(atoms);
      }
    }
 
    bool contains_atom(int id) const { return graph_.contains_node(id); }
    bool contains_atom(typename GraphType::ConstNodeRef atom) const {
      return graph_.contains_node(atom);
    }
 
    MutableAtom operator[](int idx) { return graph_[idx]; }
    Atom operator[](int idx) const { return graph_[idx]; }
 
    MutableAtom atom(int idx) { return graph_.node(idx); }
    Atom atom(int idx) const { return graph_.node(idx); }
 
    iterator find_atom(int id) { return graph_.find_node(id); }
    iterator find_atom(typename GraphType::ConstNodeRef atom) {
      return graph_.find_node(atom);
    }
 
    const_iterator find_atom(int id) const { return graph_.find_node(id); }
    const_iterator find_atom(typename GraphType::ConstNodeRef atom) const {
      return graph_.find_node(atom);
    }
 
    void erase_atom(int idx) { graph_.erase_node(idx); }
    void erase_atom(Atom atom) { graph_.erase_node(atom); }
 
    void erase_atoms(const_iterator begin, const_iterator end) {
      graph_.erase_nodes(begin, end);
    }
 
    void erase_atom_of(int id) { graph_.erase_node_of(id); }
 
    void erase_atom_of(typename GraphType::ConstNodeRef atom) {
      graph_.erase_node_of(atom);
    }
 
    template <class UnaryPred>
    void erase_atoms_if(UnaryPred &&pred) {
      graph_.erase_nodes_if(std::forward<UnaryPred>(pred));
    }
 
    iterator begin() { return graph_.begin(); }
    iterator end() { return graph_.end(); }
 
    const_iterator begin() const { return cbegin(); }
    const_iterator end() const { return cend(); }
 
    const_iterator cbegin() const { return graph_.cbegin(); }
    const_iterator cend() const { return graph_.cend(); }
 
    atom_iterator atom_begin() { return graph_.node_begin(); }
    atom_iterator atom_end() { return graph_.node_end(); }
 
    const_atom_iterator atom_begin() const { return atom_cbegin(); }
    const_atom_iterator atom_end() const { return atom_cend(); }
 
    const_atom_iterator atom_cbegin() const { return graph_.node_cbegin(); }
    const_atom_iterator atom_cend() const { return graph_.node_cend(); }
 
    const std::vector<int> &atom_ids() const { return graph_.node_ids(); }
 
    BondsWrapper bonds() { return graph_.edges(); }
    ConstBondsWrapper bonds() const { return graph_.edges(); }
 
    void clear_bonds() noexcept { graph_.clear_edges(); }
 
    void update_bonds(IndexSet &&bonds) noexcept {
      graph_.update_edges(std::move(bonds));
    }
 
    void refresh_bonds() { graph_.refresh_edges(); }
 
    void reserve_bonds(int n) { graph_.reserve_edges(n); }
 
    void add_bond(int id) { graph_.add_edge(id); }
 
    void add_bonds(const IndexSet &bonds) { graph_.add_edges(bonds); }
 
    bool contains_bond(int id) const { return graph_.contains_edge(id); }
    bool contains_bond(typename GraphType::ConstEdgeRef bond) const {
      return graph_.contains_edge(bond);
    }
 
    MutableBond bond(int idx) { return graph_.edge(idx); }
    Bond bond(int idx) const { return graph_.edge(idx); }
 
    bond_iterator find_bond(int id) { return graph_.find_edge(id); }
    bond_iterator find_bond(typename GraphType::ConstEdgeRef bond) {
      return graph_.find_edge(bond);
    }
 
    const_bond_iterator find_bond(int id) const { return graph_.find_edge(id); }
    const_bond_iterator find_bond(typename GraphType::ConstEdgeRef bond) const {
      return graph_.find_edge(bond);
    }
 
    void erase_bond(int idx) { graph_.erase_edge(idx); }
    void erase_bond(Bond bond) { graph_.erase_edge(bond); }
 
    void erase_bonds(const_bond_iterator begin, const_bond_iterator end) {
      graph_.erase_edges(begin, end);
    }
 
    void erase_bond_of(int id) { graph_.erase_edge_of(id); }
 
    void erase_bond_of(typename GraphType::ConstEdgeRef bond) {
      graph_.erase_edge_of(bond);
    }
 
    template <class UnaryPred>
    void erase_bonds_if(UnaryPred &&pred) {
      graph_.erase_edges_if(std::forward<UnaryPred>(pred));
    }
 
    bond_iterator bond_begin() { return graph_.edge_begin(); }
    bond_iterator bond_end() { return graph_.edge_end(); }
 
    const_bond_iterator bond_begin() const { return bond_cbegin(); }
    const_bond_iterator bond_end() const { return bond_cend(); }
 
    const_bond_iterator bond_cbegin() const { return graph_.edge_cbegin(); }
    const_bond_iterator bond_cend() const { return graph_.edge_cend(); }
 
    const std::vector<int> &bond_ids() const { return graph_.edge_ids(); }
 
    int degree(int id) const { return graph_.degree(id); }
 
    bond_iterator find_bond(Atom src, Atom dst) {
      return graph_.find_edge(src, dst);
    }
 
    const_bond_iterator find_bond(Atom src, Atom dst) const {
      return graph_.find_edge(src, dst);
    }
 
    bond_iterator find_bond(typename GraphType::ConstNodeRef src,
                            typename GraphType::ConstNodeRef dst) {
      return graph_.find_edge(src, dst);
    }
 
    const_bond_iterator find_bond(typename GraphType::ConstNodeRef src,
                                  typename GraphType::ConstNodeRef dst) const {
      return graph_.find_edge(src, dst);
    }
 
    neighbor_iterator find_neighbor(Atom src, Atom dst) {
      return graph_.find_adjacent(src, dst);
    }
 
    const_neighbor_iterator find_neighbor(Atom src, Atom dst) const {
      return graph_.find_adjacent(src, dst);
    }
 
    neighbor_iterator neighbor_begin(int id) { return graph_.adj_begin(id); }
    neighbor_iterator neighbor_end(int id) { return graph_.adj_end(id); }
 
    const_neighbor_iterator neighbor_begin(int id) const {
      return graph_.adj_begin(id);
    }
    const_neighbor_iterator neighbor_end(int id) const {
      return graph_.adj_end(id);
    }
 
    const_neighbor_iterator neighbor_cbegin(int id) const {
      return graph_.adj_cbegin(id);
    }
    const_neighbor_iterator neighbor_cend(int id) const {
      return graph_.adj_cend(id);
    }
 
    std::string &name() { return name_; }
 
    const std::string &name() const { return name_; }
 
    void set_id(int id) { id_ = id; }
 
    int id() const { return id_; }
 
    SubstructCategory &category() { return cat_; }
 
    SubstructCategory category() const { return cat_; }
 
    template <class KT, class VT>
    void add_prop(KT &&key, VT &&val) {
      internal::set_key(props_, std::forward<KT>(key), std::forward<VT>(val));
    }
 
    internal::PropertyMap &props() { return props_; }
 
    const internal::PropertyMap &props() const { return props_; }
 
  private:
    friend Molecule;
    friend MoleculeMutator;
 
    void rebind(typename SubgraphType::parent_type &parent) {
      graph_.rebind(parent);
    }
 
    SubgraphType graph_;
 
    std::string name_;
    int id_ = 0;
    SubstructCategory cat_;
    internal::PropertyMap props_;
  };
 
  extern std::pair<std::vector<int>, bool>
  place_trailing_hydrogens_initial(const Molecule &mol, Matrix3Xd &conf,
                                   int h_begin);
 
  extern bool optimize_free_hydrogens(const Molecule &mol, Matrix3Xd &conf,
                                      const std::vector<int> &free_hs);
}  // namespace internal
 
using Substructure = internal::Substructure<false>;
using ConstSubstructure = internal::Substructure<true>;

class Molecule {
public:
  using GraphType = Graph<AtomData, BondData>;
 
  using MutableAtom = GraphType::NodeRef;
  using Atom = GraphType::ConstNodeRef;
  using iterator = GraphType::iterator;
  using const_iterator = GraphType::const_iterator;
  using atom_iterator = iterator;
  using const_atom_iterator = const_iterator;
 
  using MutableBond = GraphType::EdgeRef;
  using Bond = GraphType::ConstEdgeRef;
  using bond_iterator = GraphType::edge_iterator;
  using const_bond_iterator = GraphType::const_edge_iterator;
 
  using MutableNeighbor = GraphType::AdjRef;
  using Neighbor = GraphType::ConstAdjRef;
  using neighbor_iterator = GraphType::adjacency_iterator;
  using const_neighbor_iterator = GraphType::const_adjacency_iterator;
 
  friend MoleculeMutator;

  Molecule() noexcept = default;
  ~Molecule() noexcept = default;
 
  Molecule(const Molecule &other) noexcept;
  Molecule(Molecule &&other) noexcept;
  Molecule &operator=(const Molecule &other) noexcept;
  Molecule &operator=(Molecule &&other) noexcept;

  template <class Iterator,
            class = internal::enable_if_compatible_iter_t<Iterator, AtomData>>
  Molecule(Iterator begin, Iterator end);
 
  std::string &name() { return name_; }
 
  const std::string &name() const { return name_; }
 
  void reserve(int num_atoms) { graph_.reserve(num_atoms); }
 
  void reserve_bonds(int num_bonds) { graph_.reserve_edges(num_bonds); }

  bool empty() const { return graph_.empty(); }

  int size() const { return graph_.size(); }

  int num_atoms() const { return graph_.num_nodes(); }

  int count_heavy_atoms() const {
    return absl::c_count_if(graph_, [](Molecule::Atom a) {
      return a.data().atomic_number() != 1;
    });
  }

  MutableAtom atom(int atom_idx) { return graph_.node(atom_idx); }

  Atom atom(int atom_idx) const { return graph_.node(atom_idx); }

  MutableAtom operator[](int atom_idx) { return graph_[atom_idx]; }

  Atom operator[](int atom_idx) const { return graph_[atom_idx]; }

  iterator begin() { return graph_.begin(); }
  const_iterator begin() const { return graph_.begin(); }

  iterator atom_begin() { return begin(); }
  const_iterator atom_begin() const { return begin(); }

  iterator end() { return graph_.end(); }
  const_iterator end() const { return graph_.end(); }

  iterator atom_end() { return end(); }
  const_iterator atom_end() const { return end(); }

  bool bond_empty() const { return graph_.edge_empty(); }

  int num_bonds() const { return graph_.num_edges(); }

  MutableBond bond(int bond_idx) { return graph_.edge(bond_idx); }

  Bond bond(int bond_idx) const { return graph_.edge(bond_idx); }

  bond_iterator find_bond(int src, int dst) {
    return graph_.find_edge(src, dst);
  }

  const_bond_iterator find_bond(int src, int dst) const {
    return graph_.find_edge(src, dst);
  }

  bond_iterator find_bond(Atom src, Atom dst) {
    return graph_.find_edge(src, dst);
  }

  const_bond_iterator find_bond(Atom src, Atom dst) const {
    return graph_.find_edge(src, dst);
  }

  auto bonds() { return graph_.edges(); }

  auto bonds() const { return graph_.edges(); }

  auto cbonds() const { return graph_.edges(); }

  bond_iterator bond_begin() { return graph_.edge_begin(); }
  const_bond_iterator bond_begin() const { return graph_.edge_begin(); }
  bond_iterator bond_end() { return graph_.edge_end(); }
  const_bond_iterator bond_end() const { return graph_.edge_end(); }

  int count_heavy_bonds() const {
    return absl::c_count_if(bonds(), [](Molecule::Bond bond) {
      return bond.src().data().atomic_number() != 1
             && bond.dst().data().atomic_number() != 1;
    });
  }

  int num_neighbors(int atom) const { return graph_.degree(atom); }

  neighbor_iterator find_neighbor(int src, int dst) {
    return graph_.find_adjacent(src, dst);
  }

  const_neighbor_iterator find_neighbor(int src, int dst) const {
    return graph_.find_adjacent(src, dst);
  }

  neighbor_iterator find_neighbor(Atom src, Atom dst) {
    return graph_.find_adjacent(src, dst);
  }

  const_neighbor_iterator find_neighbor(Atom src, Atom dst) const {
    return graph_.find_adjacent(src, dst);
  }

  neighbor_iterator neighbor_begin(int atom_idx) {
    return graph_.adj_begin(atom_idx);
  }

  const_neighbor_iterator neighbor_begin(int atom_idx) const {
    return graph_.adj_begin(atom_idx);
  }

  neighbor_iterator neighbor_end(int atom_idx) {
    return graph_.adj_end(atom_idx);
  }

  const_neighbor_iterator neighbor_end(int atom_idx) const {
    return graph_.adj_end(atom_idx);
  }

  MoleculeMutator mutator();

  void clear() noexcept;

  void clear_atoms() noexcept;

  void clear_bonds() noexcept;

  const GraphType &raw_graph() const { return graph_; }

  ABSL_MUST_USE_RESULT bool add_hydrogens(bool update_confs = true,
                                          bool optimize = true);

  void erase_hydrogens();

  bool is_3d() const { return !conformers_.empty(); }

  std::vector<Matrix3Xd> &confs() { return conformers_; }

  const std::vector<Matrix3Xd> &confs() const { return conformers_; }

  void transform(const Eigen::Isometry3d &trans) {
    for (Matrix3Xd &m: conformers_)
      m = trans * m;
  }

  void transform(int i, const Eigen::Isometry3d &trans) {
    Matrix3Xd &m = conformers_[i];
    m = trans * m;
  }

  double distsq(int src, int dst, int conf = 0) const;

  double distsq(Bond bond, int conf = 0) const {
    return distsq(bond.src().id(), bond.dst().id(), conf);
  }

  double distance(int src, int dst, int conf = 0) const {
    return std::sqrt(distsq(src, dst, conf));
  }

  double distance(Bond bond, int conf = 0) const {
    return distance(bond.src().id(), bond.dst().id(), conf);
  }

  ArrayXd bond_lengths(int conf = 0) const;

  bool rotate_bond(int ref_atom, int pivot_atom, double angle);

  bool rotate_bond(int bid, double angle);

  bool rotate_bond_conf(int i, int ref_atom, int pivot_atom, double angle);

  bool rotate_bond_conf(int i, int bid, double angle);

  Substructure
  substructure(SubstructCategory cat = SubstructCategory::kUnknown) {
    return { make_subgraph(graph_), cat };
  }

  Substructure
  substructure(internal::IndexSet &&atoms, internal::IndexSet &&bonds,
               SubstructCategory cat = SubstructCategory::kUnknown) {
    return { make_subgraph(graph_, std::move(atoms), std::move(bonds)), cat };
  }

  Substructure
  atom_substructure(internal::IndexSet &&atoms,
                    SubstructCategory cat = SubstructCategory::kUnknown) {
    return { subgraph_from_nodes(graph_, std::move(atoms)), cat };
  }

  Substructure bond_substructure(
      internal::IndexSet &&bonds,
      SubstructCategory cat = SubstructCategory::kUnknown) noexcept {
    return { subgraph_from_edges(graph_, std::move(bonds)), cat };
  }

  ConstSubstructure
  substructure(SubstructCategory cat = SubstructCategory::kUnknown) const {
    return { Subgraph(graph_), cat };
  }

  ConstSubstructure
  substructure(internal::IndexSet &&atoms, internal::IndexSet &&bonds,
               SubstructCategory cat = SubstructCategory::kUnknown) const {
    return { make_subgraph(graph_, std::move(atoms), std::move(bonds)), cat };
  }

  ConstSubstructure
  atom_substructure(internal::IndexSet &&atoms,
                    SubstructCategory cat = SubstructCategory::kUnknown) const {
    return { subgraph_from_nodes(graph_, std::move(atoms)), cat };
  }

  ConstSubstructure
  bond_substructure(internal::IndexSet &&bonds,
                    SubstructCategory cat = SubstructCategory::kUnknown) const {
    return { subgraph_from_edges(graph_, std::move(bonds)), cat };
  }

  std::vector<Substructure> &substructures() { return substructs_; }

  const std::vector<Substructure> &substructures() const { return substructs_; }

  void update_topology();

  int num_sssr() const { return num_bonds() - num_atoms() + num_fragments(); }

  const std::vector<std::vector<int>> &ring_groups() const {
    return ring_groups_;
  }

  int num_fragments() const { return num_fragments_; }

  template <class MoleculeLike>
  void merge(const MoleculeLike &other) {
    graph_.merge(other.graph_);
    update_topology();
 
    for (Matrix3Xd &conf: conformers_)
      conf.conservativeResize(Eigen::NoChange, size());
  }
 
  template <class KT, class VT>
  void add_prop(KT &&key, VT &&val) {
    internal::set_key(props_, std::forward<KT>(key), std::forward<VT>(val));
  }
 
  internal::PropertyMap &props() { return props_; }
 
  const internal::PropertyMap &props() const { return props_; }
 
private:
  Molecule(GraphType &&graph, std::vector<Matrix3Xd> &&conformers) noexcept
      : graph_(std::move(graph)), conformers_(std::move(conformers)) { }
 
  void rebind_substructs() noexcept;
 
  bool rotate_bond_common(int i, int ref_atom, int pivot_atom, double angle);
 
  GraphType graph_;
  std::vector<Matrix3Xd> conformers_;
  std::string name_;
  internal::PropertyMap props_;
 
  std::vector<Substructure> substructs_;
 
  std::vector<std::vector<int>> ring_groups_;
  int num_fragments_ = 0;
};

class MoleculeMutator {
public:
  MoleculeMutator(Molecule &mol)
      : mol_(&mol), prev_num_atoms_(mol.num_atoms()) { }
 
  MoleculeMutator() = delete;
  MoleculeMutator(const MoleculeMutator &) = delete;
  MoleculeMutator &operator=(const MoleculeMutator &) = delete;
  MoleculeMutator(MoleculeMutator &&) noexcept = default;
  MoleculeMutator &operator=(MoleculeMutator &&) noexcept = default;
 
  ~MoleculeMutator() noexcept { finalize(); }

  int add_atom(const AtomData &atom) { return mol().graph_.add_node(atom); }

  int add_atom(AtomData &&atom) noexcept {
    return mol().graph_.add_node(std::move(atom));
  }

  int add_atom(const Molecule::Atom &atom) { return add_atom(atom.data()); }

  void mark_atom_erase(int atom_idx) { erased_atoms_.push_back(atom_idx); }

  void mark_atom_erase(Molecule::Atom atom) { mark_atom_erase(atom.id()); }

  void clear_atoms() noexcept;

  std::pair<int, bool> register_bond(int src, int dst, const BondData &bond);

  std::pair<int, bool> register_bond(int src, int dst,
                                     BondData &&bond) noexcept;

  std::pair<int, bool> register_bond(int src, int dst, Molecule::Bond bond) {
    return register_bond(src, dst, bond.data());
  }

  void mark_bond_erase(int bid) { erased_bonds_.push_back(bid); }

  void mark_bond_erase(int src, int dst);

  void mark_bond_erase(Molecule::Atom src, Molecule::Atom dst) {
    mark_bond_erase(src.id(), dst.id());
  }

  void clear_bonds() noexcept;

  void clear() noexcept;

  void finalize() noexcept;
 
  // GCOV_EXCL_START
  Molecule &mol() noexcept {
    ABSL_ASSUME(mol_ != nullptr);
    return *mol_;
  }
 
  const Molecule &mol() const noexcept {
    ABSL_ASSUME(mol_ != nullptr);
    return *mol_;
  }
  // GCOV_EXCL_STOP
 
private:
  void discard_bonds() noexcept;
 
  void discard() noexcept;
 
  Molecule *mol_;
  int prev_num_atoms_;
 
  std::vector<Molecule::GraphType::StoredEdge> bond_registry_;
  absl::flat_hash_map<int, std::vector<std::pair<int, int>>> delta_adj_;
 
  std::vector<int> erased_atoms_;
  std::vector<int> erased_bonds_;
};
 
class MoleculeSanitizer {
public:
  explicit MoleculeSanitizer(Molecule &molecule);
 
  MoleculeSanitizer(MoleculeSanitizer &&) noexcept = default;
  ~MoleculeSanitizer() noexcept = default;
 
  MoleculeSanitizer() = delete;
  MoleculeSanitizer(const MoleculeSanitizer &) = delete;
  MoleculeSanitizer &operator=(const MoleculeSanitizer &) = delete;
  MoleculeSanitizer &operator=(MoleculeSanitizer &&) noexcept = delete;

  ABSL_MUST_USE_RESULT bool sanitize_conjugated();

  ABSL_MUST_USE_RESULT bool sanitize_aromaticity();

  ABSL_MUST_USE_RESULT bool sanitize_hybridization();

  ABSL_MUST_USE_RESULT bool sanitize_valence();

  ABSL_MUST_USE_RESULT bool sanitize_all() {
    return sanitize_conjugated() && sanitize_aromaticity()
           && sanitize_hybridization() && sanitize_valence();
  }
 
  Molecule &mol() noexcept {
    ABSL_ASSUME(mol_ != nullptr);
    return *mol_;
  }
 
  const Molecule &mol() const noexcept {
    ABSL_ASSUME(mol_ != nullptr);
    return *mol_;
  }
 
private:
  Molecule *mol_;
  absl::FixedArray<int> valences_;
};
 
/* Out-of-line definitions for molecule */
 
template <class Iterator, class>
Molecule::Molecule(Iterator begin, Iterator end): Molecule() {
  MoleculeMutator m = mutator();
  for (auto it = begin; it != end; ++it) {
    m.add_atom(*it);
  }
}
 
inline MoleculeMutator Molecule::mutator() {
  return MoleculeMutator(*this);
}
 
/* Utility functions */
 
namespace internal {
  inline int common_valence(const Element &effective) {
    const int val_electrons = effective.valence_electrons(),
              common_valence = val_electrons <= 4 ? val_electrons
                                                  : 8 - val_electrons;
    return common_valence;
  }
 
  extern const Element &
  effective_element_or_element(const AtomData &data) noexcept;
 
  inline const Element &
  effective_element_or_element(Molecule::Atom atom) noexcept {
    return effective_element_or_element(atom.data());
  }
 
  extern int sum_bond_order_raw(Molecule::Atom atom, int implicit_hydrogens,
                                bool aromatic_correct);
 
  inline int sum_bond_order(Molecule::Atom atom, bool aromatic_correct) {
    return sum_bond_order_raw(atom, atom.data().implicit_hydrogens(),
                              aromatic_correct);
  }
 
  extern int steric_number(int total_degree, int nb_electrons);
 
  extern constants::Hybridization from_degree(int total_degree,
                                              int nb_electrons);
 
  extern int nonbonding_electrons(const AtomData &data, int total_valence);
 
  extern int count_pi_e(Molecule::Atom atom, int total_valence);
 
  extern int aromatic_pi_e(Molecule::Atom atom, int total_valence);
}  // namespace internal

inline int all_neighbors(Molecule::Atom atom) {
  return atom.degree() + atom.data().implicit_hydrogens();
}

extern int count_heavy(Molecule::Atom atom);

extern int count_hydrogens(Molecule::Atom atom);

inline int sum_bond_order(Molecule::Atom atom) {
  return internal::sum_bond_order(atom, true);
}

inline int nonbonding_electrons(Molecule::Atom atom) {
  return internal::nonbonding_electrons(atom.data(), sum_bond_order(atom));
}

inline int steric_number(Molecule::Atom atom) {
  int nbe = nonbonding_electrons(atom);
  if (nbe < 0)
    return nbe;
 
  return internal::steric_number(all_neighbors(atom), nbe);
}

inline const Element *effective_element(const AtomData &data) {
  const int effective_z = data.atomic_number() - data.formal_charge();
  return kPt.find_element(effective_z);
}

inline const Element *effective_element(Molecule::Atom atom) {
  return effective_element(atom.data());
}

extern std::vector<std::vector<int>> fragments(const Molecule &mol);
}  // namespace nuri
 
#endif /* NURI_CORE_MOLECULE_H_ */