Add change mass/charge command. (#1048)

* Start adding change mass command. * Finish initial incarnation of change mass * Add change mass test * Add change mass from data set test * Report data set name * Add ability to change charge as well * Fix help text. Add change mass|charge to manual. * 6.20.4. Revision bump for change mass|charge * Protect test in parallel * Fix printf format
Amber-MD · Aug 23, 2023 · 21def2c · 21def2c
1 parent e520de5
commit 21def2c
Show file tree

Hide file tree

Showing 8 changed files with 293 additions and 15 deletions.
diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx
@@ -13743,7 +13743,14 @@ change [parm <name> | parmindex <#> | <#> |
 
 \begin_layout LyX-Code
       bondparm <make1> [<mask2>] {setrk|scalerk|setreq|scalereq} <value>
- }
+\end_layout
+
+\begin_layout LyX-Code
+      {mass|charge} [of <mask>] {to <value>|fromset <data set>}
+\end_layout
+
+\begin_layout LyX-Code
+    }
 \end_layout
 
 \begin_deeper
@@ -14033,6 +14040,41 @@ setreq Set bond equilibrium lengths to <value>.
 scalereq Scale bond equilibirum lengths by <value>.
 \end_layout
 
+\end_deeper
+\begin_layout Description
+mass|charge Change mass or charge in specified topology.
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+of
+\begin_inset space ~
+\end_inset
+
+<mask> Atoms to change mass/charge of.
+\end_layout
+
+\begin_layout Description
+to
+\begin_inset space ~
+\end_inset
+
+<value> Value to change mass/charge to.
+\end_layout
+
+\begin_layout Description
+fromset
+\begin_inset space ~
+\end_inset
+
+<data
+\begin_inset space ~
+\end_inset
+
+set> Use values in <data set> for mass/charge; must have the same number
+ of values as atoms selected by <mask>.
+\end_layout
+
 \end_deeper
 \end_deeper
 \begin_layout Standard

diff --git a/src/Exec_Change.cpp b/src/Exec_Change.cpp
@@ -3,6 +3,7 @@
 #include "CpptrajStdio.h"
 #include "TypeNameHolder.h"
 #include "ParameterTypes.h"
+#include "DataSet_1D.h"
 
 // Exec_Change::Help()
 void Exec_Change::Help() const
@@ -17,6 +18,8 @@ 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}\n"
           "  Change specified parts of topology or topology of a COORDS data set.\n",
           DataSetList::TopArgs);
 }
@@ -26,7 +29,8 @@ Exec::RetType Exec_Change::Execute(CpptrajState& State, ArgList& argIn)
 {
   // Change type
   enum ChangeType { UNKNOWN = 0, RESNAME, CHAINID, ORESNUMS, ICODES,
-                    ATOMNAME, ADDBOND, REMOVEBONDS, SPLITRES, BONDPARM };
+                    ATOMNAME, ADDBOND, REMOVEBONDS, SPLITRES, BONDPARM,
+                    MASS, CHARGE };
   ChangeType type = UNKNOWN;
   if (argIn.hasKey("resname"))
     type = RESNAME;
@@ -46,6 +50,10 @@ Exec::RetType Exec_Change::Execute(CpptrajState& State, ArgList& argIn)
     type = BONDPARM;
   else if (argIn.hasKey("splitres"))
     type = SPLITRES;
+  else if (argIn.hasKey("mass"))
+    type = MASS;
+  else if (argIn.hasKey("charge"))
+    type = CHARGE;
   if (type == UNKNOWN) {
     mprinterr("Error: No change type specified.\n");
     return CpptrajState::ERR;
@@ -77,6 +85,8 @@ Exec::RetType Exec_Change::Execute(CpptrajState& State, ArgList& argIn)
     case REMOVEBONDS : err = RemoveBonds(State, *parm, argIn); break;
     case BONDPARM    : err = ChangeBondParameters(*parm, argIn); break;
     case SPLITRES    : err = ChangeSplitRes(*parm, argIn); break;
+    case MASS        : err = ChangeMassOrCharge(*parm, argIn, State.DSL(), 0); break;
+    case CHARGE      : err = ChangeMassOrCharge(*parm, argIn, State.DSL(), 1); break;
     case UNKNOWN  : err = 1; // sanity check
   }
   if (err != 0) return CpptrajState::ERR;
@@ -590,5 +600,89 @@ int Exec_Change::ChangeBondParameters(Topology& topIn, ArgList& argIn) const {
     topIn.AddBond(it->A1(), it->A2(), bp);
   }
 
+  return 0;
+}
+
+/** Function to change specific value in topology. */
+static inline void changeTopVal(Topology& topIn, int atnum, int typeIn, double newVal,
+                                const char** desc)
+{
+  Atom& currentAtom = topIn.SetAtom(atnum);
+  double oldVal;
+  switch (typeIn) {
+    case 0 :
+      oldVal = currentAtom.Mass();
+      currentAtom.SetMass( newVal );
+      break;
+    case 1 :
+      oldVal = currentAtom.Charge();
+      currentAtom.SetCharge( 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.
+  * \param typeIn: 0=mass, 1=charge
+  */
+int Exec_Change::ChangeMassOrCharge(Topology& topIn, ArgList& argIn,
+                                    DataSetList const& DSL, int typeIn)
+const
+{
+  // sanity check
+  if (typeIn > 1 || typeIn < 0) {
+    mprinterr("Internal Error: typeIn is not 0 or 1.\n");
+    return 1;
+  }
+  static const char* desc[] = { "mass", "charge" };
+
+  std::string maskExpression = argIn.GetStringKey("of");
+  AtomMask atomsToChange;
+  if (atomsToChange.SetMaskString( maskExpression )) {
+    mprinterr("Error: Could not set mask expression.\n");
+    return 1;
+  }
+  if (topIn.SetupIntegerMask( atomsToChange )) {
+    mprinterr("Error: Could not set up mask.\n");
+    return 1;
+  }
+  if (atomsToChange.None()) {
+    mprintf("Warning: Mask '%s' selects no atoms.\n", atomsToChange.MaskString());
+    return 0;
+  }
+  atomsToChange.MaskInfo();
+
+  std::string fromSet = argIn.GetStringKey("fromset");
+  if (!fromSet.empty()) {
+    DataSet* ds = DSL.GetDataSet( fromSet );
+    if (ds == 0) {
+      mprinterr("Error: No set selected by '%s'\n", fromSet.c_str());
+      return 1;
+    }
+    if (ds->Group() != DataSet::SCALAR_1D) {
+      mprinterr("Error: Data set '%s' is not scalar 1D.\n", ds->legend());
+      return 1;
+    }
+    mprintf("\tUsing data from '%s' for %s.\n", ds->legend(), desc[typeIn]);
+    if (ds->Size() != (unsigned int)atomsToChange.Nselected()) {
+      mprinterr("Error: %i atoms to change %s of, but set '%s' has %zu elements.\n",
+                atomsToChange.Nselected(), desc[typeIn], ds->legend(), ds->Size());
+      return 1;
+    }
+    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);
+    }
+  } else {
+    if (!argIn.Contains("to")) {
+      mprinterr("Error: Expected either 'fromset' or 'to' 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);
+    }
+  }
+
   return 0;
 } 
diff --git a/src/Exec_Change.h b/src/Exec_Change.h
@@ -21,5 +21,6 @@ class Exec_Change : public Exec {
     int AddBond(Topology&, ArgList&) const;
     int RemoveBonds(CpptrajState&, Topology&, ArgList&) const;
     int ChangeBondParameters(Topology&, ArgList&) const;
+    int ChangeMassOrCharge(Topology&, ArgList&, DataSetList const&, int) const;
 };
 #endif
diff --git a/src/Version.h b/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.20.3"
+#define CPPTRAJ_INTERNAL_VERSION "V6.20.4"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif
diff --git a/src/cpptrajdepend b/src/cpptrajdepend
@@ -283,7 +283,7 @@ Exec_AddMissingRes.o : Exec_AddMissingRes.cpp Action.h ActionFrameCounter.h Acti
 Exec_Analyze.o : Exec_Analyze.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cmd.h CmdList.h Command.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Analyze.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.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 ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h
 Exec_Calc.o : Exec_Calc.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Calc.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h RPNcalc.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h
 Exec_CatCrd.o : Exec_CatCrd.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_CatCrd.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.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 ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h
-Exec_Change.o : Exec_Change.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Change.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.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 ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h
+Exec_Change.o : Exec_Change.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Change.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.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 ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h
 Exec_ClusterMap.o : Exec_ClusterMap.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h ClusterMap.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_2D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_MatrixFlt.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_ClusterMap.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.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 ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h
 Exec_CombineCoords.o : Exec_CombineCoords.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_CombineCoords.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.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 ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h
 Exec_Commands.o : Exec_Commands.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Commands.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.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 ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h

diff --git a/test/Test_Change/AFV.fluctMass.dat.save b/test/Test_Change/AFV.fluctMass.dat.save
@@ -0,0 +1,50 @@
+#Atom Name  #Res Name  #Mol Type   Charge     Mass GBradius El     rVDW     eVDW
+    1 N        1 ALA      1 N3     0.1414   1.3392   1.5500  N   1.8240   0.1700
+    2 H1       1 ALA      1 H      0.1997   1.8307   1.3000  H   0.6000   0.0157
+    3 H2       1 ALA      1 H      0.1997   1.7955   1.3000  H   0.6000   0.0157
+    4 H3       1 ALA      1 H      0.1997   1.8237   1.3000  H   0.6000   0.0157
+    5 CA       1 ALA      1 CX     0.0962   0.8400   1.7000  C   1.9080   0.1094
+    6 HA       1 ALA      1 HP     0.0889   1.3625   1.3000  H   1.1000   0.0157
+    7 CB       1 ALA      1 CT    -0.0597   1.5867   1.7000  C   1.9080   0.1094
+    8 HB1      1 ALA      1 HC     0.0300   2.0004   1.3000  H   1.4870   0.0157
+    9 HB2      1 ALA      1 HC     0.0300   2.2382   1.3000  H   1.4870   0.0157
+   10 HB3      1 ALA      1 HC     0.0300   2.0209   1.3000  H   1.4870   0.0157
+   11 C        1 ALA      1 C      0.6163   0.6606   1.7000  C   1.9080   0.0860
+   12 O        1 ALA      1 O     -0.5722   1.0294   1.5000  O   1.6612   0.2100
+   13 N        2 PHE      1 N     -0.4157   0.6194   1.5500  N   1.8240   0.1700
+   14 H        2 PHE      1 H      0.2719   0.9363   1.3000  H   0.6000   0.0157
+   15 CA       2 PHE      1 CX    -0.0024   0.5888   1.7000  C   1.9080   0.1094
+   16 HA       2 PHE      1 H1     0.0978   0.6938   1.3000  H   1.3870   0.0157
+   17 CB       2 PHE      1 CT    -0.0343   0.8447   1.7000  C   1.9080   0.1094
+   18 HB2      2 PHE      1 HC     0.0295   1.4507   1.3000  H   1.4870   0.0157
+   19 HB3      2 PHE      1 HC     0.0295   1.5059   1.3000  H   1.4870   0.0157
+   20 CG       2 PHE      1 CA     0.0118   0.3897   1.7000  C   1.9080   0.0860
+   21 CD1      2 PHE      1 CA    -0.1256   1.2616   1.7000  C   1.9080   0.0860
+   22 HD1      2 PHE      1 HA     0.1330   1.9513   1.3000  H   1.4590   0.0150
+   23 CE1      2 PHE      1 CA    -0.1704   1.7951   1.7000  C   1.9080   0.0860
+   24 HE1      2 PHE      1 HA     0.1430   2.5937   1.3000  H   1.4590   0.0150
+   25 CZ       2 PHE      1 CA    -0.1072   1.8400   1.7000  C   1.9080   0.0860
+   26 HZ       2 PHE      1 HA     0.1297   2.4711   1.3000  H   1.4590   0.0150
+   27 CE2      2 PHE      1 CA    -0.1704   1.7855   1.7000  C   1.9080   0.0860
+   28 HE2      2 PHE      1 HA     0.1430   2.5780   1.3000  H   1.4590   0.0150
+   29 CD2      2 PHE      1 CA    -0.1256   1.2652   1.7000  C   1.9080   0.0860
+   30 HD2      2 PHE      1 HA     0.1330   1.9666   1.3000  H   1.4590   0.0150
+   31 C        2 PHE      1 C      0.5973   0.4730   1.7000  C   1.9080   0.0860
+   32 O        2 PHE      1 O     -0.5679   0.6642   1.5000  O   1.6612   0.2100
+   33 N        3 VAL      1 N     -0.3821   0.5251   1.5500  N   1.8240   0.1700
+   34 H        3 VAL      1 H      0.2681   0.7312   1.3000  H   0.6000   0.0157
+   35 CA       3 VAL      1 CX    -0.3438   0.5870   1.7000  C   1.9080   0.1094
+   36 HA       3 VAL      1 H1     0.1438   1.0249   1.3000  H   1.3870   0.0157
+   37 CB       3 VAL      1 CT     0.1940   1.1005   1.7000  C   1.9080   0.1094
+   38 HB       3 VAL      1 HC     0.0308   1.5142   1.3000  H   1.4870   0.0157
+   39 CG1      3 VAL      1 CT    -0.3064   1.5705   1.7000  C   1.9080   0.1094
+   40 HG11     3 VAL      1 HC     0.0836   2.1092   1.3000  H   1.4870   0.0157
+   41 HG12     3 VAL      1 HC     0.0836   2.0668   1.3000  H   1.4870   0.0157
+   42 HG13     3 VAL      1 HC     0.0836   1.8651   1.3000  H   1.4870   0.0157
+   43 CG2      3 VAL      1 CT    -0.3064   1.3964   1.7000  C   1.9080   0.1094
+   44 HG21     3 VAL      1 HC     0.0836   1.9660   1.3000  H   1.4870   0.0157
+   45 HG22     3 VAL      1 HC     0.0836   1.7581   1.3000  H   1.4870   0.0157
+   46 HG23     3 VAL      1 HC     0.0836   1.6678   1.3000  H   1.4870   0.0157
+   47 C        3 VAL      1 C      0.8350   0.8819   1.7000  C   1.9080   0.0860
+   48 O        3 VAL      1 O2    -0.8173   1.3691   1.5000  O   1.6612   0.2100
+   49 OXT      3 VAL      1 O2    -0.8173   1.7206   1.5000  O   1.6612   0.2100