Add ability to read gbnsr6 output. Add new output format keywords to 'avg' command.

doc/DocumentChecksums.txt

@@ -1,3 +1,3 @@
 f6f8cb1a79951d80a9d2656fd9c30f55  CpptrajDevelopmentGuide.lyx
-f7837faf109fc9c36db4e0274ff6d92f  cpptraj.lyx
+cbbe066a99553d484a9113510c740f1a  cpptraj.lyx
 5d9b5b5ed47a3ded57b6464df99b3585  CpptrajManual.lyx

doc/cpptraj.lyx

@@ -3141,7 +3141,7 @@ status open
 \begin_layout Plain Layout
 \align center
 \begin_inset Tabular
-<lyxtabular version="3" rows="27" columns="5">
+<lyxtabular version="3" rows="28" columns="5">
 <features tabularvalignment="middle">
 <column alignment="center" valignment="top" width="20text%">
 <column alignment="center" valignment="top" width="20text%">
@@ -3685,6 +3685,54 @@ pH data only
 <cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
 \begin_inset Text
 
+\begin_layout Plain Layout
+Amber GBNSR6 output
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+.gbnsr6
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+gbnsr6
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+-
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Energy information,
+ Read Only
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
 \begin_layout Plain Layout
 Density Peaks
 \end_layout
@@ -44909,7 +44957,8 @@ surftension
 \end_layout
 
 \begin_layout Description
-<mask> Atoms used to build the instantaneous interface (typically heavy atoms of the liquid, e.g.
+<mask> Atoms used to build the instantaneous interface (typically heavy atoms of the liquid,
+ e.g.
  :WAT@O).
 \end_layout
 
@@ -44962,7 +45011,8 @@ interface
 \begin_inset space ~
 \end_inset
 
-{willard|itim} Instantaneous interface: Willard–Chandler density isosurface (default) or ITIM per-column min/max.
+{willard|itim} Instantaneous interface:
+ Willard–Chandler density isosurface (default) or ITIM per-column min/max.
 \end_layout
 
 \begin_layout Description
@@ -44996,8 +45046,7 @@ qmin / qmax
 \begin_inset space ~
 \end_inset
 
-Fit window in Ang.
-^-1.
+Fit window in Ang.^-1.
  If qmin is omitted it is 2*pi/max(L1,L2) from the first good frame.
 \end_layout
 
@@ -45015,7 +45064,8 @@ nblock / dt / blocktime
 \end_inset
 
 Block averaging.
- blocktime (ps) with dt (analyzed-frame spacing, ps) sets nblock.
+ blocktime (ps) with dt (analyzed-frame spacing,
+ ps) sets nblock.
 \end_layout
 
 \begin_layout Description
@@ -45041,7 +45091,9 @@ Optional output files.
 \begin_inset space ~
 \end_inset
 
-Force Grace or gnuplot format for spectrum, roughness, or block files.
+Force Grace or gnuplot format for spectrum,
+ roughness,
+ or block files.
 \end_layout
 
 \end_deeper
@@ -45058,7 +45110,9 @@ literal "true"
  Frames are recentered so the film sits at mid-box along the chosen normal.
  The default interface is a Willard–Chandler isosurface of a Gaussian-smoothed number-density field;
  interface itim uses per-column min/max of <mask> (probe radius to 0).
- The height field is Fourier transformed (numpy fft2 convention), shell-averaged to S(q), and gamma (mN/m) is taken from the small-q plateau of q^2 S(q).
+ The height field is Fourier transformed (numpy fft2 convention),
+ shell-averaged to S(q),
+ and gamma (mN/m) is taken from the small-q plateau of q^2 S(q).
  A Helfrich fit of 1/(q^2 S) vs q^2 on the same window gives kappa in kT.
  The trajectory is assumed NVT (fixed lateral box lengths).
  The 3-D Gaussian filter is OpenMP-parallel.
@@ -48984,6 +49038,10 @@ avg <dset0> [<dset1> ...] [torsion] [out <file>] [oversets]
 
 \end_layout
 
+\begin_layout LyX-Code
+    [ffmt {GEN|DBL|SCI}] [fwidth <width>] [fprec <precision>]
+\end_layout
+
 \begin_deeper
 \begin_layout Description
 <dsetX> Data set(s) to calculate the average for.
@@ -49020,6 +49078,32 @@ avg <dset0> [<dset1> ...] [torsion] [out <file>] [oversets]
 [nostdout] If 'nostdout' specified do not write averages to STDOUT when 'out' not specified.
 \end_layout
 
+\begin_layout Description
+[ffmt
+\begin_inset space ~
+\end_inset
+
+{GEN|DBL|SCI}] Change the default floating point output format to GEN=general,
+ DBL=double,
+ or SCI=scientific.
+\end_layout
+
+\begin_layout Description
+[fwidth
+\begin_inset space ~
+\end_inset
+
+<width>] Change the default floating point output width to <width>.
+\end_layout
+
+\begin_layout Description
+[fprec
+\begin_inset space ~
+\end_inset
+
+<precision>] Change the default floating point output precision to <precision>.
+\end_layout
+
 \begin_layout Standard
 DataSets Created (not 
 \series bold

src/AmberEterm.cpp

@@ -0,0 +1,128 @@
+#include "AmberEterm.h"
+#include "CpptrajStdio.h"
+#include "StringRoutines.h" // NoTrailingWhitespace, validDouble
+
+using namespace Cpptraj;
+
+AmberEterm::AmberEterm() {
+  // Populate the term name to index map. In some cases, multiple term names
+  // map to the same index.
+  termIdxMap_.insert(NameIdxPair("Etot", ETOT));
+  termIdxMap_.insert(NameIdxPair("EPtot", EPTOT));
+  termIdxMap_.insert(NameIdxPair("GMAX", GMAX)); // Not necessary?
+  termIdxMap_.insert(NameIdxPair("BOND", BOND));
+  termIdxMap_.insert(NameIdxPair("ANGLE", ANGLE));
+  termIdxMap_.insert(NameIdxPair("DIHED", DIHED));
+  termIdxMap_.insert(NameIdxPair("VDWAALS", VDWAALS));
+  termIdxMap_.insert(NameIdxPair("EEL", EEL));
+  termIdxMap_.insert(NameIdxPair("EELEC", EEL));
+  termIdxMap_.insert(NameIdxPair("EGB", EGB));
+  termIdxMap_.insert(NameIdxPair("EPB", EPB));
+  termIdxMap_.insert(NameIdxPair("ECAVITY", ECAVITY));
+  termIdxMap_.insert(NameIdxPair("EDISPER", EDISPER));
+  termIdxMap_.insert(NameIdxPair("1-4 VDW", VDW14));
+  termIdxMap_.insert(NameIdxPair("1-4 NB", VDW14));
+  termIdxMap_.insert(NameIdxPair("1-4 EEL", EEL14));
+  termIdxMap_.insert(NameIdxPair("RESTRAINT", RESTRAINT));
+  termIdxMap_.insert(NameIdxPair("EAMBER", EAMBER));
+  termIdxMap_.insert(NameIdxPair("Density", DENSITY));
+  termIdxMap_.insert(NameIdxPair("RMS", RMS)); // Not necessary?
+  termIdxMap_.insert(NameIdxPair("EKtot", EKTOT));
+  termIdxMap_.insert(NameIdxPair("ESURF", ESURF));
+  termIdxMap_.insert(NameIdxPair("EAMD_BOOST", EAMD_BOOST));
+  termIdxMap_.insert(NameIdxPair("VOLUME", VOLUME));
+  termIdxMap_.insert(NameIdxPair("TEMP(K)", TEMP));
+  termIdxMap_.insert(NameIdxPair("PRESS", PRESS));
+  termIdxMap_.insert(NameIdxPair("DV/DL", DVDL));
+  termIdxMap_.insert(NameIdxPair("CMAP", CMAP));
+}
+
+/** Names corresponding to FieldType. */
+const char* AmberEterm::Enames_[] = {
+  "Etot",   "EPtot",  "GMAX",  "BOND",
+  "ANGLE",  "DIHED",  "VDW",   "EELEC",      "EGB",     "EPB", "ECAVITY", "EDISPER",
+  "VDW1-4", "EEL1-4", "RST",   "EAMBER",     "Density",
+  "RMS",    "EKtot",  "ESURF", "EAMD_BOOST", "VOLUME",  "TEMP",
+  "PRESS",  "DVDL",   "CMAP",  0
+};
+
+/** \return FieldType corresponding to given term name, or N_FIELDTYPES if
+  *         not recognized.
+  */
+AmberEterm::FieldType AmberEterm::getTermIdx(std::string const& name) const {
+  NameIdxMap::const_iterator it = termIdxMap_.find( name );
+  if (it == termIdxMap_.end()) {
+    return (FieldType)N_FIELDTYPES;
+  } else {
+    return (FieldType)it->second;
+  }
+}
+
+/** Allocate an array with enough space for all energy terms. */
+AmberEterm::Darray AmberEterm::AllocEnergyArray() {
+  return Darray(N_FIELDTYPES, 0);
+}
+
+/** Allocate a boolean array to indicate whether the term was seen by GetAmberEterms. */
+std::vector<bool> AmberEterm::AllocExistsArray() {
+  return std::vector<bool>(N_FIELDTYPES, false);
+}
+
+/** Parse the given line for energy terms of format <name>=<value>. */
+int AmberEterm::GetAmberEterms(const char* ptr, Darray& Energy, std::vector<bool>& EnergyExists) const {
+  //mprintf("DBG: [%s]\n", ptr);
+  if (ptr == 0 || ptr[0] == '|') return 0;
+  const char* beg = ptr;
+  //          111111111122222222223
+  //0123456789012345678901234567890
+  // NSTEP =        0   TIME(PS) =       0.000  TEMP(K) =   435.99  PRESS =-10207.6
+  bool eol = false;
+  while (!eol) {
+    // Skip leading whitespace
+    while (*beg == ' ' && *beg != '\0') ++beg;
+    if (*beg == '\0') {
+      // Line is blank or no more terms. Bail out.
+      break;
+    }
+    //mprintf("DBG: beg= %c\n", *beg);
+    // Search for next '='
+    const char* eq = beg + 1;
+    while (*eq != '=' && *eq != '\0') ++eq;
+    if (*eq == '\0')
+      eol = true;
+    else {
+      // Search for end token. Start just after '='.
+      const char* val = eq + 1;
+      // Skip leading whitespace
+      while (*val == ' ' && *val != '\0') ++val;
+      if (*val == '\0') {
+        eol = true;
+        mprintf("Warning: EOL encountered before energy term could be read.\n");
+        return 1;
+      } else {
+        //mprintf("DBG: val= %c\n", *val);
+        // Search for next whitespace or line end.
+        const char* end = val + 1;
+        while (*end != ' ' && *end != '\0' && *end != '\n' && *end != '\r') ++end;
+        // Term is now complete. Convert.
+        std::string valstr(val, end);
+        //mprintf("DBG: valstr= '%s'\n", valstr.c_str());
+        std::string termName = NoTrailingWhitespace(std::string(beg,eq));
+        FieldType Eindex = getTermIdx(termName);
+        if (Eindex != N_FIELDTYPES) {
+          if (!validDouble(valstr)) {
+            mprintf("Warning: Invalid number detected: %s = %s\n", termName.c_str(), valstr.c_str());
+          } else {
+            //mprintf("DBG: %s = %s\n", termName.c_str(), valstr.c_str());
+            Energy[Eindex] = atof( valstr.c_str() );
+            EnergyExists[Eindex] = true;
+          }
+        }
+        beg = end;
+      }
+    }
+  } // END loop over line
+
+  return 0;
+}
+

src/AmberEterm.h

@@ -0,0 +1,35 @@
+#ifndef INC_AMBERETERM_H
+#define INC_AMBERETERM_H
+#include <vector>
+#include <string>
+#include <map>
+namespace Cpptraj {
+class AmberEterm {
+  public:
+    typedef std::vector<double> Darray;
+
+    enum FieldType { ETOT= 0, EPTOT, GMAX, BOND,
+                     ANGLE, DIHED, VDWAALS, EEL, EGB, EPB, ECAVITY, EDISPER,
+                     VDW14, EEL14, RESTRAINT, EAMBER, DENSITY,
+                     RMS, EKTOT, ESURF, EAMD_BOOST, VOLUME, TEMP,
+                     PRESS, DVDL, CMAP, N_FIELDTYPES };
+
+    AmberEterm();
+    Darray AllocEnergyArray();
+    std::vector<bool> AllocExistsArray();
+    int GetAmberEterms(const char*, Darray&, std::vector<bool>&) const;
+    static int NenergyTerms() { return (int)N_FIELDTYPES; }
+    static const char* Ename(int i) { return Enames_[i]; }
+  private:
+    static const char* Enames_[];
+
+    typedef std::map<std::string, unsigned int> NameIdxMap;
+    typedef std::pair<std::string, unsigned int> NameIdxPair;
+
+    FieldType getTermIdx(std::string const&) const;
+
+    /// Map field names to indices into energy sets.
+    NameIdxMap termIdxMap_;
+};
+}
+#endif