/*
 * Read a molecule from CIF
 *
 * (c) 2014 Schrodinger, Inc.
 */

#include <algorithm>
#include <complex>
#include <string>
#include <map>
#include <set>
#include <vector>
#include <memory>
#include <array>

#include <glm/glm.hpp>
#include <glm/gtc/type_ptr.hpp>

#ifdef TRACE
#include <iostream>
#include <iomanip>
#include <glm/geometric.hpp> // For glm::length, glm::dot, glm::normalize, glm::distance
#include <glm/gtx/string_cast.hpp> // For printing glm::vec3
#endif // TRACE

#include "os_predef.h"
#include "os_std.h"

#include "MemoryDebug.h"
#include "Err.h"

#include "AssemblyHelpers.h"
#include "AtomInfo.h"
#include "Base.h"
#include "Executive.h"
#include "P.h"
#include "Util.h"
#include "Scene.h"
#include "Rep.h"
#include "ObjectMolecule.h"
#include "ObjectMap.h"
#include "CifFile.h"
#include "CifBondDict.h"
#include "Util2.h"
#include "Vector.h"
#include "Lex.h"
#include "strcasecmp.h"
#include "pymol/zstring_view.h"
#include "Feedback.h"

#include "pocketfft_hdronly.h"

#ifdef _PYMOL_IP_PROPERTIES
#endif

using pymol::cif_array;

constexpr auto EPSILON = std::numeric_limits<float>::epsilon();

// canonical amino acid three letter codes
const char * aa_three_letter[] = {
  "ALA", // A
  "ASX", // B  for ambiguous asparagine/aspartic-acid
  "CYS", // C
  "ASP", // D
  "GLU", // E
  "PHE", // F
  "GLY", // G
  "HIS", // H
  "ILE", // I
  nullptr,  // J
  "LYS", // K
  "LEU", // L
  "MET", // M
  "ASN", // N
  "HOH", // O  for water
  "PRO", // P
  "GLN", // Q
  "ARG", // R
  "SER", // S
  "THR", // T
  nullptr,  // U
  "VAL", // V
  "TRP", // W
  nullptr,  // X  for other
  "TYR", // Y
  "GLX", // Z  for ambiguous glutamine/glutamic acid
};

// amino acid one-to-three letter code translation
static const char * aa_get_three_letter(char aa) {
  if (aa < 'A' || aa > 'Z')
    return "UNK";
  const char * three = aa_three_letter[aa - 'A'];
  return (three) ? three : "UNK";
}

// dictionary content types
enum CifDataType {
  CIF_UNKNOWN,
  CIF_CORE,       // small molecule
  CIF_MMCIF,      // macromolecular structure
  CIF_CHEM_COMP   // chemical component
};

// simple 1-indexed string storage
class seqvec_t : public std::vector<std::string> {
public:
  void set(int i, const char * mon_id) {
    if (i < 1) {
      printf("error: i(%d) < 1\n", i);
      return;
    }
    if (i > size())
      resize(i);
    (*this)[i - 1] = mon_id;
  }

  const char * get(int i) const {
    if (i < 1 || i > size())
      return nullptr;
    return (*this)[i - 1].c_str();
  }
};

// structure to collect information about a data block
struct CifContentInfo {
  PyMOLGlobals * G;
  CifDataType type;
  bool fractional;
  bool use_auth;
  std::set<lexidx_t> chains_filter;
  std::set<std::string> polypeptide_entities; // entity ids
  std::map<std::string, seqvec_t> sequences;  // entity_id -> [resn1, resn2, ...]

  bool is_excluded_chain(const char * chain) const {
    if (chains_filter.empty())
      return false;

    auto borrowed = LexBorrow(G, chain);
    if (borrowed != LEX_BORROW_NOTFOUND)
      return is_excluded_chain(borrowed);

    return false;
  }

  bool is_excluded_chain(const lexborrow_t& chain) const {
    return (!chains_filter.empty() &&
        chains_filter.count(reinterpret_cast<const lexidx_t&>(chain)) == 0);
  }

  bool is_polypeptide(const char * entity_id) const {
    return polypeptide_entities.count(entity_id);
  }

  CifContentInfo(PyMOLGlobals * G, bool use_auth=true) :
    G(G),
    type(CIF_UNKNOWN),
    fractional(false),
    use_auth(use_auth) {}
};

/**
 * Make a string key that represents the collection of alt id, asym id,
 * atom id, comp id and seq id components of the label for a macromolecular
 * atom site.
 */
static std::string make_mm_atom_site_label(PyMOLGlobals * G, const AtomInfoType * a) {
  char resi[8];
  AtomResiFromResv(resi, sizeof(resi), a);

  std::string key(LexStr(G, a->chain));
  key += '/';
  key += LexStr(G, a->resn);
  key += '/';
  key += resi;
  key += '/';
  key += LexStr(G, a->name);
  key += '/';
  key += a->alt;
  return key;
}
static std::string make_mm_atom_site_label(PyMOLGlobals * G, const char * asym_id,
    const char * comp_id, const char * seq_id, const char * ins_code,
    const char * atom_id, const char * alt_id) {
  std::string key(asym_id);
  key += '/';
  key += comp_id;
  key += '/';
  key += seq_id;
  key += ins_code;
  key += '/';
  key += atom_id;
  key += '/';
  key += alt_id;
  return key;
}

/**
 * Like strncpy, but only copy alphabetic characters.
 */
static void strncpy_alpha(char* dest, const char* src, size_t n)
{
  for (size_t i = 0; i != n; ++i) {
    if (!isalpha(src[i])) {
      memset(dest + i, 0, n - i);
      break;
    }
    dest[i] = src[i];
  }
}

/**
 * Get first non-NULL element
 */
template <typename T>
static T VLAGetFirstNonNULL(T * vla) {
  int n = VLAGetSize(vla);
  for (int i = 0; i < n; ++i)
    if (vla[i])
      return vla[i];
  return nullptr;
}

/**
 * Lookup one key in a map, return true if found and
 * assign output reference `value1`
 */
template <typename Map, typename Key, typename T>
inline bool find1(Map& dict, T& value1, const Key& key1) {
  auto it = dict.find(key1);
  if (it == dict.end())
    return false;
  value1 = it->second;
  return true;
}

/**
 * Lookup two keys in a map, return true if both found and
 * assign output references `value1` and `value2`.
 */
template <typename Map, typename Key, typename T>
inline bool find2(Map& dict,
    T& value1, const Key& key1,
    T& value2, const Key& key2) {
  if (!find1(dict, value1, key1))
    return false;
  if (!find1(dict, value2, key2))
    return false;
  return true;
}

static void AtomInfoSetEntityId(PyMOLGlobals * G, AtomInfoType * ai, const char * entity_id) {
  ai->custom = LexIdx(G, entity_id);

#ifdef _PYMOL_IP_PROPERTIES
  PropertySet(G, ai, "entity_id", entity_id);
#endif
}

/**
 * Initialize a bond. Only one of symmetry_1 or symmetry_2 must be non-default.
 * If symmetry_2 is default and symmetry_1 is non-default, then swap the
 * indices.
 */
static bool BondTypeInit3(PyMOLGlobals* G, BondType* bond, unsigned i1,
    unsigned i2, const char* symmetry_1, const char* symmetry_2, int order = 1)
{
  auto symop_1 = pymol::SymOp(symmetry_1);
  auto symop_2 = pymol::SymOp(symmetry_2);

  if (symop_1) {
    if (symop_2) {
      PRINTFB(G, FB_Executive, FB_Warnings)
      " Warning: Bonds with two symmetry operations not supported\n" ENDFB(G);
      return false;
    }

    std::swap(i1, i2);
    std::swap(symop_1, symop_2);
  }

  BondTypeInit2(bond, i1, i2, order);
  bond->symop_2 = symop_2;

  return true;
}

/**
 * Add one bond without checking if it already exists
 */
static void ObjectMoleculeAddBond2(ObjectMolecule * I, int i1, int i2, int order) {
  VLACheck(I->Bond, BondType, I->NBond);
  BondTypeInit2(I->Bond + I->NBond, i1, i2, order);
  I->NBond++;
}

/**
 * Get the distance between two atoms in ObjectMolecule
 */
static float GetDistance(ObjectMolecule * I, int i1, int i2) {
  const CoordSet *cset;
  int idx1 = -1, idx2 = -1;

  // find first coordset which contains both atoms
  if (I->DiscreteFlag) {
    cset = I->DiscreteCSet[i1];
    if (cset == I->DiscreteCSet[i2]) {
      idx1 = I->DiscreteAtmToIdx[i1];
      idx2 = I->DiscreteAtmToIdx[i2];
    }
  } else {
    for (int i = 0; i < I->NCSet; ++i) {
      if ((cset = I->CSet[i])) {
        if ((idx1 = cset->AtmToIdx[i1]) != -1 &&
            (idx2 = cset->AtmToIdx[i2]) != -1) {
          break;
        }
      }
    }
  }

  if (idx1 == -1 || idx2 == -1)
    return 999.f;

  float v[3];
  subtract3f(
      cset->coordPtr(idx1),
      cset->coordPtr(idx2), v);
  return length3f(v);
}

/**
 * Bond order string to int
 */
static int bondOrderLookup(const char * order) {
  if (p_strcasestartswith(order, "doub"))
    return 2;
  if (p_strcasestartswith(order, "trip"))
    return 3;
  if (p_strcasestartswith(order, "arom"))
    return 4;
  if (p_strcasestartswith(order, "delo"))
    return 4;
  // single
  return 1;
}

/**
 * Read bonds from CHEM_COMP_BOND in `bond_dict` dictionary
 */
static bool read_chem_comp_bond_dict(const pymol::cif_data * data, bond_dict_t &bond_dict) {
  const cif_array *arr_id_1, *arr_id_2, *arr_order, *arr_comp_id;

  if( !(arr_id_1  = data->get_arr("_chem_comp_bond.atom_id_1")) ||
      !(arr_id_2  = data->get_arr("_chem_comp_bond.atom_id_2")) ||
      !(arr_order = data->get_arr("_chem_comp_bond.value_order")) ||
      !(arr_comp_id = data->get_arr("_chem_comp_bond.comp_id"))) {

    if ((arr_comp_id = data->get_arr("_chem_comp_atom.comp_id"))) {
      // atom(s) but no bonds (e.g. metals)
      bond_dict.set_unknown(arr_comp_id->as_s());
      return true;
    }

    return false;
  }

  const char *name1, *name2, *resn;
  int order_value;
  int nrows = arr_id_1->size();

  for (int i = 0; i < nrows; i++) {
    resn = arr_comp_id->as_s(i);
    name1 = arr_id_1->as_s(i);
    name2 = arr_id_2->as_s(i);

    const char *order = arr_order->as_s(i);
    order_value = bondOrderLookup(order);

    bond_dict[resn].set(name1, name2, order_value);
  }

  // alt_atom_id -> atom_id
  if ((arr_comp_id = data->get_arr("_chem_comp_atom.comp_id")) &&
      (arr_id_1 = data->get_arr("_chem_comp_atom.atom_id")) &&
      (arr_id_2 = data->get_arr("_chem_comp_atom.alt_atom_id"))) {
    nrows = arr_id_1->size();

    // set of all non-alt ids
    std::set<pymol::zstring_view> atom_ids;
    for (int i = 0; i < nrows; ++i) {
      atom_ids.insert(arr_id_1->as_s(i));
    }

    for (int i = 0; i < nrows; ++i) {
      resn = arr_comp_id->as_s(i);
      name1 = arr_id_1->as_s(i);
      name2 = arr_id_2->as_s(i);

      // skip identity mapping
      if (strcmp(name1, name2) == 0) {
        continue;
      }

      // alt id must not also be a non-alt id (PYMOL-3470)
      if (atom_ids.count(name2)) {
        fprintf(stderr,
            "Warning: _chem_comp_atom.alt_atom_id %s/%s ignored for bonding\n",
            resn, name2);
        continue;
      }

      bond_dict[resn].add_alt_name(name1, name2);
    }
  }

  return true;
}

/**
 * parse $PYMOL_DATA/chem_comp_bond-top100.cif (subset of components.cif) into
 * a static (global) dictionary.
 */
static bond_dict_t * get_global_components_bond_dict(PyMOLGlobals * G) {
  static bond_dict_t bond_dict;

  if (bond_dict.empty()) {
    const char * pymol_data = getenv("PYMOL_DATA");
    if (!pymol_data || !pymol_data[0])
      return nullptr;

    std::string path(pymol_data);
    path.append(PATH_SEP).append("chem_comp_bond-top100.cif");
    cif_file_with_error_capture cif;
    if (!cif.parse_file(path.c_str())) {
      PRINTFB(G, FB_Executive, FB_Warnings)
        " Warning: Loading '%s' failed: %s\n", path.c_str(),
        cif.m_error_msg.c_str() ENDFB(G);
      return nullptr;
    }

    for (const auto& [code, datablock] : cif.datablocks()) {
      read_chem_comp_bond_dict(&datablock, bond_dict);
    }
  }

  return &bond_dict;
}

/**
 * True for N-H1 and N-H3, those are not in the chemical components dictionary.
 */
static bool is_N_H1_or_H3(PyMOLGlobals * G,
    const AtomInfoType * a1,
    const AtomInfoType * a2) {
  if (a2->name == G->lex_const.N) {
    a2 = a1;
  } else if (a1->name != G->lex_const.N) {
    return false;
  }

  return (a2->name == G->lex_const.H1 || a2->name == G->lex_const.H3);
}

/**
 * Add bonds for one residue, with atoms spanning from i_start to i_end-1,
 * based on components.cif
 */
static void ConnectComponent(ObjectMolecule * I, int i_start, int i_end,
    bond_dict_t * bond_dict) {

  if (i_end - i_start < 2)
    return;

  auto G = I->G;
  const AtomInfoType *a1, *a2, *ai = I->AtomInfo.data();
  int order;

  // get residue bond dictionary
  auto res_dict = bond_dict->get(G, LexStr(G, ai[i_start].resn));
  if (res_dict == nullptr)
    return;

  // for all pairs of atoms in given set
  for (int i1 = i_start + 1; i1 < i_end; i1++) {
    for (int i2 = i_start; i2 < i1; i2++) {
      a1 = ai + i1;
      a2 = ai + i2;

      // don't connect different alt codes
      if (a1->alt[0] && a2->alt[0] && strcmp(a1->alt, a2->alt) != 0) {
        continue;
      }

      // restart if we hit the next residue in bulk solvent (atoms must
      // not be sorted for this)
      // TODO artoms are sorted at this point
      if (a1->name == a2->name) {
        i_start = i1;
        break;
      }

      // lookup if atoms are bonded
      order = res_dict->get(LexStr(G, a1->name), LexStr(G, a2->name));
      if (order < 0) {
        if (!is_N_H1_or_H3(G, a1, a2) || GetDistance(I, i1, i2) > 1.2)
          continue;

        order = 1;
      }

      // make bond
      ObjectMoleculeAddBond2(I, i1, i2, order);
    }
  }
}

/**
 * Add intra residue bonds based on components.cif, and common polymer
 * connecting bonds (C->N, O3*->P)
 */
static int ObjectMoleculeConnectComponents(ObjectMolecule * I,
    bond_dict_t * bond_dict=nullptr) {

  PyMOLGlobals * G = I->G;
  int i_start = 0;
  std::vector<int> i_prev_c[2], i_prev_o3[2];

  const lexborrow_t lex_O3s = LexBorrow(G, "O3*");
  const lexborrow_t lex_O3p = LexBorrow(G, "O3'");

  if (!bond_dict) {
    // read components.cif
    if (!(bond_dict = get_global_components_bond_dict(G)))
      return false;
  }

  // reserve some memory for new bonds
  I->Bond.reserve(I->NAtom * 4);

  for (int i = 0; i < I->NAtom; ++i) {
    auto const& atom = I->AtomInfo[i];

    // intra-residue
    if(!AtomInfoSameResidue(G, I->AtomInfo + i_start, I->AtomInfo + i)) {
      ConnectComponent(I, i_start, i, bond_dict);
      i_start = i;
      i_prev_c[0] = std::move(i_prev_c[1]);
      i_prev_o3[0] = std::move(i_prev_o3[1]);
      i_prev_c[1].clear();
      i_prev_o3[1].clear();
    }

    // inter-residue polymer bonds
    if (atom.name == G->lex_const.C) {
      i_prev_c[1].push_back(i);
    } else if (atom.name == lex_O3s || atom.name == lex_O3p) {
      i_prev_o3[1].push_back(i);
    } else {
      auto const* i_prev_ptr =
        (atom.name == G->lex_const.N) ? i_prev_c :
        (atom.name == G->lex_const.P) ? i_prev_o3 : nullptr;

      if (i_prev_ptr && !i_prev_ptr->empty()) {
        for (int i_prev : *i_prev_ptr) {
          bool alt_check = !atom.alt[0] || !I->AtomInfo[i_prev].alt[0] ||
                           atom.alt[0] == I->AtomInfo[i_prev].alt[0];
          if (alt_check && GetDistance(I, i_prev, i) < 1.8) {
            // make bond
            ObjectMoleculeAddBond2(I, i_prev, i, 1);
          }
        }
      }
    }
  }

  // final residue
  ConnectComponent(I, i_start, I->NAtom, bond_dict);

  // clean up
  VLASize(I->Bond, BondType, I->NBond);

  return true;
}

/**
 * secondary structure hash
 */
class sshashkey {
public:
  lexborrow_t chain; // borrowed ref
  int resv;
  char inscode;

  void assign(const lexborrow_t& asym_id_, int resv_, char ins_code_ = '\0') {
    chain = asym_id_;
    resv = resv_;
    inscode = ins_code_;
  }

  // comparable to sshashkey and AtomInfoType
  template <typename T> int compare(const T &other) const {
    int test = resv - other.resv;
    if (test == 0) {
      test = (chain - other.chain);
      if (test == 0)
        test = inscode - other.inscode;
    }
    return test;
  }
  bool operator<(const sshashkey &other) const { return compare(other) < 0; }
  bool operator>(const sshashkey &other) const { return compare(other) > 0; }
};
class sshashvalue {
public:
  char ss;
  sshashkey end;
};
typedef std::map<sshashkey, sshashvalue> sshashmap;

// PDBX_STRUCT_OPER_LIST type
typedef std::map<std::string, std::array<float, 16> > oper_list_t;

// type for parsed PDBX_STRUCT_OPER_LIST
typedef std::vector<std::vector<std::string> > oper_collection_t;

/**
 * Parse operation expressions like (1,2)(3-6)
 */
static oper_collection_t parse_oper_expression(const std::string &expr) {
  oper_collection_t collection;

  // first step to split parenthesized chunks
  std::vector<std::string> a_vec = strsplit(expr, ')');

  // loop over chunks (still include leading '(')
  for (auto& a_item : a_vec) {
    const char * a_chunk = a_item.c_str();

    // finish chunk
    while (*a_chunk == '(')
      ++a_chunk;

    // skip empty chunks
    if (!*a_chunk)
      continue;

    collection.resize(collection.size() + 1);
    oper_collection_t::reference ids = collection.back();

    // split chunk by commas
    std::vector<std::string> b_vec = strsplit(a_chunk, ',');

    // look for ranges
    for (auto& b_item : b_vec) {
      // "c_d" will have either one (no range) or two items
      std::vector<std::string> c_d = strsplit(b_item, '-');

      ids.push_back(c_d[0]);

      if (c_d.size() == 2)
        for (int i = atoi(c_d[0].c_str()) + 1,
                 j = atoi(c_d[1].c_str()) + 1; i < j; ++i)
        {
          char i_str[16];
          snprintf(i_str, sizeof(i_str), "%d", i);
          ids.push_back(i_str);
        }
    }
  }

  return collection;
}

/**
 * Get chains which are part of the assembly
 *
 * assembly_chains: output set
 * assembly_id: ID of the assembly or nullptr to use first assembly
 */
static bool get_assembly_chains(PyMOLGlobals * G,
    const pymol::cif_data * data,
    std::set<lexidx_t> &assembly_chains,
    const char * assembly_id) {

  const cif_array *arr_id, *arr_asym_id_list;

  if ((arr_id           = data->get_arr("_pdbx_struct_assembly_gen.assembly_id")) == nullptr ||
      (arr_asym_id_list = data->get_arr("_pdbx_struct_assembly_gen.asym_id_list")) == nullptr)
    return false;

  for (unsigned i = 0, nrows = arr_id->size(); i < nrows; ++i) {
    if (strcmp(assembly_id, arr_id->as_s(i)))
      continue;

    const char * asym_id_list = arr_asym_id_list->as_s(i);
    std::vector<std::string> chains = strsplit(asym_id_list, ',');
    for (auto& chain : chains) {
      assembly_chains.insert(LexIdx(G, chain.c_str()));
    }
  }

  return !assembly_chains.empty();
}

/**
 * Read assembly
 *
 * atInfo: atom info array to use for chain check
 * cset: template coordinate set to create assembly coordsets from
 * assembly_id: assembly identifier
 *
 * return: assembly coordinates as VLA of coordinate sets
 */
static
CoordSet ** read_pdbx_struct_assembly(PyMOLGlobals * G,
    const pymol::cif_data * data,
    const AtomInfoType * atInfo,
    const CoordSet * cset,
    const char * assembly_id) {

  const cif_array *arr_id, *arr_assembly_id, *arr_oper_expr, *arr_asym_id_list;


  if ((arr_id           = data->get_arr("_pdbx_struct_oper_list.id")) == nullptr ||
      (arr_assembly_id  = data->get_arr("_pdbx_struct_assembly_gen.assembly_id")) == nullptr ||
      (arr_oper_expr    = data->get_arr("_pdbx_struct_assembly_gen.oper_expression")) == nullptr ||
      (arr_asym_id_list = data->get_arr("_pdbx_struct_assembly_gen.asym_id_list")) == nullptr)
    return nullptr;

  const cif_array * arr_matrix[] = {
    data->get_opt("_pdbx_struct_oper_list.matrix[1][1]"),
    data->get_opt("_pdbx_struct_oper_list.matrix[1][2]"),
    data->get_opt("_pdbx_struct_oper_list.matrix[1][3]"),
    data->get_opt("_pdbx_struct_oper_list.vector[1]"),
    data->get_opt("_pdbx_struct_oper_list.matrix[2][1]"),
    data->get_opt("_pdbx_struct_oper_list.matrix[2][2]"),
    data->get_opt("_pdbx_struct_oper_list.matrix[2][3]"),
    data->get_opt("_pdbx_struct_oper_list.vector[2]"),
    data->get_opt("_pdbx_struct_oper_list.matrix[3][1]"),
    data->get_opt("_pdbx_struct_oper_list.matrix[3][2]"),
    data->get_opt("_pdbx_struct_oper_list.matrix[3][3]"),
    data->get_opt("_pdbx_struct_oper_list.vector[3]")
  };

  // build oper_list from _pdbx_struct_oper_list

  oper_list_t oper_list;

  for (unsigned i = 0, nrows = arr_id->size(); i < nrows; ++i) {
    float * matrix = oper_list[arr_id->as_s(i)].data();

    identity44f(matrix);

    for (int j = 0; j < 12; ++j) {
      matrix[j] = arr_matrix[j]->as_d(i);
    }
  }

  CoordSet ** csets = nullptr;
  int csetbeginidx = 0;

  // assembly
  for (unsigned i = 0, nrows = arr_oper_expr->size(); i < nrows; ++i) {
    if (strcmp(assembly_id, arr_assembly_id->as_s(i)))
      continue;

    const char * oper_expr    = arr_oper_expr->as_s(i);
    const char * asym_id_list = arr_asym_id_list->as_s(i);

    oper_collection_t collection = parse_oper_expression(oper_expr);
    std::vector<std::string> chains = strsplit(asym_id_list, ',');
    std::set<lexborrow_t> chains_set;
    for (auto& chain : chains) {
      auto borrowed = LexBorrow(G, chain.c_str());
      if (borrowed != LEX_BORROW_NOTFOUND) {
        chains_set.insert(borrowed);
      }
    }

    // new coord set VLA
    int ncsets = 1;
    for (const auto& c_item : collection) {
      ncsets *= c_item.size();
    }

    if (!csets) {
      csets = VLACalloc(CoordSet*, ncsets);
    } else {
      csetbeginidx = VLAGetSize(csets);
      VLASize(csets, CoordSet*, csetbeginidx + ncsets);
    }

    // for cartesian product
    int c_src_len = 1;

    // coord set for subset of atoms
    CoordSet ** c_csets = csets + csetbeginidx;
    c_csets[0] = CoordSetCopyFilterChains(cset, atInfo, chains_set);

    // build new coord sets
    for (auto c_it = collection.rbegin(); c_it != collection.rend(); ++c_it) {
      // copy
      int j = c_src_len;
      while (j < c_src_len * c_it->size()) {
        // cartesian product
        for (int k = 0; k < c_src_len; ++k, ++j) {
          c_csets[j] = CoordSetCopy(c_csets[k]);
        }
      }

      // transform
      j = 0;
      for (auto& s_item : *c_it) {
        const float * matrix = oper_list[s_item].data();

        // cartesian product
        for (int k = 0; k < c_src_len; ++k, ++j) {
          CoordSetTransform44f(c_csets[j], matrix);
        }
      }

      // cartesian product
      // Note: currently, "1m4x" seems to be the only structure in the PDB
      // which uses a cartesian product expression
      c_src_len *= c_it->size();
    }
  }

  // return assembly coordsets
  return csets;
}

/**
 * Set ribbon_trace_atoms and cartoon_trace_atoms for CA/P only models
 */
static bool read_pdbx_coordinate_model(PyMOLGlobals * G, const pymol::cif_data * data, ObjectMolecule * mol) {
  const cif_array * arr_type = data->get_arr("_pdbx_coordinate_model.type");
  const cif_array * arr_asym = data->get_arr("_pdbx_coordinate_model.asym_id");

  if (!arr_type || !arr_asym)
    return false;

  // affected chains
  std::set<pymol::zstring_view> asyms;

  // collect CA/P-only chain identifiers
  for (unsigned i = 0, nrows = arr_type->size(); i < nrows; ++i) {
    const char * type = arr_type->as_s(i);
    // no need anymore to check "CA ATOMS ONLY", since nonbonded CA are
    // now (v1.8.2) detected automatically in RepCartoon and RepRibbon
    if (strcmp(type, "P ATOMS ONLY") == 0) {
      asyms.insert(arr_asym->as_s(i));
    }
  }

  if (asyms.empty())
    return false;

  // set on atom-level
  for (int i = 0, nrows = VLAGetSize(mol->AtomInfo); i < nrows; ++i) {
    AtomInfoType * ai = mol->AtomInfo + i;
    if (asyms.count(LexStr(G, ai->segi))) {
      SettingSet(G, cSetting_cartoon_trace_atoms, true, ai);
      SettingSet(G, cSetting_ribbon_trace_atoms,  true, ai);
    }
  }

  return true;
}

/**
 * Read CELL and SYMMETRY
 */
static CSymmetry * read_symmetry(PyMOLGlobals * G, const pymol::cif_data * data) {
  const cif_array * cell[6] = {
    data->get_arr("_cell?length_a"),
    data->get_arr("_cell?length_b"),
    data->get_arr("_cell?length_c"),
    data->get_arr("_cell?angle_alpha"),
    data->get_arr("_cell?angle_beta"),
    data->get_arr("_cell?angle_gamma")
  };

  for (int i = 0; i < 6; i++)
    if (cell[i] == nullptr)
      return nullptr;

  CSymmetry * symmetry = new CSymmetry(G);
  if (!symmetry)
    return nullptr;

  float cellparams[6];
  for (int i = 0; i < 6; ++i) {
    cellparams[i] = cell[i]->as_d();
  }

  symmetry->Crystal.setDims(cellparams);
  symmetry->Crystal.setAngles(cellparams + 3);
  symmetry->setSpaceGroup(
      data->get_opt("_symmetry?space_group_name_h-m",
                    "_space_group?name_h-m_alt")->as_s());

  symmetry->PDBZValue = data->get_opt("_cell.z_pdb")->as_i(0, 1);

  // register symmetry operations if given
  const cif_array * arr_as_xyz = data->get_arr(
      "_symmetry_equiv?pos_as_xyz",
      "_space_group_symop?operation_xyz");
  if (arr_as_xyz) {
    std::vector<std::string> sym_op;
    for (unsigned i = 0, n = arr_as_xyz->size(); i < n; ++i) {
      sym_op.push_back(arr_as_xyz->as_s(i));
    }
    SymmetrySpaceGroupRegister(G, symmetry->spaceGroup(), sym_op);
  }

  return symmetry;
}

/**
 * Read CHEM_COMP_ATOM
 */
static CoordSet ** read_chem_comp_atom_model(PyMOLGlobals * G, const pymol::cif_data * data,
    AtomInfoType ** atInfoPtr) {

  const cif_array *arr_x, *arr_y = nullptr, *arr_z = nullptr;

  // setting to exclude one or more coordinate columns
  unsigned mask = SettingGetGlobal_i(G, cSetting_chem_comp_cartn_use);
  const char * feedback = "";

  if (!mask) {
    mask = 0xFF;
  }

  if ((mask & 0x01)
      && (arr_x = data->get_arr("_chem_comp_atom.pdbx_model_cartn_x_ideal"))
      && !arr_x->is_missing_all()) {
    arr_y = data->get_arr("_chem_comp_atom.pdbx_model_cartn_y_ideal");
    arr_z = data->get_arr("_chem_comp_atom.pdbx_model_cartn_z_ideal");
    feedback = ".pdbx_model_Cartn_{x,y,z}_ideal";
  } else if ((mask & 0x02)
      && (arr_x = data->get_arr("_chem_comp_atom.model_cartn_x"))) {
    arr_y = data->get_arr("_chem_comp_atom.model_cartn_y");
    arr_z = data->get_arr("_chem_comp_atom.model_cartn_z");
    feedback = ".model_Cartn_{x,y,z}";
  } else if ((mask & 0x04)
      && (arr_x = data->get_arr("_chem_comp_atom.x"))
      && !arr_x->is_missing_all()) {
    arr_y = data->get_arr("_chem_comp_atom.y");
    arr_z = data->get_arr("_chem_comp_atom.z");
    feedback = ".{x,y,z}";
  }

  if (!arr_x || !arr_y || !arr_z) {
    return nullptr;
  }

  PRINTFB(G, FB_Executive, FB_Details)
    " ExecutiveLoad-Detail: Detected chem_comp CIF (%s)\n", feedback
    ENDFB(G);

  const cif_array * arr_name            = data->get_opt("_chem_comp_atom.atom_id");
  const cif_array * arr_symbol          = data->get_opt("_chem_comp_atom.type_symbol");
  const cif_array * arr_resn            = data->get_opt("_chem_comp_atom.comp_id");
  const cif_array * arr_partial_charge  = data->get_opt("_chem_comp_atom.partial_charge");
  const cif_array * arr_formal_charge   = data->get_opt("_chem_comp_atom.charge");
  const cif_array * arr_stereo          = data->get_opt("_chem_comp_atom.pdbx_stereo_config");

  int nrows = arr_x->size();
  AtomInfoType *ai;
  int atomCount = 0, nAtom = nrows;
  float * coord = VLAlloc(float, 3 * nAtom);
  int auto_show = RepGetAutoShowMask(G);

  for (int i = 0; i < nrows; i++) {
    if (arr_x->is_missing(i))
      continue;

    VLACheck(*atInfoPtr, AtomInfoType, atomCount);
    ai = *atInfoPtr + atomCount;
    memset((void*) ai, 0, sizeof(AtomInfoType));

    ai->rank = atomCount;
    ai->id = atomCount + 1;

    LexAssign(G, ai->name, arr_name->as_s(i));
    LexAssign(G, ai->resn, arr_resn->as_s(i));
    strncpy(ai->elem, arr_symbol->as_s(i), cElemNameLen);

    ai->partialCharge = arr_partial_charge->as_d(i);
    ai->formalCharge = arr_formal_charge->as_i(i);

    ai->hetatm = true;

    ai->visRep = auto_show;
    AtomInfoSetStereo(ai, arr_stereo->as_s(i));

    AtomInfoAssignParameters(G, ai);
    AtomInfoAssignColors(G, ai);

    coord[atomCount * 3 + 0] = arr_x->as_d(i);
    coord[atomCount * 3 + 1] = arr_y->as_d(i);
    coord[atomCount * 3 + 2] = arr_z->as_d(i);

    atomCount++;
  }

  VLASize(coord, float, 3 * atomCount);
  VLASize(*atInfoPtr, AtomInfoType, atomCount);

  CoordSet ** csets = VLACalloc(CoordSet*, 1);
  csets[0] = CoordSetNew(G);
  csets[0]->NIndex = atomCount;
  csets[0]->Coord= pymol::vla_take_ownership(coord);

  return csets;
}

/**
 * Map model number to state (1-based)
 */
class ModelStateMapper {
  bool remap;
  std::map<int, int> mapping;
public:
  ModelStateMapper(bool remap) : remap(remap) {}

  int operator()(int model) {
    if (!remap)
      return model;

    int state = mapping[model];

    if (!state) {
      state = mapping.size();
      mapping[model] = state;
    }

    return state;
  }
};

/**
 * Read ATOM_SITE
 *
 * atInfoPtr: atom info array to fill
 * info: data content configuration to populate with collected information
 *
 * return: models as VLA of coordinate sets
 */
static CoordSet** read_atom_site(PyMOLGlobals* G, const pymol::cif_data* data,
    AtomInfoType** atInfoPtr, CifContentInfo& info, bool discrete, bool quiet)
{

  const cif_array *arr_x, *arr_y, *arr_z;
  const cif_array *arr_name = nullptr, *arr_resn = nullptr, *arr_resi = nullptr,
            *arr_chain = nullptr, *arr_symbol,
            *arr_group_pdb, *arr_alt, *arr_ins_code = nullptr, *arr_b, *arr_u,
            *arr_q, *arr_ID, *arr_mod_num, *arr_entity_id, *arr_segi;

  if ((arr_x = data->get_arr("_atom_site?cartn_x")) &&
      (arr_y = data->get_arr("_atom_site?cartn_y")) &&
      (arr_z = data->get_arr("_atom_site?cartn_z"))) {
  } else if (
      (arr_x = data->get_arr("_atom_site?fract_x")) &&
      (arr_y = data->get_arr("_atom_site?fract_y")) &&
      (arr_z = data->get_arr("_atom_site?fract_z"))) {
    info.fractional = true;
  } else {
    return nullptr;
  }

  if (info.use_auth) {
    arr_name      = data->get_arr("_atom_site.auth_atom_id");
    arr_resn      = data->get_arr("_atom_site.auth_comp_id");
    arr_resi      = data->get_arr("_atom_site.auth_seq_id");
    arr_chain     = data->get_arr("_atom_site.auth_asym_id");
    arr_ins_code  = data->get_arr("_atom_site.pdbx_pdb_ins_code");
  }

  if (!arr_name) arr_name = data->get_arr("_atom_site.label_atom_id");
  if (!arr_resn) arr_resn = data->get_opt("_atom_site.label_comp_id");

  const cif_array *arr_label_seq_id = data->get_opt("_atom_site.label_seq_id");

  // PDBe provides unique seq_ids for bulk het groups
  if (!arr_resi) arr_resi = data->get_arr("_atom_site.pdbe_label_seq_id");
  if (!arr_resi) arr_resi = arr_label_seq_id;

  if (arr_name) {
    info.type = CIF_MMCIF;
    if (!quiet) {
      PRINTFB(G, FB_Executive, FB_Details)
        " ExecutiveLoad-Detail: Detected mmCIF\n" ENDFB(G);
    }
  } else {
    arr_name      = data->get_opt("_atom_site_label");
    info.type = CIF_CORE;
    if (!quiet) {
      PRINTFB(G, FB_Executive, FB_Details)
        " ExecutiveLoad-Detail: Detected small molecule CIF\n" ENDFB(G);
    }
  }

  arr_segi        = data->get_opt("_atom_site.label_asym_id");
  arr_symbol      = data->get_opt("_atom_site?type_symbol", "_atom_site_label");
  arr_group_pdb   = data->get_opt("_atom_site.group_pdb");
  arr_alt         = data->get_opt("_atom_site.label_alt_id");
  arr_b           = data->get_opt("_atom_site?b_iso_or_equiv");
  arr_u           = data->get_arr("_atom_site?u_iso_or_equiv"); // nullptr
  arr_q           = data->get_opt("_atom_site?occupancy");
  arr_ID          = data->get_opt("_atom_site.id",
                                  "_atom_site_label");
  arr_mod_num     = data->get_opt("_atom_site.pdbx_pdb_model_num");
  arr_entity_id   = data->get_arr("_atom_site.label_entity_id"); // nullptr

  const cif_array * arr_color = data->get_arr("_atom_site.pymol_color");
  const cif_array * arr_reps  = data->get_arr("_atom_site.pymol_reps");
  const cif_array * arr_ss    = data->get_opt("_atom_site.pymol_ss");
  const cif_array * arr_label = data->get_opt("_atom_site.pymol_label");
  const cif_array * arr_vdw   = data->get_opt("_atom_site.pymol_vdw");
  const cif_array * arr_elec_radius = data->get_opt("_atom_site.pymol_elec_radius");
  const cif_array * arr_partial_charge = data->get_opt("_atom_site.pymol_partial_charge");
  const cif_array * arr_formal_charge = data->get_opt("_atom_site.pdbx_formal_charge");

  if (!arr_chain)
    arr_chain = arr_segi;

  ModelStateMapper model_to_state(!SettingGetGlobal_i(G, cSetting_pdb_honor_model_number));
  int nrows = arr_x->size();
  AtomInfoType *ai;
  int atomCount = 0;
  int auto_show = RepGetAutoShowMask(G);
  int first_model_num = model_to_state(arr_mod_num->as_i(0, 1));
  CoordSet * cset;
  int mod_num, ncsets = 0;

  // collect number of atoms per model and number of coord sets
  std::map<int, int> atoms_per_model;
  for (int i = 0, n = nrows; i < n; i++) {
    mod_num = model_to_state(arr_mod_num->as_i(i, 1));

    if (mod_num < 1) {
      PRINTFB(G, FB_ObjectMolecule, FB_Errors)
        " Error: model numbers < 1 not supported: %d\n", mod_num ENDFB(G);
      return nullptr;
    }

    atoms_per_model[mod_num - 1] += 1;

    if (ncsets < mod_num)
      ncsets = mod_num;
  }

  // set up coordinate sets
  CoordSet ** csets = VLACalloc(CoordSet*, ncsets);
  for (auto it = atoms_per_model.begin(); it != atoms_per_model.end(); ++it) {
    csets[it->first] = cset = CoordSetNew(G);
    cset->Coord.resize(3 * it->second);
    cset->IdxToAtm.resize(it->second);
  }

  // mm_atom_site_label -> atom index (1-indexed)
  std::map<std::string, int> name_dict;

  for (int i = 0, n = nrows; i < n; i++) {
    lexidx_t segi = LexIdx(G, arr_segi->as_s(i));

    if (info.is_excluded_chain(segi)) {
      LexDec(G, segi);
      continue;
    }

    mod_num = model_to_state(arr_mod_num->as_i(i, 1));

    // copy coordinates into coord set
    cset = csets[mod_num - 1];
    int idx = cset->NIndex++;
    float * coord = cset->coordPtr(idx);
    coord[0] = arr_x->as_d(i);
    coord[1] = arr_y->as_d(i);
    coord[2] = arr_z->as_d(i);

    if (!discrete && ncsets > 1) {
      // mm_atom_site_label aggregate
      std::string key = make_mm_atom_site_label(G,
          arr_chain->as_s(i),
          arr_resn->as_s(i),
          arr_resi->as_s(i),
          arr_ins_code ? arr_ins_code->as_s(i) : "",
          arr_name->as_s(i),
          arr_alt->as_s(i));

      // check if this is not a new atom
      if (mod_num != first_model_num) {
        int atm = name_dict[key] - 1;
        if (atm >= 0) {
          cset->IdxToAtm[idx] = atm;
          continue;
        }
      }

      name_dict[key] = atomCount + 1;
    }

    cset->IdxToAtm[idx] = atomCount;

    VLACheck(*atInfoPtr, AtomInfoType, atomCount);
    ai = *atInfoPtr + atomCount;

    ai->rank = atomCount;
    ai->alt[0] = arr_alt->as_s(i)[0];

    ai->id = arr_ID->as_i(i);
    ai->b = (arr_u != nullptr) ?
             arr_u->as_d(i) * 78.95683520871486 : // B = U * 8 * pi^2
             arr_b->as_d(i);
    ai->q = arr_q->as_d(i, 1.0);

    strncpy_alpha(ai->elem, arr_symbol->as_s(i), cElemNameLen);

    ai->chain = LexIdx(G, arr_chain->as_s(i));
    ai->name = LexIdx(G, arr_name->as_s(i));
    ai->resn = LexIdx(G, arr_resn->as_s(i));
    ai->segi = std::move(segi); // steal reference

    if ('H' == arr_group_pdb->as_s(i)[0]) {
      ai->hetatm = true;
      ai->flags = cAtomFlag_ignore;
    }

    ai->resv = arr_resi->as_i(i);
    ai->temp1 = arr_label_seq_id->as_i(i); // for add_missing_ca

    if (arr_ins_code) {
      ai->setInscode(arr_ins_code->as_s(i)[0]);
    }

    if (arr_reps) {
      ai->visRep = arr_reps->as_i(i, auto_show);
      ai->flags |= cAtomFlag_inorganic; // suppress auto_show_classified
    } else {
      ai->visRep = auto_show;
    }

    ai->ssType[0] = arr_ss->as_s(i)[0];
    ai->formalCharge = arr_formal_charge->as_i(i);
    ai->partialCharge = arr_partial_charge->as_d(i);
    ai->elec_radius = arr_elec_radius->as_d(i);
    ai->vdw = arr_vdw->as_d(i);
    ai->label = LexIdx(G, arr_label->as_s(i));

    AtomInfoAssignParameters(G, ai);

    if (arr_color) {
      ai->color = arr_color->as_i(i);
    } else {
      AtomInfoAssignColors(G, ai);
    }

    if (arr_entity_id != nullptr) {
      AtomInfoSetEntityId(G, ai, arr_entity_id->as_s(i));
    }

    atomCount++;
  }

  VLASize(*atInfoPtr, AtomInfoType, atomCount);

  return csets;
}

/**
 * Update `info` with entity polymer information
 */
static bool read_entity_poly(PyMOLGlobals * G, const pymol::cif_data * data, CifContentInfo &info) {
  const cif_array *arr_entity_id = nullptr, *arr_type = nullptr,
        *arr_num = nullptr, *arr_mon_id = nullptr;

  if (!(arr_entity_id     = data->get_arr("_entity_poly.entity_id")) ||
      !(arr_type          = data->get_arr("_entity_poly.type")))
    return false;

  const cif_array * arr_seq_one_letter = data->get_arr("_entity_poly.pdbx_seq_one_letter_code");

  // polypeptides
  for (unsigned i = 0, n = arr_entity_id->size(); i < n; i++) {
    if (!strncasecmp("polypeptide", arr_type->as_s(i), 11)) {
      const char * entity_id = arr_entity_id->as_s(i);
      info.polypeptide_entities.insert(entity_id);

      if (arr_seq_one_letter) {
        // sequences
        auto& entity_sequence = info.sequences[entity_id];
        const char * one = arr_seq_one_letter->as_s(i);
        for (int i = 0; *one; ++one) {
          if (strchr(" \t\r\n", *one)) // skip whitespace
            continue;

          if (*one == '(') {
            const char * end = strchr(one, ')');
            if (!end)
              break;

            std::string three(one + 1, end - one - 1);
            entity_sequence.set(++i, three.c_str());
            one = end;
          } else {
            entity_sequence.set(++i, aa_get_three_letter(*one));
          }
        }
      }
    }
  }

  if (!arr_seq_one_letter) {
    // sequences
    if ((arr_entity_id     = data->get_arr("_entity_poly_seq.entity_id")) &&
        (arr_num           = data->get_arr("_entity_poly_seq.num")) &&
        (arr_mon_id        = data->get_arr("_entity_poly_seq.mon_id"))) {
      for (unsigned i = 0, n = arr_entity_id->size(); i < n; i++) {
        info.sequences[arr_entity_id->as_s(i)].set(
            arr_num->as_i(i),
            arr_mon_id->as_s(i));
      }
    }
  }

  return true;
}

/**
 * Sub-routine for `add_missing_ca`
 *
 * @param i_ref Atom index of the next observed residue if `!at_terminus`,
 * otherwise of the last observed residue in this chain.
 * @param at_terminus True if adding residues beyond the last observed residue
 * in this chain.
 */
static void add_missing_ca_sub(PyMOLGlobals * G,
    pymol::vla<AtomInfoType>& atInfo,
    int& current_resv,
    int& atomCount,
    const int i_ref, int resv,
    const seqvec_t * current_seq,
    const char * entity_id,
    bool at_terminus = true)
{
  if (!atInfo[i_ref].temp1)
    return;

  if (current_resv == 0) {
    at_terminus = true;
  }

  for (++current_resv; current_resv < resv; ++current_resv) {
    const char * resn = current_seq->get(current_resv);
    if (!resn)
      continue;

    int added_resv = current_resv + (atInfo[i_ref].resv - atInfo[i_ref].temp1);

    if (!at_terminus && ((i_ref > 0 && added_resv <= atInfo[i_ref - 1].resv) ||
                            added_resv >= atInfo[i_ref].resv)) {
      // don't use insertion codes
      continue;
    }

    AtomInfoType *ai = atInfo.check(atomCount);

    ai->rank = atomCount;
    ai->id = -1;

    ai->elem[0] = 'C';
    LexAssign(G, ai->name, "CA");
    LexAssign(G, ai->resn, resn);
    LexAssign(G, ai->segi, atInfo[i_ref].segi);
    LexAssign(G, ai->chain, atInfo[i_ref].chain);

    ai->temp1 = current_resv;
    ai->resv = added_resv;

    AtomInfoAssignParameters(G, ai);
    AtomInfoAssignColors(G, ai);
    AtomInfoSetEntityId(G, ai, entity_id);

    ++atomCount;
  }
}

/**
 * Read missing residues / full sequence
 *
 * This function relies on the label_seq_id numbering which must be available
 * in the `temp1` kludge field.
 *
 * Use the _entity_poly and _entity_poly_seq information to identify
 * missing residues in partially present chains. Add CA atoms for those
 * to present complete sequences in the sequence viewer.
 */
static bool add_missing_ca(PyMOLGlobals * G,
    pymol::vla<AtomInfoType>& atInfo, CifContentInfo &info) {

  int oldAtomCount = atInfo.size();
  int atomCount = oldAtomCount;
  int current_resv = 0;
  const seqvec_t * current_seq = nullptr;
  const char * current_entity_id = "";

  for (int i = 0; i < oldAtomCount; ++i) {
    const char * entity_id = LexStr(G, atInfo[i].custom);

    if (i == 0
        || atInfo[i].chain != atInfo[i - 1].chain
        || strcmp(entity_id, current_entity_id)) {
      // finish prev seq
      if (current_seq && i > 0) {
         add_missing_ca_sub(G,
             atInfo, current_resv, atomCount,
             i - 1, current_seq->size() + 1,
             current_seq, current_entity_id);
      }

      current_resv = 0;
      current_seq = nullptr;
      current_entity_id = entity_id;

      if (info.is_polypeptide(entity_id) && !info.is_excluded_chain(atInfo[i].segi)) {
        // get new sequence
        auto it = info.sequences.find(entity_id);
        if (it != info.sequences.end()) {
          current_seq = &it->second;
        }
      }

    } else if (i > 0 && atInfo[i].temp1 == atInfo[i - 1].temp1) {
      continue;
    }

    if (current_seq) {
      add_missing_ca_sub(G,
          atInfo, current_resv, atomCount,
          i, atInfo[i].temp1,
          current_seq, entity_id, false);
    }
  }

  // finish last seq
  if (current_seq) {
    add_missing_ca_sub(G,
        atInfo, current_resv, atomCount,
        oldAtomCount - 1, current_seq->size() + 1,
        current_seq, current_entity_id);
  }

  atInfo.resize(atomCount);

  return true;
}

/**
 * Read secondary structure from STRUCT_CONF or STRUCT_SHEET_RANGE
 */
static bool read_ss_(PyMOLGlobals * G, const pymol::cif_data * data, char ss,
    sshashmap &ssrecords, CifContentInfo &info)
{
  const cif_array *arr_beg_chain = nullptr, *arr_beg_resi = nullptr,
                  *arr_end_chain = nullptr, *arr_end_resi = nullptr,
                  *arr_beg_ins_code = nullptr, *arr_end_ins_code = nullptr;

  std::string prefix = "_struct_conf.";
  if (ss == 'S')
    prefix = "_struct_sheet_range.";

  if (info.use_auth &&
      (arr_beg_chain = data->get_arr((prefix + "beg_auth_asym_id").c_str())) &&
      (arr_beg_resi  = data->get_arr((prefix + "beg_auth_seq_id").c_str())) &&
      (arr_end_chain = data->get_arr((prefix + "end_auth_asym_id").c_str())) &&
      (arr_end_resi  = data->get_arr((prefix + "end_auth_seq_id").c_str()))) {
    // auth only
    arr_beg_ins_code = data->get_arr((prefix + "pdbx_beg_pdb_ins_code").c_str());
    arr_end_ins_code = data->get_arr((prefix + "pdbx_end_pdb_ins_code").c_str());
  } else if (
      !(arr_beg_chain = data->get_arr((prefix + "beg_label_asym_id").c_str())) ||
      !(arr_beg_resi  = data->get_arr((prefix + "beg_label_seq_id").c_str())) ||
      !(arr_end_chain = data->get_arr((prefix + "end_label_asym_id").c_str())) ||
      !(arr_end_resi  = data->get_arr((prefix + "end_label_seq_id").c_str()))) {
    return false;
  }

  const cif_array *arr_conf_type_id = (ss == 'S') ? nullptr :
    data->get_arr("_struct_conf.conf_type_id");

  int nrows = arr_beg_chain->size();
  sshashkey key;

  for (int i = 0; i < nrows; i++) {
    // first character of conf_type_id (one of H, S, T)
    char ss_i = arr_conf_type_id ? arr_conf_type_id->as_s(i)[0] : ss;

    // exclude TURN_* (include HELX_* and STRN)
    if (ss_i == 'T')
      continue;

    key.assign(
      LexBorrow(G, arr_beg_chain->as_s(i)),
      arr_beg_resi->as_i(i),
      arr_beg_ins_code ? arr_beg_ins_code->as_s(i)[0] : '\0');

    sshashvalue &value = ssrecords[key];
    value.ss = ss_i;
    value.end.assign(
        LexBorrow(G, arr_end_chain->as_s(i)),
        arr_end_resi->as_i(i),
        arr_end_ins_code ? arr_end_ins_code->as_s(i)[0] : '\0');
  }

  return true;
}

/**
 * Read secondary structure
 */
static bool read_ss(PyMOLGlobals * G, const pymol::cif_data * datablock,
    pymol::vla<AtomInfoType>& atInfo, CifContentInfo &info)
{
  sshashmap ssrecords;

  read_ss_(G, datablock, 'H', ssrecords, info);
  read_ss_(G, datablock, 'S', ssrecords, info);

  if (ssrecords.empty())
    return false;

  AtomInfoType *aj, *ai, *atoms_end = atInfo + VLAGetSize(atInfo);
  sshashkey key;

  for (ai = atInfo.data(); ai < atoms_end;) {
    // advance to the next residue
    aj = ai;
    while (++ai < atoms_end &&
        AtomInfoSameResidue(G, aj, ai)) {}

    // check if residue is the beginning of a secondary structure element
    key.assign(aj->chain, aj->resv, aj->inscode);
    sshashmap::iterator it = ssrecords.find(key);

    if (it == ssrecords.end())
      continue;

    sshashvalue &value = it->second;

    // assign ss type to all atoms in the segment
    bool hit_end_residue = false;
    for (; aj < atoms_end; aj++) {
      if (value.end.compare(*aj) == 0) {
        hit_end_residue = true;
      } else if (hit_end_residue) {
        break;
      }
      aj->ssType[0] = value.ss;
    }
  }

  return true;
}

/**
 * Read the SCALEn matrix into 4x4 `matrix`
 */
static bool read_atom_site_fract_transf(PyMOLGlobals * G, const pymol::cif_data * data, float * matrix) {
  const cif_array *arr_transf[12];

  if (!(arr_transf[0] = data->get_arr("_atom_sites.fract_transf_matrix[1][1]", "_atom_sites_fract_tran_matrix_11")))
    return false;

  arr_transf[1]  = data->get_opt("_atom_sites.fract_transf_matrix[1][2]", "_atom_sites_fract_tran_matrix_12");
  arr_transf[2]  = data->get_opt("_atom_sites.fract_transf_matrix[1][3]", "_atom_sites_fract_tran_matrix_13");
  arr_transf[3]  = data->get_opt("_atom_sites.fract_transf_vector[1]", "_atom_sites_fract_tran_vector_1");
  arr_transf[4]  = data->get_opt("_atom_sites.fract_transf_matrix[2][1]", "_atom_sites_fract_tran_matrix_21");
  arr_transf[5]  = data->get_opt("_atom_sites.fract_transf_matrix[2][2]", "_atom_sites_fract_tran_matrix_22");
  arr_transf[6]  = data->get_opt("_atom_sites.fract_transf_matrix[2][3]", "_atom_sites_fract_tran_matrix_23");
  arr_transf[7]  = data->get_opt("_atom_sites.fract_transf_vector[2]", "_atom_sites_fract_tran_vector_2");
  arr_transf[8]  = data->get_opt("_atom_sites.fract_transf_matrix[3][1]", "_atom_sites_fract_tran_matrix_31");
  arr_transf[9]  = data->get_opt("_atom_sites.fract_transf_matrix[3][2]", "_atom_sites_fract_tran_matrix_32");
  arr_transf[10] = data->get_opt("_atom_sites.fract_transf_matrix[3][3]", "_atom_sites_fract_tran_matrix_33");
  arr_transf[11] = data->get_opt("_atom_sites.fract_transf_vector[3]", "_atom_sites_fract_tran_vector_3");

  for (int i = 0; i < 12; ++i)
    matrix[i] = arr_transf[i]->as_d(0);

  zero3f(matrix + 12);
  matrix[15] = 1.f;

  return true;
}

/**
 * Read anisotropic temperature factors from ATOM_SITE or ATOM_SITE_ANISOTROP
 */
static bool read_atom_site_aniso(PyMOLGlobals * G, const pymol::cif_data * data,
    pymol::vla<AtomInfoType>& atInfo) {

  const cif_array *arr_label, *arr_u11, *arr_u22, *arr_u33, *arr_u12, *arr_u13, *arr_u23;
  bool mmcif = true;
  float factor = 1.0;

  if ((arr_label = data->get_arr("_atom_site_anisotrop.id", "_atom_site.id"))) {
    // mmCIF, assume _atom_site_id is numeric and look up by atom ID
    // Warning: according to mmCIF spec, id can be any alphanumeric string
  } else if ((arr_label = data->get_arr("_atom_site_aniso_label"))) {
    // small molecule CIF, lookup by atom name
    mmcif = false;
  } else {
    return false;
  }

  if ((arr_u11 = data->get_arr("_atom_site_anisotrop.u[1][1]", "_atom_site_aniso_u_11", "_atom_site.aniso_u[1][1]"))) {
    // U
    arr_u22 = data->get_opt("_atom_site_anisotrop.u[2][2]", "_atom_site_aniso_u_22", "_atom_site.aniso_u[2][2]");
    arr_u33 = data->get_opt("_atom_site_anisotrop.u[3][3]", "_atom_site_aniso_u_33", "_atom_site.aniso_u[3][3]");
    arr_u12 = data->get_opt("_atom_site_anisotrop.u[1][2]", "_atom_site_aniso_u_12", "_atom_site.aniso_u[1][2]");
    arr_u13 = data->get_opt("_atom_site_anisotrop.u[1][3]", "_atom_site_aniso_u_13", "_atom_site.aniso_u[1][3]");
    arr_u23 = data->get_opt("_atom_site_anisotrop.u[2][3]", "_atom_site_aniso_u_23", "_atom_site.aniso_u[2][3]");
  } else if (
      (arr_u11 = data->get_arr("_atom_site_anisotrop.b[1][1]", "_atom_site_aniso_b_11", "_atom_site.aniso_b[1][1]"))) {
    // B
    factor = 0.012665147955292222; // U = B / (8 * pi^2)
    arr_u22 = data->get_opt("_atom_site_anisotrop.b[2][2]", "_atom_site_aniso_b_22", "_atom_site.aniso_b[2][2]");
    arr_u33 = data->get_opt("_atom_site_anisotrop.b[3][3]", "_atom_site_aniso_b_33", "_atom_site.aniso_b[3][3]");
    arr_u12 = data->get_opt("_atom_site_anisotrop.b[1][2]", "_atom_site_aniso_b_12", "_atom_site.aniso_b[1][2]");
    arr_u13 = data->get_opt("_atom_site_anisotrop.b[1][3]", "_atom_site_aniso_b_13", "_atom_site.aniso_b[1][3]");
    arr_u23 = data->get_opt("_atom_site_anisotrop.b[2][3]", "_atom_site_aniso_b_23", "_atom_site.aniso_b[2][3]");
  } else {
    return false;
  }

  AtomInfoType *ai;
  int nAtom = VLAGetSize(atInfo);

  std::map<int, AtomInfoType*> id_dict;
  std::map<std::string, AtomInfoType*> name_dict;

  // build dictionary
  for (int i = 0; i < nAtom; i++) {
    ai = atInfo + i;
    if (mmcif) {
      id_dict[ai->id] = ai;
    } else {
      std::string key(LexStr(G, ai->name));
      name_dict[key] = ai;
    }
  }

  // read aniso table
  for (unsigned i = 0; i < arr_u11->size(); i++) {
    ai = nullptr;

    if (mmcif) {
      find1(id_dict, ai, arr_label->as_i(i));
    } else {
      find1(name_dict, ai, arr_label->as_s(i));
    }

    if (!ai) {
      // expected for multi-models
      continue;
    }

    float * anisou = ai->get_anisou();
    anisou[0] = arr_u11->as_d(i) * factor;
    anisou[1] = arr_u22->as_d(i) * factor;
    anisou[2] = arr_u33->as_d(i) * factor;
    anisou[3] = arr_u12->as_d(i) * factor;
    anisou[4] = arr_u13->as_d(i) * factor;
    anisou[5] = arr_u23->as_d(i) * factor;
  }

  return true;
}

/**
 * Read GEOM_BOND
 *
 * return: BondType VLA
 */
static pymol::vla<BondType> read_geom_bond(PyMOLGlobals * G, const pymol::cif_data * data,
    const pymol::vla<AtomInfoType>& atInfo) {

  const cif_array *arr_ID_1, *arr_ID_2;
  if ((arr_ID_1 = data->get_arr("_geom_bond.atom_site_id_1",
                                "_geom_bond_atom_site_label_1")) == nullptr ||
      (arr_ID_2 = data->get_arr("_geom_bond.atom_site_id_2",
                                "_geom_bond_atom_site_label_2")) == nullptr)
    return {};

  const cif_array *arr_symm_1 = data->get_opt("_geom_bond?site_symmetry_1");
  const cif_array *arr_symm_2 = data->get_opt("_geom_bond?site_symmetry_2");

  int nrows = arr_ID_1->size();
  int nAtom = VLAGetSize(atInfo);
  int nBond = 0;

  auto bondvla = pymol::vla<BondType>(6 * nAtom);

  // name -> atom index
  std::map<std::string, int> name_dict;

  // build dictionary
  for (int i = 0; i < nAtom; i++) {
    std::string key(LexStr(G, atInfo[i].name));
    name_dict[key] = i;
  }

  // read table
  for (int i = 0; i < nrows; i++) {
    std::string key1(arr_ID_1->as_s(i));
    std::string key2(arr_ID_2->as_s(i));

    int i1, i2;
    if (find2(name_dict, i1, key1, i2, key2)) {
      auto const bond = bondvla.check(nBond);
      if (BondTypeInit3(G, bond, i1, i2, //
              arr_symm_1->as_s(i),       //
              arr_symm_2->as_s(i))) {
        ++nBond;
      }
    } else {
      PRINTFB(G, FB_Executive, FB_Details)
        " Executive-Detail: _geom_bond name lookup failed: %s %s\n",
        key1.c_str(), key2.c_str() ENDFB(G);
    }
  }

  if (nBond) {
    VLASize(bondvla, BondType, nBond);
  } else {
    VLAFreeP(bondvla);
  }

  return bondvla;
}

/**
 * Read CHEMICAL_CONN_BOND
 *
 * return: BondType VLA
 */
static pymol::vla<BondType> read_chemical_conn_bond(PyMOLGlobals * G, const pymol::cif_data * data) {

  const cif_array *arr_number, *arr_atom_1, *arr_atom_2, *arr_type;

  if ((arr_number = data->get_arr("_atom_site?chemical_conn_number")) == nullptr ||
      (arr_atom_1 = data->get_arr("_chemical_conn_bond?atom_1")) == nullptr ||
      (arr_atom_2 = data->get_arr("_chemical_conn_bond?atom_2")) == nullptr ||
      (arr_type   = data->get_arr("_chemical_conn_bond?type")) == nullptr)
    return {};

  int nAtom = arr_number->size();
  int nBond = arr_atom_1->size();

  auto bondvla = pymol::vla<BondType>(nBond);
  auto bond = bondvla.data();

  // chemical_conn_number -> atom index
  std::map<int, int> number_dict;

  // build dictionary
  for (int i = 0; i < nAtom; i++) {
    number_dict[arr_number->as_i(i)] = i;
  }

  // read table
  int i1, i2;
  for (int i = 0; i < nBond; i++) {
    if (find2(number_dict,
          i1, arr_atom_1->as_i(i),
          i2, arr_atom_2->as_i(i))) {
      BondTypeInit2(bond++, i1, i2,
          bondOrderLookup(arr_type->as_s(i)));
    } else {
      PRINTFB(G, FB_Executive, FB_Details)
        " Executive-Detail: _chemical_conn_bond name lookup failed\n" ENDFB(G);
    }
  }

  return bondvla;
}

/**
 * Read bonds from STRUCT_CONN
 *
 * Output:
 *   cset->TmpBond
 *   cset->NTmpBond
 */
static bool read_struct_conn_(PyMOLGlobals * G, const pymol::cif_data * data,
    const pymol::vla<AtomInfoType>& atInfo, CoordSet * cset,
    CifContentInfo &info) {

  const cif_array *col_type_id = data->get_arr("_struct_conn.conn_type_id");

  if (!col_type_id)
    return false;

  const cif_array
    *col_asym_id[2] = {nullptr, nullptr},
    *col_comp_id[2] = {nullptr, nullptr},
    *col_seq_id[2] = {nullptr, nullptr},
    *col_atom_id[2] = {nullptr, nullptr},
    *col_alt_id[2] = {nullptr, nullptr},
    *col_ins_code[2] = {nullptr, nullptr},
    *col_symm[2] = {nullptr, nullptr};

  if (info.use_auth) {
    col_asym_id[0] = data->get_arr("_struct_conn.ptnr1_auth_asym_id");
    col_comp_id[0] = data->get_arr("_struct_conn.ptnr1_auth_comp_id");
    col_seq_id[0]  = data->get_arr("_struct_conn.ptnr1_auth_seq_id");
    col_atom_id[0] = data->get_arr("_struct_conn.ptnr1_auth_atom_id");
    col_asym_id[1] = data->get_arr("_struct_conn.ptnr2_auth_asym_id");
    col_comp_id[1] = data->get_arr("_struct_conn.ptnr2_auth_comp_id");
    col_seq_id[1]  = data->get_arr("_struct_conn.ptnr2_auth_seq_id");
    col_atom_id[1] = data->get_arr("_struct_conn.ptnr2_auth_atom_id");

    col_alt_id[0]  = data->get_arr("_struct_conn.pdbx_ptnr1_auth_alt_id");
    col_alt_id[1]  = data->get_arr("_struct_conn.pdbx_ptnr2_auth_alt_id");

    // auth only
    col_ins_code[0] = data->get_arr("_struct_conn.pdbx_ptnr1_pdb_ins_code");
    col_ins_code[1] = data->get_arr("_struct_conn.pdbx_ptnr2_pdb_ins_code");
  }

  // for assembly chain filtering
  const cif_array *col_label_asym_id[2] = {
    data->get_arr("_struct_conn.ptnr1_label_asym_id"),
    data->get_arr("_struct_conn.ptnr2_label_asym_id")
  };

  if ((!col_asym_id[0] && !(col_asym_id[0] = col_label_asym_id[0])) ||
      (!col_comp_id[0] && !(col_comp_id[0] = data->get_arr("_struct_conn.ptnr1_label_comp_id"))) ||
      (!col_seq_id[0]  && !(col_seq_id[0]  = data->get_arr("_struct_conn.ptnr1_label_seq_id"))) ||
      (!col_atom_id[0] && !(col_atom_id[0] = data->get_arr("_struct_conn.ptnr1_label_atom_id"))) ||
      (!col_asym_id[1] && !(col_asym_id[1] = col_label_asym_id[1])) ||
      (!col_comp_id[1] && !(col_comp_id[1] = data->get_arr("_struct_conn.ptnr2_label_comp_id"))) ||
      (!col_seq_id[1]  && !(col_seq_id[1]  = data->get_arr("_struct_conn.ptnr2_label_seq_id"))) ||
      (!col_atom_id[1] && !(col_atom_id[1] = data->get_arr("_struct_conn.ptnr2_label_atom_id"))))
    return false;

  if (!col_alt_id[0]) col_alt_id[0] = data->get_opt("_struct_conn.pdbx_ptnr1_label_alt_id");
  if (!col_alt_id[1]) col_alt_id[1] = data->get_opt("_struct_conn.pdbx_ptnr2_label_alt_id");

  col_symm[0]     = data->get_opt("_struct_conn.ptnr1_symmetry");
  col_symm[1]     = data->get_opt("_struct_conn.ptnr2_symmetry");

  const cif_array *col_order = data->get_opt("_struct_conn.pdbx_value_order");

  int nrows = col_type_id->size();
  int nAtom = VLAGetSize(atInfo);
  int nBond = 0;

  cset->TmpBond = pymol::vla<BondType>(6 * nAtom);

  // identifiers -> coord set index
  std::map<std::string, int> name_dict;

  for (int i = 0; i < nAtom; i++) {
    int idx = cset->atmToIdx(i);
    if (idx != -1)
      name_dict[make_mm_atom_site_label(G, atInfo + i)] = idx;
  }

#ifdef _PYMOL_IP_EXTRAS
  bool metalc_as_zero = SettingGetGlobal_b(G, cSetting_cif_metalc_as_zero_order_bonds);
#endif

  for (int i = 0; i < nrows; i++) {
    const char * type_id = col_type_id->as_s(i);
    if (strncasecmp(type_id, "covale", 6) &&
        strcasecmp(type_id, "modres") &&
#ifdef _PYMOL_IP_EXTRAS
        !(metalc_as_zero && strcasecmp(type_id, "metalc") == 0) &&
#endif
        strcasecmp(type_id, "disulf"))
      // ignore non-covalent bonds (saltbr, hydrog)
      continue;

    std::string key[2];
    for (int j = 0; j < 2; j++) {
      const char * asym_id = col_asym_id[j]->as_s(i);

      if (col_label_asym_id[j] &&
          info.is_excluded_chain(col_label_asym_id[j]->as_s(i)))
        goto next_row;

      // doen't work with label_seq_id and bulk solvent
      const char * seq_id = col_seq_id[j]->as_s(i);
      if (!seq_id[0])
        goto next_row;

      key[j] = make_mm_atom_site_label(G,
          asym_id,
          col_comp_id[j]->as_s(i),
          seq_id,
          col_ins_code[j] ? col_ins_code[j]->as_s(i) : "",
          col_atom_id[j]->as_s(i),
          col_alt_id[j]->as_s(i));
    }

    int i1, i2;
    if (find2(name_dict, i1, key[0], i2, key[1])) {
      // zero-order bond for metal coordination
      int order = strcasecmp(type_id, "metalc") ? 1 : 0;

      if (order) {
        order = bondOrderLookup(col_order->as_s(i));
      }

      auto const bond = cset->TmpBond.check(nBond);
      if (BondTypeInit3(G, bond, i1, i2, //
              col_symm[0]->as_s(i),      //
              col_symm[1]->as_s(i), order)) {
        ++nBond;
      }
    } else {
      PRINTFB(G, FB_Executive, FB_Details)
        " Executive-Detail: _struct_conn name lookup failed: %s %s\n",
        key[0].c_str(), key[1].c_str() ENDFB(G);
    }

    // label to "continue" from inner for-loop
next_row:;
  }

  if (nBond) {
    VLASize(cset->TmpBond, BondType, nBond);
    cset->NTmpBond = nBond;
  } else {
    VLAFreeP(cset->TmpBond);
  }

  return true;
}

/**
 * Read bonds from CHEM_COMP_BOND
 *
 * return: BondType VLA
 */
static pymol::vla<BondType> read_chem_comp_bond(PyMOLGlobals * G, const pymol::cif_data * data,
    const pymol::vla<AtomInfoType>& atInfo) {

  const cif_array *col_ID_1, *col_ID_2, *col_comp_id;

  if ((col_ID_1    = data->get_arr("_chem_comp_bond.atom_id_1")) == nullptr ||
      (col_ID_2    = data->get_arr("_chem_comp_bond.atom_id_2")) == nullptr ||
      (col_comp_id = data->get_arr("_chem_comp_bond.comp_id")) == nullptr)
    return {};

  // "_chem_comp_bond.type" seems to be non-standard here. It's found in the
  // wild with values like "double" and "aromatic". mmcif_nmr-star.dic defines
  // it, but with different vocabulary (e.g. "amide", "ether", etc.).

  const cif_array *col_order = data->get_opt(
      "_chem_comp_bond.value_order",
      "_chem_comp_bond.type");

  int nrows = col_ID_1->size();
  int nAtom = VLAGetSize(atInfo);
  int nBond = 0;

  auto bondvla = pymol::vla<BondType>(6 * nAtom);
  auto bond = bondvla.data();

  // name -> atom index
  std::map<std::string, int> name_dict;

  for (int i = 0; i < nAtom; i++) {
    std::string key(LexStr(G, atInfo[i].name));
    name_dict[key] = i;
  }

  for (int i = 0; i < nrows; i++) {
    std::string key1(col_ID_1->as_s(i));
    std::string key2(col_ID_2->as_s(i));
    const char * order = col_order->as_s(i);

    int i1, i2;
    if (find2(name_dict, i1, key1, i2, key2)) {
      int order_value = bondOrderLookup(order);

      nBond++;
      BondTypeInit2(bond++, i1, i2, order_value);

    } else {
      PRINTFB(G, FB_Executive, FB_Details)
        " Executive-Detail: _chem_comp_bond name lookup failed: %s %s\n",
        key1.c_str(), key2.c_str() ENDFB(G);
    }
  }

  if (nBond) {
    VLASize(bondvla, BondType, nBond);
  } else {
    VLAFreeP(bondvla);
  }

  return bondvla;
}

/**
 * Read bonds from _pymol_bond (non-standard extension)
 *
 * return: BondType VLA
 */
static pymol::vla<BondType> read_pymol_bond(PyMOLGlobals * G, const pymol::cif_data * data,
    const pymol::vla<AtomInfoType>& atInfo) {

  const cif_array *col_ID_1, *col_ID_2, *col_order;

  if ((col_ID_1    = data->get_arr("_pymol_bond.atom_site_id_1")) == nullptr ||
      (col_ID_2    = data->get_arr("_pymol_bond.atom_site_id_2")) == nullptr ||
      (col_order   = data->get_arr("_pymol_bond.order")) == nullptr)
    return {};

  int nrows = col_ID_1->size();
  int nAtom = VLAGetSize(atInfo);

  auto bondvla = pymol::vla<BondType>(nrows);
  auto bond = bondvla.data();

  // ID -> atom index
  std::map<int, int> id_dict;

  for (int atm = 0; atm < nAtom; ++atm) {
    id_dict[atInfo[atm].id] = atm;
  }

  for (int i = 0; i < nrows; i++) {
    auto key1 = col_ID_1->as_i(i);
    auto key2 = col_ID_2->as_i(i);
    auto order_value = col_order->as_i(i);

    int i1, i2;
    if (find2(id_dict, i1, key1, i2, key2)) {
      BondTypeInit2(bond++, i1, i2, order_value);
    } else {
      PRINTFB(G, FB_Executive, FB_Details)
        " Executive-Detail: _pymol_bond name lookup failed: %d %d\n",
        key1, key2 ENDFB(G);
    }
  }

  return bondvla;
}

/**
 * Create a new (multi-state) object-molecule from datablock
 */
ObjectMolecule *ObjectMoleculeReadCifData(PyMOLGlobals * G,
    const pymol::cif_data * datablock, int discrete, bool quiet)
{
  CoordSet ** csets = nullptr;
  int ncsets;
  CifContentInfo info(G, SettingGetGlobal_b(G, cSetting_cif_use_auth));
  const char * assembly_id = SettingGetGlobal_s(G, cSetting_assembly);

  // title "echo tag"
  const char * title = datablock->get_opt("_struct.title")->as_s();
  if (!quiet && title[0] &&
      strstr(SettingGetGlobal_s(G, cSetting_pdb_echo_tags), "TITLE")) {
    PRINTFB(G, FB_ObjectMolecule, FB_Details)
      "TITLE     %s\n", title ENDFB(G);
  }

  if (assembly_id && assembly_id[0]) {
    if (!get_assembly_chains(G, datablock, info.chains_filter, assembly_id))
      PRINTFB(G, FB_Executive, FB_Details)
        " ExecutiveLoad-Detail: No such assembly: '%s'\n", assembly_id ENDFB(G);
  }

  // allocate ObjectMolecule
  ObjectMolecule * I = new ObjectMolecule(G, (discrete > 0));
  I->Color = AtomInfoUpdateAutoColor(G);

  // read coordsets from datablock
  if ((csets = read_atom_site(
           G, datablock, &I->AtomInfo, info, I->DiscreteFlag, quiet))) {
    // anisou
    read_atom_site_aniso(G, datablock, I->AtomInfo);

    // secondary structure
    read_ss(G, datablock, I->AtomInfo, info);

    // trace atoms
    read_pdbx_coordinate_model(G, datablock, I);

    // polymer information
    read_entity_poly(G, datablock, info);

    // missing residues
    if (!I->DiscreteFlag && !SettingGetGlobal_i(G, cSetting_retain_order)) {
      add_missing_ca(G, I->AtomInfo, info);
    }
  } else if ((csets = read_chem_comp_atom_model(G, datablock, &I->AtomInfo))) {
    info.type = CIF_CHEM_COMP;
  } else {
    DeleteP(I);
    return nullptr;
  }

  // get number of atoms and coordinate sets
  I->NAtom = VLAGetSize(I->AtomInfo);
  ncsets = VLAGetSize(csets);

  // initialize the new coordsets (not data, but indices, etc.)
  for (int i = 0; i < ncsets; i++) {
    if (csets[i]) {
      csets[i]->Obj = I;
      if (csets[i]->IdxToAtm.empty())
        csets[i]->enumIndices();
    }
  }

  // get coordinate sets into ObjectMolecule
  VLAFreeP(I->CSet);
  I->CSet = pymol::vla_take_ownership(csets);
  I->NCSet = ncsets;
  I->updateAtmToIdx();

  // handle symmetry and update fractional -> cartesian
  I->Symmetry.reset(read_symmetry(G, datablock));
  if (I->Symmetry) {
    float sca[16];

    if (info.fractional) {
      for (int i = 0; i < ncsets; i++) {
        if (csets[i])
          CoordSetFracToReal(csets[i], &I->Symmetry->Crystal);
      }
      } else if (info.chains_filter.empty() &&
          read_atom_site_fract_transf(G, datablock, sca)) {
        // don't do this for assemblies
        for (int i = 0; i < ncsets; i++) {
          if (csets[i])
            CoordSetInsureOrthogonal(G, csets[i], sca, &I->Symmetry->Crystal);
        }
      }
  }

  // coord set to use for distance based bonding and for attaching TmpBond
  CoordSet * cset = VLAGetFirstNonNULL(csets);

  // create bonds
  switch (info.type) {
    case CIF_CHEM_COMP:
      I->Bond = read_chem_comp_bond(G, datablock, I->AtomInfo);
      break;
    case CIF_CORE:
      I->Bond = read_geom_bond(G, datablock, I->AtomInfo);

      if (!I->Bond)
        I->Bond = read_chemical_conn_bond(G, datablock);

      break;
    case CIF_MMCIF:
      I->Bond = read_pymol_bond(G, datablock, I->AtomInfo);

      if (cset && !I->Bond) {
        // sort atoms internally
        ObjectMoleculeSort(I);

        // bonds from file, goes to cset->TmpBond
        read_struct_conn_(G, datablock, I->AtomInfo, cset, info);

        // macromolecular bonding
        bond_dict_t bond_dict_local;
        if (read_chem_comp_bond_dict(datablock, bond_dict_local)) {
          ObjectMoleculeConnectComponents(I, &bond_dict_local);
        } else if(SettingGetGlobal_i(G, cSetting_connect_mode) == 4) {
          // read components.cif
          ObjectMoleculeConnectComponents(I);
        }
      }
      break;
    case CIF_UNKNOWN:
      printf("coding error...\n");
  }

  // if non of the above created I->Bond, then do distance based bonding
  if (!I->Bond) {
    if (I->DiscreteFlag) {
      ObjectMoleculeConnectDiscrete(I, true, 3);
    } else if (cset) {
      ObjectMoleculeConnect(I, cset, true, 3);
    }

    // guess valences for distance based bonding
    if (SettingGetGlobal_b(G, cSetting_pdb_hetatm_guess_valences)) {
      ObjectMoleculeGuessValences(I, 0, nullptr, nullptr, false);
    }
  } else {
    if (!I->NBond)
      I->NBond = VLAGetSize(I->Bond);

    // bonds from coordset
    if (cset && cset->TmpBond && cset->NTmpBond) {
      for (int i = 0; i < cset->NTmpBond; ++i) {
        ObjectMoleculeAddBond2(I,
            cset->IdxToAtm[cset->TmpBond[i].index[0]],
            cset->IdxToAtm[cset->TmpBond[i].index[1]],
            cset->TmpBond[i].order);
      }
      VLASize(I->Bond, BondType, I->NBond);
      VLAFreeP(cset->TmpBond);
    }
  }

  // assemblies
  if (cset && !info.chains_filter.empty()) {
    PRINTFB(G, FB_Executive, FB_Details)
      " ExecutiveLoad-Detail: Creating assembly '%s'\n", assembly_id ENDFB(G);

    CoordSet **assembly_csets = read_pdbx_struct_assembly(G, datablock,
        I->AtomInfo, cset, assembly_id);

    ObjectMoleculeSetAssemblyCSets(I, assembly_csets);
  }

  // computationally intense update tasks
  SceneCountFrames(G);
  I->invalidate(cRepAll, cRepInvAll, -1);
  ObjectMoleculeUpdateIDNumbers(I);
  ObjectMoleculeUpdateNonbonded(I);
  ObjectMoleculeAutoDisableAtomNameWildcard(I);

  // hetatm classification if `group_PDB` record missing
  if (info.type == CIF_MMCIF && !datablock->get_arr("_atom_site.group_pdb")) {
    I->need_hetatm_classification = true;
  }

  return I;
}

/**
 * @brief Structure to hold reciprocal space parameters
 */
struct ReciprocalSpaceParams {
  glm::vec3 abc{};
  glm::vec3 angles{};
  float volume{};
};

/**
 * @brief Calculate the reciprocal space parameters for a given crystal
 * @param crystal The crystal object
 * @return The reciprocal space parameters
 * @note Fundamentals of Crystallography, 3rd edition (Table 2.1)
 */
ReciprocalSpaceParams CalculateReciprocalSpaceParam(const CCrystal& crystal) {
  auto dims = crystal.dims();
  auto a = dims[0];
  auto b = dims[1];
  auto c = dims[2];
  auto angles = crystal.angles();
  auto alpha = glm::radians(angles[0]);
  auto beta = glm::radians(angles[1]);
  auto gamma = glm::radians(angles[2]);
  auto cosAlpha = std::cos(alpha);
  auto cosBeta = std::cos(beta);
  auto cosGamma = std::cos(gamma);
  auto sinAlpha = std::sin(alpha);
  auto sinBeta = std::sin(beta);
  auto sinGamma = std::sin(gamma);
  auto volume =
      a * b * c *
      std::sqrt(1.0 - cosAlpha * cosAlpha - cosBeta * cosBeta -
                cosGamma * cosGamma + 2 * cosAlpha * cosBeta * cosGamma);
  auto inv_volume = 1.0 / volume;
  auto aStar = b * c * sinAlpha * inv_volume;
  auto bStar = a * c * sinBeta * inv_volume;
  auto cStar = a * b * sinGamma * inv_volume;
  auto cosAlphaStar = (cosBeta * cosGamma - cosAlpha) / (sinBeta * sinGamma);
  auto cosBetaStar = (cosAlpha * cosGamma - cosBeta) / (sinAlpha * sinGamma);
  auto cosGammaStar = (cosAlpha * cosBeta - cosGamma) / (sinAlpha * sinBeta);
  ReciprocalSpaceParams params{};
  params.abc = {aStar, bStar, cStar};
  params.angles = {cosAlphaStar, cosBetaStar, cosGammaStar};
  params.volume = static_cast<float>(volume);
  return params;
}

/**
 * @brief Calculate the squared distance in reciprocal space for a given hkl
 * @param param The reciprocal space parameters
 * @param hkl The Miller indices
 * @return The squared distance in reciprocal space
 * @note Fundamentals of Crystallography, 3rd edition (Eqn 2.17b)
 *      triclinic system is the most general case
 */
float CalculateReciprocalSpaceDistanceSquared(
    const ReciprocalSpaceParams& param, const glm::vec3& hkl)
{
  auto a_star = param.abc[0];
  auto b_star = param.abc[1];
  auto c_star = param.abc[2];
  auto cos_alpha_star = param.angles[0];
  auto cos_beta_star = param.angles[1];
  auto cos_gamma_star = param.angles[2];
  auto h = hkl.x;
  auto k = hkl.y;
  auto l = hkl.z;
  return (h * h * a_star * a_star +                      //
          k * k * b_star * b_star +                      //
          l * l * c_star * c_star +                      //
          2 * h * k * a_star * b_star * cos_gamma_star + //
          2 * h * l * a_star * c_star * cos_beta_star +  //
          2 * k * l * b_star * c_star * cos_alpha_star);
}

/**
 * @brief Structure to hold reflection information
 */
struct ReflectionInfo {
  glm::ivec3 hkl{};
  float fwt{};
  float phwt{};
  float fom{};
};

/**
 * @brief Class to handle the reflection data from a CIF data block
 */
class ReflnBlock
{
public:
  /**
   * @brief Constructor
   * @param datablock The CIF data block containing reflection data
   * @pre the CIF data block must contain the required reflection arrays
   */
  ReflnBlock(const pymol::cif_data& datablock)
      : m_datablock(&datablock)
      , m_size(datablock.get_arr("_refln.index_h")->size())
  {
  }

  /**
   * @brief Get the reflection information for a given index
   * @param index The index of the reflection
   * @return The reflection information
   */
  pymol::Result<ReflectionInfo> getReflectionInfo(int index) const
  {
    if (index < 0 || index >= m_size) {
      return pymol::make_error("Index out of bounds");
    }
    auto* millerH = m_datablock->get_arr("_refln.index_h");
    auto* millerK = m_datablock->get_arr("_refln.index_k");
    auto* millerL = m_datablock->get_arr("_refln.index_l");
    auto* fwt = m_datablock->get_arr("_refln.pdbx_fwt");
    auto* phwt = m_datablock->get_arr("_refln.pdbx_phwt");
    auto* fom = m_datablock->get_arr("_refln.fom");

    return ReflectionInfo{glm::ivec3(millerH->as_i(index), millerK->as_i(index),
                              millerL->as_i(index)),
        fwt->as_f(index), phwt->as_f(index), fom->as_f(index)};
  }

  /**
   * @return The number of reflections in the block
   */
  auto size() const { return m_size; }

private:
  const pymol::cif_data* m_datablock{};
  const std::size_t m_size{};
};


/**
 * @brief Convert HKL to flat index
 * @param hkl The Miller indices
 * @param shape The shape of the map
 * @return The flat index corresponding to the HKL
 */
static std::size_t HKLToFlatIndex(
    const glm::ivec3& hkl, const glm::ivec3& shape)
{
  return static_cast<std::size_t>(
      hkl.z + hkl.y * shape[2] + hkl.x * shape[1] * shape[2]);
}

/**
 * @brief Convert flat index to HKL
 * @param flat_idx The flat index
 * @param shape The shape of the map
 * @return The HKL corresponding to the flat index
 */
static glm::ivec3 FlatIndexToHKL(
    std::size_t flat_idx, const glm::ivec3& shape)
{
  int iz = flat_idx % shape[2];
  int iy = (flat_idx / shape[2]) % shape[1];
  int ix = flat_idx / (shape[1] * shape[2]);
  return glm::ivec3(ix, iy, iz);
}

#ifdef TRACE
/**
 * @brief Print statistics for the reciprocal grid
 * @param recip_map The reciprocal map to analyze
 */
static void PrintReciprocalGridStatistics(
    const CFieldTyped<std::complex<float>>& recip_map)
{
  double sum_abs_real_for_scale = 0.0;
  double max_abs_imag = 0.0;
  auto* map_flat_data = recip_map.ptr(0, 0, 0);
  auto total_map_elements = recip_map.size() / sizeof(std::complex<float>);
  for (std::size_t i = 0; i < total_map_elements; ++i) {
    auto val = map_flat_data[i];
    sum_abs_real_for_scale += std::abs(val.real());
    max_abs_imag = std::max(max_abs_imag, std::abs((double)val.imag()));
  }

  std::cout << "Reciprocal Map - Max absolute imaginary part: "
            << max_abs_imag << "\n";
  std::cout << "Average absolute real part (for scale): "
            << sum_abs_real_for_scale / total_map_elements << "\n";
  if (sum_abs_real_for_scale > EPSILON &&
      (max_abs_imag / sum_abs_real_for_scale > 1e-5)) {
    std::cerr << "Warning: Significant imaginary components found in c2c FFT "
                 "output! Max abs imag: "
              << max_abs_imag << "\n";
  } else {
    std::cout << "c2c FFT output seems valid. Max abs imag: " << max_abs_imag
              << "\n";
  }
}

/**
 * @brief Print statistics for the raw map
 * @param refln_block The reflection block containing reflection data
 * @param map The map to analyze
 */
static void PrintRawMapStatistics(
    const ReflnBlock& refln_block, const CFieldTyped<float>& map)
{
  double sum_scaled_rho = 0.0;
  double sum_sq_scaled_rho = 0.0;
  auto* map_flat_data = map.ptr(0, 0, 0);
  float min_scaled_rho = map_flat_data[0];
  float max_scaled_rho = map_flat_data[0];
  auto total_map_elements = map.size() / sizeof(float);
  for (std::size_t i = 0; i < total_map_elements; ++i) {
    float val = map_flat_data[i];
    sum_scaled_rho += val;
    sum_sq_scaled_rho += static_cast<double>(val) * val;
    if (val < min_scaled_rho)
      min_scaled_rho = val;
    if (val > max_scaled_rho)
      max_scaled_rho = val;
  }

  float mean_scaled_rho = 0.0f;
  float sd_scaled_rho = 0.0f;

  if (total_map_elements > 0) {
    mean_scaled_rho = static_cast<float>(sum_scaled_rho / total_map_elements);
    double variance = (sum_sq_scaled_rho / total_map_elements) -
                      (static_cast<double>(mean_scaled_rho) * mean_scaled_rho);
    sd_scaled_rho =
        (variance > 0) ? static_cast<float>(std::sqrt(variance)) : 0.0f;
  }

  std::cout << "SCALED Map - Min: " << min_scaled_rho
            << ", Max: " << max_scaled_rho << ", Mean: " << mean_scaled_rho
            << ", SD: " << sd_scaled_rho << "\n";

  double sum_signed_rho = 0.0;
  for (std::size_t i = 0; i < total_map_elements; ++i) {
    sum_signed_rho += map_flat_data[i];
  }
  double mean_signed_rho = sum_signed_rho / total_map_elements;
  std::cout << "Mean of signed map density: " << mean_signed_rho << "\n";

  float max_fwt = 0.0;
  glm::ivec3 hkl_strong{};
  for (std::size_t i = 0; i < refln_block.size(); ++i) {
    auto refl_info = refln_block.getReflectionInfo(i);
    if (refl_info && refl_info->fwt > max_fwt) {
      max_fwt = refl_info->fwt;
      hkl_strong = refl_info->hkl;
    }
  }
  std::cout << "Strongest reflection: F(" << hkl_strong.x << "," << hkl_strong.y
            << "," << hkl_strong.z << ") with amplitude " << max_fwt
            << "\n";

  float max_density_peak = -std::numeric_limits<float>::infinity();
  std::size_t peak_flat_idx = 0;

  for (std::size_t i = 0; i < total_map_elements; ++i) {
    if (map_flat_data[i] > max_density_peak) {
      max_density_peak = map_flat_data[i];
      peak_flat_idx = i;
    }
  }

  glm::ivec3 shape(map.dim[0], map.dim[1], map.dim[2]);
  auto i_peak = FlatIndexToHKL(peak_flat_idx, shape);

  std::cout << "Highest density peak in map: " << max_density_peak
            << " at grid index (" << i_peak.x << ", " << i_peak.y << ", "
            << i_peak.z << ")\n";
}

/**
 * @brief Get the Cartesian point from the map state
 * @param ms The object map state
 * @param i, j, k The grid indices
 * @return The Cartesian coordinates of the point
 */
static glm::vec3 GetCartesianPointFromMapState(
    const ObjectMapState& ms, int i, int j, int k)
{
  auto* points_cfield = ms.Field->points.get();
  return glm::vec3(points_cfield->get<float>(i, j, k, 0),
      points_cfield->get<float>(i, j, k, 1),
      points_cfield->get<float>(i, j, k, 2));
}

/**
 * @brief Validate the map voxel geometry
 * @param ms The object map state
 */
void ValidateMapVoxelGeometry(const ObjectMapState& ms)
{
  std::cout << "\n--- Validating Map Voxel Geometry ---\n";

  const float* crystal_dims_arr = ms.Symmetry->Crystal.dims();
  float L_a = crystal_dims_arr[0];
  float L_b = crystal_dims_arr[1];
  float L_c = crystal_dims_arr[2];
  const float* angles_deg = ms.Symmetry->Crystal.angles();

  int N_g0 = ms.FDim[0];
  int N_g1 = ms.FDim[1];
  int N_g2 = ms.FDim[2];

  if (N_g0 <= 1 || N_g1 <= 1 || N_g2 <= 1) {
    std::cout << "Grid dimensions are too small for step validation (<=1 along "
                 "an axis)."
              << "\n";
    return;
  }

  std::cout << std::fixed << std::setprecision(6);
  std::cout << "Cell (from ms.Symmetry->Crystal): a=" << L_a << " b=" << L_b
            << " c=" << L_c << " alpha=" << angles_deg[0]
            << " beta=" << angles_deg[1] << " gamma=" << angles_deg[2]
            << "\n";
  std::cout << "Grid (ms.FDim): N_g0=" << N_g0 << " N_g1=" << N_g1
            << " N_g2=" << N_g2 << "\n";
  std::cout << "Grid Divisions (ms.Div): Div0=" << ms.Div[0]
            << " Div1=" << ms.Div[1] << " Div2=" << ms.Div[2] << "\n";

  float expected_step_g0 =
      (ms.Div[0] > 0) ? L_a / static_cast<float>(ms.Div[0]) : 0.f;
  float expected_step_g1 =
      (ms.Div[1] > 0) ? L_b / static_cast<float>(ms.Div[1]) : 0.f;
  float expected_step_g2 =
      (ms.Div[2] > 0) ? L_c / static_cast<float>(ms.Div[2]) : 0.f;
  std::cout << "Expected Cartesian step length (using Div for sampling along "
               "cell axes): "
            << "dx=" << expected_step_g0 << ", dy=" << expected_step_g1
            << ", dz=" << expected_step_g2 << "\n";

  glm::vec3 P000 = GetCartesianPointFromMapState(ms, 0, 0, 0);
  glm::vec3 P100 = GetCartesianPointFromMapState(ms, 1, 0, 0);
  glm::vec3 P010 = GetCartesianPointFromMapState(ms, 0, 1, 0);
  glm::vec3 P001 = GetCartesianPointFromMapState(ms, 0, 0, 1);

  std::cout << "P(grid 0,0,0): " << glm::to_string(P000) << "\n";
  std::cout << "P(grid 1,0,0): " << glm::to_string(P100) << "\n";
  std::cout << "P(grid 0,1,0): " << glm::to_string(P010) << "\n";
  std::cout << "P(grid 0,0,1): " << glm::to_string(P001) << "\n";

  glm::vec3 step_G0 = P100 - P000;
  glm::vec3 step_G1 = P010 - P000;
  glm::vec3 step_G2 = P001 - P000;

  std::cout << "Cartesian Step vector along grid dim 0 (ms.FDim[0]-axis): "
            << glm::to_string(step_G0) << "\n";
  std::cout << "Cartesian Step vector along grid dim 1 (ms.FDim[1]-axis): "
            << glm::to_string(step_G1) << "\n";
  std::cout << "Cartesian Step vector along grid dim 2 (ms.FDim[2]-axis): "
            << glm::to_string(step_G2) << "\n";

  float len_G0 = glm::length(step_G0);
  float len_G1 = glm::length(step_G1);
  float len_G2 = glm::length(step_G2);

  std::cout << "Magnitude of step_G0: " << len_G0 << "\n";
  std::cout << "Magnitude of step_G1: " << len_G1 << "\n";
  std::cout << "Magnitude of step_G2: " << len_G2 << "\n";

  bool is_cubic_cell =
      (std::abs(angles_deg[0] - 90.0f) < 1e-3f &&
          std::abs(angles_deg[1] - 90.0f) < 1e-3f &&
          std::abs(angles_deg[2] - 90.0f) < 1e-3f &&
          std::abs(L_a - L_b) < 1e-3f && std::abs(L_a - L_c) < 1e-3f);
  bool grid_dims_iso = (N_g0 == N_g1 && N_g0 == N_g2);

  if (is_cubic_cell && grid_dims_iso) {
    if (std::abs(len_G0 - len_G1) > 1e-3f ||
        std::abs(len_G0 - len_G2) > 1e-3f ||
        std::abs(len_G1 - len_G2) > 1e-3f) {
      std::cout << "Warning: Step magnitudes are unequal for a cubic cell with "
                   "isotropic grid sampling! Potential stretching source."
                << "\n";
    } else {
      std::cout << "Step magnitudes are (nearly) equal, as expected for cubic "
                   "cell with isotropic grid."
                << "\n";
    }
  } else {
    std::cout << "Note: Unequal step magnitudes may be expected for non-cubic "
                 "cells or anisotropic grid sampling."
              << "\n";
  }

  glm::vec3 norm_G0 =
      (len_G0 > EPSILON) ? glm::normalize(step_G0) : glm::vec3(0.f);
  glm::vec3 norm_G1 =
      (len_G1 > EPSILON) ? glm::normalize(step_G1) : glm::vec3(0.f);
  glm::vec3 norm_G2 =
      (len_G2 > EPSILON) ? glm::normalize(step_G2) : glm::vec3(0.f);

  float dot_G0G1 = glm::dot(norm_G0, norm_G1);
  float dot_G0G2 = glm::dot(norm_G0, norm_G2);
  float dot_G1G2 = glm::dot(norm_G1, norm_G2);

  std::cout << "Dot product (norm_G0 . norm_G1): " << dot_G0G1
            << " (cos angle between grid axis 0 and 1)\n";
  std::cout << "Dot product (norm_G0 . norm_G2): " << dot_G0G2
            << " (cos angle between grid axis 0 and 2)\n";
  std::cout << "Dot product (norm_G1 . norm_G2): " << dot_G1G2
            << " (cos angle between grid axis 1 and 2)\n";

  float cos_gamma_cell = std::cos(glm::radians(angles_deg[2]));
  float cos_beta_cell = std::cos(glm::radians(angles_deg[1]));
  float cos_alpha_cell = std::cos(glm::radians(angles_deg[0]));

  // approx cos(89 deg) or cos(91 deg) - allow ~1 degree deviation
  float angle_check_threshold = 0.0175f;

  bool angles_match_cell = true;
  if (std::abs(dot_G0G1 - cos_gamma_cell) > angle_check_threshold) {
    std::cout << "Warning: Angle between grid axes 0 & 1 (from dot_G0G1: "
              << std::acos(std::clamp(dot_G0G1, -1.0f, 1.0f)) * 180.f /
                     glm::pi<float>()
              << " deg) inconsistent with crystal gamma (" << angles_deg[2]
              << " deg).\n";
    angles_match_cell = false;
  }
  if (std::abs(dot_G0G2 - cos_beta_cell) > angle_check_threshold) {
    std::cout << "Warning: Angle between grid axes 0 & 2 (from dot_G0G2: "
              << std::acos(std::clamp(dot_G0G2, -1.0f, 1.0f)) * 180.f /
                     glm::pi<float>()
              << " deg) inconsistent with crystal beta (" << angles_deg[1]
              << " deg).\n";
    angles_match_cell = false;
  }
  if (std::abs(dot_G1G2 - cos_alpha_cell) > angle_check_threshold) {
    std::cout << "Warning: Angle between grid axes 1 & 2 (from dot_G1G2: "
              << std::acos(std::clamp(dot_G1G2, -1.0f, 1.0f)) * 180.f /
                     glm::pi<float>()
              << " deg) inconsistent with crystal alpha (" << angles_deg[0]
              << " deg).\n";
    angles_match_cell = false;
  }

  if (!angles_match_cell) {
    std::cout << "Warning: Grid step vectors' angles do not accurately reflect "
                 "cell angles! Potential shear/slant source.\n";
  } else {
    std::cout
        << "Grid step vectors' angles appear consistent with cell angles.\n";
  }

  // Manual FracToReal Check vs. GetPointComponent
  std::cout << "\n--- Manual FracToReal Check vs. GetPointComponent ---\n"
  const float* f2r_matrix = ms.Symmetry->Crystal.fracToReal();
  float frac_v_manual[3];
  float manual_cart_vr[3];

  // Point (grid 1,0,0)
  // Grid indices for fractional calculation: (1+Min[0]), (0+Min[1]), (0+Min[2])
  frac_v_manual[0] =
      (static_cast<float>(1) + ms.Min[0]) / static_cast<float>(ms.Div[0]);
  frac_v_manual[1] =
      (static_cast<float>(0) + ms.Min[1]) / static_cast<float>(ms.Div[1]);
  frac_v_manual[2] =
      (static_cast<float>(0) + ms.Min[2]) / static_cast<float>(ms.Div[2]);
  transform33f3f(f2r_matrix, frac_v_manual, manual_cart_vr);
  std::cout << "Grid (1,0,0): ManualFracToReal: (" << manual_cart_vr[0] << ","
            << manual_cart_vr[1] << "," << manual_cart_vr[2] << ")"
            << " GetPoint: " << glm::to_string(P100) << "\n";
  if (glm::distance(glm::make_vec3(manual_cart_vr), P100) > 1e-4f) {
    std::cout << "MISMATCH for P(1,0,0)!\n";
  }

  // Point (grid 0,1,0)
  frac_v_manual[0] =
      (static_cast<float>(0) + ms.Min[0]) / static_cast<float>(ms.Div[0]);
  frac_v_manual[1] =
      (static_cast<float>(1) + ms.Min[1]) / static_cast<float>(ms.Div[1]);
  frac_v_manual[2] =
      (static_cast<float>(0) + ms.Min[2]) / static_cast<float>(ms.Div[2]);
  transform33f3f(f2r_matrix, frac_v_manual, manual_cart_vr);
  std::cout << "Grid (0,1,0): ManualFracToReal: (" << manual_cart_vr[0] << ","
            << manual_cart_vr[1] << "," << manual_cart_vr[2] << ")"
            << " GetPoint: " << glm::to_string(P010) << "\n";
  if (glm::distance(glm::make_vec3(manual_cart_vr), P010) > 1e-4f) {
    std::cout << "MISMATCH for P(0,1,0)!\n";
  }

  // Point (grid 0,0,1)
  frac_v_manual[0] =
      (static_cast<float>(0) + ms.Min[0]) / static_cast<float>(ms.Div[0]);
  frac_v_manual[1] =
      (static_cast<float>(0) + ms.Min[1]) / static_cast<float>(ms.Div[1]);
  frac_v_manual[2] =
      (static_cast<float>(1) + ms.Min[2]) / static_cast<float>(ms.Div[2]);
  transform33f3f(f2r_matrix, frac_v_manual, manual_cart_vr);
  std::cout << "Grid (0,0,1): ManualFracToReal: (" << manual_cart_vr[0] << ","
            << manual_cart_vr[1] << "," << manual_cart_vr[2] << ")"
            << " GetPoint: " << glm::to_string(P001) << "\n";
  if (glm::distance(glm::make_vec3(manual_cart_vr), P001) > 1e-4f) {
    std::cout << "MISMATCH for P(0,0,1)!\n";
  }
  std::cout << "--- End of Manual FracToReal Check ---\n";
  std::cout << "--- End of Voxel Geometry Validation ---\n";
}
#endif // TRACE

/**
 * @brief Create an ObjectMapState from a CFieldTyped map and symmetry
 * @param map The map to convert
 * @param sym The symmetry to use
 * @return The created ObjectMapState
 */
ObjectMapState ObjectMapStateFromField(PyMOLGlobals* G,
    const CFieldTyped<float>& map, std::unique_ptr<CSymmetry> sym)
{
  ObjectMapState ms(G);
  glm::ivec3 grid_size(map.dim[0], map.dim[1], map.dim[2]);
  ms.MapSource = cMapSourceCCP4;
  ms.Active = true;
  ms.Symmetry.reset(sym.release());
  ms.FDim[0] = grid_size[0];
  ms.FDim[1] = grid_size[1];
  ms.FDim[2] = grid_size[2];
  ms.FDim[3] = 3;

  ms.Min[0] = 0;
  ms.Min[1] = 0;
  ms.Min[2] = 0;
  ms.Max[0] = grid_size.x - 1;
  ms.Max[1] = grid_size.y - 1;
  ms.Max[2] = grid_size.z - 1;

  ms.Div[0] = grid_size.x;
  ms.Div[1] = grid_size.y;
  ms.Div[2] = grid_size.z;

  ms.Field.reset(new Isofield(G, ms.FDim));

  float v[3];
  v[2] = (ms.Min[2]) / ((float) ms.Div[2]);
  v[1] = (ms.Min[1]) / ((float) ms.Div[1]);
  v[0] = (ms.Min[0]) / ((float) ms.Div[0]);

  transform33f3f(ms.Symmetry->Crystal.fracToReal(), v, ms.ExtentMin);

  v[2] = ((ms.FDim[2] - 1) + ms.Min[2]) / ((float) ms.Div[2]);
  v[1] = ((ms.FDim[1] - 1) + ms.Min[1]) / ((float) ms.Div[1]);
  v[0] = ((ms.FDim[0] - 1) + ms.Min[0]) / ((float) ms.Div[0]);

  transform33f3f(ms.Symmetry->Crystal.fracToReal(), v, ms.ExtentMax);

  auto& field = *ms.Field;
  std::copy(map.data.begin(), map.data.end(), field.data->data.begin());

  ms.Field->save_points = false;

  auto& crystal = ms.Symmetry->Crystal;

  ObjectMapStateRegeneratePoints(&ms);
  float frac_coords[3];
  float cart_coords[3];
  int corner_write_idx = 0;

  for (int z_loop = 0; z_loop < 2; ++z_loop) {
    int grid_idx_z = (z_loop == 0) ? ms.Min[2] : (ms.Min[2] + ms.FDim[2] - 1);
    frac_coords[2] =
        static_cast<float>(grid_idx_z) / static_cast<float>(ms.Div[2]);

    for (int y_loop = 0; y_loop < 2; ++y_loop) {
      int grid_idx_y = (y_loop == 0) ? ms.Min[1] : (ms.Min[1] + ms.FDim[1] - 1);
      frac_coords[1] =
          static_cast<float>(grid_idx_y) / static_cast<float>(ms.Div[1]);

      for (int x_loop = 0; x_loop < 2; ++x_loop) {
        int grid_idx_x =
            (x_loop == 0) ? ms.Min[0] : (ms.Min[0] + ms.FDim[0] - 1);
        frac_coords[0] =
            static_cast<float>(grid_idx_x) / static_cast<float>(ms.Div[0]);

        transform33f3f(
            ms.Symmetry->Crystal.fracToReal(), frac_coords, cart_coords);

        ms.Corner[corner_write_idx * 3 + 0] = cart_coords[0];
        ms.Corner[corner_write_idx * 3 + 1] = cart_coords[1];
        ms.Corner[corner_write_idx * 3 + 2] = cart_coords[2];
        corner_write_idx++;
      }
    }
  }
#ifdef TRACE
  ValidateMapVoxelGeometry(ms);
#endif // TRACE
  return ms;
}

/**
 * @brief Find the optimal grid size for the FFT
 * @param refln_block The reflection block containing reflection data
 * @param reciprocal_space_params The parameters for reciprocal space
 * @return The optimal grid size as a glm::ivec3
 */
glm::ivec3 FindOptimalGridSize(const ReflnBlock& refln_block,
    const ReciprocalSpaceParams& reciprocal_space_params)
{
  std::array<int, 3> dimensions{};
  float inv_dist_sq{};
  for (std::size_t i = 0; i < refln_block.size(); ++i) {
    auto refl_info = refln_block.getReflectionInfo(i);
    if (!refl_info) {
      continue;
    }

    auto cand_inv_dist_sq = CalculateReciprocalSpaceDistanceSquared(
        reciprocal_space_params, refl_info->hkl);
    inv_dist_sq = std::max(inv_dist_sq, cand_inv_dist_sq);
  }
  double inv_dist = std::sqrt(inv_dist_sq);
  // Nyquist
  double nyquist_sampling_rate = 2.0; // Perhaps make into a setting?
  glm::dvec3 cell_size{};
  for (int i = 0; i < 3; ++i) {
    cell_size[i] = std::max(cell_size[i],
        nyquist_sampling_rate * inv_dist / reciprocal_space_params.abc[i]);
  }

  glm::ivec3 grid_size{};
  for (int i = 0; i < 3; ++i) {
    auto round_up = static_cast<int>(std::ceil(cell_size[i]));
    grid_size[i] = pocketfft::detail::util::good_size_real(round_up);
  }
  return grid_size;
}

/**
 * @brief Calculate the phase shift for a given HKL and symmetry operation
 * @param hkl_in The input HKL indices
 * @param sym_op_T_3x1 The symmetry operation translation vector
 * @return The calculated phase shift (2 * PI * (h.trans))
 */
double CalculatePhaseShift(
    const glm::ivec3& hkl, const glm::vec3& trans)
{
  double dot_product = glm::dot(glm::dvec3(hkl), glm::dvec3(trans));
  return 2.0 * glm::pi<double>() * dot_product;
}

/**
 * @brief Normalizes the complex map to float density map
 * @param reciprocal_space_params The parameters for reciprocal space
 * @param recip_grid The reciprocal grid data
 * @return Density Map as Float Field
 */
static CFieldTyped<float> NormalizeComplexMapToFloat(PyMOLGlobals* G,
    const ReciprocalSpaceParams& reciprocal_space_params,
    const CFieldTyped<std::complex<float>>& recip_grid)
{
  CFieldTyped<float> map(reinterpret_cast<const int*>(recip_grid.dim.data()), 3);
  auto total_map_elements = map.size() / sizeof(float);

  float inv_Vcell = 1.0f;
  if (std::abs(reciprocal_space_params.volume) > EPSILON) {
    inv_Vcell = 1.0f / static_cast<float>(reciprocal_space_params.volume);
  }

  auto* flat_recip_data = recip_grid.ptr(0, 0, 0);
  auto* map_flat_data = map.ptr(0, 0, 0);

  if (!SettingGet<bool>(G, cSetting_normalize_ccp4_maps)) {
    for (std::size_t i = 0; i < total_map_elements; ++i) {
      map_flat_data[i] = flat_recip_data[i].real() * inv_Vcell;
    }
    return map;
  }

  double sum_rho_scaled = 0.0;
  double sum_sq_rho_scaled = 0.0;
  for (std::size_t i = 0; i < total_map_elements; ++i) {
    map_flat_data[i] = flat_recip_data[i].real() * inv_Vcell;
    sum_rho_scaled += map_flat_data[i];
    sum_sq_rho_scaled +=
        static_cast<double>(map_flat_data[i]) * map_flat_data[i];
  }

  // Normalize the map to mean ~0 and standard deviation ~1
  float current_mean = 0.0f;
  float current_sd = 0.0f;

  if (total_map_elements > 0) {
    current_mean = static_cast<float>(sum_rho_scaled / total_map_elements);
    double variance = (sum_sq_rho_scaled / total_map_elements) -
                      (static_cast<double>(current_mean) * current_mean);
    current_sd =
        (variance > 0) ? static_cast<float>(std::sqrt(variance)) : 0.0f;
  }

  if (std::abs(current_sd) > EPSILON) { // Avoid division by zero if SD is tiny
    for (std::size_t i = 0; i < total_map_elements; ++i) {
      map_flat_data[i] = (map_flat_data[i] - current_mean) / current_sd;
    }
  } else {
    // Shift everything to force mean to zero
    if (std::abs(current_mean) > EPSILON) {
      for (std::size_t i = 0; i < total_map_elements; ++i) {
        map_flat_data[i] = map_flat_data[i] - current_mean;
      }
    }
  }
  return map;
}

/**
 * @brief Structure to hold symmetry matrices in GLM format
 */
struct SymMatsGLM {
  glm::mat3 rot;
  glm::vec3 trans;
};

/**
 * @brief Calculate the symmetry matrices in GLM format.
 * @param symmetry The symmetry object containing the symmetry matrices.
 * @return A vector of SymMatsGLM structures containing the rotation and
 * translation components of the symmetry matrices.
 */
static std::vector<SymMatsGLM> CalculateSymMatricesGLM(
    const CSymmetry& symmetry)
{
  std::vector<SymMatsGLM> sym_matrices(symmetry.getNSymMat());
  for (int i = 0; i < symmetry.getNSymMat(); ++i) {
    auto* sym_matrix_4x4 = symmetry.getSymMat(i);
    const float R_j[9] = {sym_matrix_4x4[0], sym_matrix_4x4[1],
        sym_matrix_4x4[2], sym_matrix_4x4[4], sym_matrix_4x4[5],
        sym_matrix_4x4[6], sym_matrix_4x4[8], sym_matrix_4x4[9],
        sym_matrix_4x4[10]};
    sym_matrices[i].rot = glm::transpose(glm::make_mat3(R_j));
    sym_matrices[i].trans =
        glm::vec3(sym_matrix_4x4[3], sym_matrix_4x4[7], sym_matrix_4x4[11]);
  }
  return sym_matrices;
}

/**
 * @brief Checks if HKL indices are within the bounds of a grid.
 * @param idx The HKL indices to check.
 * @param shape The shape of the grid.
 * @return True if the indices are within bounds, false otherwise.
 */
static bool IsHKLWithinShape(const glm::ivec3& idx, const glm::ivec3& shape)
{
  return (idx.x >= 0 && idx.x < shape[0] && //
          idx.y >= 0 && idx.y < shape[1] && //
          idx.z >= 0 && idx.z < shape[2]);
}

/**
 * @brief Apply symmetry matrices to the reflections and store the results in
 * the reciprocal space data.
 * @param sym_matrices The symmetry matrices to apply.
 * @param hkl_orig The original Miller indices.
 * @param shape The shape of the reciprocal space grid.
 * @param f_amp_weighted The weighted amplitude of the reflection.
 * @param phase_rad_orig The original phase in radians.
 * @param recip_1D_total_elements The total number of elements in the 1D
 * reciprocal space data.
 * @param[out] flat_recip_data The flat reciprocal space data array.
 * @param[out] is_grid_point_set A vector indicating whether a grid point has been
 * set.
 */
void ApplySymMatricesToReflections(
    const std::vector<SymMatsGLM>& sym_matrices, const glm::ivec3& hkl_orig,
    const glm::ivec3& shape,
    float f_amp_weighted, float phase_rad_orig, int recip_1D_total_elements,
    std::complex<float>* flat_recip_data,
    std::vector<bool>& is_grid_point_set)
{
  for (auto& sym_mat : sym_matrices) {
    glm::ivec3 hkl_symm = glm::round(sym_mat.rot * hkl_orig);
    auto phase_shift_rad = CalculatePhaseShift(hkl_orig, sym_mat.trans);

    auto F_symm_complex = std::polar(
        f_amp_weighted, phase_rad_orig + static_cast<float>(phase_shift_rad));

    int idx_h_s = (hkl_symm.x >= 0) ? hkl_symm.x : hkl_symm.x + shape[0];
    int idx_k_s = (hkl_symm.y >= 0) ? hkl_symm.y : hkl_symm.y + shape[1];
    int idx_l_s = (hkl_symm.z >= 0) ? hkl_symm.z : hkl_symm.z + shape[2];
    auto idx_hkl_s = glm::ivec3(idx_h_s, idx_k_s, idx_l_s);

    if (IsHKLWithinShape(idx_hkl_s, shape)) {
      auto flat_idx_s =
          HKLToFlatIndex(idx_hkl_s, shape);
      if (flat_idx_s < recip_1D_total_elements &&
          !is_grid_point_set[flat_idx_s]) {
        flat_recip_data[flat_idx_s] = F_symm_complex;
        is_grid_point_set[flat_idx_s] = true;
      }
    }
  }
}

/**
 * @brief Calculate inverse Friedel mate indices for a given HKL and shape.
 * @param hkl_orig The original HKL indices.
 * @param shape The shape of the reciprocal space grid.
 * @return The Friedel mate indices as a glm::ivec3.
 */
glm::ivec3 CalculateFriedelMateIndex(
    const glm::ivec3& hkl_orig, const glm::ivec3& shape)
{
  auto h_idx = hkl_orig.x;
  auto k_idx = hkl_orig.y;
  auto l_idx = hkl_orig.z;

  int h =
      (h_idx < shape[0] / 2 + shape[0] % 2) ? h_idx : h_idx - shape[0];
  int k =
      (k_idx < shape[1] / 2 + shape[1] % 2) ? k_idx : k_idx - shape[1];
  int l =
      (l_idx < shape[2] / 2 + shape[2] % 2) ? l_idx : l_idx - shape[2];

  // Friedel mate indices
  int h_bar = -h;
  int k_bar = -k;
  int l_bar = -l;

  // Map (-h,-k,-l) to Friedel grid indices
  int idx_h_bar = (h_bar >= 0) ? h_bar : h_bar + shape[0];
  int idx_k_bar = (k_bar >= 0) ? k_bar : k_bar + shape[1];
  int idx_l_bar = (l_bar >= 0) ? l_bar : l_bar + shape[2];
  return glm::ivec3(idx_h_bar, idx_k_bar, idx_l_bar);
}

/**
 * @brief Add friedel pair to the reciprocal space data.
 * @param symmetry The symmetry object.
 * @param refln_block The reflection block containing reflection data.
 * @param reciprocal_space_params The parameters for reciprocal space.
 * @return A CFieldTyped<float> object representing the Fourier transformed map.
 */
void AddFriedelMate(const glm::ivec3& hkl_orig,
    std::size_t recip_1D_total_elements, const glm::ivec3& shape,
    std::complex<float>* flat_recip_data, std::vector<bool>& is_grid_point_set)
{
  auto flat_idx_hkl = HKLToFlatIndex(hkl_orig, shape);

  auto idx_bar = CalculateFriedelMateIndex(hkl_orig, shape);

  if (IsHKLWithinShape(idx_bar, shape)) {
    auto flat_idx_bar = HKLToFlatIndex(idx_bar, shape);

    if (flat_idx_bar < recip_1D_total_elements &&
        !is_grid_point_set[flat_idx_bar]) {
      flat_recip_data[flat_idx_bar] =
          std::conj(flat_recip_data[flat_idx_hkl]);
      is_grid_point_set[flat_idx_bar] = true;
    }
    // If flat_idx_hkl == flat_idx_bar (centrosymmetric FFT point),
    // ensure it's real
    bool centrosymmetric = flat_idx_hkl == flat_idx_bar;
    if (centrosymmetric) {
      auto& recip_data = flat_recip_data[flat_idx_hkl];
      bool around_zero_and_much_less_than_real =
          std::abs(recip_data.imag()) > 1e-5f * std::abs(recip_data.real()) &&
          std::abs(recip_data.imag()) > 1e-5f;
      if (around_zero_and_much_less_than_real) {
        recip_data.imag(0.0f);
      }
    }
  }
}

/**
 * @brief Perform a Fast Fourier Transform in place on the reciprocal grid.
 * @param grid The reciprocal grid data.
 */
static void PerformFastFourierTransformInPlace(
    CFieldTyped<std::complex<float>>& grid)
{
  auto* flat_grid = grid.ptr(0, 0, 0);
  pocketfft::shape_t pocket_shape{static_cast<std::size_t>(grid.dim[0]),
      static_cast<std::size_t>(grid.dim[1]),
      static_cast<std::size_t>(grid.dim[2])};
  auto base_size = grid.base_size;

  // CField is internally a 1D array, so we need to calculate the strides
  // for the 3D data. Z is contiguous; thus the axis order is ZYX.
  glm::uvec3 u_strides = {base_size * pocket_shape[2] * pocket_shape[1],
      base_size * pocket_shape[2], base_size};
  pocketfft::stride_t stride_in(3);

  for (int i = 0; i < 3; ++i) {
    stride_in[i] = static_cast<pocketfft::stride_t::value_type>(u_strides[i]);
  }
  auto stride_out = stride_in;
  pocketfft::shape_t axes{2, 1, 0}; // ZYX
  bool direction = pocketfft::BACKWARD;
  float norm = 1.0f;

  pocketfft::c2c<float>(pocket_shape, stride_in, stride_out, axes, direction,
      flat_grid, flat_grid, norm);
}

/**
 * @brief Fourier transform structure factors to a map.
 * @param symmetry The symmetry object.
 * @param refln_block The reflection block containing reflection data.
 * @param reciprocal_space_params The parameters for reciprocal space.
 * @return A CFieldTyped<float> object representing the Fourier transformed map.
 */
CFieldTyped<float> FourierTransformStructureFactorsToMap(PyMOLGlobals* G,
    const CSymmetry& symmetry,
    const ReflnBlock& refln_block,
    const ReciprocalSpaceParams& reciprocal_space_params)
{
  auto shape = FindOptimalGridSize(refln_block, reciprocal_space_params);
  CFieldTyped<std::complex<float>> recip_grid(glm::value_ptr(shape), 3);

  auto recip_1D_total_elements = shape[0] * shape[1] * shape[2];
  std::vector<bool> is_grid_point_set(recip_1D_total_elements, false);

  auto* flat_recip_data = recip_grid.ptr(0, 0, 0);

  auto sym_matrices = CalculateSymMatricesGLM(symmetry);

  for (std::size_t i = 0; i < refln_block.size(); ++i) {
    auto refl_info = refln_block.getReflectionInfo(i);
    if (!refl_info) {
      continue;
    }

    auto hkl_orig = refl_info->hkl;
    auto f_amp_raw = refl_info->fwt;
    auto ph_deg_orig = refl_info->phwt;
    auto fom_val = refl_info->fom;

    auto f_amp_weighted = f_amp_raw * fom_val;
    if (f_amp_weighted == 0.0f &&
        !(hkl_orig.x == 0 && hkl_orig.y == 0 &&
            hkl_orig.z == 0)) { // Skip zero amplitude reflections unless F000
      continue;
    }

    auto phase_rad_orig = glm::radians(ph_deg_orig);
    auto F_orig_complex = std::polar(f_amp_weighted, phase_rad_orig);

    ApplySymMatricesToReflections(
        sym_matrices, hkl_orig, shape, f_amp_weighted, phase_rad_orig,
        recip_1D_total_elements, flat_recip_data, is_grid_point_set);
  } // End loop over reflections from CIF

  for (std::size_t h_idx = 0; h_idx < shape[0]; ++h_idx) {
    for (std::size_t k_idx = 0; k_idx < shape[1]; ++k_idx) {
      for (std::size_t l_idx = 0; l_idx < shape[2]; ++l_idx) {
        auto flat_idx_hkl =
            HKLToFlatIndex(glm::ivec3(h_idx, k_idx, l_idx), shape);
        if (is_grid_point_set[flat_idx_hkl]) {
          AddFriedelMate(glm::ivec3(h_idx, k_idx, l_idx), recip_1D_total_elements,
              shape, flat_recip_data, is_grid_point_set);
        }
      }
    }
  }

  PerformFastFourierTransformInPlace(recip_grid);
  auto map = NormalizeComplexMapToFloat(G, reciprocal_space_params, recip_grid);
#ifdef TRACE
  PrintReciprocalGridStatistics(recip_grid);
  PrintRawMapStatistics(refln_block, map);
#endif // TRACE
  return map;
}

pymol::Result<ObjectMap*> ObjectMapReadCifStr(
    PyMOLGlobals* G, const pymol::cif_data& datablock)
{
  auto I = new ObjectMap(G);
  initializeTTT44f(I->TTT);
  I->TTTFlag = false;

  ReflnBlock refln_block(datablock);
  std::unique_ptr<CSymmetry> sym(read_symmetry(G, &datablock));

  auto reciprocal_space_params = CalculateReciprocalSpaceParam(sym->Crystal);
  auto map = FourierTransformStructureFactorsToMap(
      G, *sym, refln_block, reciprocal_space_params);

  I->State.push_back(ObjectMapStateFromField(G, map, std::move(sym)));
  ObjectMapUpdateExtents(I);

  return I;
}

/**
 * Read one or multiple object-molecules from a CIF file. If there is only one
 * or multiplex=0, then return the object-molecule. Otherwise, create each
 * object - named by its data block name - and return nullptr.
 */
pymol::Result<ObjectMolecule*> ObjectMoleculeReadCifStr(PyMOLGlobals * G, ObjectMolecule * I,
                                      const char *st, int frame,
                                      int discrete, int quiet, int multiplex,
                                      int zoom)
{
  if (I) {
    return pymol::make_error("loading mmCIF into existing object not supported, "
                        "please use 'create' to append to an existing object.");
  }

  if (multiplex > 0) {
    return pymol::make_error("loading mmCIF with multiplex=1 not supported, please "
                        "use 'split_states' after loading the object.");
  }

  auto cif = std::make_shared<cif_file_with_error_capture>();
  if (!cif->parse_string(st)) {
    return pymol::make_error("Parsing CIF file failed: ", cif->m_error_msg);
  }

  for (const auto& [code, datablock] : cif->datablocks()) {
    ObjectMolecule * obj = ObjectMoleculeReadCifData(G, &datablock, discrete, quiet);

    if (!obj) {
      PRINTFB(G, FB_ObjectMolecule, FB_Warnings)
        " mmCIF-Warning: no coordinates found in data_%s\n", datablock.code() ENDFB(G);
      continue;
    }

#ifndef _PYMOL_NOPY
    // we only provide access from the Python API so far
    if (SettingGetGlobal_b(G, cSetting_cif_keepinmemory)) {
      obj->m_cifdata = &datablock;
      obj->m_ciffile = cif;
    }
#endif

    if (cif->datablocks().size() == 1 || multiplex == 0)
      return obj;

    // multiplexing
    ObjectSetName(obj, datablock.code());
    ExecutiveDelete(G, obj->Name);
    ExecutiveManageObject(G, obj, zoom, true);
  }

  return nullptr;
}

/**
 * Bond dictionary getter, with on-demand download of residue dictionaries
 */
const bond_dict_t::mapped_type * bond_dict_t::get(PyMOLGlobals * G, const char * resn, bool try_download) {
  auto key = make_key(resn);
  auto it = m_map.find(key);

  if (it != m_map.end())
    return &it->second;

  if (unknown_resn.count(key))
    return nullptr;

#ifndef _PYMOL_NOPY
  if (try_download) {
    pymol::GIL_Ensure gil;
    bool downloaded = false;

    // call into Python
    unique_PyObject_ptr pyfilename(
        PyObject_CallMethod(G->P_inst->cmd, "download_chem_comp", "siO", resn,
            !Feedback(G, FB_Executive, FB_Details), G->P_inst->cmd));

    if (pyfilename) {
      const char* filename = PyString_AsString(pyfilename.get());

      // update
      if ((downloaded = (filename && filename[0]))) {
        cif_file_with_error_capture cif;
        if (!cif.parse_file(filename)) {
          PRINTFB(G, FB_Executive, FB_Warnings)
            " Warning: Loading _chem_comp_bond CIF data for residue '%s' "
            "failed: %s\n",
            resn, cif.m_error_msg.c_str() ENDFB(G);
          return nullptr;
        }

        for (auto& [code, item] : cif.datablocks())
          read_chem_comp_bond_dict(&item, *this);
      }
    }

    if (downloaded) {
      // second attempt to look up, from eventually updated dictionary
      return get(G, resn, false);
    }
  }
#endif

  PRINTFB(G, FB_Executive, FB_Warnings)
    " ExecutiveLoad-Warning: No _chem_comp_bond data for residue '%s'\n", resn
    ENDFB(G);

  // don't try downloading again
  unknown_resn.insert(key);

  return nullptr;
}


///////////////////////////////////////

pymol::Result<ObjectMolecule*> ObjectMoleculeReadBCif(PyMOLGlobals* G,
    ObjectMolecule* I, const char* bytes, std::size_t size, int frame,
    int discrete, int quiet, int multiplex, int zoom)
{
#ifdef _PYMOL_NO_MSGPACKC
  PRINTFB(G, FB_ObjectMolecule, FB_Errors)
    " Error: This build has no BinaryCIF support.\n"
    " Please install/enable msgpack-c.\n"
  ENDFB(G);
  return nullptr;
#endif

  if (I) {
    return pymol::make_error("loading BCIF into existing object not supported, "
                        "please use 'create' to append to an existing object.");
  }

  if (multiplex > 0) {
    return pymol::make_error("loading BCIF with multiplex=1 not supported, please "
                        "use 'split_states' after loading the object.");
  }

  auto cif = std::make_shared<pymol::cif_file>();
  cif->parse_bcif(bytes, size);
  
  for (const auto& [code, datablock] : cif->datablocks()) {
    auto obj = ObjectMoleculeReadCifData(G, &datablock, discrete, quiet);
    if (!obj) {
      PRINTFB(G, FB_ObjectMolecule, FB_Warnings)
        " BCIF-Warning: no coordinates found in data_%s\n", datablock.code() ENDFB(G);
      continue;
    }

#ifndef _PYMOL_NOPY
    // we only provide access from the Python API so far
    if (SettingGet<bool>(G, cSetting_cif_keepinmemory)) {
      obj->m_cifdata = &datablock;
      obj->m_ciffile = cif;
    }
#endif

    if (cif->datablocks().size() == 1 || multiplex == 0)
      return obj;
  }
  return nullptr;
}

// vi:sw=2:ts=2:expandtab