mdlmol.h

// -*-c++-*-
#ifndef STORMM_MDL_FILE_H
#define STORMM_MDL_FILE_H
 
#include <fstream>
#include <string>
#include <vector>
#include "copyright.h"
#include "Accelerator/hybrid.h"
#include "Chemistry/znumber.h"
#include "Chemistry/chemical_features.h"
#include "Constants/behavior.h"
#include "DataTypes/common_types.h"
#include "DataTypes/stormm_vector_types.h"
#include "FileManagement/file_util.h"
#include "Parsing/ascii_numbers.h"
#include "Parsing/parse.h"
#include "Potential/energy_enumerators.h"
#include "Restraints/restraint_apparatus.h"
#include "Synthesis/condensate.h"
#include "Synthesis/phasespace_synthesis.h"
#include "Topology/atomgraph.h"
#include "Topology/atomgraph_abstracts.h"
#include "Trajectory/coordinateframe.h"
#include "Trajectory/coordinate_series.h"
#include "Trajectory/phasespace.h"
#include "Trajectory/trajectory_enumerators.h"
#include "molecule_file_io.h"
#include "molecule_format_enumerators.h"
#include "mdlmol_atomlist.h"
#include "mdlmol_bond.h"
#include "mdlmol_dataitem.h"
#include "mdlmol_property.h"
 
namespace stormm {
namespace structure {
 
using card::HybridTargetLevel;
using chemistry::symbolToZNumber;
using chemistry::ChemicalFeatures;
using constants::CartesianDimension;
using constants::CaseSensitivity;
using constants::ExceptionResponse;
using data_types::isFloatingPointScalarType;
using diskutil::PrintSituation;
using energy::StateVariable;
using parse::TextFile;
using restraints::RestraintApparatus;
using synthesis::Condensate;
using synthesis::CondensateReader;
using synthesis::PhaseSpaceSynthesis;
using synthesis::PsSynthesisReader;
using topology::AtomGraph;
using topology::ChemicalDetailsKit;
using topology::NonbondedKit;
using trajectory::CoordinateCycle;
using trajectory::CoordinateFrame;
using trajectory::CoordinateFrameReader;
using trajectory::CoordinateFrameWriter;
using trajectory::CoordinateSeries;
using trajectory::CoordinateSeriesReader;
using trajectory::CoordinateSeriesWriter;
using trajectory::PhaseSpace;
using trajectory::PhaseSpaceReader;
using trajectory::PhaseSpaceWriter;

constexpr char4 default_mdl_atomic_symbol = { ' ', ' ', ' ', ' ' };
constexpr int default_mdl_atomic_number = -1;
constexpr int default_mdl_formal_charge = 0;
constexpr RadicalState default_mdl_radical_state = RadicalState::NONE;
constexpr int default_mdl_isotopic_shift = 0;
constexpr MolObjAtomStereo default_mdl_stereo_parity = MolObjAtomStereo::NOT_STEREO;
constexpr HydrogenAssignment default_hydrogenation = HydrogenAssignment::DO_NOT_HYDROGENATE;
constexpr int default_mdl_implicit_hydrogen = 0;
constexpr int default_mdl_valence_connections = 0;
constexpr int default_mdl_map_count = 0;
constexpr bool default_mdl_stereo_considerations = false;
constexpr bool default_mdl_exact_change = false;
constexpr StereoRetention default_mdl_stereo_retention = StereoRetention::NOT_APPLIED;

class MolObjSTextGroup {
public:
private:
};
  
  
class MdlMol {
public:

  MdlMol(ExceptionResponse policy_in = ExceptionResponse::WARN);
 
  MdlMol(const std::string &filename, ExceptionResponse policy_in = ExceptionResponse::WARN,
         ModificationPolicy dimod_policy = ModificationPolicy::DO_NOT_MODIFY,
         ExceptionResponse dimod_notify = ExceptionResponse::WARN);
 
  MdlMol(const char* filename, ExceptionResponse policy_in = ExceptionResponse::WARN,
         ModificationPolicy dimod_policy = ModificationPolicy::DO_NOT_MODIFY,
         ExceptionResponse dimod_notify = ExceptionResponse::WARN);
 
  MdlMol(const TextFile &tf, int line_start = 0, int line_end = -1,
         CaseSensitivity capitalization = CaseSensitivity::YES,
         ExceptionResponse policy_in = ExceptionResponse::WARN,
         ModificationPolicy dimod_policy = ModificationPolicy::DO_NOT_MODIFY,
         ExceptionResponse dimod_notify = ExceptionResponse::WARN);
 
  template <typename T>
  MdlMol(const ChemicalFeatures *chemfe, const T* xcrd, const T* ycrd, const T* zcrd,
         double inv_scale, int molecule_index = 0);
 
  template <typename T>
  MdlMol(const ChemicalFeatures *chemfe, const T* xcrd, const T* ycrd, const T* zcrd,
         int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures *chemfe, const CoordinateFrameReader cfr, int molecule_index = 0);
  
  MdlMol(const ChemicalFeatures *chemfe, const CoordinateFrame *cf, int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures &chemfe, const CoordinateFrame &cf, int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures *chemfe, const PhaseSpaceReader psr, int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures *chemfe, const PhaseSpace *ps, int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures &chemfe, const PhaseSpace &ps, int molecule_index = 0);
 
  template <typename T>
  MdlMol(const ChemicalFeatures *chemfe, const CoordinateSeriesReader<T> csr, int frame_index,
         int molecule_index = 0);
 
  template <typename T>
  MdlMol(const ChemicalFeatures *chemfe, const CoordinateSeries<T> *cs, int frame_index,
         int molecule_index = 0);
 
  template <typename T>
  MdlMol(const ChemicalFeatures &chemfe, const CoordinateSeries<T> &cs, int frame_index,
         int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures *chemfe, const PsSynthesisReader poly_psr, int system_index,
         int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures *chemfe, const PhaseSpaceSynthesis *poly_ps, int system_index,
         int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures &chemfe, const PhaseSpaceSynthesis &poly_ps, int system_index,
         int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures *chemfe, const CondensateReader *cdnsr, int system_index,
         int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures *chemfe, const Condensate *cdns, int system_index,
         int molecule_index = 0);
 
  MdlMol(const ChemicalFeatures &chemfe, const Condensate &cdns, int system_index,
         int molecule_index = 0);

  MdlMol(const MdlMol &original) = default;
  MdlMol(MdlMol &&original) = default;
  MdlMol& operator=(const MdlMol &original) = default;
  MdlMol& operator=(MdlMol &&original) = default;

  const std::string& getTitle() const;
  
  int getAtomCount() const;
  
  int getBondCount() const;

  int getPropertiesCount() const;

  int getDataItemCount() const;

  double3 getCoordinates(int index) const;
  const std::vector<double3>& getCoordinates() const;
  std::vector<double> getCoordinates(CartesianDimension dim) const;

  const MdlMol* getSelfPointer() const;
  
  PhaseSpace exportPhaseSpace() const ;

  CoordinateFrame exportCoordinateFrame() const;

  char4 getAtomSymbol(int index) const;

  const std::vector<char4>& getAtomSymbols() const;
  
  int getAtomicNumber(int index) const;

  const std::vector<int>& getAtomicNumbers() const;

  int getFormalCharge(int index) const;

  const std::vector<int>& getFormalCharges() const;
  
  template <typename T> void impartCoordinates(const T* xcrd, const T* ycrd, const T* zcrd,
                                               double scale_factor);
 
  template <typename T> void impartCoordinates(const T* xcrd, const T* ycrd, const T* zcrd,
                                               double scale_factor, const ChemicalDetailsKit &cdk,
                                               int molecule_index);
  
  void impartCoordinates(const llint* xcrd, const llint* ycrd, const llint* zcrd,
                         const int* xcrd_ovrf, const int* ycrd_ovrf, const int* zcrd_ovrf,
                         double scale_factor);
 
  void impartCoordinates(const llint* xcrd, const llint* ycrd, const llint* zcrd,
                         const int* xcrd_ovrf, const int* ycrd_ovrf, const int* zcrd_ovrf,
                         double scale_factor, const ChemicalDetailsKit &cdk, int molecule_index);
  
  void impartCoordinates(const PhaseSpaceReader &psr);
 
  void impartCoordinates(const PhaseSpaceWriter &psw);
 
  void impartCoordinates(const PhaseSpace *ps, CoordinateCycle orientation,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  void impartCoordinates(const PhaseSpace *ps,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  void impartCoordinates(const PhaseSpace &ps, CoordinateCycle orientation,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  void impartCoordinates(const PhaseSpace &ps,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  void impartCoordinates(const CoordinateFrameReader &cfr);
 
  void impartCoordinates(const CoordinateFrameWriter &cfw);
 
  void impartCoordinates(const CoordinateFrame *cf,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  void impartCoordinates(const CoordinateFrame &cf,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  template <typename T>
  void impartCoordinates(const CoordinateSeriesReader<T> &cs, int frame_index,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  template <typename T>
  void impartCoordinates(const CoordinateSeriesWriter<T> &cs, int frame_index,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  template <typename T>
  void impartCoordinates(const CoordinateSeries<T> *cs, int frame_index,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  template <typename T>
  void impartCoordinates(const CoordinateSeries<T> &cs, int frame_index,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  void impartCoordinates(const PhaseSpaceSynthesis *poly_ps, int system_index,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);
 
  void impartCoordinates(const PhaseSpaceSynthesis &poly_ps, int system_index,
                         HybridTargetLevel tier = HybridTargetLevel::HOST);

  void addProperty(const MdlMolProperty *mmprop);
  void addProperty(const MdlMolProperty &mmprop);

  void addDataItem(const MdlMolDataRequest &ask, const AtomGraph &ag,
                   const RestraintApparatus &ra);

  MdlMolDataItemKind getDataItemKind(int item_index) const;

  StateVariable getTrackedState(int item_index) const;

  void addLineToDataItem(const std::string &text, int item_index);

  const std::string& getDataItemName(int item_index) const;

  const std::string& getDataItemOutputName(int item_index) const;

  int getDataItemIndex(const std::string &item_name) const;
  
  const std::vector<std::string>& getDataItemContent(int item_index) const;
  std::vector<std::string> getDataItemContent(const std::string &item_name) const;

  void writeMdl(std::ofstream *foutp, MdlMolVersion vformat) const;
 
  void writeMdl(const std::string &fname, MdlMolVersion vformat, PrintSituation expectation) const;
 
  std::string writeMdl(MdlMolVersion vformat = MdlMolVersion::V3000) const;

  void writeDataItems(std::ofstream *foutp, int mol_index = 0) const;
 
  void writeDataItems(const std::string &fname, PrintSituation expectation,
                      int mol_index = 0) const;
 
  std::string writeDataItems(int mol_index = 0) const;
  
private:
 
  // Items describing quantities of information (most of them from the counts line)
  ExceptionResponse policy;   
  MdlMolVersion version_no;   
  int atom_count;             
  int bond_count;             
  int list_count;             
  int stext_entry_count;      
  int properties_count;       
  int sgroup_count;           
  int constraint_count;       
  MolObjChirality chirality;  
  int registry_number;        
  int data_item_count;        
 
  // Atomic properties
  bool property_formal_charges;             
  bool property_radicals;                   
  bool property_isotopes;                   
  bool property_element_lists;              
  std::vector<double3> coordinates;         
  std::vector<char4> atomic_symbols;        
  std::vector<int> atomic_numbers;          
  std::vector<int> formal_charges;          
  std::vector<RadicalState> radicals;       
  std::vector<int> isotopic_shifts;         
  std::vector<MolObjAtomStereo> parities;   
  std::vector<int> implicit_hydrogens;      
  std::vector<bool> stereo_considerations;  
  std::vector<int> valence_connections;     
  std::vector<int> atom_atom_mapping_count; 
  std::vector<bool> exact_change_enforced;  

  std::vector<HydrogenAssignment> hydrogenation_protocol;

  std::vector<StereoRetention> orientation_stability;

  std::vector<MdlMolBond> bonds;

  std::vector<MdlMolAtomList> element_lists;

  std::vector<MolObjSTextGroup> stext_entries;
  
  std::vector<MdlMolProperty> properties;

  std::vector<MdlMolDataItem> data_items;
  
  std::string title;

  std::string software_details;

  std::string general_comment;

  std::string external_regno;

  void validateAtomIndex(int index, const char* caller) const;
  
  void allocate();
  void allocate(const ChemicalDetailsKit &cdk, const NonbondedKit<double> &nbk, int mol_idx = 0);
  
  int getIsotopicShiftCode(int atom_index) const;
  
  void interpretFormalCharge(int charge_in, int atom_index);

  int getFormalChargeCode(int atom_index) const;
  
  MolObjAtomStereo interpretStereoParity(int setting_in);

  void interpretImplicitHydrogenContent(int nh_in, int atom_index);

  int getImplicitHydrogenCode(const int atom_index) const;
  
  bool interpretBooleanValue(int value_in, const std::string &desc);

  int interpretValenceNumber(int count_in);

  StereoRetention interpretStereoStability(int code_in);

  void updateV2kAtomAttributes();
  
  void hydrogenate();

  void compareExternalRegistryNumbers(const std::string &regno_in);

  void checkAtomCount(int ext_atom_count) const;

  void checkDataItemIndex(int item_index, const char* caller) const;

  void addAtomIntProperty(const std::vector<int2> &notable_data, const char4 prcode);
  
  void transferTopologicalDetails(const ChemicalFeatures *chemfe, int molecule_index = 0);
};

std::vector<MdlMol> readStructureDataFile(const TextFile &tf, int low_frame_limit,
                                          int high_frame_limit,
                                          CaseSensitivity capitalization = CaseSensitivity::YES,
                                          ExceptionResponse policy = ExceptionResponse::WARN);
std::vector<MdlMol> readStructureDataFile(const std::string &file_name,
                                          CaseSensitivity capitalization = CaseSensitivity::YES,
                                          ExceptionResponse policy = ExceptionResponse::WARN);
 
std::vector<MdlMol> readStructureDataFile(const TextFile &tf,
                                          CaseSensitivity capitalization = CaseSensitivity::YES,
                                          ExceptionResponse policy = ExceptionResponse::WARN);
  
} // namespace structure
} // namespace stormm
 
#include "mdlmol.tpp"
 
#endif