Add hmassrepartition command

doc/cpptraj.lyx

@@ -11726,7 +11726,7 @@ The following topology related commands are available:
 \begin_layout Standard
 \align center
 \begin_inset Tabular
-<lyxtabular version="3" rows="21" columns="2">
+<lyxtabular version="3" rows="22" columns="2">
 <features tabularvalignment="middle">
 <column alignment="center" valignment="top">
 <column alignment="center" valignment="top">
@@ -11898,6 +11898,26 @@ Print dihedral info for selected atoms.
 <cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
 \begin_inset Text
 
+\begin_layout Plain Layout
+hmassrepartition
+\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
+Perform hydrogen mass repartitioning.
+\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
 improperinfo, impropers, printimpropers
 \end_layout
@@ -13043,6 +13063,86 @@ If 4 masks are given instead of 1, print info for dihedrals with first atom
  in <mask4>.
 \end_layout
 
+\begin_layout Subsection
+hmassrepartition
+\end_layout
+
+\begin_layout LyX-Code
+hmassrepartition [parm <name> | crdset <set> | parmindex <#> | <#>]
+\end_layout
+
+\begin_layout LyX-Code
+             	[<mask>] [hmass <hydrogen new mass>] [dowater]
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+parm
+\begin_inset space ~
+\end_inset
+
+<name> Modify topology selected by name.
+\end_layout
+
+\begin_layout Description
+crdset
+\begin_inset space ~
+\end_inset
+
+<set> Modify topology of COORDS set.
+\end_layout
+
+\begin_layout Description
+parmindex
+\begin_inset space ~
+\end_inset
+
+<#>
+\begin_inset space ~
+\end_inset
+
+|
+\begin_inset space ~
+\end_inset
+
+<#> Modify topology selected by index <#> (starting from 0).
+\end_layout
+
+\begin_layout Description
+<mask> Atoms to modify (all solute atoms by default).
+\end_layout
+
+\begin_layout Description
+hmass
+\begin_inset space ~
+\end_inset
+
+<hydrogen
+\begin_inset space ~
+\end_inset
+
+new
+\begin_inset space ~
+\end_inset
+
+mass> Mass to change hydrogens to (3.024 u by default).
+\end_layout
+
+\begin_layout Description
+dowater If specified, modify water hydrogen mass as well.
+\end_layout
+
+\end_deeper
+\begin_layout Standard
+Perform hydrogen mass repartitioning on the specified topology.
+ Hydrogen mass repartitioning means that for a given heavy atom, the mass
+ of all bonded hydrogens are increased (to 3.024 u by default) and the mass
+ of that heavy atom is decreased so as to maintain the same overall mass.
+ The main use case is to allow longer time steps for molecular dynamics
+ integration due to reduced frequency of vibration of bonds to hydrogen
+ atoms.
+\end_layout
+
 \begin_layout Subsection
 improperinfo | impropers | printimpropers
 \end_layout

src/Command.cpp

@@ -57,6 +57,7 @@
 // ----- TOPOLOGY --------------------------------------------------------------
 #include "Exec_Change.h"
 #include "Exec_CompareTop.h"
+#include "Exec_HmassRepartition.h"
 #include "Exec_ParmBox.h"
 #include "Exec_ParmSolvent.h"
 #include "Exec_ParmStrip.h"
@@ -280,6 +281,7 @@ void Command::Init() {
   Command::AddCmd( new Exec_ChargeInfo(),    Cmd::EXE, 1, "charge" );
   Command::AddCmd( new Exec_CompareTop(),    Cmd::EXE, 1, "comparetop" );
   Command::AddCmd( new Exec_DihedralInfo(),Cmd::EXE, 3,"dihedrals","dihedralinfo","printdihedrals");
+  Command::AddCmd( new Exec_HmassRepartition(),Cmd::EXE, 1, "hmassrepartition" );
   Command::AddCmd( new Exec_ImproperInfo(),Cmd::EXE, 3,"impropers","improperinfo","printimpropers");
   Command::AddCmd( new Exec_MassInfo(),      Cmd::EXE, 1, "mass" );
   Command::AddCmd( new Exec_MolInfo(),       Cmd::EXE, 1, "molinfo" );

src/Exec_HmassRepartition.cpp

@@ -0,0 +1,112 @@
+#include "Exec_HmassRepartition.h"
+#include "CharMask.h"
+#include "CpptrajStdio.h"
+
+// Exec_HmassRepartition::Help()
+void Exec_HmassRepartition::Help() const
+{
+  mprintf("\t[%s]\n"
+          "\t[<mask>] [hmass <hydrogen new mass>] [dowater]\n", DataSetList::TopIdxArgs);
+  mprintf("  Perform hydrogen mass repartitioning for atoms selected by <mask>\n"
+          "  (all solute atoms by default).\n");
+}
+
+// Exec_HmassRepartition::Execute()
+Exec::RetType Exec_HmassRepartition::Execute(CpptrajState& State, ArgList& argIn)
+{
+  double hmass = argIn.getKeyDouble("hmass", 3.024);
+  bool do_water = argIn.hasKey("dowater");
+
+  std::string maskArg = argIn.GetMaskNext();
+  CharMask mask;
+  if (mask.SetMaskString(maskArg)) {
+    mprinterr("Error: Could not process mask string.\n");
+    return CpptrajState::ERR;
+  }
+  // Get Topology
+  Topology* topIn = State.DSL().GetTopByIndex( argIn );
+  if (topIn == 0) return CpptrajState::ERR;
+
+  mprintf("\tPerforming hydrogen mask repartitioning for atoms selected by mask '%s'\n",
+          mask.MaskString());
+  mprintf("\tTopology is '%s'\n", topIn->c_str());
+  mprintf("\tHydrogen masses will be changed to %f amu\n", hmass);
+  if (do_water)
+    mprintf("\tRepartitioning solvent hydrogen masses as well.\n");
+  else
+    mprintf("\tSkipping solvent.\n");
+
+  if (topIn->SetupCharMask( mask )) {
+    mprinterr("Error: Could not set up atom mask.\n");
+    return CpptrajState::ERR;
+  }
+  mask.MaskInfo();
+  if (mask.None()) {
+    mprintf("Warning: No atoms selected.\n");
+    return CpptrajState::OK;
+  }
+
+  if (repartition(*topIn, hmass, mask, do_water, State.Debug()))
+    return CpptrajState::ERR; 
+
+  return CpptrajState::OK;
+}
+
+/** Do hydrogen mass repartitioning. Change mass of all hydrogens to 
+  * the given mass. Adjust mass of bonded heavy atom to maintain the
+  * same overall mass.
+  */
+int Exec_HmassRepartition::repartition(Topology& topIn, double hmass, CharMask const& cmaskIn,
+                                       bool do_water, int debugIn)
+{
+  AtomMask amask( cmaskIn.ConvertToIntMask(), cmaskIn.Natom() );
+
+  unsigned int n_h_changed = 0;
+  unsigned int n_heavy_changed = 0;
+  for (AtomMask::const_iterator at = amask.begin(); at != amask.end(); ++at)
+  {
+    // Skip solvent
+    if (!do_water) {
+      int molnum = topIn[*at].MolNum();
+      if (topIn.Mol(molnum).IsSolvent()) continue;
+    }
+
+    if (topIn[*at].Element() != Atom::HYDROGEN) {
+      Atom& heavyAtom = topIn.SetAtom(*at);
+      double delta = 0;
+      int n_selected_h_atoms = 0;
+      for (Atom::bond_iterator bat = heavyAtom.bondbegin();
+                               bat != heavyAtom.bondend(); ++bat)
+      {
+        if (cmaskIn.AtomInCharMask(*bat) && topIn[*bat].Element() == Atom::HYDROGEN)
+        {
+          n_selected_h_atoms++;
+          Atom& hAtom = topIn.SetAtom(*bat);
+          double diff = hmass - hAtom.Mass();
+          hAtom.SetMass( hmass );
+          n_h_changed++;
+          delta += diff;
+        }
+      }
+      if (n_selected_h_atoms > 0) {
+        if (debugIn > 0)
+          mprintf("\tHeavy atom %s bonded to %i hydrogens, subtract %f amu.\n",
+                  topIn.TruncResAtomNameNum(*at).c_str(),
+                  n_selected_h_atoms, delta);
+        double newMass = heavyAtom.Mass() - delta;
+        // Sanity check
+        if (newMass < 0) {
+          mprinterr("Error: New mass for '%s' would be less than 0.\n",
+                    topIn.TruncResAtomNameNum(*at).c_str());
+          return 1;
+        }
+        heavyAtom.SetMass( newMass );
+        n_heavy_changed++;
+      }
+    }
+  }
+  mprintf("\t%u hydrogen masses modified.\n", n_h_changed);
+  mprintf("\t%u heavy atom masses modified.\n", n_heavy_changed);
+
+  return 0;
+}

src/Exec_HmassRepartition.h

@@ -0,0 +1,14 @@
+#ifndef INC_EXEC_HMASSREPARTITION_H
+#define INC_EXEC_HMASSREPARTITION_H
+#include "Exec.h"
+/// Do hydrogen mass repartitioning on a Topology 
+class Exec_HmassRepartition : public Exec {
+  public:
+    Exec_HmassRepartition() : Exec(PARM) {}
+    void Help() const;
+    DispatchObject* Alloc() const { return (DispatchObject*)new Exec_HmassRepartition(); }
+    RetType Execute(CpptrajState&, ArgList&);
+  private:
+    static int repartition(Topology&, double, CharMask const&, bool, int);
+};
+#endif

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.2.3"
+#define CPPTRAJ_INTERNAL_VERSION "V6.2.4"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif