Amber-MD · drroe · Jul 29, 2022 · Apr 5, 2018 · Apr 6, 2018 · Apr 6, 2018
diff --git a/README.md b/README.md
@@ -249,3 +249,5 @@ External code/libraries bundled with CPPTRAJ
 * Support for reading DTR trajectories uses the VMD DTR plugin.
 
 * CPPTRAJ uses code for the [permuted congruent pseudo-random number generator](https://www.pcg-random.org/index.html) PCG32 by Melissa O'Neill and the [Xoshiro 128++ pseudo-random number generator](http://prng.di.unimi.it) by David Blackman and Sebastino Vigna.
+
+* The code for quaternion RMSD calculation was adapted from code in [qcprot.c](https://theobald.brandeis.edu/qcp/qcprot.c) originally written by Douglas L. Theobald and Pu Lio (Brandeis University).
diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx
@@ -42042,7 +42042,8 @@ cluster
 \end_layout
 
 \begin_layout LyX-Code
-	[{dme|rms|srmsd} [mass] [nofit] [<mask>]] [{euclid|manhattan}] [wgt <list>]
+	[{dme|rms|srmsd|qrmsd} [mass] [nofit] [<mask>]] [{euclid|manhattan}] [wgt
+ <list>]
 \end_layout
 
 \begin_layout LyX-Code
@@ -42580,6 +42581,26 @@ reference "subsec:cpptraj_symmrmsd"
 [nofit] Do not fit structures onto each other prior to calculating RMSD.
 \end_layout
 
+\end_deeper
+\begin_layout Description
+qrmsd
+\begin_inset space ~
+\end_inset
+
+[<mask>] For COORDS data, distance between coordinate frames calculated
+ using best-fit quaternion RMSD (can be 15-20% faster than regular RMSD)
+ using atoms in 
+\series bold
+<mask>
+\series default
+.
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+[mass] Mass-weight the RMSD.
+\end_layout
+
 \end_deeper
 \begin_layout Description
 dme
@@ -48322,7 +48343,7 @@ rms2d
 \end_layout
 
 \begin_layout LyX-Code
-      [dme | nofit | srmsd] [mass]
+      [{dme | nofit | srmsd | qrmsd}] [mass]
 \end_layout
 
 \begin_layout LyX-Code
@@ -48384,6 +48405,10 @@ reference "subsec:cpptraj_symmrmsd"
 ).
 \end_layout
 
+\begin_layout Description
+[qrmsd] Use quaternion RMSD calculation (can be 15-20% faster).
+\end_layout
+
 \begin_layout Description
 [mass] Mass-weight RMSD.
 \end_layout

diff --git a/src/Analysis_Rms2d.cpp b/src/Analysis_Rms2d.cpp
@@ -3,6 +3,7 @@
 #include "CpptrajStdio.h"
 #include "ProgressBar.h"
 #include "DataSet_Coords_TRJ.h"
+#include "QuaternionRMSD.h"
 #ifdef _OPENMP
 #  include <omp.h>
 #endif
@@ -21,7 +22,7 @@ Analysis_Rms2d::Analysis_Rms2d() :
 
 void Analysis_Rms2d::Help() const {
   mprintf("\t[crdset <crd set>] [<name>] [<mask>] [out <filename>]\n"
-          "\t[dme | nofit | srmsd] [mass]\n"
+          "\t[{dme | nofit | srmsd | qrmsd}] [mass]\n"
           "\t[reftraj <traj> [parm <parmname> | parmindex <#>] [<refmask>]]\n"
           "\t[corr <corrfilename>]\n"
           "  Calculate RMSD between all frames in <crd set>, or between frames in\n"
@@ -30,7 +31,7 @@ void Analysis_Rms2d::Help() const {
 }
 
 const char* Analysis_Rms2d::ModeStrings_[] = {
-  "RMSD (fit)", "RMSD (no fitting)", "DME", "symmetry-corrected RMSD"
+  "RMSD (fit)", "RMSD (no fitting)", "DME", "symmetry-corrected RMSD", "quaternion RMSD"
 };
 
 // Analysis_Rms2d::Setup()
@@ -52,6 +53,8 @@ Analysis::RetType Analysis_Rms2d::Setup(ArgList& analyzeArgs, AnalysisSetup& set
     mode_ = DME;
   else if (analyzeArgs.hasKey("srmsd"))
     mode_ = SRMSD;
+  else if (analyzeArgs.hasKey("qrmsd"))
+    mode_ = QUAT;
   else
     mode_ = RMS_FIT;
   useMass_ = analyzeArgs.hasKey("mass");
@@ -138,6 +141,8 @@ Analysis::RetType Analysis_Rms2d::Setup(ArgList& analyzeArgs, AnalysisSetup& set
   if (corrfile != 0)
     mprintf("\tRMSD auto-correlation will be calculated and output to '%s'\n",
             corrfile->DataFilename().full());
+  if (mode_ == QUAT)
+    PrintQrmsdCitation();
 
   return Analysis::OK;
 }
@@ -257,7 +262,7 @@ int Analysis_Rms2d::Calculate_2D() {
 #     endif
         RefTraj_->GetFrame( nref, SelectedRef, RefMask_ );
       // Select and pre-center reference atoms (if fitting)
-      if (mode_ == RMS_FIT || mode_ == SRMSD)
+      if (mode_ == QUAT || mode_ == RMS_FIT || mode_ == SRMSD)
         SelectedRef.CenterOnOrigin(useMass_);
       // LOOP OVER TARGET FRAMES
       if (calculateFullMatrix)
@@ -272,6 +277,8 @@ int Analysis_Rms2d::Calculate_2D() {
           case RMS_NOFIT: R = (float)SelectedTgt.RMSD_NoFit(SelectedRef, useMass_); break;
           case DME:       R = (float)SelectedTgt.DISTRMSD(SelectedRef); break;
           case SRMSD:     R = (float)SRMSD_.SymmRMSD_CenteredRef(SelectedTgt, SelectedRef); break;
+          case QUAT:
+            R = (float)QuaternionRMSD_CenteredRef(SelectedRef, SelectedTgt,  useMass_); break;
         }
         rmsdataset_->SetElement(nref, ntgt, R);
         // DEBUG

diff --git a/src/Analysis_Rms2d.h b/src/Analysis_Rms2d.h
@@ -22,7 +22,7 @@ class Analysis_Rms2d: public Analysis {
     void CalcAutoCorr();
 
     static const char* ModeStrings_[];
-    enum ModeType { RMS_FIT = 0, RMS_NOFIT, DME, SRMSD };
+    enum ModeType { RMS_FIT = 0, RMS_NOFIT, DME, SRMSD, QUAT };
     ModeType mode_;
     DataSet_Coords* TgtTraj_;  ///< Hold coords from input frames.
     bool useReferenceTraj_;    ///< If true read from reference trajectory.

diff --git a/src/Cluster/CMakeLists.txt b/src/Cluster/CMakeLists.txt
@@ -18,6 +18,7 @@ target_sources(cpptraj_common_obj PRIVATE
   ${CMAKE_CURRENT_LIST_DIR}/List.cpp
   ${CMAKE_CURRENT_LIST_DIR}/MetricArray.cpp
   ${CMAKE_CURRENT_LIST_DIR}/Metric_DME.cpp
+  ${CMAKE_CURRENT_LIST_DIR}/Metric_QuatRMSD.cpp
   ${CMAKE_CURRENT_LIST_DIR}/Metric_RMS.cpp
   ${CMAKE_CURRENT_LIST_DIR}/Metric_Scalar.cpp
   ${CMAKE_CURRENT_LIST_DIR}/Metric_SRMSD.cpp

diff --git a/src/Cluster/Metric.h b/src/Cluster/Metric.h
@@ -14,7 +14,7 @@ const int UNCLASSIFIED = -2;
 /// Abstract base class for calculating distance between points or determining centroid.
 class Metric {
   public:
-    enum Type { RMS=0, DME, SRMSD, SCALAR, TORSION, MATRIX2D, UNKNOWN_METRIC };
+    enum Type { RMS=0, DME, SRMSD, SCALAR, TORSION, MATRIX2D, QRMSD, UNKNOWN_METRIC };
     enum CentOpType { ADDFRAME=0, SUBTRACTFRAME };
 
     /// CONSTRUCTOR

diff --git a/src/Cluster/MetricArray.cpp b/src/Cluster/MetricArray.cpp
@@ -17,12 +17,14 @@
 #include "../OnlineVarT.h" // for metric stats average
 #include "../ProgressBar.h"
 #include "../StringRoutines.h"
+#include "../QuaternionRMSD.h" // for printing qrmsd citation
 // Metric classes
 #include "Metric_RMS.h"
 #include "Metric_DME.h"
 #include "Metric_Scalar.h"
 #include "Metric_SRMSD.h"
 #include "Metric_Torsion.h"
+#include "Metric_QuatRMSD.h"
 
 /** CONSTRUCTOR */
 Cpptraj::Cluster::MetricArray::MetricArray() :
@@ -277,7 +279,8 @@ Cpptraj::Cluster::Metric const*
   {
     if ( (*it)->MetricType() == Metric::RMS ||
          (*it)->MetricType() == Metric::DME ||
-         (*it)->MetricType() == Metric::SRMSD )
+         (*it)->MetricType() == Metric::SRMSD ||
+         (*it)->MetricType() == Metric::QRMSD )
       return *it;
   }
   return 0;
@@ -291,6 +294,7 @@ Cpptraj::Cluster::Metric* Cpptraj::Cluster::MetricArray::AllocateMetric(Metric::
     case Metric::RMS       : met = new Metric_RMS(); break;
     case Metric::DME       : met = new Metric_DME(); break;
     case Metric::SRMSD     : met = new Metric_SRMSD(); break;
+    case Metric::QRMSD     : met = new Metric_QuatRMSD(); break;
     case Metric::SCALAR    : met = new Metric_Scalar(); break;
     case Metric::TORSION   : met = new Metric_Torsion(); break;
     default: mprinterr("Error: Unhandled Metric in AllocateMetric.\n");
@@ -300,7 +304,7 @@ Cpptraj::Cluster::Metric* Cpptraj::Cluster::MetricArray::AllocateMetric(Metric::
 
 /** Recognized metric args. */
 const char* Cpptraj::Cluster::MetricArray::MetricArgs_ =
-  "[{dme|rms|srmsd} [mass] [nofit] [<mask>]] [{euclid|manhattan}] [wgt <list>]";
+  "[{dme|rms|srmsd|qrmsd} [mass] [nofit] [<mask>]] [{euclid|manhattan}] [wgt <list>]";
 
 /** Initialize with given sets and arguments. */
 int Cpptraj::Cluster::MetricArray::initMetricArray(DataSetList const& setsToCluster, ArgList& analyzeArgs)
@@ -311,6 +315,7 @@ int Cpptraj::Cluster::MetricArray::initMetricArray(DataSetList const& setsToClus
   int usedme = (int)analyzeArgs.hasKey("dme");
   int userms = (int)analyzeArgs.hasKey("rms");
   int usesrms = (int)analyzeArgs.hasKey("srmsd");
+  int useqrms = (int)analyzeArgs.hasKey("qrmsd");
   bool useMass = analyzeArgs.hasKey("mass");
   bool nofit   = analyzeArgs.hasKey("nofit");
   std::string maskExpr = analyzeArgs.GetMaskNext();
@@ -329,15 +334,18 @@ int Cpptraj::Cluster::MetricArray::initMetricArray(DataSetList const& setsToClus
   std::string wgtArgStr = analyzeArgs.GetStringKey("wgt");
 
   // Check args
-  if (usedme + userms + usesrms > 1) {
-    mprinterr("Error: Specify either 'dme', 'rms', or 'srmsd'.\n");
+  if (usedme + userms + usesrms + useqrms > 1) {
+    mprinterr("Error: Specify either 'dme', 'rms', 'srmsd', or 'qrmsd'.\n");
     return 1;
   }
   Metric::Type coordsMetricType;
   if      (usedme)  coordsMetricType = Metric::DME;
   else if (userms)  coordsMetricType = Metric::RMS;
   else if (usesrms) coordsMetricType = Metric::SRMSD;
-  else coordsMetricType = Metric::RMS; // default
+  else if (useqrms) {
+    coordsMetricType = Metric::QRMSD;
+    PrintQrmsdCitation();
+  } else coordsMetricType = Metric::RMS; // default
 
   // For each input set, set up the appropriate metric
   for (DataSetList::const_iterator ds = setsToCluster.begin(); ds != setsToCluster.end(); ++ds)
@@ -373,6 +381,8 @@ int Cpptraj::Cluster::MetricArray::initMetricArray(DataSetList const& setsToClus
         err = ((Metric_DME*)met)->Init((DataSet_Coords*)*ds, AtomMask(maskExpr)); break;
       case Metric::SRMSD :
         err = ((Metric_SRMSD*)met)->Init((DataSet_Coords*)*ds, AtomMask(maskExpr), nofit, useMass, debug_); break;
+      case Metric::QRMSD :
+        err = ((Metric_QuatRMSD*)met)->Init((DataSet_Coords*)*ds, AtomMask(maskExpr), useMass); break;
       case Metric::SCALAR :
         err = ((Metric_Scalar*)met)->Init((DataSet_1D*)*ds); break;
       case Metric::TORSION :
Original file line number	Diff line number	Diff line change
Expand Up		@@ -249,3 +249,5 @@ External code/libraries bundled with CPPTRAJ
		* Support for reading DTR trajectories uses the VMD DTR plugin.

		* CPPTRAJ uses code for the [permuted congruent pseudo-random number generator](https://www.pcg-random.org/index.html) PCG32 by Melissa O'Neill and the [Xoshiro 128++ pseudo-random number generator](http://prng.di.unimi.it) by David Blackman and Sebastino Vigna.

		* The code for quaternion RMSD calculation was adapted from code in [qcprot.c](https://theobald.brandeis.edu/qcp/qcprot.c) originally written by Douglas L. Theobald and Pu Lio (Brandeis University).