Implement ability to calculate vector magnitude in 'vectormath' comma…

…nd. (#1105) * Start code for reading vector magnitudes * Actually read the vector magnitudes * Add magnitude mode to vectormath for calculating vector magnitude * Add vectormath magnitude test * Update help text * Update manual text * Update documentation files * 6.29.3. Revision bump for adding 'magnitude' keyword to 'readdata vector' and 'vectormath'.
Amber-MD · Oct 2, 2024 · 8d1a163 · 8d1a163
1 parent db9e1ae
commit 8d1a163
Show file tree

Hide file tree

Showing 10 changed files with 222 additions and 70 deletions.
diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf
diff --git a/doc/DocumentChecksums.txt b/doc/DocumentChecksums.txt
@@ -1,3 +1,3 @@
 b37726e7a841f6fc695ecd7fb040ffbf  CpptrajDevelopmentGuide.lyx
-c9b19c2bade122b63363f364f86a0a69  cpptraj.lyx
+736ee725ca9cc67d349f67d714f41c20  cpptraj.lyx
 07c4039e732fc2eb1df0c1e0863cb949  CpptrajManual.lyx
diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx
@@ -5640,7 +5640,7 @@ Read
 \end_layout
 
 \begin_layout LyX-Code
- [vector] |
+ [vector [magnitude]] |
 \end_layout
 
 \begin_layout LyX-Code
@@ -5737,14 +5737,22 @@ nosquare2d Read data as XYZ matrix (i.e.
 \begin_layout Description
 vector Read data as vector.
  If indices are present they will be skipped.
- Assume first 3 columns after the index volumn are vector X,
+ Assume first 3 columns after the index column are vector X,
  Y,
  and Z,
  and (if present) the next 3 columns contain vector origin X,
  Y,
  and Z.
 \end_layout
 
+\begin_deeper
+\begin_layout Description
+magnitude If specified,
+ assume vector file last column contains vector magnitudes;
+ read these into a set with Aspect 'Mag'.
+\end_layout
+
+\end_deeper
 \begin_layout Description
 mat3x3 Read data as 3x3 matrix.
  If indices are present they will be skipped.
@@ -55687,7 +55695,7 @@ vectormath vec1 <vecname1> vec2 <vecname2> [out <filename>] [norm] [name <setnam
 \end_layout
 
 \begin_layout LyX-Code
-           [ dotproduct | dotangle | crossproduct ]
+           [ dotproduct | dotangle | crossproduct | magnitude ]
 \begin_inset Separator latexpar
 \end_inset
 
@@ -55740,12 +55748,24 @@ vec2
 [crossproduct] Calculate cross-product of the two vectors.
 \end_layout
 
+\begin_layout Description
+[magnitude] Calculate the magnitude of the vectors selected by 
+\series bold
+vec1
+\series default
+ (no need to specify 
+\series bold
+vec2
+\series default
+).
+\end_layout
+
 \begin_layout Description
 [norm] Normalize the vectors;
  this will affect any subsequent calculations with the vectors.
  This is turned on automatically if 
 \series bold
-dotangle
+dotangle/magnitude
 \series default
  specified.
 \end_layout

diff --git a/src/Analysis_VectorMath.cpp b/src/Analysis_VectorMath.cpp
@@ -1,4 +1,5 @@
 #include <algorithm> // std::min, std::max
+#include <cmath> // sqrt
 #include "Analysis_VectorMath.h"
 #include "CpptrajStdio.h"
 #include "Constants.h"
@@ -7,7 +8,7 @@
 
 /// Strings corresponding to modes, used in output.
 const char* Analysis_VectorMath::ModeString[] = {
-  "Dot product", "Angle from dot product", "Cross product" };
+  "Dot product", "Angle from dot product", "Cross product", "Magnitude" };
 
 // CONSTRUCTOR
 Analysis_VectorMath::Analysis_VectorMath() :
@@ -17,46 +18,24 @@ Analysis_VectorMath::Analysis_VectorMath() :
 
 void Analysis_VectorMath::Help() const {
   mprintf("\tvec1 <vecname1> vec2 <vecname2> [out <filename>] [norm] [name <setname>]\n"
-          "\t[ dotproduct | dotangle | crossproduct ]\n"
+          "\t[ dotproduct | dotangle | crossproduct | magnitude ]\n"
           "  Calculate dot product, angle from dot product (degrees), or cross product\n"
           "  for specified vectors. Either vec1 or vec2 can be size 1, otherwise they\n"
-          "  must both be the same size.\n");
+          "  must both be the same size. If 'magnitude' is specified, just calculate\n"
+          "  the magnitudes of the vectors selected by 'vec1' (no need to specify\n"
+          "  'vec2').\n");
 }
 
 // Analysis_VectorMath::Setup()
 Analysis::RetType Analysis_VectorMath::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn)
 {
-  // Get Vectors
-  DataSetList vsets1 = setup.DSL().SelectGroupSets( analyzeArgs.GetStringKey("vec1"),
-                                                  DataSet::VECTOR_1D );
-  if (vsets1.empty()) {
-    mprinterr("Error: 'vec1' not found.\n");
-    return Analysis::ERR;
-  }
-  DataSetList vsets2 = setup.DSL().SelectGroupSets( analyzeArgs.GetStringKey("vec2"),
-                                                  DataSet::VECTOR_1D );
-  if (vsets2.empty()) {
-    mprinterr("Error: 'vec2' not found.\n");
-    return Analysis::ERR;
-  }
-
-  if (vsets1.size() != vsets2.size()) {
-    mprinterr("Error: 'vec1' (%zu) and 'vec2' (%zu) do not select the same number of sets.\n",
-              vsets1.size(), vsets2.size());
-    return Analysis::ERR;
-  }
-
-  for (DataSetList::const_iterator it = vsets1.begin(); it != vsets1.end(); ++it)
-    vinfo1_.push_back( static_cast<DataSet_Vector*>( *it ) );
-  for (DataSetList::const_iterator it = vsets2.begin(); it != vsets2.end(); ++it)
-    vinfo2_.push_back( static_cast<DataSet_Vector*>( *it ) );
-
   std::string setname = analyzeArgs.GetStringKey("name");
   norm_ = analyzeArgs.hasKey("norm");
   // Check for dotproduct/crossproduct keywords. Default is dotproduct.
   mode_ = DOTPRODUCT;
   DataSet::DataType dtype = DataSet::DOUBLE;
   const char* dname = "Dot";
+  bool requires_two_vecs = true;
   if (analyzeArgs.hasKey("dotproduct")) {
     mode_ = DOTPRODUCT;
   } else if (analyzeArgs.hasKey("dotangle")) {
@@ -67,7 +46,47 @@ Analysis::RetType Analysis_VectorMath::Setup(ArgList& analyzeArgs, AnalysisSetup
     mode_ = CROSSPRODUCT;
     dtype = DataSet::VECTOR;
     dname = "Cross";
+  } else if (analyzeArgs.hasKey("magnitude")) {
+    mode_ = MAGNITUDE;
+    dname = "Mag";
+    requires_two_vecs = false;
+    if (norm_) {
+      mprintf("Warning: 'norm' does not make sense with 'magnitude', ignoring.\n");
+      norm_ = false;
+    }
   }
+
+  // Get Vectors
+  DataSetList vsets1 = setup.DSL().SelectGroupSets( analyzeArgs.GetStringKey("vec1"),
+                                                    DataSet::VECTOR_1D );
+  if (vsets1.empty()) {
+    mprinterr("Error: 'vec1' not found.\n");
+    return Analysis::ERR;
+  }
+
+  DataSetList vsets2;
+  if (requires_two_vecs) {
+    vsets2 = setup.DSL().SelectGroupSets( analyzeArgs.GetStringKey("vec2"),
+                                          DataSet::VECTOR_1D );
+    if (vsets2.empty()) {
+      mprinterr("Error: 'vec2' not found.\n");
+      return Analysis::ERR;
+    }
+
+    if (vsets1.size() != vsets2.size()) {
+      mprinterr("Error: 'vec1' (%zu) and 'vec2' (%zu) do not select the same number of sets.\n",
+                vsets1.size(), vsets2.size());
+      return Analysis::ERR;
+    }
+  }
+
+  for (DataSetList::const_iterator it = vsets1.begin(); it != vsets1.end(); ++it)
+    vinfo1_.push_back( static_cast<DataSet_Vector*>( *it ) );
+  if (requires_two_vecs) {
+    for (DataSetList::const_iterator it = vsets2.begin(); it != vsets2.end(); ++it)
+      vinfo2_.push_back( static_cast<DataSet_Vector*>( *it ) );
+  }
+
   // Set up output file in DataFileList if necessary
   DataFile* outfile = setup.DFL().AddDataFile( analyzeArgs.GetStringKey("out"), analyzeArgs );
   // Set up output data sets based on mode
@@ -89,12 +108,22 @@ Analysis::RetType Analysis_VectorMath::Setup(ArgList& analyzeArgs, AnalysisSetup
 
   // Print Status
   mprintf("    VECTORMATH:");
-  if (vinfo1_.size() == 1)
-    mprintf(" Calculating %s of vectors %s and %s\n", ModeString[mode_], vinfo1_[0]->legend(), vinfo2_[0]->legend());
-  else {
-    mprintf(" Calculating %s of:\n", ModeString[mode_]);
-    for (unsigned int ii = 0; ii < vinfo1_.size(); ii++)
-      mprintf("\t  %s and %s\n", vinfo1_[ii]->legend(), vinfo2_[ii]->legend());
+  if (requires_two_vecs) {
+    if (vinfo1_.size() == 1)
+      mprintf(" Calculating %s of vectors %s and %s\n", ModeString[mode_], vinfo1_[0]->legend(), vinfo2_[0]->legend());
+    else {
+      mprintf(" Calculating %s of:\n", ModeString[mode_]);
+      for (unsigned int ii = 0; ii < vinfo1_.size(); ii++)
+        mprintf("\t  %s and %s\n", vinfo1_[ii]->legend(), vinfo2_[ii]->legend());
+    }
+  } else {
+    if (vinfo1_.size() == 1)
+      mprintf(" Calculating %s of vector %s\n", ModeString[mode_], vinfo1_[0]->legend());
+    else {
+      mprintf(" Calculating %s of:\n", ModeString[mode_]);
+      for (DVarray::const_iterator it = vinfo1_.begin(); it != vinfo1_.end(); ++it)
+        mprintf("\t  %s\n", (*it)->legend());
+    }
   }
   if (norm_) mprintf("\tVectors will be normalized.\n");
   if (outfile != 0)
@@ -103,6 +132,19 @@ Analysis::RetType Analysis_VectorMath::Setup(ArgList& analyzeArgs, AnalysisSetup
   return Analysis::OK;
 }
 
+/** Calculate the magnitude of each vector. */
+int Analysis_VectorMath::Magnitude(DataSet* Dout, DataSet_Vector const& V1)
+const
+{
+  DataSet_double& Out = static_cast<DataSet_double&>( *Dout );
+  Out.Resize( V1.Size() );
+  for (unsigned int idx = 0; idx < V1.Size(); idx++)
+  {
+    Out[idx] = sqrt( V1[idx].Magnitude2() );
+  }
+  return 0;
+}
+
 // Analysis_VectorMath::DotProduct()
 int Analysis_VectorMath::DotProduct(DataSet* Dout, DataSet_Vector& V1, DataSet_Vector& V2,
                                     unsigned int vmax, unsigned int v1inc, unsigned int v2inc)
@@ -182,10 +224,30 @@ const
 
 // Analysis_VectorMath::Analyze()
 Analysis::RetType Analysis_VectorMath::Analyze() {
-  for (unsigned int ii = 0; ii < vinfo1_.size(); ii++)
-  {
-    int err = DoMath( DataOut_[ii], *(vinfo1_[ii]), *(vinfo2_[ii]) );
-    if (err != 0) return Analysis::ERR;
+  if (vinfo2_.empty()) {
+    // Single set ops
+    for (unsigned int ii = 0; ii < vinfo1_.size(); ii++)
+    {
+      if (vinfo1_[ii]->Size() < 1) {
+        mprintf("Warning: Vector set '%s' is empty.\n", vinfo1_[ii]->legend());
+        continue;
+      }
+      int err = 1;
+      if (mode_ == MAGNITUDE)
+        err = Magnitude( DataOut_[ii], *(vinfo1_[ii]) );
+      if (err != 0) {
+        mprinterr("Error: A problem occurred when performing vector math for set '%s'\n",
+                  vinfo1_[ii]->legend());
+        return Analysis::ERR;
+      }
+    }
+  } else {
+    // Two set ops
+    for (unsigned int ii = 0; ii < vinfo1_.size(); ii++)
+    {
+      int err = DoMath( DataOut_[ii], *(vinfo1_[ii]), *(vinfo2_[ii]) );
+      if (err != 0) return Analysis::ERR;
+    }
   }
   return Analysis::OK;
 }
diff --git a/src/Analysis_VectorMath.h b/src/Analysis_VectorMath.h
@@ -11,9 +11,10 @@ class Analysis_VectorMath : public Analysis {
     Analysis::RetType Setup(ArgList&, AnalysisSetup&, int);
     Analysis::RetType Analyze();
   private:
-    enum ModeType {  DOTPRODUCT = 0, DOTANGLE, CROSSPRODUCT };
+    enum ModeType {  DOTPRODUCT = 0, DOTANGLE, CROSSPRODUCT, MAGNITUDE };
     static const char* ModeString[];
 
+    int Magnitude(DataSet*, DataSet_Vector const&) const;
     int DotProduct(DataSet*, DataSet_Vector&, DataSet_Vector&,
                    unsigned int, unsigned int, unsigned int) const;
     int CrossProduct(DataSet*, DataSet_Vector&, DataSet_Vector&,