diff --git a/README.md b/README.md index 2c6a812fbf..404d00b831 100644 --- 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 index 394dc72e2e..6ce3580f4a 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -42042,7 +42042,8 @@ cluster \end_layout \begin_layout LyX-Code - [{dme|rms|srmsd} [mass] [nofit] []] [{euclid|manhattan}] [wgt ] + [{dme|rms|srmsd|qrmsd} [mass] [nofit] []] [{euclid|manhattan}] [wgt + ] \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 + +[] 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 + +\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 index f425bda036..1a4ae98ff1 100644 --- 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 #endif @@ -21,7 +22,7 @@ Analysis_Rms2d::Analysis_Rms2d() : void Analysis_Rms2d::Help() const { mprintf("\t[crdset ] [] [] [out ]\n" - "\t[dme | nofit | srmsd] [mass]\n" + "\t[{dme | nofit | srmsd | qrmsd}] [mass]\n" "\t[reftraj [parm | parmindex <#>] []]\n" "\t[corr ]\n" " Calculate RMSD between all frames in , 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 index 592cb7b199..e9932f11b8 100644 --- 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 index b1a686e9ce..e1aba18ac7 100644 --- 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 index bfa561c35e..9a3432fdbd 100644 --- 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 index 5a2a6b3dc8..f667ce56c4 100644 --- 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] []] [{euclid|manhattan}] [wgt ]"; + "[{dme|rms|srmsd|qrmsd} [mass] [nofit] []] [{euclid|manhattan}] [wgt ]"; /** 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 : diff --git a/src/Cluster/Metric_QuatRMSD.cpp b/src/Cluster/Metric_QuatRMSD.cpp new file mode 100644 index 0000000000..924ec81d85 --- /dev/null +++ b/src/Cluster/Metric_QuatRMSD.cpp @@ -0,0 +1,170 @@ +#include "Metric_QuatRMSD.h" +#include "Centroid_Coord.h" +#include "Cframes.h" +#include "../CpptrajStdio.h" +#include "../QuaternionRMSD.h" + +/** Initialize the metric. */ +int Cpptraj::Cluster::Metric_QuatRMSD::Init(DataSet_Coords* dIn, AtomMask const& maskIn, + bool useMass) +{ + // TODO better error handles + if (dIn == 0) { + mprinterr("Internal Error: Metric_QuatRMSD::Init() called with null data set.\n"); + return 1; + } +# ifdef DEBUG_CLUSTER + mprintf("DEBUG: Init QRMSD metric for '%s', mask '%s', nofit=%i, usemass=%i\n", + dIn->legend(), maskIn.MaskString(), 0, (int)useMass); +# endif + coords_ = dIn; + mask_ = maskIn; +// nofit_ = nofit; + useMass_ = useMass; + + return 0; +} + +/** Set up the metric. */ +int Cpptraj::Cluster::Metric_QuatRMSD::Setup() { + if (coords_->Top().SetupIntegerMask( mask_ )) return 1; + mask_.MaskInfo(); +# ifdef DEBUG_CLUSTER + mprintf("DEBUG: QRMSD metric topology: %s %s %i\n", coords_->legend(), + coords_->Top().c_str(), coords_->Top().Natom()); +# endif + if (frm1_.SetupFrameFromMask(mask_, coords_->Top().Atoms())) return 1; + frm2_ = frm1_; +# ifdef DEBUG_CLUSTER + mprintf("DEBUG: Setup QRMSD metric for %i atoms, %zu frames.\n", frm1_.Natom(), coords_->Size()); +# endif + return 0; +} + +/** \return RMSD between two given frames. */ +double Cpptraj::Cluster::Metric_QuatRMSD::FrameDist(int f1, int f2) { + coords_->GetFrame( f1, frm1_, mask_ ); + coords_->GetFrame( f2, frm2_, mask_ ); +# ifdef DEBUG_CLUSTER + double rms; + frm1_.printAtomCoord(0); + frm2_.printAtomCoord(0); + //if (nofit_) + // rms = frm1_.RMSD_NoFit( frm2_, useMass_ ); + //else + frm2_.CenterOnOrigin(useMass_); + rms = QuaternionRMSD_CenteredRef( frm2_, frm1_, useMass_ ); + mprintf("\tMetric_QuatRMSD::FrameDist(%i, %i)= %g\n", f1, f2, rms); + return rms; +# else + //if (nofit_) + // return frm1_.RMSD_NoFit( frm2_, useMass_ ); + //else + frm2_.CenterOnOrigin(useMass_); + return QuaternionRMSD_CenteredRef( frm2_, frm1_, useMass_ ); +# endif +} + +/** \return RMSD between two given centroids. */ +double Cpptraj::Cluster::Metric_QuatRMSD::CentroidDist(Centroid* c1, Centroid* c2) { + //if (nofit_) + // return ((Centroid_Coord*)c1)->Cframe().RMSD_NoFit( ((Centroid_Coord*)c2)->Cframe(), useMass_ ); + //else // Centroid is already at origin. + return QuaternionRMSD_CenteredRef( ((Centroid_Coord*)c2)->Cframe(), + ((Centroid_Coord*)c1)->Cframe(), useMass_ ); +} + +/** \return RMSD between given frame and centroid. */ +double Cpptraj::Cluster::Metric_QuatRMSD::FrameCentroidDist(int f1, Centroid* c1) { + coords_->GetFrame( f1, frm1_, mask_ ); + //if (nofit_) + // return frm1_.RMSD_NoFit( ((Centroid_Coord*)c1)->Cframe(), useMass_ ); + //else // Centroid is already at origin. + return QuaternionRMSD_CenteredRef( ((Centroid_Coord*)c1)->Cframe(), frm1_, useMass_ ); +} + +/** Compute the centroid (avg) coords for each atom from all frames in this + * cluster. If fitting, RMS fit to centroid as it is being built. + */ +void Cpptraj::Cluster::Metric_QuatRMSD::CalculateCentroid(Centroid* centIn, Cframes const& cframesIn) { + Matrix_3x3 Rot; + Vec3 Trans; + Centroid_Coord* cent = (Centroid_Coord*)centIn; + // Reset atom count for centroid. + cent->Cframe().ClearAtoms(); + for (Cframes_it frm = cframesIn.begin(); frm != cframesIn.end(); ++frm) + { + coords_->GetFrame( *frm, frm1_, mask_ ); + if (cent->Cframe().empty()) { + cent->Cframe() = frm1_; + //if (!nofit_) + cent->Cframe().CenterOnOrigin(useMass_); + } else { + //if (!nofit_) { + QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, Trans, useMass_ ); + //Rot.Print("CalculateCentroid"); // DEBUG + frm1_.InverseRotate( Rot ); + //} + cent->Cframe() += frm1_; + } + } + //mprintf("DEBUG: Metric_QuatRMSD::CalculateCentroid divide by %zu\n", cframesIn.size()); + cent->Cframe().Divide( (double)cframesIn.size() ); + //mprintf("\t\tFirst 3 centroid coords (of %i): %f %f %f\n", cent->Cframe().Natom(), + // cent->cent->Cframe()[0], cent->Cframe()[1],cent->Cframe()[2]); +} + +/** \return Average structure of given frames. */ +Cpptraj::Cluster::Centroid* Cpptraj::Cluster::Metric_QuatRMSD::NewCentroid( Cframes const& cframesIn ) { + // TODO: Incorporate mass? + Centroid_Coord* cent = new Centroid_Coord( mask_.Nselected() ); + CalculateCentroid( cent, cframesIn ); + return cent; +} + +// Subtract Notes +// TODO: Handle single frame +// TODO: Check if frame is in cluster? +void Cpptraj::Cluster::Metric_QuatRMSD::FrameOpCentroid(int frame, Centroid* centIn, double oldSize, + CentOpType OP) +{ + Matrix_3x3 Rot; + Vec3 Trans; + Centroid_Coord* cent = (Centroid_Coord*)centIn; + coords_->GetFrame( frame, frm1_, mask_ ); + //if (!nofit_) { + QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, Trans, useMass_ ); + frm1_.InverseRotate( Rot ); + //} + cent->Cframe().Multiply( oldSize ); + if (OP == ADDFRAME) { + cent->Cframe() += frm1_; + cent->Cframe().Divide( oldSize + 1 ); + } else { // SUBTRACTFRAME + cent->Cframe() -= frm1_; + cent->Cframe().Divide( oldSize - 1 ); + } +} + +/** \return Description of RMS calc. */ +std::string Cpptraj::Cluster::Metric_QuatRMSD::Description() const { + std::string description("qrmsd " + mask_.MaskExpression()); + //if (nofit_) description.append(" nofit"); + if (useMass_) description.append(" mass"); + return description; +} + +void Cpptraj::Cluster::Metric_QuatRMSD::Info() const { + mprintf("\tMetric: Quaternion RMSD"); + if (mask_.MaskExpression() == "*") + mprintf(" (all atoms)"); + else + mprintf(" (mask '%s')", mask_.MaskString()); + if (useMass_) + mprintf(", mass-weighted"); + //if (nofit_) + // mprintf(", no fitting"); + //else + mprintf(" best-fit"); + mprintf("\n"); +} diff --git a/src/Cluster/Metric_QuatRMSD.h b/src/Cluster/Metric_QuatRMSD.h new file mode 100644 index 0000000000..cf959d08e4 --- /dev/null +++ b/src/Cluster/Metric_QuatRMSD.h @@ -0,0 +1,42 @@ +#ifndef INC_CLUSTER_METRIC_QUATRMSD_H +#define INC_CLUSTER_METRIC_QUATRMSD_H +#include "Metric.h" +#include "../AtomMask.h" +#include "../DataSet_Coords.h" +#include "../Frame.h" +namespace Cpptraj { +namespace Cluster { + +/// RMS cluster distance calc for Coords DataSet +class Metric_QuatRMSD : public Metric { + public: + Metric_QuatRMSD() : Metric(QRMSD), coords_(0), useMass_(false) {} + int Setup(); + double FrameDist(int, int); + double CentroidDist( Centroid*, Centroid* ); + double FrameCentroidDist(int, Centroid*); + void CalculateCentroid(Centroid*, Cframes const&); + Centroid* NewCentroid(Cframes const&); + Metric* Copy() { return new Metric_QuatRMSD( *this ); } + void FrameOpCentroid(int, Centroid*, double, CentOpType); + std::string Description() const; + void Info() const; + unsigned int Ntotal() const { return (unsigned int)coords_->Size(); } + // ------------------------------------------- + int Init(DataSet_Coords*, AtomMask const&, bool); + /// \return whether RMS is mass-weighted + bool UseMass() const { return useMass_; } + /// \return Atom mask + AtomMask const& Mask() const { return mask_; } + private: + DataSet_Coords* coords_; + AtomMask mask_; + //bool nofit_; + bool useMass_; + Frame frm1_; + Frame frm2_; +}; + +} +} +#endif diff --git a/src/Cluster/Metric_RMS.cpp b/src/Cluster/Metric_RMS.cpp index 9569248a37..bbf02110e0 100644 --- a/src/Cluster/Metric_RMS.cpp +++ b/src/Cluster/Metric_RMS.cpp @@ -99,6 +99,7 @@ void Cpptraj::Cluster::Metric_RMS::CalculateCentroid(Centroid* centIn, Cframes } else { if (!nofit_) { frm1_.RMSD_CenteredRef( cent->Cframe(), Rot, Trans, useMass_ ); + //Rot.Print("CalculateCentroid"); // DEBUG frm1_.Rotate( Rot ); } cent->Cframe() += frm1_; diff --git a/src/Cluster/clusterdepend b/src/Cluster/clusterdepend index ae2b0552ac..d731f3ce12 100644 --- a/src/Cluster/clusterdepend +++ b/src/Cluster/clusterdepend @@ -14,8 +14,9 @@ DBI.o : DBI.cpp CentroidArray.h Cframes.h DBI.h List.h Metric.h MetricArray.h No DrawGraph.o : DrawGraph.cpp ../AssociatedData.h ../Atom.h ../Constants.h ../CpptrajFile.h ../CpptrajStdio.h ../DataSet.h ../DataSet_1D.h ../Dimension.h ../FileIO.h ../FileName.h ../MetaData.h ../NameType.h ../PDBfile.h ../Parallel.h ../Range.h ../Residue.h ../SymbolExporting.h ../TextFormat.h ../Vec3.h Cframes.h DrawGraph.h Metric.h MetricArray.h DynamicMatrix.o : DynamicMatrix.cpp ../ArrayIterator.h ../CpptrajStdio.h ../Matrix.h DynamicMatrix.h List.o : List.cpp ../AssociatedData.h ../CpptrajFile.h ../CpptrajStdio.h ../DataSet.h ../DataSet_1D.h ../DataSet_integer.h ../Dimension.h ../FileIO.h ../FileName.h ../MetaData.h ../Parallel.h ../ProgressBar.h ../Range.h ../TextFormat.h CentroidArray.h Cframes.h List.h Metric.h MetricArray.h Node.h -MetricArray.o : MetricArray.cpp ../ArgList.h ../ArrayIterator.h ../AssociatedData.h ../Atom.h ../AtomMap.h ../AtomMask.h ../AtomType.h ../BaseIOtype.h ../Box.h ../Cluster/Cframes.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DataFile.h ../DataFileList.h ../DataSet.h ../DataSetList.h ../DataSet_Coords.h ../DataSet_Coords_REF.h ../DataSet_PairwiseCache.h ../Dimension.h ../FileIO.h ../FileName.h ../FileTypes.h ../Frame.h ../Hungarian.h ../MapAtom.h ../MaskToken.h ../Matrix.h ../Matrix_3x3.h ../MetaData.h ../Molecule.h ../NameType.h ../OnlineVarT.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 ../SymmetricRmsdCalc.h ../TextFormat.h ../Timer.h ../Topology.h ../TypeNameHolder.h ../Unit.h ../Vec3.h CentroidArray.h Metric.h MetricArray.h Metric_DME.h Metric_RMS.h Metric_SRMSD.h Metric_Scalar.h Metric_Torsion.h +MetricArray.o : MetricArray.cpp ../ArgList.h ../ArrayIterator.h ../AssociatedData.h ../Atom.h ../AtomMap.h ../AtomMask.h ../AtomType.h ../BaseIOtype.h ../Box.h ../Cluster/Cframes.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DataFile.h ../DataFileList.h ../DataSet.h ../DataSetList.h ../DataSet_Coords.h ../DataSet_Coords_REF.h ../DataSet_PairwiseCache.h ../Dimension.h ../FileIO.h ../FileName.h ../FileTypes.h ../Frame.h ../Hungarian.h ../MapAtom.h ../MaskToken.h ../Matrix.h ../Matrix_3x3.h ../MetaData.h ../Molecule.h ../NameType.h ../OnlineVarT.h ../Parallel.h ../ParameterHolders.h ../ParameterSet.h ../ParameterTypes.h ../ProgressBar.h ../QuaternionRMSD.h ../Range.h ../ReferenceFrame.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../StringRoutines.h ../SymbolExporting.h ../SymmetricRmsdCalc.h ../TextFormat.h ../Timer.h ../Topology.h ../TypeNameHolder.h ../Unit.h ../Vec3.h CentroidArray.h Metric.h MetricArray.h Metric_DME.h Metric_QuatRMSD.h Metric_RMS.h Metric_SRMSD.h Metric_Scalar.h Metric_Torsion.h Metric_DME.o : Metric_DME.cpp ../AssociatedData.h ../Atom.h ../AtomMask.h ../AtomType.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DataSet.h ../DataSet_Coords.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 ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../TextFormat.h ../Topology.h ../TypeNameHolder.h ../Unit.h ../Vec3.h Centroid.h Centroid_Coord.h Cframes.h Metric.h Metric_DME.h +Metric_QuatRMSD.o : Metric_QuatRMSD.cpp ../AssociatedData.h ../Atom.h ../AtomMask.h ../AtomType.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DataSet.h ../DataSet_Coords.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 ../QuaternionRMSD.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../TextFormat.h ../Topology.h ../TypeNameHolder.h ../Unit.h ../Vec3.h Centroid.h Centroid_Coord.h Cframes.h Metric.h Metric_QuatRMSD.h Metric_RMS.o : Metric_RMS.cpp ../AssociatedData.h ../Atom.h ../AtomMask.h ../AtomType.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DataSet.h ../DataSet_Coords.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 ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../TextFormat.h ../Topology.h ../TypeNameHolder.h ../Unit.h ../Vec3.h Centroid.h Centroid_Coord.h Cframes.h Metric.h Metric_RMS.h Metric_SRMSD.o : Metric_SRMSD.cpp ../ArrayIterator.h ../AssociatedData.h ../Atom.h ../AtomMap.h ../AtomMask.h ../AtomType.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DataSet.h ../DataSet_Coords.h ../Dimension.h ../FileIO.h ../FileName.h ../Frame.h ../Hungarian.h ../MapAtom.h ../MaskToken.h ../Matrix.h ../Matrix_3x3.h ../MetaData.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterHolders.h ../ParameterSet.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../SymmetricRmsdCalc.h ../TextFormat.h ../Topology.h ../TypeNameHolder.h ../Unit.h ../Vec3.h Centroid.h Centroid_Coord.h Cframes.h Metric.h Metric_SRMSD.h Metric_Scalar.o : Metric_Scalar.cpp ../AssociatedData.h ../CpptrajFile.h ../CpptrajStdio.h ../DataSet.h ../DataSet_1D.h ../Dimension.h ../FileIO.h ../FileName.h ../MetaData.h ../Parallel.h ../Range.h ../TextFormat.h Centroid.h Centroid_Num.h Cframes.h Metric.h Metric_Scalar.h diff --git a/src/Cluster/clusterfiles b/src/Cluster/clusterfiles index 293bad8a4c..f4188b5dcc 100644 --- a/src/Cluster/clusterfiles +++ b/src/Cluster/clusterfiles @@ -18,6 +18,7 @@ CLUSTER_SOURCES= \ List.cpp \ MetricArray.cpp \ Metric_DME.cpp \ + Metric_QuatRMSD.cpp \ Metric_RMS.cpp \ Metric_Scalar.cpp \ Metric_SRMSD.cpp \ diff --git a/src/Frame.h b/src/Frame.h index 84d6fc6880..6e3a8d8836 100644 --- a/src/Frame.h +++ b/src/Frame.h @@ -239,6 +239,8 @@ class Frame { inline void Rotate(Matrix_3x3 const&, AtomMask const&); /// Apply inverse of rotation defined by matrix to all atoms in mask inline void InverseRotate(Matrix_3x3 const&, AtomMask const&); + /// Apply invese of rotation defined by matrix to all coords. + inline void InverseRotate(Matrix_3x3 const&); /// Apply translation followed by rotation followed by second translation inline void Trans_Rot_Trans(Vec3 const&, Matrix_3x3 const&, Vec3 const&); /// Apply translation, rotation, 2nd translation for atoms in mask @@ -523,6 +525,17 @@ void Frame::InverseRotate(Matrix_3x3 const& RotMatrix, AtomMask const& mask) { } } +void Frame::InverseRotate(Matrix_3x3 const& T) { + for (int i = 0; i < ncoord_; i += 3) { + double x = X_[i ]; + double y = X_[i+1]; + double z = X_[i+2]; + X_[i ] = (x*T[0]) + (y*T[3]) + (z*T[6]); + X_[i+1] = (x*T[1]) + (y*T[4]) + (z*T[7]); + X_[i+2] = (x*T[2]) + (y*T[5]) + (z*T[8]); + } +} + void Frame::Trans_Rot_Trans(Vec3 const& t1, Matrix_3x3 const& R, Vec3 const& t2) { for (int i = 0; i < ncoord_; i+=3) { double x = X_[i ] + t1[0]; diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp new file mode 100644 index 0000000000..2d0a28b84b --- /dev/null +++ b/src/QuaternionRMSD.cpp @@ -0,0 +1,135 @@ +#include "QuaternionRMSD.h" +#include "qcprot.h" +#include "Constants.h" +#include "CpptrajStdio.h" + +/** Print citation to stdout. */ +void PrintQrmsdCitation() { + mprintf("# Citation: Theobald, D. L.; \"Rapid calculation of RMSD using a\n" + "# quaternion-based characteristic polynomial.\"\n" + "# Acta Crystallographica A 61(4):478-480.\n" + "# Citation: Liu, P.; Agrafiotis, D. K.; Theobald, D. L.;\n" + "# \"Fast determination of the optimal rotational matrix\n" + "# for macromolecular superpositions.\"\n" + "# Journal of Computational Chemistry 31(7):1561-1563.\n"); +} + +/** Calculate the coordinate covariance matrix and perform the + * quaternion RMSD calculation. + * \param Ref Reference coordinates centered at the origin; will not be modified. + * \param Tgt Target coordinates; will be centered at the origin. + * \param U If specified, will hold rotation vectors in columns. + * \param Trans Will be set to translation from Target to origin. + * \param useMass If true, weight by Target atom masses. + */ +double do_quaternion_rmsd(Frame const& Ref, Frame& Tgt, + double* U, double* Trans, + bool useMass) +{ + // Center Tgt on the origin + Trans[0] = 0; + Trans[1] = 0; + Trans[2] = 0; + int ncoord = Ref.size(); + double total_mass; + if (useMass) { + total_mass = 0.0; + int im = 0; + for (int ix = 0; ix < ncoord; ix += 3, im++) { + double mass = Tgt.Mass(im); + total_mass += mass; + Trans[0] += (Tgt[ix ] * mass); + Trans[1] += (Tgt[ix+1] * mass); + Trans[2] += (Tgt[ix+2] * mass); + } + } else { + total_mass = (double)Ref.Natom(); + int im = 0; + for (int ix = 0; ix < ncoord; ix += 3, im++) { + Trans[0] += Tgt[ix ]; + Trans[1] += Tgt[ix+1]; + Trans[2] += Tgt[ix+2]; + } + } + if (total_mass < Constants::SMALL) { + mprinterr("Error: RMSD: Divide by zero.\n"); + return -1; + } + Trans[0] /= total_mass; + Trans[1] /= total_mass; + Trans[2] /= total_mass; + Trans[0] = -Trans[0]; + Trans[1] = -Trans[1]; + Trans[2] = -Trans[2]; + Tgt.Translate(Trans); + + // Calculate covariance matrix of Coords and Reference (R = Xt * Ref) + // Calculate the Kabsch matrix: R = (rij) = Sum(yni*xnj) + double mwss = 0.0; + Matrix_3x3 rot(0.0); + if (useMass) { + int im = 0; + for (int i = 0; i < ncoord; i += 3, im++) + { + double xt = Tgt[i ]; + double yt = Tgt[i+1]; + double zt = Tgt[i+2]; + double xr = Ref[i ]; + double yr = Ref[i+1]; + double zr = Ref[i+2]; + double atom_mass = Tgt.Mass(im); + mwss += atom_mass * ( (xt*xt)+(yt*yt)+(zt*zt)+(xr*xr)+(yr*yr)+(zr*zr) ); + rot[0] += atom_mass*xt*xr; + rot[1] += atom_mass*xt*yr; + rot[2] += atom_mass*xt*zr; + rot[3] += atom_mass*yt*xr; + rot[4] += atom_mass*yt*yr; + rot[5] += atom_mass*yt*zr; + rot[6] += atom_mass*zt*xr; + rot[7] += atom_mass*zt*yr; + rot[8] += atom_mass*zt*zr; + } + } else { + for (int i = 0; i < ncoord; i += 3) + { + double xt = Tgt[i ]; + double yt = Tgt[i+1]; + double zt = Tgt[i+2]; + double xr = Ref[i ]; + double yr = Ref[i+1]; + double zr = Ref[i+2]; + mwss += ( (xt*xt)+(yt*yt)+(zt*zt)+(xr*xr)+(yr*yr)+(zr*zr) ); + rot[0] += xt*xr; + rot[1] += xt*yr; + rot[2] += xt*zr; + rot[3] += yt*xr; + rot[4] += yt*yr; + rot[5] += yt*zr; + rot[6] += zt*xr; + rot[7] += zt*yr; + rot[8] += zt*zr; + } + } + mwss *= 0.5; // E0 = 0.5*Sum(xn^2+yn^2) + + double rmsd; + //int err = + Cpptraj::QCPRot::FastCalcRMSDAndRotation(U, rot.Dptr(), &rmsd, mwss, total_mass); + //mprintf("DEBUG: qcprot returned %i\n", err); + return rmsd; +} + +/** Calculate quaternion RMSD with translation vector and rotation matrix. */ +double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, + Matrix_3x3& U, Vec3& Trans, + bool useMass) +{ + return do_quaternion_rmsd(Ref, Tgt, U.Dptr(), Trans.Dptr(), useMass); +} + +/** Calculate quaternion RMSD only. */ +double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, bool useMass) +{ + Vec3 Trans; + return do_quaternion_rmsd(Ref, Tgt, 0, Trans.Dptr(), useMass); +} diff --git a/src/QuaternionRMSD.h b/src/QuaternionRMSD.h new file mode 100644 index 0000000000..83170b3e91 --- /dev/null +++ b/src/QuaternionRMSD.h @@ -0,0 +1,10 @@ +#ifndef INC_QUATERNIONRMSD_H +#define INC_QUATERNIONRMSD_H +#include "Frame.h" +/// Print quaternion RMSD citation to stdout +void PrintQrmsdCitation(); +/// \return Quaternion RMSD, set rotation matrix (in cols) and target translation vector. +double QuaternionRMSD_CenteredRef(Frame const&, Frame&, Matrix_3x3&, Vec3&, bool); +/// \return Quaternion RMSD +double QuaternionRMSD_CenteredRef(Frame const&, Frame&, bool); +#endif diff --git a/src/Version.h b/src/Version.h index 814a916b93..93052a137a 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V6.12.1" +#define CPPTRAJ_INTERNAL_VERSION "V6.12.2" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index b5309e59f2..8e704c1e10 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -119,7 +119,7 @@ Analysis_Overlap.o : Analysis_Overlap.cpp ActionState.h Analysis.h AnalysisState Analysis_PhiPsi.o : Analysis_PhiPsi.cpp ActionState.h Analysis.h AnalysisState.h Analysis_PhiPsi.h ArgList.h Array1D.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.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 FileIO.h FileName.h FileTypes.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 Analysis_Regression.o : Analysis_Regression.cpp ActionState.h Analysis.h AnalysisState.h Analysis_Regression.h ArgList.h Array1D.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mesh.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.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 Spline.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Analysis_RemLog.o : Analysis_RemLog.cpp ActionState.h Analysis.h AnalysisState.h Analysis_Lifetime.h Analysis_RemLog.h ArgList.h Array1D.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mesh.h DataSet_RemLog.h DataSet_integer.h DataSet_integer_mem.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.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 Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h -Analysis_Rms2d.o : Analysis_Rms2d.cpp ActionState.h Analysis.h AnalysisState.h Analysis_Rms2d.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_2D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_MatrixFlt.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h Hungarian.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.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 SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TypeNameHolder.h Unit.h Vec3.h +Analysis_Rms2d.o : Analysis_Rms2d.cpp ActionState.h Analysis.h AnalysisState.h Analysis_Rms2d.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_2D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_MatrixFlt.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h Hungarian.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h QuaternionRMSD.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TypeNameHolder.h Unit.h Vec3.h Analysis_RmsAvgCorr.o : Analysis_RmsAvgCorr.cpp ActionState.h Analysis.h AnalysisState.h Analysis_RmsAvgCorr.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.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 SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Analysis_Rotdif.o : Analysis_Rotdif.cpp ActionState.h Analysis.h AnalysisState.h Analysis_Rotdif.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h ComplexArray.h Constants.h CoordinateInfo.h Corr.h CpptrajFile.h CpptrajStdio.h CurveFit.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mat3x3.h DataSet_Mesh.h DataSet_Vector.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h PubFFT.h Random.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SimplexMin.h Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Analysis_RunningAvg.o : Analysis_RunningAvg.cpp ActionState.h Analysis.h AnalysisState.h Analysis_RunningAvg.h ArgList.h Array1D.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mesh.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.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 Spline.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h @@ -164,8 +164,9 @@ Cluster/DBI.o : Cluster/DBI.cpp Cluster/CentroidArray.h Cluster/Cframes.h Cluste Cluster/DrawGraph.o : Cluster/DrawGraph.cpp AssociatedData.h Atom.h Cluster/Cframes.h Cluster/DrawGraph.h Cluster/Metric.h Cluster/MetricArray.h Constants.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h Dimension.h FileIO.h FileName.h MetaData.h NameType.h PDBfile.h Parallel.h Range.h Residue.h SymbolExporting.h TextFormat.h Vec3.h Cluster/DynamicMatrix.o : Cluster/DynamicMatrix.cpp ArrayIterator.h Cluster/DynamicMatrix.h CpptrajStdio.h Matrix.h Cluster/List.o : Cluster/List.cpp AssociatedData.h Cluster/CentroidArray.h Cluster/Cframes.h Cluster/List.h Cluster/Metric.h Cluster/MetricArray.h Cluster/Node.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h DataSet_integer.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h ProgressBar.h Range.h TextFormat.h -Cluster/MetricArray.o : Cluster/MetricArray.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cluster/CentroidArray.h Cluster/Cframes.h Cluster/Metric.h Cluster/MetricArray.h Cluster/Metric_DME.h Cluster/Metric_RMS.h Cluster/Metric_SRMSD.h Cluster/Metric_Scalar.h Cluster/Metric_Torsion.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_PairwiseCache.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h Hungarian.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OnlineVarT.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 SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h +Cluster/MetricArray.o : Cluster/MetricArray.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cluster/CentroidArray.h Cluster/Cframes.h Cluster/Metric.h Cluster/MetricArray.h Cluster/Metric_DME.h Cluster/Metric_QuatRMSD.h Cluster/Metric_RMS.h Cluster/Metric_SRMSD.h Cluster/Metric_Scalar.h Cluster/Metric_Torsion.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_PairwiseCache.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h Hungarian.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OnlineVarT.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h QuaternionRMSD.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Cluster/Metric_DME.o : Cluster/Metric_DME.cpp AssociatedData.h Atom.h AtomMask.h AtomType.h Box.h Cluster/Centroid.h Cluster/Centroid_Coord.h Cluster/Cframes.h Cluster/Metric.h Cluster/Metric_DME.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Coords.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 ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Topology.h TypeNameHolder.h Unit.h Vec3.h +Cluster/Metric_QuatRMSD.o : Cluster/Metric_QuatRMSD.cpp AssociatedData.h Atom.h AtomMask.h AtomType.h Box.h Cluster/Centroid.h Cluster/Centroid_Coord.h Cluster/Cframes.h Cluster/Metric.h Cluster/Metric_QuatRMSD.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Coords.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 QuaternionRMSD.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Topology.h TypeNameHolder.h Unit.h Vec3.h Cluster/Metric_RMS.o : Cluster/Metric_RMS.cpp AssociatedData.h Atom.h AtomMask.h AtomType.h Box.h Cluster/Centroid.h Cluster/Centroid_Coord.h Cluster/Cframes.h Cluster/Metric.h Cluster/Metric_RMS.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Coords.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 ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Topology.h TypeNameHolder.h Unit.h Vec3.h Cluster/Metric_SRMSD.o : Cluster/Metric_SRMSD.cpp ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h Box.h Cluster/Centroid.h Cluster/Centroid_Coord.h Cluster/Cframes.h Cluster/Metric.h Cluster/Metric_SRMSD.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Coords.h Dimension.h FileIO.h FileName.h Frame.h Hungarian.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Topology.h TypeNameHolder.h Unit.h Vec3.h Cluster/Metric_Scalar.o : Cluster/Metric_Scalar.cpp AssociatedData.h Cluster/Centroid.h Cluster/Centroid_Num.h Cluster/Cframes.h Cluster/Metric.h Cluster/Metric_Scalar.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h @@ -391,6 +392,7 @@ PubFFT.o : PubFFT.cpp ArrayIterator.h ComplexArray.h CpptrajStdio.h PubFFT.h Pucker_PuckerMask.o : Pucker_PuckerMask.cpp Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Pucker_PuckerMask.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h Pucker_PuckerSearch.o : Pucker_PuckerSearch.cpp ArgList.h Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Pucker.h Pucker_PuckerMask.h Pucker_PuckerSearch.h Pucker_PuckerToken.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h Pucker_PuckerToken.o : Pucker_PuckerToken.cpp Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Pucker_PuckerMask.h Pucker_PuckerToken.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h +QuaternionRMSD.o : QuaternionRMSD.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h QuaternionRMSD.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Unit.h Vec3.h qcprot.h RNG.o : RNG.cpp CpptrajStdio.h RNG.h RNG_Marsaglia.o : RNG_Marsaglia.cpp CpptrajStdio.h RNG.h RNG_Marsaglia.h RNG_MersenneTwister.o : RNG_MersenneTwister.cpp CpptrajStdio.h RNG.h RNG_MersenneTwister.h @@ -456,5 +458,6 @@ Vec3.o : Vec3.cpp Constants.h CpptrajStdio.h Vec3.h ViewRst.o : ViewRst.cpp ActionFrameCounter.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.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 ViewRst.h main.o : main.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 Cpptraj.h CpptrajFile.h CpptrajState.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.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 molsurf.o : molsurf.c molsurf.h +qcprot.o : qcprot.cpp CpptrajStdio.h qcprot.h vmdplugin/dtrplugin.o : vmdplugin/dtrplugin.cpp ByteRoutines.h vmdplugin/dtrplugin.hxx vmdplugin/vmddir.h xoshiro128plusplus.o : xoshiro128plusplus.cpp diff --git a/src/cpptrajfiles b/src/cpptrajfiles index 2ce3ad56b7..8f58c5b8ae 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -359,6 +359,8 @@ COMMON_SOURCES= \ ProgressBar.cpp \ ProgressTimer.cpp \ PubFFT.cpp \ + qcprot.cpp \ + QuaternionRMSD.cpp \ Pucker_PuckerMask.cpp \ Pucker_PuckerSearch.cpp \ Pucker_PuckerToken.cpp \ @@ -427,7 +429,7 @@ COMMON_SOURCES= \ vmdplugin/dtrplugin.cpp \ xoshiro128plusplus.cpp -CSOURCES= molsurf.c +CSOURCES= molsurf.c # The files below are used by cpptraj and libcpptraj but must be compiled # differently for libcpptraj. diff --git a/src/qcprot.cpp b/src/qcprot.cpp new file mode 100644 index 0000000000..72404c5349 --- /dev/null +++ b/src/qcprot.cpp @@ -0,0 +1,427 @@ +/******************************************************************************* + * -/_|:|_|_\- + * + * File: qcprot.c + * Version: 1.5 + * + * Function: Rapid calculation of the least-squares rotation using a + * quaternion-based characteristic polynomial and + * a cofactor matrix + * + * Author(s): Douglas L. Theobald + * Department of Biochemistry + * MS 009 + * Brandeis University + * 415 South St + * Waltham, MA 02453 + * USA + * + * dtheobald@brandeis.edu + * + * Pu Liu + * Johnson & Johnson Pharmaceutical Research and Development, L.L.C. + * 665 Stockton Drive + * Exton, PA 19341 + * USA + * + * pliu24@its.jnj.com + * + * + * If you use this QCP rotation calculation method in a publication, please + * reference: + * + * Douglas L. Theobald (2005) + * "Rapid calculation of RMSD using a quaternion-based characteristic + * polynomial." + * Acta Crystallographica A 61(4):478-480. + * + * Pu Liu, Dmitris K. Agrafiotis, and Douglas L. Theobald (2009) + * "Fast determination of the optimal rotational matrix for macromolecular + * superpositions." + * Journal of Computational Chemistry 31(7):1561-1563. + * + * Copyright (c) 2009-2016 Pu Liu and Douglas L. Theobald + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, are permitted + * provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, this list of + * conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, this list + * of conditions and the following disclaimer in the documentation and/or other materials + * provided with the distribution. + * * Neither the name of the nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific prior written + * permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * Source: started anew. + * + * Change History: + * 2009/04/13 Started source + * 2010/03/28 Modified FastCalcRMSDAndRotation() to handle tiny qsqr + * If trying all rows of the adjoint still gives too small + * qsqr, then just return identity matrix. (DLT) + * 2010/06/30 Fixed prob in assigning A[9] = 0 in InnerProduct() + * invalid mem access + * 2011/02/21 Made CenterCoords use weights + * 2011/05/02 Finally changed CenterCoords declaration in qcprot.h + * Also changed some functions to static + * 2011/07/08 Put in fabs() to fix taking sqrt of small neg numbers, fp error + * 2012/07/26 Minor changes to comments and main.c, more info (v.1.4) + * 2016/07/13 Fixed normalization of RMSD in FastCalcRMSDAndRotation(), should divide by + * sum of weights (thanks to Geoff Skillman) + * + ******************************************************************************/ +#include "qcprot.h" +#include +#include "CpptrajStdio.h" +/* +static double +InnerProduct(double *A, double **coords1, double **coords2, const int len, const double *weight) +{ + double x1, x2, y1, y2, z1, z2; + int i; + const double *fx1 = coords1[0], *fy1 = coords1[1], *fz1 = coords1[2]; + const double *fx2 = coords2[0], *fy2 = coords2[1], *fz2 = coords2[2]; + double G1 = 0.0, G2 = 0.0; + + A[0] = A[1] = A[2] = A[3] = A[4] = A[5] = A[6] = A[7] = A[8] = 0.0; + + if (weight != NULL) + { + for (i = 0; i < len; ++i) + { + x1 = weight[i] * fx1[i]; + y1 = weight[i] * fy1[i]; + z1 = weight[i] * fz1[i]; + + G1 += x1 * fx1[i] + y1 * fy1[i] + z1 * fz1[i]; + + x2 = fx2[i]; + y2 = fy2[i]; + z2 = fz2[i]; + + G2 += weight[i] * (x2 * x2 + y2 * y2 + z2 * z2); + + A[0] += (x1 * x2); + A[1] += (x1 * y2); + A[2] += (x1 * z2); + + A[3] += (y1 * x2); + A[4] += (y1 * y2); + A[5] += (y1 * z2); + + A[6] += (z1 * x2); + A[7] += (z1 * y2); + A[8] += (z1 * z2); + } + } + else + { + for (i = 0; i < len; ++i) + { + x1 = fx1[i]; + y1 = fy1[i]; + z1 = fz1[i]; + + G1 += x1 * x1 + y1 * y1 + z1 * z1; + + x2 = fx2[i]; + y2 = fy2[i]; + z2 = fz2[i]; + + G2 += (x2 * x2 + y2 * y2 + z2 * z2); + + A[0] += (x1 * x2); + A[1] += (x1 * y2); + A[2] += (x1 * z2); + + A[3] += (y1 * x2); + A[4] += (y1 * y2); + A[5] += (y1 * z2); + + A[6] += (z1 * x2); + A[7] += (z1 * y2); + A[8] += (z1 * z2); + } + } + + return (G1 + G2) * 0.5; +} +*/ + +/// Calculate quaternion RMSD. +int Cpptraj::QCPRot::FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double len) +{ + double Sxx, Sxy, Sxz, Syx, Syy, Syz, Szx, Szy, Szz; + double Szz2, Syy2, Sxx2, Sxy2, Syz2, Sxz2, Syx2, Szy2, Szx2, + SyzSzymSyySzz2, Sxx2Syy2Szz2Syz2Szy2, Sxy2Sxz2Syx2Szx2, + SxzpSzx, SyzpSzy, SxypSyx, SyzmSzy, + SxzmSzx, SxymSyx, SxxpSyy, SxxmSyy; + double C[4]; + int i; + double mxEigenV; + double oldg = 0.0; + double b, a, delta, rms, qsqr; + double q1, q2, q3, q4, normq; + double a11, a12, a13, a14, a21, a22, a23, a24; + double a31, a32, a33, a34, a41, a42, a43, a44; + double a2, x2, y2, z2; + double xy, az, zx, ay, yz, ax; + double a3344_4334, a3244_4234, a3243_4233, a3143_4133,a3144_4134, a3142_4132; + double evecprec = 1e-6; + double evalprec = 1e-11; + + Sxx = A[0]; Sxy = A[1]; Sxz = A[2]; + Syx = A[3]; Syy = A[4]; Syz = A[5]; + Szx = A[6]; Szy = A[7]; Szz = A[8]; + + Sxx2 = Sxx * Sxx; + Syy2 = Syy * Syy; + Szz2 = Szz * Szz; + + Sxy2 = Sxy * Sxy; + Syz2 = Syz * Syz; + Sxz2 = Sxz * Sxz; + + Syx2 = Syx * Syx; + Szy2 = Szy * Szy; + Szx2 = Szx * Szx; + + SyzSzymSyySzz2 = 2.0*(Syz*Szy - Syy*Szz); + Sxx2Syy2Szz2Syz2Szy2 = Syy2 + Szz2 - Sxx2 + Syz2 + Szy2; + + C[2] = -2.0 * (Sxx2 + Syy2 + Szz2 + Sxy2 + Syx2 + Sxz2 + Szx2 + Syz2 + Szy2); + C[1] = 8.0 * (Sxx*Syz*Szy + Syy*Szx*Sxz + Szz*Sxy*Syx - Sxx*Syy*Szz - Syz*Szx*Sxy - Szy*Syx*Sxz); + + SxzpSzx = Sxz + Szx; + SyzpSzy = Syz + Szy; + SxypSyx = Sxy + Syx; + SyzmSzy = Syz - Szy; + SxzmSzx = Sxz - Szx; + SxymSyx = Sxy - Syx; + SxxpSyy = Sxx + Syy; + SxxmSyy = Sxx - Syy; + Sxy2Sxz2Syx2Szx2 = Sxy2 + Sxz2 - Syx2 - Szx2; + + C[0] = Sxy2Sxz2Syx2Szx2 * Sxy2Sxz2Syx2Szx2 + + (Sxx2Syy2Szz2Syz2Szy2 + SyzSzymSyySzz2) * (Sxx2Syy2Szz2Syz2Szy2 - SyzSzymSyySzz2) + + (-(SxzpSzx)*(SyzmSzy)+(SxymSyx)*(SxxmSyy-Szz)) * (-(SxzmSzx)*(SyzpSzy)+(SxymSyx)*(SxxmSyy+Szz)) + + (-(SxzpSzx)*(SyzpSzy)-(SxypSyx)*(SxxpSyy-Szz)) * (-(SxzmSzx)*(SyzmSzy)-(SxypSyx)*(SxxpSyy+Szz)) + + (+(SxypSyx)*(SyzpSzy)+(SxzpSzx)*(SxxmSyy+Szz)) * (-(SxymSyx)*(SyzmSzy)+(SxzpSzx)*(SxxpSyy+Szz)) + + (+(SxypSyx)*(SyzmSzy)+(SxzmSzx)*(SxxmSyy-Szz)) * (-(SxymSyx)*(SyzpSzy)+(SxzmSzx)*(SxxpSyy-Szz)); + + /* Newton-Raphson */ + mxEigenV = E0; + for (i = 0; i < 50; ++i) + { + oldg = mxEigenV; + x2 = mxEigenV*mxEigenV; + b = (x2 + C[2])*mxEigenV; + a = b + C[1]; + delta = ((a*mxEigenV + C[0])/(2.0*x2*mxEigenV + b + a)); + mxEigenV -= delta; + /* printf("\n diff[%3d]: %16g %16g %16g", i, mxEigenV - oldg, evalprec*mxEigenV, mxEigenV); */ + if (fabs(mxEigenV - oldg) < fabs(evalprec*mxEigenV)) + break; + } + + if (i == 50) + mprinterr("\nMore than %d iterations needed!\n", i); + + /* the fabs() is to guard against extremely small, but *negative* numbers due to floating point error */ + rms = sqrt(fabs(2.0 * (E0 - mxEigenV)/len)); + (*rmsd) = rms; + /* printf("\n\n %16g %16g %16g \n", rms, E0, 2.0 * (E0 - mxEigenV)/len); */ + if (rot == 0) return (-1); + +// if (minScore > 0) +// if (rms < minScore) +// return (-1); // Don't bother with rotation. + + a11 = SxxpSyy + Szz-mxEigenV; a12 = SyzmSzy; a13 = - SxzmSzx; a14 = SxymSyx; + a21 = SyzmSzy; a22 = SxxmSyy - Szz-mxEigenV; a23 = SxypSyx; a24= SxzpSzx; + a31 = a13; a32 = a23; a33 = Syy-Sxx-Szz - mxEigenV; a34 = SyzpSzy; + a41 = a14; a42 = a24; a43 = a34; a44 = Szz - SxxpSyy - mxEigenV; + a3344_4334 = a33 * a44 - a43 * a34; a3244_4234 = a32 * a44-a42*a34; + a3243_4233 = a32 * a43 - a42 * a33; a3143_4133 = a31 * a43-a41*a33; + a3144_4134 = a31 * a44 - a41 * a34; a3142_4132 = a31 * a42-a41*a32; + q1 = a22*a3344_4334-a23*a3244_4234+a24*a3243_4233; + q2 = -a21*a3344_4334+a23*a3144_4134-a24*a3143_4133; + q3 = a21*a3244_4234-a22*a3144_4134+a24*a3142_4132; + q4 = -a21*a3243_4233+a22*a3143_4133-a23*a3142_4132; + + qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4; + +/* The following code tries to calculate another column in the adjoint matrix when the norm of the + current column is too small. + Usually this block will never be activated. To be absolutely safe this should be + uncommented, but it is most likely unnecessary. +*/ + if (qsqr < evecprec) + { + q1 = a12*a3344_4334 - a13*a3244_4234 + a14*a3243_4233; + q2 = -a11*a3344_4334 + a13*a3144_4134 - a14*a3143_4133; + q3 = a11*a3244_4234 - a12*a3144_4134 + a14*a3142_4132; + q4 = -a11*a3243_4233 + a12*a3143_4133 - a13*a3142_4132; + qsqr = q1*q1 + q2 *q2 + q3*q3+q4*q4; + + if (qsqr < evecprec) + { + double a1324_1423 = a13 * a24 - a14 * a23, a1224_1422 = a12 * a24 - a14 * a22; + double a1223_1322 = a12 * a23 - a13 * a22, a1124_1421 = a11 * a24 - a14 * a21; + double a1123_1321 = a11 * a23 - a13 * a21, a1122_1221 = a11 * a22 - a12 * a21; + + q1 = a42 * a1324_1423 - a43 * a1224_1422 + a44 * a1223_1322; + q2 = -a41 * a1324_1423 + a43 * a1124_1421 - a44 * a1123_1321; + q3 = a41 * a1224_1422 - a42 * a1124_1421 + a44 * a1122_1221; + q4 = -a41 * a1223_1322 + a42 * a1123_1321 - a43 * a1122_1221; + qsqr = q1*q1 + q2 *q2 + q3*q3+q4*q4; + + if (qsqr < evecprec) + { + q1 = a32 * a1324_1423 - a33 * a1224_1422 + a34 * a1223_1322; + q2 = -a31 * a1324_1423 + a33 * a1124_1421 - a34 * a1123_1321; + q3 = a31 * a1224_1422 - a32 * a1124_1421 + a34 * a1122_1221; + q4 = -a31 * a1223_1322 + a32 * a1123_1321 - a33 * a1122_1221; + qsqr = q1*q1 + q2 *q2 + q3*q3 + q4*q4; + + if (qsqr < evecprec) + { + /* if qsqr is still too small, return the identity matrix. */ + rot[0] = rot[4] = rot[8] = 1.0; + rot[1] = rot[2] = rot[3] = rot[5] = rot[6] = rot[7] = 0.0; + + return(0); + } + } + } + } + + normq = sqrt(qsqr); + q1 /= normq; + q2 /= normq; + q3 /= normq; + q4 /= normq; + + a2 = q1 * q1; + x2 = q2 * q2; + y2 = q3 * q3; + z2 = q4 * q4; + + xy = q2 * q3; + az = q1 * q4; + zx = q4 * q2; + ay = q1 * q3; + yz = q3 * q4; + ax = q1 * q2; + + rot[0] = a2 + x2 - y2 - z2; + rot[1] = 2 * (xy + az); + rot[2] = 2 * (zx - ay); + rot[3] = 2 * (xy - az); + rot[4] = a2 - x2 + y2 - z2; + rot[5] = 2 * (yz + ax); + rot[6] = 2 * (zx + ay); + rot[7] = 2 * (yz - ax); + rot[8] = a2 - x2 - y2 + z2; + + return (1); +} + +/* +static void +CenterCoords(double **coords, const int len, const double *weight) +{ + int i; + double xsum, ysum, zsum, wsum; + double *x = coords[0], *y = coords[1], *z = coords[2]; + + xsum = ysum = zsum = 0.0; + + if (weight != NULL) + { + wsum = 0.0; + for (i = 0; i < len; ++i) + { + xsum += weight[i] * x[i]; + ysum += weight[i] * y[i]; + zsum += weight[i] * z[i]; + + wsum += weight[i]; + } + + xsum /= wsum; + ysum /= wsum; + zsum /= wsum; + } + else + { + for (i = 0; i < len; ++i) + { + xsum += x[i]; + ysum += y[i]; + zsum += z[i]; + } + + xsum /= len; + ysum /= len; + zsum /= len; + } + + for (i = 0; i < len; ++i) + { + x[i] -= xsum; + y[i] -= ysum; + z[i] -= zsum; + } +} +*/ + +/* Superposition coords2 onto coords1 -- in other words, coords2 is rotated, coords1 is held fixed */ +/* +double +CalcRMSDRotationalMatrix(double **coords1, double **coords2, const int len, double *rot, const double *weight) +{ + int i; + double A[9], rmsd, wsum; + + // center the structures -- if precentered you can omit this step + CenterCoords(coords1, len, weight); + CenterCoords(coords2, len, weight); + + if (weight == NULL) + { + wsum = len; + } + else + { + wsum = 0.0; + for (i = 0; i < len; ++i) + { + wsum += weight[i]; + } + } + + // calculate the (weighted) inner product of two structures + double E0 = InnerProduct(A, coords1, coords2, len, weight); + + // calculate the RMSD & rotational matrix + FastCalcRMSDAndRotation(rot, A, &rmsd, E0, wsum, -1); + + return rmsd; +} +*/ diff --git a/src/qcprot.h b/src/qcprot.h new file mode 100644 index 0000000000..ed435c6be0 --- /dev/null +++ b/src/qcprot.h @@ -0,0 +1,151 @@ +#ifndef INC_QCPROT_H +#define INC_QCPROT_H +namespace Cpptraj { +namespace QCPRot { +/******************************************************************************* + * -/_|:|_|_\- + * + * File: qcprot.h + * Version: 1.5 + * + * Function: Rapid calculation of the least-squares rotation using a + * quaternion-based characteristic polynomial and + * a cofactor matrix + * + * Author(s): Douglas L. Theobald + * Department of Biochemistry + * MS 009 + * Brandeis University + * 415 South St + * Waltham, MA 02453 + * USA + * + * dtheobald@brandeis.edu + * + * Pu Liu + * Johnson & Johnson Pharmaceutical Research and Development, L.L.C. + * 665 Stockton Drive + * Exton, PA 19341 + * USA + * + * pliu24@its.jnj.com + * + * + * If you use this QCP rotation calculation method in a publication, please + * reference: + * + * Douglas L. Theobald (2005) + * "Rapid calculation of RMSD using a quaternion-based characteristic + * polynomial." + * Acta Crystallographica A 61(4):478-480. + * + * Pu Liu, Dmitris K. Agrafiotis, and Douglas L. Theobald (2009) + * "Fast determination of the optimal rotational matrix for macromolecular + * superpositions." + * Journal of Computational Chemistry 31(7):1561-1563. + * + * Copyright (c) 2009-2016 Pu Liu and Douglas L. Theobald + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, are permitted + * provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, this list of + * conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, this list + * of conditions and the following disclaimer in the documentation and/or other materials + * provided with the distribution. + * * Neither the name of the nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific prior written + * permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * Source: started anew. + * + * Change History: + * 2009/04/13 Started source + * + ******************************************************************************/ + +/* Calculate the RMSD & rotational matrix. + + Input: + coords1 -- reference structure + coords2 -- candidate structure + len -- the size of the system + weight -- the weight array of size len; set to NULL if not needed + Output: + rot[9] -- rotation matrix + Return: + RMSD value + +*/ +//double CalcRMSDRotationalMatrix(double **coords1, double **coords2, const int len, double *rot, const double *weight); + +/* Calculate the inner product of two structures. + If weight array is not NULL, calculate the weighted inner product. + + Input: + coords1 -- reference structure + coords2 -- candidate structure + len -- the size of the system + weight -- the weight array of size len: set to NULL if not needed + Output: + A[9] -- the inner product matrix + Return: + (G1 + G2) * 0.5; used as E0 in function 'FastCalcRMSDAndRotation' + + Warning: + 1. You MUST center the structures, coords1 and coords2, before calling this function. + + 2. Please note how the structure coordinates are stored in the double **coords + arrays. They are 3xN arrays, not Nx3 arrays as is also commonly + used (where the x, y, z axes are interleaved). The difference is + something like this for storage of a structure with 8 atoms: + + Nx3: xyzxyzxyzxyzxyzxyzxyzxyz + 3xN: xxxxxxxxyyyyyyyyzzzzzzzz + + The functions can be easily modified, however, to accomodate any + data format preference. I chose this format because it is readily + used in vectorized functions (SIMD, Altivec, MMX, SSE2, etc.). +*/ +//double InnerProduct(double *A, double **coords1, double **coords2, const int len, const double *weight); + +/** Calculate the RMSD, and/or the optimal rotation matrix. + * \param rot -- Output rotation matrix (order of xx, xy, xz, yx, yy, yz, zx, zy, zz). + * If null, do not calculate rotation matrix. + * \param A[9] -- the inner product (coordinate covariance matrix) of two structures. + * \param rmsd -- Output RMSD value. + * \param E0 -- (G1 + G2) * 0.5 + * \param len -- The size of the system (total mass if mass-weighted, # atoms otherwise). + * \param minScore -- if( minScore > 0 && rmsd < minScore) then calculate only the rmsd; + otherwise, calculate both the RMSD & the rotation matrix. + * \return -1 if only the RMSD was calculated. + * \return >= 0 If rotation matrix and RMSD were calculated. + */ +int FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double len); + +/* Center the coordinates. + + Warning: + If you are doing a full superposition (the usual least squares way), + you MUST center each structure first. That is, you must translate + each structure so that its centroid is at the origin. + You can use CenterCoords() for this. +*/ +//void CenterCoords(double **coords, const int len); +} +} +#endif diff --git a/test/Makefile b/test/Makefile index 1a2c3b3de4..ba6266e4ef 100644 --- a/test/Makefile +++ b/test/Makefile @@ -113,6 +113,7 @@ test.cluster: @-cd Test_Cluster_Cache && ./RunTest.sh $(OPT) @-cd Test_Cluster_ReadInfo && ./RunTest.sh $(OPT) @-cd Test_Cluster_CoordsAndData && ./RunTest.sh $(OPT) + @-cd Test_Cluster_QuatRMSD && ./RunTest.sh $(OPT) test.ired: @-cd Test_IRED && ./RunTest.sh $(OPT) diff --git a/test/Test_2DRMS/RunTest.sh b/test/Test_2DRMS/RunTest.sh index b1744de3b7..69c8c6373d 100755 --- a/test/Test_2DRMS/RunTest.sh +++ b/test/Test_2DRMS/RunTest.sh @@ -4,7 +4,7 @@ # Clean CleanFiles rms.in rmsd1.dat rmsd2.dat ref.nc rmsd.mass.dat dme.dat trp.dat \ - nofit.dat rgacc.dat + nofit.dat rgacc.dat quat.rmsd2.dat TESTNAME='2D RMSD tests' Requires netcdf maxthreads 10 @@ -87,6 +87,16 @@ EOF RunCpptraj "2D RMSD test, read coords with velocity info" DoTest rgacc.dat.save rgacc.dat +# Test 8 - 2drms with quaternion +cat > rms.in < cluster.in < cluster.in <