Add ability to offset charge/mass, save charge/mass to variables (#1077)

* Add test of changing charge

* Add 'by' keyword for change charge/mass

* Add change charge/mass by tests

* 6.25.1. Revision bump for adding 'change {mass|charge} by'

* Add charge and mass to 'set inmask'

* Allow string vars to be written to dat files for easier debug/check

* Add set charge/mass test

* Fix help text

* Update manual

doc/cpptraj.lyx

@@ -2652,7 +2652,11 @@ set { <variable> <OP> <value> |
 \end_layout
 
 \begin_layout LyX-Code
-	                   resnums|oresnums|molnums} inmask <mask>
+	                   resnums|oresnums|molnums|
+\end_layout
+
+\begin_layout LyX-Code
+                       charge|mass} inmask <mask>
 \end_layout
 
 \begin_layout LyX-Code
@@ -2702,8 +2706,9 @@ inmask
 \end_inset
 
 <mask> Set/append a script variable to/by the total number of atoms/residues/mol
-ecules or a range expression of selected atom #s/residue #s/original residue
- #s/molecule #s in the given mask expression.
+ecules in, a range expression of selected atom #s/residue #s/original residue
+ #s/molecule #s in, or the total charge/mass of atoms selected by the given
+ mask expression.
 \end_layout
 
 \begin_deeper
@@ -14837,7 +14842,8 @@ change [parm <name> | parmindex <#> | <#> |
 \end_layout
 
 \begin_layout LyX-Code
-      {mass|charge} [of <mask>] {to <value>|fromset <data set>} |
+      {mass|charge} [of <mask>] {to <value>|by <offset>|fromset <data set>}
+ |
 \end_layout
 
 \begin_layout LyX-Code
@@ -15157,6 +15163,14 @@ to
 <value> Value to change mass/charge to.
 \end_layout
 
+\begin_layout Description
+by
+\begin_inset space ~
+\end_inset
+
+<offset> Value to offset masses/charges by.
+\end_layout
+
 \begin_layout Description
 fromset
 \begin_inset space ~

src/DataIO_Std.cpp

@@ -14,6 +14,7 @@
 #include "DataSet_double.h" // For reading TODO remove dependency?
 #include "DataSet_float.h" // For reading TODO remove dependency?
 #include "DataSet_string.h" // For reading TODO remove dependency?
+#include "DataSet_StringVar.h"
 #include "DataSet_Vector.h" // For reading TODO remove dependency?
 #include "DataSet_Mat3x3.h" // For reading TODO remove dependency?
 #include "DataSet_PairwiseCache_MEM.h" // For reading
@@ -44,6 +45,7 @@ DataIO_Std::DataIO_Std() :
   dims_[0] = 0;
   dims_[1] = 0;
   dims_[2] = 0;
+  SetValid( DataSet::STRINGVAR );
 }
 
 static void PrintColumnError(int idx) {
@@ -1161,6 +1163,9 @@ int DataIO_Std::WriteData(FileName const& fname, DataSetList const& SetList)
     if (SetList[0]->Group() == DataSet::PWCACHE) {
       // Special case of 2D - may have sieved frames.
       err = WriteCmatrix(file, SetList);
+    } else if (SetList[0]->Type() == DataSet::STRINGVAR) {
+      // For debugging string variables
+      err = WriteStringVars(file, SetList);
     } else if (SetList[0]->Ndim() == 1) {
       if (group_ == NO_TYPE) {
         if (isInverted_)
@@ -1178,6 +1183,24 @@ int DataIO_Std::WriteData(FileName const& fname, DataSetList const& SetList)
   return err;
 }
 
+/** For debugging string variables, write out to a file. */
+int DataIO_Std::WriteStringVars(CpptrajFile& file, DataSetList const& Sets) const {
+  for (DataSetList::const_iterator ds = Sets.begin(); ds != Sets.end(); ++ds)
+  {
+    if ( (*ds)->Type() != DataSet::STRINGVAR ) {
+      mprintf("Warning: First variable written to %s was a string variable, skipping '%s'\n",
+              file.Filename().full(), (*ds)->legend());
+    } else {
+      DataSet_StringVar const& var = static_cast<DataSet_StringVar const&>( *(*ds) );
+      file.Write( var.Meta().Name().c_str(), var.Meta().Name().size() );
+      file.Printf(" = ");
+      file.Write( var.Value().c_str(), var.Value().size() );
+      file.Printf("\n");
+    }
+  }
+  return 0;
+}
+
 // DataIO_Std::WriteCmatrix()
 int DataIO_Std::WriteCmatrix(CpptrajFile& file, DataSetList const& Sets) {
   for (DataSetList::const_iterator ds = Sets.begin(); ds != Sets.end(); ++ds)

src/DataIO_Std.h

@@ -31,6 +31,7 @@ class DataIO_Std : public DataIO {
     int Read_Mat3x3(std::string const&,DataSetList&,std::string const&);
     static void WriteNameToBuffer(CpptrajFile&, std::string const&, int,  bool);
     int WriteByGroup(CpptrajFile&, DataSetList const&, GroupType);
+    int WriteStringVars(CpptrajFile&, DataSetList const&) const;
     int WriteCmatrix(CpptrajFile&, DataSetList const&);
     int WriteDataNormal(CpptrajFile&,DataSetList const&);
     int WriteDataInverted(CpptrajFile&,DataSetList const&);

src/Exec_Change.cpp

@@ -18,7 +18,7 @@ void Exec_Change::Help() const
           "\t  addbond <mask1> <mask2> [req <length> <rk> <force constant>] |\n"
           "\t  removebonds <mask1> [<mask2>] [out <file>] |\n"
           "\t  bondparm <mask1> [<mask2>] {setrk|scalerk|setreq|scalereq} <value> |\n"
-          "\t  {mass|charge} [of <mask>] {to <value>|fromset <data set>} |\n"
+          "\t  {mass|charge} [of <mask>] {to <value>|by <offset>|fromset <data set>} |\n"
           "\t  mergeres firstres <start res#> lastres <stop res#>\n"
           "\t}\n"
           "  Change specified parts of topology or topology of a COORDS data set.\n",
@@ -608,22 +608,32 @@ int Exec_Change::ChangeBondParameters(Topology& topIn, ArgList& argIn) const {
 }
 
 /** Function to change specific value in topology. */
-static inline void changeTopVal(Topology& topIn, int atnum, int typeIn, double newVal,
-                                const char** desc)
+void Exec_Change::changeTopVal(Topology& topIn, int atnum, ChangeType typeIn, double newVal)
 {
   Atom& currentAtom = topIn.SetAtom(atnum);
   double oldVal;
   switch (typeIn) {
-    case 0 :
+    case MASS_TO :
       oldVal = currentAtom.Mass();
+      mprintf("\tChanging mass of atom '%s' from %g to %g\n", topIn.AtomMaskName(atnum).c_str(), oldVal, newVal);
       currentAtom.SetMass( newVal );
       break;
-    case 1 :
+    case CHARGE_TO :
       oldVal = currentAtom.Charge();
-      currentAtom.SetCharge( newVal ); break;
+      mprintf("\tChanging charge of atom '%s' from %g to %g\n", topIn.AtomMaskName(atnum).c_str(), oldVal, newVal);
+      currentAtom.SetCharge( newVal );
+      break;
+    case MASS_BY :
+      oldVal = currentAtom.Mass();
+      mprintf("\tChanging mass of atom '%s' by %g to %g\n", topIn.AtomMaskName(atnum).c_str(), newVal, oldVal+newVal);
+      currentAtom.SetMass( oldVal + newVal );
+      break;
+    case CHARGE_BY :
+      oldVal = currentAtom.Charge();
+      mprintf("\tChanging charge of atom '%s' by %g to %g\n", topIn.AtomMaskName(atnum).c_str(), newVal, oldVal+newVal);
+      currentAtom.SetCharge( oldVal + newVal );
+      break;
   }
-  mprintf("\tChanging %s of atom '%s' from %g to %g\n", desc[typeIn],
-          topIn.AtomMaskName(atnum).c_str(), oldVal, newVal);
 }
 
 /** Change mass/charge in topology.
@@ -658,6 +668,7 @@ const
 
   std::string fromSet = argIn.GetStringKey("fromset");
   if (!fromSet.empty()) {
+    // Get charges from a data set
     DataSet* ds = DSL.GetDataSet( fromSet );
     if (ds == 0) {
       mprinterr("Error: No set selected by '%s'\n", fromSet.c_str());
@@ -673,18 +684,46 @@ const
                 atomsToChange.Nselected(), desc[typeIn], ds->legend(), ds->Size());
       return 1;
     }
+    ChangeType ctype;
+    if (typeIn == 0)
+      ctype = MASS_TO;
+    else
+      ctype = CHARGE_TO;
     DataSet_1D const& dset = static_cast<DataSet_1D const&>( *ds );
     for (int idx = 0; idx != atomsToChange.Nselected(); idx++) {
-      changeTopVal(topIn, atomsToChange[idx], typeIn, dset.Dval(idx), desc);
+      changeTopVal(topIn, atomsToChange[idx], ctype, dset.Dval(idx));
     }
   } else {
-    if (!argIn.Contains("to")) {
-      mprinterr("Error: Expected either 'fromset' or 'to' for 'change %s'.\n", desc[typeIn]);
+    // Get charge/offset from command line
+    ChangeType ctype;
+    bool change_to = argIn.Contains("to");
+    bool change_by = argIn.Contains("by");
+    if (!change_to && !change_by) {
+      mprinterr("Error: Expected either 'fromset', 'to', or 'by' for 'change %s'.\n", desc[typeIn]);
       return 1;
     }
-    double newVal = argIn.getKeyDouble("to", 0.0);
-    for (AtomMask::const_iterator at = atomsToChange.begin(); at != atomsToChange.end(); ++at) {
-      changeTopVal(topIn, *at, typeIn, newVal, desc);
+    if (change_to && change_by) {
+      mprinterr("Error: Specify either 'to' or 'by', not both.\n");
+      return 1;
+    }
+    if (change_to) {
+      if (typeIn == 0)
+        ctype = MASS_TO;
+      else
+        ctype = CHARGE_TO;
+      double newVal = argIn.getKeyDouble("to", 0.0);
+      for (AtomMask::const_iterator at = atomsToChange.begin(); at != atomsToChange.end(); ++at) {
+        changeTopVal(topIn, *at, ctype, newVal);
+      }
+    } else if (change_by) {
+      if (typeIn == 0)
+        ctype = MASS_BY;
+      else
+        ctype = CHARGE_BY;
+      double offset = argIn.getKeyDouble("by", 0.0);
+      for (AtomMask::const_iterator at = atomsToChange.begin(); at != atomsToChange.end(); ++at) {
+        changeTopVal(topIn, *at, ctype, offset);
+      }
     }
   }
 

src/Exec_Change.h

@@ -10,6 +10,10 @@ class Exec_Change : public Exec {
     DispatchObject* Alloc() const { return (DispatchObject*)new Exec_Change(); }
     RetType Execute(CpptrajState&, ArgList&);
   private:
+    enum ChangeType { MASS_TO = 0, CHARGE_TO, MASS_BY, CHARGE_BY };
+    /// Change mass or charge
+    static inline void changeTopVal(Topology&, int, ChangeType, double);
+
     int ChangeSplitRes(Topology&, ArgList&) const;
     int ChangeResidueName(Topology&, ArgList&) const;
     int ChangeOresNums(Topology&, ArgList&) const;

src/Exec_Set.cpp

@@ -6,15 +6,16 @@
 void Exec_Set::Help() const {
   mprintf("\t{ <variable> <OP> <value> |\n"
           "\t  <variable> <OP> {atoms|residues|molecules|atomnums|\n"
-          "\t                   resnums|oresnums|molnums} inmask <mask>\n"
+          "\t                   resnums|oresnums|molnums|\n"
+          "\t                   charge|mass} inmask <mask>\n"
           "\t    [%s]\n"
           "\t  <variable> <OP> trajinframes }\n",
           DataSetList::TopIdxArgs);
   mprintf("  Set (<OP> = '=') or append (<OP> = '+=') a script variable.\n"
           "  - Set script variable <variable> to value <value>.\n"
-          "  - Set script variable to the number of atoms/residues/molecules\n"
-          "    or selected atom #s/residue #s/original residue #s/molecule #s\n"
-          "    in the given mask.\n"
+          "  - Set script variable to the number of atoms/residues/molecules in,\n"
+          "    selected atom #s/residue #s/original residue #s/molecule #s in,\n"
+          "    or total charge/mass of atoms selected by the given mask.\n"
           "  - Set script variable to the current number of frames that will\n"
           "     be read from all previous 'trajin' statements.\n");
 }
@@ -72,6 +73,16 @@ Exec::RetType Exec_Set::Execute(CpptrajState& State, ArgList& argIn)
       for (std::vector<int>::const_iterator it = resnums.begin(); it != resnums.end(); ++it)
         oresnums.push_back( top->Res(*it).OriginalResNum() );
       value = ArrayToRangeExpression( oresnums, 0 );
+    } else if (equals.hasKey("charge")) {
+      double sum_q = 0;
+      for (AtomMask::const_iterator it = mask.begin(); it != mask.end(); ++it)
+        sum_q += (*top)[*it].Charge();
+      value = doubleToString( sum_q );
+    } else if (equals.hasKey("mass")) {
+      double sum_m = 0;
+      for (AtomMask::const_iterator it = mask.begin(); it != mask.end(); ++it)
+        sum_m += (*top)[*it].Mass();
+      value = doubleToString( sum_m );
     } else {
       mprinterr("Error: Expected one of: 'atoms', 'residues', 'molecules',\n"
                 "Error:   'atomnums', 'resnums', 'oresnums', or 'molnums'.\n");

src/Version.h

@@ -12,7 +12,7 @@
  * Whenever a number that precedes <revision> is incremented, all subsequent
  * numbers should be reset to 0.
  */
-#define CPPTRAJ_INTERNAL_VERSION "V6.25.0"
+#define CPPTRAJ_INTERNAL_VERSION "V6.25.1"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif

src/cpptrajdepend

@@ -228,7 +228,7 @@ DataIO_Numpy.o : DataIO_Numpy.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h A
 DataIO_OpenDx.o : DataIO_OpenDx.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_OpenDx.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridDbl.h DataSet_GridFlt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 DataIO_Peaks.o : DataIO_Peaks.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Peaks.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Vector_Scalar.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 DataIO_RemLog.o : DataIO_RemLog.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_RemLog.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_RemLog.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
-DataIO_Std.o : DataIO_Std.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Cluster/Cframes.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Std.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_PairwiseCache.h DataSet_PairwiseCache_MEM.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_string.h Dimension.h FileIO.h FileName.h Frame.h GridBin.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
+DataIO_Std.o : DataIO_Std.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Cluster/Cframes.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Std.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_PairwiseCache.h DataSet_PairwiseCache_MEM.h DataSet_StringVar.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_string.h Dimension.h FileIO.h FileName.h Frame.h GridBin.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 DataIO_VecTraj.o : DataIO_VecTraj.cpp ActionFrameCounter.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_VecTraj.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Vector.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ParmFile.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajectoryFile.h Trajout_Single.h TypeNameHolder.h Unit.h Vec3.h
 DataIO_XVG.o : DataIO_XVG.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_XVG.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 DataIO_Xplor.o : DataIO_Xplor.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Xplor.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h

test/Test_Change/AFV.1.offset.dat.save

@@ -0,0 +1,13 @@
+#Atom Name  #Res Name  #Mol Type   Charge     Mass GBradius El     rVDW     eVDW
+    1 N        1 ALA      1 N3    -0.8586  15.0100   1.5500  N   1.8240   0.1700
+    2 H1       1 ALA      1 H     -0.8003   2.0080   1.3000  H   0.6000   0.0157
+    3 H2       1 ALA      1 H     -0.8003   2.0080   1.3000  H   0.6000   0.0157
+    4 H3       1 ALA      1 H     -0.8003   2.0080   1.3000  H   0.6000   0.0157
+    5 CA       1 ALA      1 CX    -0.9038  13.0100   1.7000  C   1.9080   0.1094
+    6 HA       1 ALA      1 HP    -0.9111   2.0080   1.3000  H   1.1000   0.0157
+    7 CB       1 ALA      1 CT    -1.0597  13.0100   1.7000  C   1.9080   0.1094
+    8 HB1      1 ALA      1 HC    -0.9700   2.0080   1.3000  H   1.4870   0.0157
+    9 HB2      1 ALA      1 HC    -0.9700   2.0080   1.3000  H   1.4870   0.0157
+   10 HB3      1 ALA      1 HC    -0.9700   2.0080   1.3000  H   1.4870   0.0157
+   11 C        1 ALA      1 C     -0.3837  13.0100   1.7000  C   1.9080   0.0860
+   12 O        1 ALA      1 O     -1.5722  17.0000   1.5000  O   1.6612   0.2100