atomgraph.h

// -*-c++-*-
#ifndef STORMM_TOPOLOGY_H
#define STORMM_TOPOLOGY_H
 
#include <cmath>
#include <vector>
#include <string>
#include <time.h>
#include "copyright.h"
#include "Constants/behavior.h"
#include "Constants/hpc_bounds.h"
#include "Constants/symbol_values.h"
#include "Accelerator/hybrid.h"
#include "DataTypes/common_types.h"
#include "DataTypes/stormm_vector_types.h"
#include "FileManagement/file_enumerators.h"
#include "Parsing/citation.h"
#include "Parsing/parsing_enumerators.h"
#include "Structure/structure_enumerators.h"
#include "atomgraph_abstracts.h"
#include "atomgraph_constants.h"
#include "atomgraph_enumerators.h"
#include "atomgraph_refinement.h"
#include "topology_limits.h"
 
namespace stormm {
namespace topology {
 
using constants::ExceptionResponse;
using constants::PrecisionModel;
using card::Hybrid;
using card::HybridTargetLevel;
using diskutil::PrintSituation;
using parse::Citation;
using parse::NumberFormat;
using parse::PolyNumeric;
using structure::ApplyConstraints;
using symbols::amber_ancient_bioq;

class AtomGraph {
public:

  AtomGraph();

  AtomGraph(const std::string &file_name, ExceptionResponse policy = ExceptionResponse::WARN,
            TopologyKind engine_format = TopologyKind::AMBER);
 
  AtomGraph(const std::string &file_name, ExceptionResponse policy, TopologyKind engine_format,
            double coulomb_constant_in, double default_elec14_screening,
            double default_vdw14_screening, double charge_rounding_tol,
            double charge_discretization,
            ApplyConstraints use_bond_constraints_in = ApplyConstraints::NO,
            ApplyConstraints use_settle_in = ApplyConstraints::NO);
 
  AtomGraph(const std::vector<AtomGraph*> &agv, const std::vector<int> &counts,
            MoleculeOrdering arrangement = MoleculeOrdering::REORDER_ALL,
            ExceptionResponse policy = ExceptionResponse::DIE);
 
  AtomGraph(const AtomGraph &ag_a, int n_a, const AtomGraph &ag_b, int n_b,
            MoleculeOrdering arrangement = MoleculeOrdering::REORDER_ALL,
            ExceptionResponse policy = ExceptionResponse::DIE);
 
  AtomGraph(const AtomGraph &ag_a, const AtomGraph &ag_b, int n_b,
            MoleculeOrdering arrangement = MoleculeOrdering::REORDER_ALL,
            ExceptionResponse policy = ExceptionResponse::DIE);
 
  AtomGraph(const AtomGraph &ag_a, const AtomGraph &ag_b,
            MoleculeOrdering arrangement = MoleculeOrdering::REORDER_ALL,
            ExceptionResponse policy = ExceptionResponse::DIE);
 
  AtomGraph(const AtomGraph &original, const std::vector<int> &atom_subset,
            ExceptionResponse policy = ExceptionResponse::DIE,
            double charge_rounding_tol = default_charge_rounding_tol,
            double charge_discretization = default_charge_precision_inc,
            ApplyConstraints use_bond_constraints_in = ApplyConstraints::NO,
            ApplyConstraints use_settle_in = ApplyConstraints::NO);
 
  AtomGraph(const AtomGraph &original, const std::vector<bool> &mask,
            ExceptionResponse policy = ExceptionResponse::DIE,
            double charge_rounding_tol = default_charge_rounding_tol,
            double charge_discretization = default_charge_precision_inc,
            ApplyConstraints use_bond_constraints_in = ApplyConstraints::NO,
            ApplyConstraints use_settle_in = ApplyConstraints::NO);

  ~AtomGraph() = default;

  AtomGraph(const AtomGraph &original);

  AtomGraph& operator=(const AtomGraph &other);

  AtomGraph(AtomGraph &&original);

  AtomGraph& operator=(AtomGraph &&other);

  void buildFromPrmtop(const std::string &file_name,
                       const ExceptionResponse policy = ExceptionResponse::WARN,
                       double coulomb_constant_in = amber_ancient_bioq,
                       double default_elec14_screening = amber_default_elec14_screen,
                       double default_vdw14_screening = amber_default_vdw14_screen,
                       double charge_rounding_tol = default_charge_rounding_tol,
                       double charge_discretization = default_charge_precision_inc);
 
  // Most getter functions (getThisStuff) will return the corresponding member variable.  Usually
  // a variable with a single value is an int, but it could also be an enumeration or a string.
  // For Hybrid object member variables, getThisHybridObjectStuff() will return the entirety of
  // the POINTER-kind object's host_data vector, up to its stated length.

  std::string getTitle() const;
  
  std::string getFileName() const;

  int getAtomCount() const;

  int getResidueCount() const;

  int getMoleculeCount() const;

  MoleculeKind getMoleculeKind(int mol_index) const;
  
  int getOrganicMoleculeCount() const;

  int getOrganicAtomCount() const;

  int getWaterResidueSize() const;

  int getWaterResidueCount() const;
  
  int getLargestResidueSize() const;

  int getLastSoluteResidue() const;

  int getLastSoluteAtom() const;

  int getFirstSolventMolecule() const;

  int getFirstSolventAtom() const;

  int getLargestMoleculeSize() const;

  double getTotalMass() const;

  int getDegreesOfFreedom() const;

  int getConstrainedDegreesOfFreedom() const;
  
  int getDescriptor(TopologyDescriptor choice) const;
  int getDescriptor(SanderDescriptor choice) const;

  const Hybrid<int>& getResidueLimits() const;
  int2 getResidueLimits(int index) const;

  std::vector<int> getResidueIndex() const;
  int getResidueIndex(int atom_index) const;

  std::vector<int> getStructuralAtomNumber() const;
  std::vector<int> getStructuralAtomNumber(int low_index, int high_index) const;
  int getStructuralAtomNumber(int index) const;

  std::vector<int> getResidueNumber() const;
  std::vector<int> getResidueNumber(int low_index, int high_index) const;
  int getResidueNumber(int index) const;

  std::vector<int> getMoleculeLimits() const;
  int2 getMoleculeLimits(int index) const;

  std::vector<int> getParticlesPerMolecule() const;
  int getParticlesPerMolecule(int index) const;

  std::vector<int> getAtomicNumber() const;
  std::vector<int> getAtomicNumber(int low_index, int high_index) const;
  int getAtomicNumber(int index) const;

  std::vector<bool> getAtomMobility() const;
  std::vector<bool> getAtomMobility(int low_index, int high_index) const;
  bool getAtomMobility(int index) const;

  std::vector<uint> getAtomMobilityMask() const;
  std::vector<uint> getAtomMobilityMask(int low_index, int high_index) const;
  
  void modifyAtomMobility(MobilitySetting movement);
  void modifyAtomMobility(int low_index, int high_index, MobilitySetting movement);
  void modifyAtomMobility(int index, MobilitySetting movement);
  void modifyAtomMobility(const std::vector<int> &mask, MobilitySetting movement);

  std::vector<int> getMoleculeMembership() const;
  std::vector<int> getMoleculeMembership(int low_index, int high_index) const;
  int getMoleculeMembership(int index) const;

  std::vector<int> getMoleculeContents() const;
  std::vector<int> getMoleculeContents(int index) const;

  template <typename T> std::vector<T> getPartialCharge() const;
  template <typename T> std::vector<T> getPartialCharge(int low_index, int high_index) const;
  template <typename T> T getPartialCharge(int index) const;

  template <typename T> std::vector<T> getAtomicMass(MassForm rep = MassForm::ORDINARY) const;
  template <typename T> std::vector<T> getAtomicMass(int low_index, int high_index,
                                                     MassForm rep = MassForm::ORDINARY) const;
  template <typename T> T getAtomicMass(int index, MassForm rep = MassForm::ORDINARY) const;

  const Hybrid<char4>& getAtomName() const;
  std::vector<char4> getAtomName(int low_index, int high_index) const;
  char4 getAtomName(int index) const;

  const Hybrid<char4>& getAtomType() const;
  std::vector<char4> getAtomType(int low_index, int high_index) const;
  char4 getAtomType(int index) const;

  std::string getFullAtomName(int index) const;
  
  std::vector<std::vector<char4>> getAtomTypeNameTable() const;
  
  const Hybrid<char4>& getResidueName() const;
  char4 getResidueName(int index) const;

  int getUreyBradleyTermCount() const;

  int getCharmmImprTermCount() const;

  int getCmapTermCount() const;

  int getUreyBradleyParameterCount() const;

  int getCharmmImprParameterCount() const;

  int getCmapSurfaceCount() const;

  int getCmapDimension(int index) const;

  template <typename T> UreyBradleyTerm<T> getUreyBradleyTerm(int index) const;

  template <typename T> CharmmImprTerm<T> getCharmmImprTerm(int index) const; 

  template <typename T> CmapTerm<T> getCmapTerm(int index) const;

  int getBondTermCount() const;

  int getAngleTermCount() const;

  int getDihedralTermCount() const;

  int getBondParameterCount() const;

  int getAngleParameterCount() const;

  int getDihedralParameterCount() const;

  template <typename T> BondTerm<T> getBondTerm(int index) const;

  template <typename T> AngleTerm<T> getAngleTerm(int index) const;

  template <typename T> DihedralTerm<T> getDihedralTerm(int index) const;

  int getVirtualSiteCount() const;

  int getVirtualSiteParameterSetCount() const;
  
  std::vector<int> findVirtualSites() const;
  std::vector<int2> findVirtualSites(int low_index, int high_index) const;
  int findVirtualSites(int index) const;

  int getVirtualSiteIndex(int atom_index, ExceptionResponse policy = ExceptionResponse::DIE) const;
  
  VirtualSiteKind getVirtualSiteFrameType(int index) const;

  int getVirtualSiteFrameAtom(int index, int nfrm) const;

  template <typename T> T getVirtualSiteFrameDimension(int index, int ndim);

  int getChargeTypeCount() const;
  
  int getLJTypeCount() const;

  int getTotalExclusions() const;

  UnitCellType getUnitCellType() const;

  ImplicitSolventModel getImplicitSolventModel() const;

  double getDielectricConstant() const;

  double getSaltConcentration() const;

  double getCoulombConstant() const;

  double getElectrostatic14Screening() const;

  double getVanDerWaals14Screening() const;

  std::string getPBRadiiSet() const;

  bool usesShake() const;

  bool usesSettle() const;

  std::string getWaterResidueName() const;
  
  int getRigidWaterCount() const;

  int getRigidWaterAtomCount() const;
  
  int getSettleGroupCount() const;

  int getBondConstraintCount() const;

  int getConstraintGroupCount() const;

  std::vector<int> getConstraintGroupAtoms(int index) const;
  
  int getConstraintGroupTotalSize() const;
  
  int getNonrigidParticleCount() const;
  
  std::vector<int> getChargeIndex() const;
  std::vector<int> getChargeIndex(int low_index, int high_index) const;
  int getChargeIndex(int index) const;

  std::vector<int> getLennardJonesIndex() const;
  std::vector<int> getLennardJonesIndex(int low_index, int high_index) const;
  int getLennardJonesIndex(int index) const;

  std::vector<int> getAtomExclusions(int index) const;

  std::vector<int> getNonbonded11Exclusions(int index) const;

  std::vector<int> getNonbonded12Exclusions(int index) const;

  std::vector<int> getNonbonded13Exclusions(int index) const;

  std::vector<int> getNonbonded14Exclusions(int index) const;

  template <typename T> T getChargeParameter(int index) const;

  template <typename T> T getLennardJonesSigma(int index_a, int index_b) const;
  template <typename T> T getLennardJonesSigma(int index_a) const;
  template <typename T> std::vector<T> getLennardJonesSigma() const;

  template <typename T> T getLennardJonesEpsilon(int index_a, int index_b) const;
  template <typename T> T getLennardJonesEpsilon(int index_a) const;
  template <typename T> std::vector<T> getLennardJonesEpsilon() const;

  template <typename T> std::vector<T> getAtomPBRadius() const;
  template <typename T> std::vector<T> getAtomPBRadius(int low_index, int high_index) const;
  template <typename T> T getAtomPBRadius(int index) const;

  template <typename T> std::vector<T> getGBScreeningFactor() const;
  template <typename T> std::vector<T> getGBScreeningFactor(int low_index, int high_index) const;
  template <typename T> T getGBScreeningFactor(int index) const;

  std::vector<char4> matchOverflowAtomName(const std::string &query) const;

  std::vector<char4> matchOverflowAtomType(const std::string &query) const;

  std::vector<char4> matchOverflowResidueName(const std::string &query) const;

  ValenceKit<double>
  getDoublePrecisionValenceKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;
  
  ValenceKit<float>
  getSinglePrecisionValenceKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  NonbondedKit<double>
  getDoublePrecisionNonbondedKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  NonbondedKit<float>
  getSinglePrecisionNonbondedKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  ImplicitSolventKit<double>
  getDoublePrecisionImplicitSolventKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  ImplicitSolventKit<float>
  getSinglePrecisionImplicitSolventKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  ChemicalDetailsKit
  getChemicalDetailsKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  VirtualSiteKit<double>
  getDoublePrecisionVirtualSiteKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  VirtualSiteKit<float>
  getSinglePrecisionVirtualSiteKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  ConstraintKit<double>
  getDoublePrecisionConstraintKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;

  ConstraintKit<float>
  getSinglePrecisionConstraintKit(HybridTargetLevel tier = HybridTargetLevel::HOST) const;
 
#ifdef STORMM_USE_HPC
#  ifdef STORMM_USE_CUDA
  ValenceKit<double> getDeviceViewToHostDPValenceKit() const;

  ValenceKit<float> getDeviceViewToHostSPValenceKit() const;

  NonbondedKit<double> getDeviceViewToHostDPNonbondedKit() const;

  NonbondedKit<float> getDeviceViewToHostSPNonbondedKit() const;

  ImplicitSolventKit<double> getDeviceViewToHostDPImplicitSolventKit() const;

  ImplicitSolventKit<float> getDeviceViewToHostSPImplicitSolventKit() const;

  ChemicalDetailsKit getDeviceViewToHostChemicalDetailsKit() const;

  VirtualSiteKit<double> getDeviceViewToHostDPVirtualSiteKit() const;

  VirtualSiteKit<float> getDeviceViewToHostSPVirtualSiteKit() const;

  ConstraintKit<double> getDeviceViewToHostDPConstraintKit() const;

  ConstraintKit<float> getDeviceViewToHostSPConstraintKit() const;
#  endif
#endif

  const AtomGraph* getSelfPointer() const;
  AtomGraph* getSelfPointer();
  
#ifdef STORMM_USE_HPC
  void upload();

  void download();
#endif

  void setSource(const std::string &new_source);

  void setBondParameters(double new_keq, double new_leq, bool set_keq, bool set_leq, int parm_idx,
                         ExceptionResponse policy = ExceptionResponse::SILENT);
 
  void setBondParameters(double new_keq, double new_leq, bool set_keq, bool set_leq, char4 type_a,
                         char4 type_b, ExceptionResponse policy = ExceptionResponse::SILENT);

  void setAngleParameters(double new_keq, double new_teq, bool set_keq, bool set_teq, int parm_idx,
                          ExceptionResponse policy = ExceptionResponse::SILENT);
 
  void setAngleParameters(double new_keq, double new_teq, bool set_keq, bool set_teq, char4 type_a,
                          char4 type_b, char4 type_c,
                          ExceptionResponse policy = ExceptionResponse::SILENT);

  void setDihedralParameters(double new_amplitude, double new_phase_angle, bool set_amplitude,
                             bool set_phase_angle, int parm_idx,
                             ExceptionResponse policy = ExceptionResponse::SILENT);
 
  void setDihedralParameters(double new_amplitude, double new_phase_angle, bool set_amplitude,
                             bool set_phase_angle, char4 type_a, char4 type_b, char4 type_c,
                             char4 type_d, double periodicity,
                             ExceptionResponse policy = ExceptionResponse::SILENT);

  void setUreyBradleyParameters(double new_keq, double new_leq, bool set_keq, bool set_leq,
                                int parm_idx,
                                ExceptionResponse policy = ExceptionResponse::SILENT);
 
  void setUreyBradleyParameters(double new_keq, double new_leq, bool set_keq, bool set_leq,
                                char4 type_a, char4 type_b, char4 type_c,
                                ExceptionResponse policy = ExceptionResponse::SILENT);

  void setCharmmImprParameters(double new_stiffness, double new_phase_angle, bool set_stiffness,
                               bool set_phase_angle, int parm_idx,
                               ExceptionResponse policy = ExceptionResponse::SILENT);
 
  void setCharmmImprParameters(double new_stiffness, double new_phase_angle, bool set_stiffness,
                               bool set_phase_angle, char4 type_a, char4 type_b, char4 type_c,
                               char4 type_d, ExceptionResponse policy = ExceptionResponse::SILENT);
  
  void setImplicitSolventModel(ImplicitSolventModel igb_in, double dielectric_in = 80.0,
                               double saltcon_in = 0.0,
                               AtomicRadiusSet radii_set = AtomicRadiusSet::NONE,
                               ExceptionResponse policy = ExceptionResponse::WARN);

  void setWaterResidueName(const char4 new_name);
  void setWaterResidueName(const std::string &new_name);

  void printToFile(const std::string &output_file = {},
                   TopologyKind output_style = TopologyKind::AMBER,
                   PrintSituation expectation = PrintSituation::OPEN_NEW,
                   ExceptionResponse pr_policy = ExceptionResponse::DIE) const;

  void printAmberFormat(const std::string &output_file = {},
                        PrintSituation expectation = PrintSituation::OPEN_NEW,
                        ExceptionResponse pr_policy = ExceptionResponse::DIE) const;
  
private:
  
  // Title, version, and date stamps found in this topology
  char version_stamp[16];  
  tm date;                 
  std::string title;       
  std::string source;      

  std::vector<Citation> force_fields;
 
  // Counts of atoms, residues, and other parts of the system
  int atom_count;                 
  int residue_count;              
  int molecule_count;             
  int largest_residue_size;       
  int last_solute_residue;        
  int last_solute_atom;           
  int first_solvent_molecule;     
  int last_atom_before_cap;       
  int implicit_copy_count;        
  int largest_molecule_size;      
  int water_residue_size;         
  int water_residue_count;        
  int unconstrained_dof;          
  int constrained_dof;            
  Hybrid<int> descriptors;        
  Hybrid<int> residue_limits;     
  Hybrid<int> atom_struc_numbers; 
  Hybrid<int> residue_numbers;    
  Hybrid<int> molecule_limits;    
 
  // Atom and residue details
  Hybrid<int> atomic_numbers;             
  Hybrid<int> mobile_atoms;               
  Hybrid<int> molecule_membership;        
  Hybrid<int> molecule_contents;          
  Hybrid<double> atomic_charges;          
  Hybrid<double> atomic_masses;           
  Hybrid<double> inverse_atomic_masses;   
  Hybrid<float> sp_atomic_charges;        
  Hybrid<float> sp_atomic_masses;         
  Hybrid<float> sp_inverse_atomic_masses; 
  Hybrid<char4> atom_names;               
  Hybrid<char4> atom_types;               
  Hybrid<char4> residue_names;            
 
  // Relevant information for CHARMM valence terms (the CHARMM force field has many terms that
  // function in the same ways as the parameters that make up all of Amber, but since Amber works
  // exclusively in these terms (harmonic angles in the theta bending dimension, cosine-based
  // dihedral terms, cosine-based impropers using the same formula) they are referred to as
  // "AMBER" force field family parameters (see ForceFieldFamily enum class above).  The "CHARMM"
  // family of parameters include Urey-Bradley angle terms, CHARMM impropers (yet another take on
  // harmonic potentials), and CHARMM CMAP terms.  The Amber ff19SB force field uses CHARMM CMAPs,
  // but such terms always belong to the "CHARMM" family of parameters.
  int urey_bradley_term_count;                
  int charmm_impr_term_count;                 
  int cmap_term_count;                        
  int urey_bradley_parameter_count;           
  int charmm_impr_parameter_count;            
  int cmap_surface_count;                     
  int urey_bradley_pert_term_count;           
  int charmm_impr_pert_term_count;            
  int cmap_pert_term_count;                   
  int urey_bradleys_in_perturbed_group;       
  int charmm_imprs_in_perturbed_group;        
  int cmaps_in_perturbed_group;               
  Hybrid<int> urey_bradley_i_atoms;           
  Hybrid<int> urey_bradley_k_atoms;           
  Hybrid<int> urey_bradley_parameter_indices; 
  Hybrid<int> urey_bradley_assigned_atoms;    
  Hybrid<int> urey_bradley_assigned_index;    
  Hybrid<int> urey_bradley_assigned_terms;    
  Hybrid<int> urey_bradley_assigned_bounds;   
  Hybrid<int> charmm_impr_i_atoms;            
  Hybrid<int> charmm_impr_j_atoms;            
  Hybrid<int> charmm_impr_k_atoms;            
  Hybrid<int> charmm_impr_l_atoms;            
  Hybrid<int> charmm_impr_parameter_indices;  
  Hybrid<int> charmm_impr_assigned_atoms;     
  Hybrid<int> charmm_impr_assigned_index;     
  Hybrid<int> charmm_impr_assigned_terms;     
  Hybrid<int> charmm_impr_assigned_bounds;    
  Hybrid<int> cmap_i_atoms;                   
  Hybrid<int> cmap_j_atoms;                   
  Hybrid<int> cmap_k_atoms;                   
  Hybrid<int> cmap_l_atoms;                   
  Hybrid<int> cmap_m_atoms;                   
  Hybrid<int> cmap_surface_dimensions;        
  Hybrid<int> cmap_surface_bounds;            
  Hybrid<int> cmap_patch_bounds;              
  Hybrid<int> cmap_surface_indices;           
  Hybrid<int> cmap_assigned_atoms;            
  Hybrid<int> cmap_assigned_index;            
  Hybrid<int> cmap_assigned_terms;            
  Hybrid<int> cmap_assigned_bounds;           
  Hybrid<double> urey_bradley_stiffnesses;    
  Hybrid<double> urey_bradley_equilibria;     
  Hybrid<double> charmm_impr_stiffnesses;     
  Hybrid<double> charmm_impr_phase_angles;    
  Hybrid<double> cmap_surfaces;               
  Hybrid<double> cmap_phi_derivatives;        
  Hybrid<double> cmap_psi_derivatives;        
  Hybrid<double> cmap_phi_psi_derivatives;    

  Hybrid<double> cmap_patches;
  Hybrid<float> sp_urey_bradley_stiffnesses;  
  Hybrid<float> sp_urey_bradley_equilibria;   
  Hybrid<float> sp_charmm_impr_stiffnesses;   
  Hybrid<float> sp_charmm_impr_phase_angles;  
  Hybrid<float> sp_cmap_surfaces;             
  Hybrid<float> sp_cmap_phi_derivatives;      
  Hybrid<float> sp_cmap_psi_derivatives;      
  Hybrid<float> sp_cmap_phi_psi_derivatives;  
  Hybrid<float> sp_cmap_patches;              
 
  // Relevant information for the standard ("Amber") bonded calculation
  int bond_term_with_hydrogen;           
  int angl_term_with_hydrogen;           
  int dihe_term_with_hydrogen;           
  int bond_term_without_hydrogen;        
  int angl_term_without_hydrogen;        
  int dihe_term_without_hydrogen;        
  int bond_term_count;                   
  int angl_term_count;                   
  int dihe_term_count;                   
  int bond_parameter_count;              
  int angl_parameter_count;              
  int dihe_parameter_count;              
  int bond_perturbation_term_count;      
  int angl_perturbation_term_count;      
  int dihe_perturbation_term_count;      
  int bonds_in_perturbed_group;          
  int angls_in_perturbed_group;          
  int dihes_in_perturbed_group;          
  int bonded_group_count;                
  Hybrid<double> bond_stiffnesses;       
  Hybrid<double> bond_equilibria;        
  Hybrid<double> angl_stiffnesses;       
  Hybrid<double> angl_equilibria;        
  Hybrid<double> dihe_amplitudes;        
  Hybrid<double> dihe_periodicities;     
  Hybrid<double> dihe_phase_angles;      
  Hybrid<float> sp_bond_stiffnesses;     
  Hybrid<float> sp_bond_equilibria;      
  Hybrid<float> sp_angl_stiffnesses;     
  Hybrid<float> sp_angl_equilibria;      
  Hybrid<float> sp_dihe_amplitudes;      
  Hybrid<float> sp_dihe_periodicities;   
  Hybrid<float> sp_dihe_phase_angles;    
  Hybrid<int> bond_i_atoms;              
  Hybrid<int> bond_j_atoms;              
  Hybrid<int> bond_parameter_indices;    
  Hybrid<int> bond_assigned_atoms;       
  Hybrid<int> bond_assigned_index;       
  Hybrid<int> bond_assigned_terms;       
  Hybrid<int> bond_assigned_bounds;      
  Hybrid<int> angl_i_atoms;              
  Hybrid<int> angl_j_atoms;              
  Hybrid<int> angl_k_atoms;              
  Hybrid<int> angl_parameter_indices;    
  Hybrid<int> angl_assigned_atoms;       
  Hybrid<int> angl_assigned_index;       
  Hybrid<int> angl_assigned_terms;       
  Hybrid<int> angl_assigned_bounds;      
  Hybrid<int> dihe_i_atoms;              
  Hybrid<int> dihe_j_atoms;              
  Hybrid<int> dihe_k_atoms;              
  Hybrid<int> dihe_l_atoms;              
  Hybrid<int> dihe_parameter_indices;    
  Hybrid<int> dihe14_parameter_indices;  
  Hybrid<int> dihe_assigned_atoms;       
  Hybrid<int> dihe_assigned_index;       
  Hybrid<int> dihe_assigned_terms;       
  Hybrid<int> dihe_assigned_bounds;      
  Hybrid<int> bond_evaluation_mask;      
  Hybrid<int> angl_evaluation_mask;      
  Hybrid<char4> bond_modifiers;          
  Hybrid<char4> angl_modifiers;          
  Hybrid<char4> dihe_modifiers;          
  Hybrid<char4> bond_assigned_mods;      
  Hybrid<char4> angl_assigned_mods;      
  Hybrid<char4> dihe_assigned_mods;      
 
  // Information relevant to virtual site placement
  int virtual_site_count;                      
  int virtual_site_parameter_set_count;        
  Hybrid<int> virtual_site_atoms;              
  Hybrid<int> virtual_site_frame_types;        
  Hybrid<int> virtual_site_frame1_atoms;       
  Hybrid<int> virtual_site_frame2_atoms;       
  Hybrid<int> virtual_site_frame3_atoms;       
  Hybrid<int> virtual_site_frame4_atoms;       
  Hybrid<int> virtual_site_parameter_indices;  
  Hybrid<double> virtual_site_frame_dim1;      
  Hybrid<double> virtual_site_frame_dim2;      
  Hybrid<double> virtual_site_frame_dim3;      
  Hybrid<float> sp_virtual_site_frame_dim1;    
  Hybrid<float> sp_virtual_site_frame_dim2;    
  Hybrid<float> sp_virtual_site_frame_dim3;    
 
  // Relevant information for the non-bonded calculation
  int charge_type_count;                 
  int lj_type_count;                     
  bool has_14_lennard_jones_data;        
  int total_exclusions;                  
  int attenuated_14_type_count;          
  int inferred_14_attenuations;          
  UnitCellType periodic_box_class;       
  ImplicitSolventModel gb_style;         
  double dielectric_constant;            
  double salt_concentration;             
  double coulomb_constant;               
  double elec14_screening_factor;        
  double vdw14_screening_factor;         
  std::string pb_radii_set;              
  Hybrid<int> charge_indices;            
  Hybrid<int> lennard_jones_indices;     
  Hybrid<int> atom_exclusion_bounds;     
  Hybrid<int> atom_exclusion_list;       
  Hybrid<int> nb11_exclusion_bounds;     
  Hybrid<int> nb11_exclusion_list;       
  Hybrid<int> nb12_exclusion_bounds;     
  Hybrid<int> nb12_exclusion_list;       
  Hybrid<int> nb13_exclusion_bounds;     
  Hybrid<int> nb13_exclusion_list;       
  Hybrid<int> nb14_exclusion_bounds;     
  Hybrid<int> nb14_exclusion_list;       
  Hybrid<int> infr14_i_atoms;            
  Hybrid<int> infr14_l_atoms;            
  Hybrid<int> infr14_parameter_indices;  
  Hybrid<int> neck_gb_indices;           
  Hybrid<double> charge_parameters;      
  Hybrid<double> lj_a_values;            
  Hybrid<double> lj_b_values;            
  Hybrid<double> lj_c_values;            
  Hybrid<double> lj_14_a_values;         
  Hybrid<double> lj_14_b_values;         
  Hybrid<double> lj_14_c_values;         
  Hybrid<double> lj_sigma_values;        
  Hybrid<double> lj_14_sigma_values;     
  Hybrid<double> lj_type_corrections;    
  Hybrid<double> attn14_elec_factors;    
  Hybrid<double> attn14_vdw_factors;     
  Hybrid<double> atomic_pb_radii;        
  Hybrid<double> gb_screening_factors;   
  Hybrid<double> gb_alpha_parameters;    
  Hybrid<double> gb_beta_parameters;     
  Hybrid<double> gb_gamma_parameters;    
  Hybrid<float> sp_charge_parameters;    
  Hybrid<float> sp_lj_a_values;          
  Hybrid<float> sp_lj_b_values;          
  Hybrid<float> sp_lj_c_values;          
  Hybrid<float> sp_lj_14_a_values;       
  Hybrid<float> sp_lj_14_b_values;       
  Hybrid<float> sp_lj_14_c_values;       
  Hybrid<float> sp_lj_sigma_values;      
  Hybrid<float> sp_lj_14_sigma_values;   
  Hybrid<float> sp_lj_type_corrections;  
  Hybrid<float> sp_attn14_elec_factors;  
  Hybrid<float> sp_attn14_vdw_factors;   
  Hybrid<float> sp_atomic_pb_radii;      
  Hybrid<float> sp_gb_screening_factors; 
  Hybrid<float> sp_gb_alpha_parameters;  
  Hybrid<float> sp_gb_beta_parameters;   
  Hybrid<float> sp_gb_gamma_parameters;  
 
  // MD propagation algorithm directives
  ApplyConstraints use_shake;                  
  ApplyConstraints use_settle;                 
  PerturbationSetting use_perturbation_info;   
  SolventCapSetting use_solvent_cap_option;    
  PolarizationSetting use_polarization;        
  char4 water_residue_name;                    
  std::string bond_constraint_mask;            
  std::string bond_constraint_omit_mask;       
  int bond_constraint_count;                   
  int nonrigid_particle_count;                 
  int settle_group_count;                      
  int settle_parameter_count;                  
  int constraint_group_count;                  
  int constraint_parameter_count;              
  Hybrid<int> settle_oxygen_atoms;             
  Hybrid<int> settle_hydro1_atoms;             
  Hybrid<int> settle_hydro2_atoms;             
  Hybrid<int> settle_parameter_indices;        
  Hybrid<int> constraint_group_atoms;          
  Hybrid<int> constraint_group_bounds;         
  Hybrid<int> constraint_parameter_indices;    
  Hybrid<int> constraint_parameter_bounds;     
  Hybrid<double> settle_mo;                    
  Hybrid<double> settle_mh;                    
  Hybrid<double> settle_moh;                   
  Hybrid<double> settle_mormt;                 
  Hybrid<double> settle_mhrmt;                 
  Hybrid<double> settle_ra;                    
  Hybrid<double> settle_rb;                    
  Hybrid<double> settle_rc;                    
  Hybrid<double> constraint_inverse_masses;    
  Hybrid<double> constraint_squared_lengths;   
  Hybrid<float> sp_settle_mo;                  
  Hybrid<float> sp_settle_mh;                  
  Hybrid<float> sp_settle_moh;                 
  Hybrid<float> sp_settle_mormt;               
  Hybrid<float> sp_settle_mhrmt;               
  Hybrid<float> sp_settle_ra;                  
  Hybrid<float> sp_settle_rb;                  
  Hybrid<float> sp_settle_rc;                  
  Hybrid<float> sp_constraint_inverse_masses;  
  Hybrid<float> sp_constraint_squared_lengths; 
 
  // Atom, atom type, and residue name overflow keys
  Hybrid<char4> atom_overflow_names;    
  Hybrid<char4> atom_overflow_types;    
  Hybrid<char4> residue_overflow_names; 
 
  // Information currently unused
  int unused_nhparm;                       
  int unused_nparm;                        
  int unused_natyp;                        
  int hbond_10_12_parameter_count;         
  int heavy_bonds_plus_constraints;        
  int heavy_angls_plus_constraints;        
  int heavy_dihes_plus_constraints;        
  Hybrid<int> tree_joining_info;           
  Hybrid<int> last_rotator_info;           
  Hybrid<double> solty_info;               
  Hybrid<double> hbond_a_values;           
  Hybrid<double> hbond_b_values;           
  Hybrid<double> hbond_cutoffs;            
  Hybrid<double> atomic_polarizabilities;  
  Hybrid<char4> tree_symbols;              

  Hybrid<int> int_data;
  Hybrid<double> double_data;
  Hybrid<float> float_data;
  Hybrid<char4> char4_data;
  
  void loadHybridArrays(const std::vector<int> &tmp_desc,
                        const std::vector<int> &tmp_residue_limits,
                        const std::vector<int> &tmp_atom_struc_numbers,
                        const std::vector<int> &tmp_residue_numbers,
                        const std::vector<int> &tmp_molecule_limits,
                        const std::vector<int> &tmp_atomic_numbers,
                        const std::vector<int> &tmp_molecule_membership,
                        const std::vector<int> &tmp_mobile_atoms,
                        const std::vector<int> &tmp_molecule_contents,
                        const std::vector<int> &tmp_cmap_surf_dims,
                        const std::vector<int> &tmp_cmap_surf_bounds,
                        const std::vector<int> &tmp_charge_indices,
                        const std::vector<int> &tmp_lennard_jones_indices,
                        const std::vector<int> &tmp_inferred_14_i_atoms,
                        const std::vector<int> &tmp_inferred_14_j_atoms,
                        const std::vector<int> &tmp_inferred_14_param_idx,
                        const std::vector<int> &tmp_neck_gb_indices,
                        const std::vector<int> &tmp_tree_joining_info,
                        const std::vector<int> &tmp_last_rotator_info,
                        const std::vector<double> &tmp_charges,
                        const std::vector<double> &tmp_masses,
                        const std::vector<double> &tmp_ub_stiffnesses,
                        const std::vector<double> &tmp_ub_equilibria,
                        const std::vector<double> &tmp_charmm_impr_stiffnesses,
                        const std::vector<double> &tmp_charmm_impr_phase_angles,
                        const std::vector<double> &tmp_cmap_surfaces,
                        const std::vector<double> &tmp_bond_stiffnesses,
                        const std::vector<double> &tmp_bond_equilibria,
                        const std::vector<double> &tmp_angl_stiffnesses,
                        const std::vector<double> &tmp_angl_equilibria,
                        const std::vector<double> &tmp_dihe_amplitudes,
                        const std::vector<double> &tmp_dihe_periodicities,
                        const std::vector<double> &tmp_dihe_phase_angles,
                        const std::vector<double> &tmp_charge_parameters,
                        const std::vector<double> &tmp_lj_a_values,
                        const std::vector<double> &tmp_lj_b_values,
                        const std::vector<double> &tmp_lj_c_values,
                        const std::vector<double> &tmp_lj_14_a_values,
                        const std::vector<double> &tmp_lj_14_b_values,
                        const std::vector<double> &tmp_lj_14_c_values,
                        const std::vector<double> &tmp_atomic_pb_radii,
                        const std::vector<double> &tmp_gb_screening_factors,
                        const std::vector<double> &tmp_gb_coef,
                        const std::vector<double> &tmp_solty_info,
                        const std::vector<double> &tmp_hbond_a_values,
                        const std::vector<double> &tmp_hbond_b_values,
                        const std::vector<double> &tmp_hbond_cutoffs,
                        const std::vector<double> &tmp_atpol,
                        const std::vector<char4> &tmp_atom_names,
                        const std::vector<char4> &tmp_atom_types,
                        const std::vector<char4> &tmp_residue_names,
                        const std::vector<char4> &tmp_tree_symbols,
                        const CmapAccessories &cmap_table,
                        const CondensedExclusions &cond_excl,
                        const BasicValenceTable &basic_vtable,
                        const CharmmValenceTable &charmm_vtable,
                        const AttenuationParameterSet &attn_parm,
                        const VirtualSiteTable &vsite_table,
                        const Map1234 &all_nb_excl, const ConstraintTable &cnst_table);
  
  void rebasePointers();

  void characterizeWaterResidue();

  void printAmberSection(std::ofstream *foutp, const std::string &section_name,
                         const std::vector<PolyNumeric> &values, int values_per_line,
                         int value_width,  int decimal_places, NumberFormat fmt,
                         const std::string &task) const;

  friend class AtomGraphStage;
};

std::vector<int2> findMoleculeOrder(const std::vector<AtomGraph*> &agv,
                                    const std::vector<int> &counts,
                                    std::vector<int2> *origins = nullptr,
                                    MoleculeOrdering arrangement = MoleculeOrdering::REORDER_ALL);

bool scanForExclusion(const std::vector<int> &list, int llim, int hlim, int partner);
  
void cullPriorExclusions(std::vector<int3> *proposals, const std::vector<int3> &priors);
 
} // namespace topology
} // namespace stormm
 
#include "atomgraph.tpp"
 
#endif