CIF fixes, preserve CHARMM PSF segment ID (#1086)

* Start reading the struct_conn block

* Since in CIF files the chain ID can be larger than 1 character, start
preparing to store multicharacter chain IDs

* Fixes for new string chainID

* Fix uninitialized warnings

* More changes for chainID -> string

* Finish converting chainid to string

* Get rid of the segmentID constructor for Residue, which could
ambiguously stand in for the old chainID as char constructor. Might be
able to eventually deprecate the segid in residue.

* Store serial blocks in the same way as loop blocks for consistency

* Fix handling of chainid keyword for PDB output

* If a blank chain ID is set, ensure the chain id string is cleared so that HasChainID functions properly

* Chain variable in pres is a string now, treat accordingly

* Do not attempt to set solvent info if we dont have molecules or didnt want to look for molecules

* Print out more structure info

* Comment out unused function

* Add function to get unit cell info from cell data block

* Change return so that no box can be indicated

* Use new box function

* Test function to add bonds. Not sure its working properly yet though.

* Hide the struct_conn stuff for now, not ready for prime time. Fix class
variable name formatting

* Use chain ID for segment ID if it is set

* Original segid is now preserved

* Segid no longer used in residue

* 6.27.0. Minor version bump for handling multichar chain IDs, detection
of bad CIF boxes, preserving charmm psf segment ID

* Try to fix appveyor configure-mingw build by going back to VS2019 image.

---------

Co-authored-by: Daniel R. Roe <daniel.roe@nih.gov>

appveyor.yml

@@ -1,7 +1,7 @@
 platform:
   - x64
 
-image: Visual Studio 2015
+image: Visual Studio 2019
 
 clone_folder: c:\projects\cpptraj
 

src/Action_AtomicFluct.cpp

@@ -232,7 +232,7 @@ void Action_AtomicFluct::Print() {
           PDBfile& adpout = static_cast<PDBfile&>( *adpoutfile_ );
           adpout.WriteANISOU(
             atom+1, (*fluctParm_)[atom].c_str(), fluctParm_->Res(resnum).c_str(),
-            fluctParm_->Res(resnum).ChainId(), fluctParm_->Res(resnum).OriginalResNum(),
+            fluctParm_->Res(resnum).ChainID_1char(), fluctParm_->Res(resnum).OriginalResNum(),
             anisou, (*fluctParm_)[atom].ElementName(), 0 );
         }
       }

src/Action_FixAtomOrder.cpp

@@ -116,7 +116,7 @@ Action::RetType Action_FixAtomOrder::PdbOrder(ActionSetup& setup) {
     else
       pt = AtomTopType::ATOM;
     atoms.push_back( AtomTopType(pt, idx, res.OriginalResNum(),
-                                 res.Icode(), res.ChainId()) );
+                                 res.Icode(), res.ChainID()) );
   }
   // Sort by PDB info
   std::sort(atoms.begin(), atoms.end());

src/Analysis_Hist.cpp

@@ -471,7 +471,7 @@ Analysis::RetType Analysis_Hist::Analyze() {
       // Create pseudo-topology/trajectory
       if (!traj3dName_.empty()) {
         Topology pseudo;
-        pseudo.AddTopAtom(Atom("H3D", 0), Residue("H3D", 1, ' ', ' '));
+        pseudo.AddTopAtom(Atom("H3D", 0), Residue("H3D", 1, ' ', ""));
         pseudo.CommonSetup();
         if (!parmoutName_.empty()) {
           ParmFile pfile;

src/AtomTopType.h

@@ -1,11 +1,12 @@
 #ifndef INC_ATOMTOPTYPE_H
 #define INC_ATOMTOPTYPE_H
+#include <string>
 /// Hold atom topology information which can be used when re-ordering
 class AtomTopType {
   public:
     enum PdbType { ATOM = 0, HETATM };
     /// Type, index, resnum, icode, chain
-    AtomTopType(PdbType t, int i, int r, char e, char c) :
+    AtomTopType(PdbType t, int i, int r, char e, std::string const& c) :
       type_(t), index_(i), resnum_(r), icode_(e), chainid_(c) {}
 
     /// \return true if less than (Chain, resnum, icode, index)
@@ -53,6 +54,6 @@ class AtomTopType {
     int index_;    ///< The atom index
     int resnum_;   ///< The residue number
     char icode_;   ///< The residue insertion code.
-    char chainid_; ///< The chain
+    std::string chainid_; ///< The chain
 };
 #endif

src/CIFfile.cpp

@@ -2,6 +2,7 @@
 #include "ArgList.h"
 #include "CIFfile.h"
 #include "CpptrajStdio.h"
+#include "StringRoutines.h" // convertToDouble
 
 static inline int LineError(const char* msg, int num, const char* ptr) {
   mprinterr("Error: CIF line %i: %s\n", num, msg);
@@ -65,7 +66,7 @@ int CIFfile::DataBlock::GetColumnData(int NexpectedCols, BufferedLine& infile, b
 {
   const char* SEP = " \t";
   // Allocate for a line of data
-  columnData_.push_back( Sarray() );
+  //columnData_.push_back( Sarray() );
   // Tokenize the initial line
   int nReadCols = 0;
   int Ncols = infile.TokenizeLine(SEP);
@@ -149,6 +150,12 @@ int CIFfile::DataBlock::AddSerialDataRecord( const char* ptr, BufferedLine& infi
   return 0;
 }
 
+/** Start a serial data block. */
+void CIFfile::DataBlock::StartSerialDataBlock() {
+  // Allocate for a line of data
+  columnData_.push_back( Sarray() );
+}
+
 /** Add column label from loop section. */
 int CIFfile::DataBlock::AddLoopColumn( const char* ptr, BufferedLine& infile ) {
   if (ptr == 0) return 1;
@@ -170,6 +177,8 @@ int CIFfile::DataBlock::AddLoopColumn( const char* ptr, BufferedLine& infile ) {
 
 /** Add loop data. */
 int CIFfile::DataBlock::AddLoopData( const char* ptr, BufferedLine& infile ) {
+  // Allocate for a line of data
+  columnData_.push_back( Sarray() );
   // Should be as much data as there are column headers
   if (GetColumnData( columnHeaders_.size(), infile, false )) return 1;
   return 0;
@@ -224,7 +233,8 @@ std::string CIFfile::DataBlock::Data(std::string const& idIn) const {
   if (columnHeaders_.empty() || columnData_.empty()) return std::string("");
   int colnum = ColumnIndex( idIn );
   if (colnum == -1) return std::string("");
-  return columnData_[colnum].front();
+  //return columnData_[colnum].front();
+  return columnData_[0][colnum];
 }
 
 // -----------------------------------------------------------------------------
@@ -258,6 +268,15 @@ int CIFfile::CIFdata::AddDataBlock( DataBlock const& block ) {
   return 0;
 }
 
+/** Print CIF data blocks */
+void CIFfile::CIFdata::PrintDataBlocks() const {
+  for (CIF_DataType::const_iterator it = cifdata_.begin(); it != cifdata_.end(); ++it)
+  {
+    mprintf("\tData block: %s\n", it->first.c_str());
+    it->second.ListData();
+  }
+}
+
 // -----------------------------------------------------------------------------
 /// Used to return empty block for GetDataBlock
 const CIFfile::DataBlock CIFfile::emptyblock = DataBlock();
@@ -317,6 +336,7 @@ int CIFfile::Read(FileName const& fnameIn, int debugIn) {
     } else if ( currentMode == SERIAL ) {
       // SERIAL data block
       DataBlock serial;
+      serial.StartSerialDataBlock();
       while ( ptr != 0 && ptr[0] == '_' ) {
         serial.AddSerialDataRecord(ptr, file_);
         ptr = file_.Line();
@@ -353,7 +373,7 @@ int CIFfile::Read(FileName const& fnameIn, int debugIn) {
 }
 
 /** Vector with DataBlocks corresponding to given header and value. */
-CIFfile::DataBlock const& CIFfile::GetBlockWithColValue(
+/*CIFfile::DataBlock const& CIFfile::GetBlockWithColValue(
                                                    std::string const& header,
                                                    std::string const& col,
                                                    std::string const& value)
@@ -371,4 +391,51 @@ const
     }
   }
   return emptyblock;
+}*/
+
+/** List all data currently in the CIFfile. */
+void CIFfile::ListAllData() const {
+  for (std::vector<CIFdata>::const_iterator it = data_.begin();
+                                            it != data_.end(); ++it)
+  {
+    mprintf("CIF data: %s\n", it->DataName().c_str());
+    it->PrintDataBlocks();
+  }
+}
+
+/** Get box info from _cell block.
+  * \return 1 if box seems invalid, -1 if not box, 0 otherwise.
+  */
+int CIFfile::cif_Box_verbose(double* cif_box) const {
+  if (cif_box == 0) {
+    mprinterr("Internal Error: CIFfile::cif_Box_verbose: Null box passed in.\n");
+    return 1;
+  }
+  int box_stat = 0;
+  DataBlock const& cellblock = GetDataBlock("_cell");
+  if (cellblock.empty()) {
+    cif_box[0] = 0;
+    cif_box[1] = 0;
+    cif_box[2] = 0;
+    cif_box[3] = 0;
+    cif_box[4] = 0;
+    cif_box[5] = 0;
+    box_stat = -1;
+  } else {
+    cif_box[0] = convertToDouble( cellblock.Data("length_a") );
+    cif_box[1] = convertToDouble( cellblock.Data("length_b") );
+    cif_box[2] = convertToDouble( cellblock.Data("length_c") );
+    if (cif_box[0] == 1.0 && cif_box[1] == 1.0 && cif_box[2] == 1.0) {
+      mprintf("Warning: CIF cell lengths are all 1.0 Ang.;"
+              " this usually indicates an invalid box.\n");
+      box_stat = 1;
+    }
+    cif_box[3] = convertToDouble( cellblock.Data("angle_alpha") );
+    cif_box[4] = convertToDouble( cellblock.Data("angle_beta" ) );
+    cif_box[5] = convertToDouble( cellblock.Data("angle_gamma") );
+    mprintf("\tRead cell info from CIF: a=%g b=%g c=%g alpha=%g beta=%g gamma=%g\n",
+              cif_box[0], cif_box[1], cif_box[2], cif_box[3], cif_box[4], cif_box[5]);
+  }
+  return box_stat;
 }
+

src/CIFfile.h

@@ -17,6 +17,7 @@ class CIFfile {
         bool empty()                const { return dataHeader_.empty(); }
         int AddHeader(std::string const&);
         int AddSerialDataRecord(const char*, BufferedLine&);
+        void StartSerialDataBlock();
         int AddLoopColumn(const char*, BufferedLine&);
         int AddLoopData(const char*, BufferedLine&);
         void Append(DataBlock const&);
@@ -48,6 +49,10 @@ class CIFfile {
         int AddDataBlock(DataBlock const&);
         /// Get data block with specified header
         DataBlock const& GetDataBlock(std::string const&) const;
+        /// Print CIF data blocks
+        void PrintDataBlocks() const;
+        /// \return CIF data name
+        std::string const& DataName() const { return dataName_; }
       private:
         typedef std::map<std::string, DataBlock> CIF_DataType;
         /// Map block names to DataBlocks
@@ -64,9 +69,13 @@ class CIFfile {
     FileName const& CIFname() const { return file_.Filename(); }
     /// Get data from most recently added data set
     DataBlock const& GetDataBlock(std::string const&h) const { return data_.back().GetDataBlock(h); }
+    /// Set box info from _cell data block if present.
+    int cif_Box_verbose(double*) const;
+    //DataBlock const& GetBlockWithColValue(std::string const&, std::string const&,
+    //                                      std::string const&) const;
 
-    DataBlock const& GetBlockWithColValue(std::string const&, std::string const&,
-                                          std::string const&) const;
+    /// For debugging, list all data currently in the CIF file.
+    void ListAllData() const;
 
   private:
     //int AddDataBlock(DataBlock const&);

src/DataIO_AmberLib.cpp

@@ -147,7 +147,7 @@ int DataIO_AmberLib::read_atoms(Topology& topOut, std::string const& line, std::
   atm.DetermineElement( elt );
   atm.SetMassFromElement();
   atm.SetCharge( charge );
-  Residue res( unitName, resx, ' ', ' ' );
+  Residue res( unitName, resx, ' ', "" );
   topOut.AddTopAtom( atm, res );
   return 0;
 }

src/DataIO_AmberPrep.cpp

@@ -195,7 +195,7 @@ const
   Topology top;
   top.SetParmName( resName, infile.Filename() );
   // Residue
-  Residue res(resName, 1, ' ', ' ');
+  Residue res(resName, 1, ' ', "");
   // Loop over entries
   line = infile.Line();
   if (CheckLine(line)) return 1;

src/DataIO_Numpy.cpp

@@ -67,7 +67,7 @@ const
   // Create a pseudo topology
   Topology top;
   for (unsigned long iat = 0; iat != natoms; iat++)
-    top.AddTopAtom( Atom("CA","C"), Residue("XXX", 1, ' ', ' ') );
+    top.AddTopAtom( Atom("CA","C"), Residue("XXX", 1, ' ', "") );
   top.SetSingleMolecule();
   top.CommonSetup(false, false);
   top.Summary();

src/DataIO_VecTraj.cpp

@@ -67,7 +67,7 @@ int DataIO_VecTraj::WriteData(FileName const& fname, DataSetList const& SetList)
   pseudo.AddBondParm( BondParmType(0.0, 1.0) );
   int natom = 0;
   for (unsigned int nres = 1; nres <= VecSets.size(); nres++) {
-    Residue vec_res("VEC", nres, ' ', ' ');
+    Residue vec_res("VEC", nres, ' ', "");
     if (VecNeedsBond[nres-1]) {
       pseudo.AddTopAtom(Atom("OXYZ", 0), vec_res);
       pseudo.AddTopAtom(Atom("VXYZ", 0), vec_res);