From 94acf04d31e4b3cccd13d38dd89ae67c2611fbaa Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 5 Apr 2018 15:21:28 -0400 Subject: [PATCH 01/38] DRR - Cpptraj: Add some debug info --- src/cuda_kernels/kernel_wrappers.cu | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/cuda_kernels/kernel_wrappers.cu b/src/cuda_kernels/kernel_wrappers.cu index 7e2a9fc156..5ef307a0d6 100644 --- a/src/cuda_kernels/kernel_wrappers.cu +++ b/src/cuda_kernels/kernel_wrappers.cu @@ -75,6 +75,7 @@ void Action_Closest_Center(const double *SolventMols_, double *D_, const double* printf("NMols = %d, NAtoms = %d\n", NMols, NAtoms); printf("active_size = %d\n", active_size); printf("NBlocks = %d\n", NBlocks); + printf("BLOCKDIM = %d\n", BLOCKDIM); printf("sizeof(double) = %d\n", sizeof(double)); printf("About to launch kernel.\n"); @@ -176,6 +177,7 @@ void Action_Closest_NoCenter(const double *SolventMols_, double *D_, const doubl printf("NMols = %d, NAtoms = %d\n", NMols, NAtoms); printf("active_size = %d\n", active_size); printf("NBlocks = %d\n", NBlocks); + printf("BLOCKDIM = %d\n", BLOCKDIM); printf("sizeof(double) = %d\n", sizeof(double)); printf("About to launch kernel.\n"); From c8ca79e57179222f3928c4abda0b49ddb3c7e79e Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 6 Apr 2018 12:26:25 -0400 Subject: [PATCH 02/38] DRR - Cpptraj: Initial test of Quaternion rmsd. Quaternion driver code from DL Theobald (https://theobald.brandeis.edu/qcp/) --- src/QuaternionRMSD.cpp | 94 +++++++++ src/QuaternionRMSD.h | 17 ++ src/cpptrajdepend | 4 +- src/cpptrajfiles | 3 +- src/qcprot.c | 424 +++++++++++++++++++++++++++++++++++++++++ src/qcprot.h | 160 ++++++++++++++++ 6 files changed, 700 insertions(+), 2 deletions(-) create mode 100644 src/QuaternionRMSD.cpp create mode 100644 src/QuaternionRMSD.h create mode 100644 src/qcprot.c create mode 100644 src/qcprot.h diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp new file mode 100644 index 0000000000..6ba0b80adf --- /dev/null +++ b/src/QuaternionRMSD.cpp @@ -0,0 +1,94 @@ +#include "QuaternionRMSD.h" +#include "qcprot.h" +#include "Constants.h" +#include "CpptrajStdio.h" + +QuaternionRMSD::~QuaternionRMSD() { + Clear(); +} + +void QuaternionRMSD::Clear() { + if (M_ != 0) { + delete[] M_; + M_ = 0; + } + if (Xtgt_ != 0) { + for (int i = 0; i < 3; i++) { + delete[] Xtgt_[i]; + delete[] Xref_[i]; + } + delete[] Xtgt_; + delete[] Xref_; + Xtgt_ = 0; + Xref_ = 0; + } + len_ = 0; +} + +int QuaternionRMSD::Init(int natom, std::vector const& mass) +{ + Clear(); + len_ = natom; + Xtgt_ = new double*[ 3 ]; + Xref_ = new double*[ 3 ]; + for (int i = 0; i < 3; i++) { + Xtgt_[i] = new double[ len_ ]; + Xref_[i] = new double[ len_ ]; + } + if (!mass.empty()) { + M_ = new double[ mass.size() ]; + std::copy(mass.begin(), mass.end(), M_); + } + return 0; +} + +double QuaternionRMSD::RMSD_CenteredRef(Frame const& Ref, Frame& Tgt, + Matrix_3x3& U, Vec3& Trans) +{ + Trans.Zero(); + double* X_ = Tgt.xAddress(); + const double* R_ = Ref.xAddress(); + int ncoord = len_ * 3; + double total_mass; + if (M_ != 0) { + total_mass = 0.0; + int im = 0; + for (int ix = 0; ix < ncoord; ix += 3, im++) { + double mass = M_[im]; + total_mass += mass; + Trans[0] += (X_[ix ] * mass); + Trans[1] += (X_[ix+1] * mass); + Trans[2] += (X_[ix+2] * mass); + } + } else { + total_mass = (double)len_; + int im = 0; + for (int ix = 0; ix < ncoord; ix += 3, im++) { + Trans[0] += X_[ix ]; + Trans[1] += X_[ix+1]; + Trans[2] += X_[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.Neg(); + Tgt.Translate(Trans); + // Save coordinates + int ix = 0; + for (int im = 0; im < len_; im++, ix += 3) + { + Xtgt_[0][im] = X_[ix ]; + Xtgt_[1][im] = X_[ix+1]; + Xtgt_[2][im] = X_[ix+2]; + Xref_[0][im] = R_[ix ]; + Xref_[1][im] = R_[ix+1]; + Xref_[2][im] = R_[ix+2]; + } + + return CalcRMSDRotationalMatrix( Xref_, Xtgt_, len_, U.Dptr(), M_ ); +} diff --git a/src/QuaternionRMSD.h b/src/QuaternionRMSD.h new file mode 100644 index 0000000000..5eee03074f --- /dev/null +++ b/src/QuaternionRMSD.h @@ -0,0 +1,17 @@ +#ifndef INC_QUATERNIONRMSD_H +#define INC_QUATERNIONRMSD_H +#include "Frame.h" +class QuaternionRMSD { + public: + QuaternionRMSD() : Xtgt_(0), Xref_(0), M_(0) {} + ~QuaternionRMSD(); + void Clear(); + int Init(int, std::vector const&); + double RMSD_CenteredRef(Frame const&, Frame&, Matrix_3x3&, Vec3&); + private: + double** Xtgt_; + double** Xref_; + double* M_; + int len_; +}; +#endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 718ff8eb34..2ebe00a3a9 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -63,7 +63,7 @@ Action_Radial.o : Action_Radial.cpp Action.h ActionState.h Action_Radial.h ArgLi Action_RandomizeIons.o : Action_RandomizeIons.cpp Action.h ActionState.h Action_RandomizeIons.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Random.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_Remap.o : Action_Remap.cpp Action.h ActionState.h Action_Remap.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h ParmFile.h ParmIO.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_ReplicateCell.o : Action_ReplicateCell.cpp Action.h ActionFrameCounter.h ActionState.h Action_ReplicateCell.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h ParmFile.h ParmIO.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h TrajectoryFile.h TrajectoryIO.h Trajout_Single.h Vec3.h -Action_Rmsd.o : Action_Rmsd.cpp Action.h ActionState.h Action_Rmsd.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h ComplexArray.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.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 ParameterTypes.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h TextFormat.h Timer.h Topology.h Vec3.h +Action_Rmsd.o : Action_Rmsd.cpp Action.h ActionState.h Action_Rmsd.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h ComplexArray.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.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 ParameterTypes.h QuaternionRMSD.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h TextFormat.h Timer.h Topology.h Vec3.h Action_Rotate.o : Action_Rotate.cpp Action.h ActionState.h Action_Rotate.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mat3x3.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_RunningAvg.o : Action_RunningAvg.cpp Action.h ActionState.h Action_RunningAvg.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_STFC_Diffusion.o : Action_STFC_Diffusion.cpp Action.h ActionState.h Action_STFC_Diffusion.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h @@ -306,6 +306,7 @@ Parm_Tinker.o : Parm_Tinker.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOty ProgressBar.o : ProgressBar.cpp CpptrajStdio.h ProgressBar.h ProgressTimer.o : ProgressTimer.cpp CpptrajStdio.h ProgressTimer.h Timer.h PubFFT.o : PubFFT.cpp ArrayIterator.h ComplexArray.h CpptrajStdio.h PubFFT.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 Vec3.h qcprot.h Random.o : Random.cpp CpptrajStdio.h Random.h Range.o : Range.cpp ArgList.h CpptrajStdio.h Range.h RPNcalc.o : RPNcalc.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h Box.h CharMask.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_MatrixDbl.h DataSet_Vector.h DataSet_double.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h RPNcalc.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h @@ -363,4 +364,5 @@ Energy_Sander.o : Energy_Sander.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h Base ReadLine.o : ReadLine.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Cmd.h CmdInput.h CmdList.h Command.h Control.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReadLine.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h VariableArray.h Vec3.h main.o : main.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h Cpptraj.h CpptrajFile.h CpptrajState.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h molsurf.o : molsurf.c molsurf.h +qcprot.o : qcprot.c qcprot.h AmbPDB.o : AmbPDB.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h BufferedFrame.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h ParmFile.h ParmIO.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h TextFormat.h Topology.h TrajFrameCounter.h Traj_AmberRestart.h TrajectoryFile.h TrajectoryIO.h Trajin.h Trajin_Single.h Trajout_Single.h Vec3.h diff --git a/src/cpptrajfiles b/src/cpptrajfiles index e7068479e3..b19c8dea84 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -308,6 +308,7 @@ COMMON_SOURCES=ActionFrameCounter.cpp \ ProgressBar.cpp \ ProgressTimer.cpp \ PubFFT.cpp \ + QuaternionRMSD.cpp \ Random.cpp \ Range.cpp \ RPNcalc.cpp \ @@ -359,7 +360,7 @@ COMMON_SOURCES=ActionFrameCounter.cpp \ Vec3.cpp \ ViewRst.cpp -CSOURCES= molsurf.c +CSOURCES= molsurf.c qcprot.c # The files below are used by cpptraj and libcpptraj but must be compiled # differently for libcpptraj. diff --git a/src/qcprot.c b/src/qcprot.c new file mode 100644 index 0000000000..24687dfaea --- /dev/null +++ b/src/qcprot.c @@ -0,0 +1,424 @@ +/******************************************************************************* + * -/_|:|_|_\- + * + * 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" + + +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; +} + + +int +FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double len, double minScore) +{ + 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) + fprintf(stderr,"\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 (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..2caf976ec7 --- /dev/null +++ b/src/qcprot.h @@ -0,0 +1,160 @@ +#ifndef INC_QCPROT_H +#define INC_QCPROT_H +#ifdef __cplusplus +extern "C" { +#endif +/******************************************************************************* + * -/_|:|_|_\- + * + * 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 + * + ******************************************************************************/ + +#include +#include +#include + +/* 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. + + Input: + A[9] -- the inner product of two structures + E0 -- (G1 + G2) * 0.5 + len -- the size of the system + minScore-- if( minScore > 0 && rmsd < minScore) then calculate only the rmsd; + otherwise, calculate both the RMSD & the rotation matrix + Output: + rot[9] -- the rotation matrix in the order of xx, xy, xz, yx, yy, yz, zx, zy, zz + rmsd -- the RMSD value + Return: + only the rmsd was calculated if < 0 + both the RMSD & rotational matrix calculated if > 0 +*/ +int FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double wsum, double minScore); + +/* 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); +#ifdef __cplusplus +} +#endif +#endif From a11f386cacd96e5dd20e347fc3536aeaf29e882f Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 6 Apr 2018 12:52:19 -0400 Subject: [PATCH 03/38] DRR - Cpptraj: Do the inner product using ordering preferred by cpptraj; eliminates need for re-ordering coordinates. --- src/QuaternionRMSD.cpp | 120 +++++++++++++++++++++++------------------ src/QuaternionRMSD.h | 14 +---- 2 files changed, 68 insertions(+), 66 deletions(-) diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp index 6ba0b80adf..2c6ea2f41a 100644 --- a/src/QuaternionRMSD.cpp +++ b/src/QuaternionRMSD.cpp @@ -3,70 +3,30 @@ #include "Constants.h" #include "CpptrajStdio.h" -QuaternionRMSD::~QuaternionRMSD() { - Clear(); -} - -void QuaternionRMSD::Clear() { - if (M_ != 0) { - delete[] M_; - M_ = 0; - } - if (Xtgt_ != 0) { - for (int i = 0; i < 3; i++) { - delete[] Xtgt_[i]; - delete[] Xref_[i]; - } - delete[] Xtgt_; - delete[] Xref_; - Xtgt_ = 0; - Xref_ = 0; - } - len_ = 0; -} - -int QuaternionRMSD::Init(int natom, std::vector const& mass) -{ - Clear(); - len_ = natom; - Xtgt_ = new double*[ 3 ]; - Xref_ = new double*[ 3 ]; - for (int i = 0; i < 3; i++) { - Xtgt_[i] = new double[ len_ ]; - Xref_[i] = new double[ len_ ]; - } - if (!mass.empty()) { - M_ = new double[ mass.size() ]; - std::copy(mass.begin(), mass.end(), M_); - } - return 0; -} - -double QuaternionRMSD::RMSD_CenteredRef(Frame const& Ref, Frame& Tgt, - Matrix_3x3& U, Vec3& Trans) +double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, + Matrix_3x3& U, Vec3& Trans, + bool useMass) { Trans.Zero(); - double* X_ = Tgt.xAddress(); - const double* R_ = Ref.xAddress(); - int ncoord = len_ * 3; + int ncoord = Ref.size(); double total_mass; - if (M_ != 0) { + if (useMass) { total_mass = 0.0; int im = 0; for (int ix = 0; ix < ncoord; ix += 3, im++) { - double mass = M_[im]; + double mass = Tgt.Mass(im); total_mass += mass; - Trans[0] += (X_[ix ] * mass); - Trans[1] += (X_[ix+1] * mass); - Trans[2] += (X_[ix+2] * mass); + Trans[0] += (Tgt[ix ] * mass); + Trans[1] += (Tgt[ix+1] * mass); + Trans[2] += (Tgt[ix+2] * mass); } } else { - total_mass = (double)len_; + total_mass = (double)Ref.Natom(); int im = 0; for (int ix = 0; ix < ncoord; ix += 3, im++) { - Trans[0] += X_[ix ]; - Trans[1] += X_[ix+1]; - Trans[2] += X_[ix+2]; + Trans[0] += Tgt[ix ]; + Trans[1] += Tgt[ix+1]; + Trans[2] += Tgt[ix+2]; } } if (total_mass < Constants::SMALL) { @@ -78,6 +38,56 @@ double QuaternionRMSD::RMSD_CenteredRef(Frame const& Ref, Frame& Tgt, Trans[2] /= total_mass; Trans.Neg(); 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) +/* // Save coordinates int ix = 0; for (int im = 0; im < len_; im++, ix += 3) @@ -91,4 +101,8 @@ double QuaternionRMSD::RMSD_CenteredRef(Frame const& Ref, Frame& Tgt, } return CalcRMSDRotationalMatrix( Xref_, Xtgt_, len_, U.Dptr(), M_ ); +*/ + double rmsd; + FastCalcRMSDAndRotation(U.Dptr(), rot.Dptr(), &rmsd, mwss, total_mass, -1); + return rmsd; } diff --git a/src/QuaternionRMSD.h b/src/QuaternionRMSD.h index 5eee03074f..7345d77586 100644 --- a/src/QuaternionRMSD.h +++ b/src/QuaternionRMSD.h @@ -1,17 +1,5 @@ #ifndef INC_QUATERNIONRMSD_H #define INC_QUATERNIONRMSD_H #include "Frame.h" -class QuaternionRMSD { - public: - QuaternionRMSD() : Xtgt_(0), Xref_(0), M_(0) {} - ~QuaternionRMSD(); - void Clear(); - int Init(int, std::vector const&); - double RMSD_CenteredRef(Frame const&, Frame&, Matrix_3x3&, Vec3&); - private: - double** Xtgt_; - double** Xref_; - double* M_; - int len_; -}; +double QuaternionRMSD_CenteredRef(Frame const&, Frame&, Matrix_3x3&, Vec3&, bool); #endif From 5629be49cf458e3bcd2a501622ccc70485a65348 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 6 Apr 2018 13:09:15 -0400 Subject: [PATCH 04/38] DRR - Cpptraj: Add version that does not care about rotation matrix --- src/QuaternionRMSD.cpp | 20 +++++++++++++++++--- src/QuaternionRMSD.h | 2 ++ 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp index 2c6ea2f41a..ae9ec207cc 100644 --- a/src/QuaternionRMSD.cpp +++ b/src/QuaternionRMSD.cpp @@ -4,8 +4,8 @@ #include "CpptrajStdio.h" double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, - Matrix_3x3& U, Vec3& Trans, - bool useMass) + Matrix_3x3& U, Vec3& Trans, + bool useMass, double minScore) { Trans.Zero(); int ncoord = Ref.size(); @@ -103,6 +103,20 @@ double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, return CalcRMSDRotationalMatrix( Xref_, Xtgt_, len_, U.Dptr(), M_ ); */ double rmsd; - FastCalcRMSDAndRotation(U.Dptr(), rot.Dptr(), &rmsd, mwss, total_mass, -1); + FastCalcRMSDAndRotation(U.Dptr(), rot.Dptr(), &rmsd, mwss, total_mass, minScore); return rmsd; } + +double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, + Matrix_3x3& U, Vec3& Trans, + bool useMass) +{ + return QuaternionRMSD_CenteredRef(Ref, Tgt, U, Trans, useMass, -1); +} + +double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, bool useMass) +{ + Matrix_3x3 U; + Vec3 Trans; + return QuaternionRMSD_CenteredRef(Ref, Tgt, U, Trans, useMass, 99999); +} diff --git a/src/QuaternionRMSD.h b/src/QuaternionRMSD.h index 7345d77586..02f884b126 100644 --- a/src/QuaternionRMSD.h +++ b/src/QuaternionRMSD.h @@ -1,5 +1,7 @@ #ifndef INC_QUATERNIONRMSD_H #define INC_QUATERNIONRMSD_H #include "Frame.h" +double QuaternionRMSD_CenteredRef(Frame const&, Frame&, Matrix_3x3&, Vec3&, bool, double); double QuaternionRMSD_CenteredRef(Frame const&, Frame&, Matrix_3x3&, Vec3&, bool); +double QuaternionRMSD_CenteredRef(Frame const&, Frame&, bool); #endif From e3f9b9e16dfb001a95514832fc92f9eed09792fd Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 6 Apr 2018 13:12:13 -0400 Subject: [PATCH 05/38] DRR - Cpptraj: add some debug info --- src/QuaternionRMSD.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp index ae9ec207cc..554bc9ac0d 100644 --- a/src/QuaternionRMSD.cpp +++ b/src/QuaternionRMSD.cpp @@ -103,7 +103,9 @@ double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, return CalcRMSDRotationalMatrix( Xref_, Xtgt_, len_, U.Dptr(), M_ ); */ double rmsd; + //int err = FastCalcRMSDAndRotation(U.Dptr(), rot.Dptr(), &rmsd, mwss, total_mass, minScore); + //mprintf("DEBUG: qcprot returned %i\n", err); return rmsd; } From aeb0de70fa17d0b243f582c039b2d062cb9bd76a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 6 Apr 2018 13:21:54 -0400 Subject: [PATCH 06/38] DRR - Cpptraj: Comment out unused functions --- src/qcprot.c | 21 +++++++++++---------- src/qcprot.h | 2 +- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/src/qcprot.c b/src/qcprot.c index 24687dfaea..077ad4ecfb 100644 --- a/src/qcprot.c +++ b/src/qcprot.c @@ -87,7 +87,7 @@ ******************************************************************************/ #include "qcprot.h" - +/* static double InnerProduct(double *A, double **coords1, double **coords2, const int len, const double *weight) { @@ -160,7 +160,7 @@ InnerProduct(double *A, double **coords1, double **coords2, const int len, const return (G1 + G2) * 0.5; } - +*/ int FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double len, double minScore) @@ -340,7 +340,7 @@ FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double return (1); } - +/* static void CenterCoords(double **coords, const int len, const double *weight) { @@ -387,18 +387,19 @@ CenterCoords(double **coords, const int len, const double *weight) 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); + // center the structures -- if precentered you can omit this step + CenterCoords(coords1, len, weight); + CenterCoords(coords2, len, weight); if (weight == NULL) { @@ -413,12 +414,12 @@ CalcRMSDRotationalMatrix(double **coords1, double **coords2, const int len, doub } } - /* calculate the (weighted) inner product of two structures */ + // calculate the (weighted) inner product of two structures double E0 = InnerProduct(A, coords1, coords2, len, weight); - /* calculate the RMSD & rotational matrix */ + // calculate the RMSD & rotational matrix FastCalcRMSDAndRotation(rot, A, &rmsd, E0, wsum, -1); return rmsd; } - +*/ diff --git a/src/qcprot.h b/src/qcprot.h index 2caf976ec7..19e371f4f1 100644 --- a/src/qcprot.h +++ b/src/qcprot.h @@ -96,7 +96,7 @@ extern "C" { RMSD value */ -double CalcRMSDRotationalMatrix(double **coords1, double **coords2, const int len, double *rot, const double *weight); +//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. From 715c2af646fd955a2aa9e893eec8b75a4875b704 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 9 Apr 2018 09:32:13 -0400 Subject: [PATCH 07/38] DRR - Cpptraj: Add quaternion option to 2dRMS. Only about 1.04x faster than "normal" RMSD, so does not seem worth pursuing at this time. --- src/Analysis_Rms2d.cpp | 7 ++++++- src/Analysis_Rms2d.h | 2 +- src/cpptrajdepend | 4 ++-- 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/Analysis_Rms2d.cpp b/src/Analysis_Rms2d.cpp index c68bd777f7..cc039aed2e 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 @@ -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("quat")) + mode_ = QUAT; else mode_ = RMS_FIT; useMass_ = analyzeArgs.hasKey("mass"); @@ -256,7 +259,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) @@ -271,6 +274,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/cpptrajdepend b/src/cpptrajdepend index 2ebe00a3a9..5c949903b1 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -63,7 +63,7 @@ Action_Radial.o : Action_Radial.cpp Action.h ActionState.h Action_Radial.h ArgLi Action_RandomizeIons.o : Action_RandomizeIons.cpp Action.h ActionState.h Action_RandomizeIons.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Random.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_Remap.o : Action_Remap.cpp Action.h ActionState.h Action_Remap.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h ParmFile.h ParmIO.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_ReplicateCell.o : Action_ReplicateCell.cpp Action.h ActionFrameCounter.h ActionState.h Action_ReplicateCell.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h ParmFile.h ParmIO.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h TrajectoryFile.h TrajectoryIO.h Trajout_Single.h Vec3.h -Action_Rmsd.o : Action_Rmsd.cpp Action.h ActionState.h Action_Rmsd.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h ComplexArray.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.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 ParameterTypes.h QuaternionRMSD.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h TextFormat.h Timer.h Topology.h Vec3.h +Action_Rmsd.o : Action_Rmsd.cpp Action.h ActionState.h Action_Rmsd.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h ComplexArray.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.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 ParameterTypes.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h TextFormat.h Timer.h Topology.h Vec3.h Action_Rotate.o : Action_Rotate.cpp Action.h ActionState.h Action_Rotate.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mat3x3.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_RunningAvg.o : Action_RunningAvg.cpp Action.h ActionState.h Action_RunningAvg.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_STFC_Diffusion.o : Action_STFC_Diffusion.cpp Action.h ActionState.h Action_STFC_Diffusion.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h @@ -110,7 +110,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 AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.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 AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h TextFormat.h Timer.h Topology.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 AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mesh.h DataSet_RemLog.h DataSet_integer.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 ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h TextFormat.h Timer.h Topology.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 AtomExtra.h AtomMap.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.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 AtomExtra.h AtomMap.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h ProgressBar.h QuaternionRMSD.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h Vec3.h Analysis_RmsAvgCorr.o : Analysis_RmsAvgCorr.cpp ActionState.h Analysis.h AnalysisState.h Analysis_RmsAvgCorr.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.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 AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h ComplexArray.h Constants.h CoordinateInfo.h Corr.h CpptrajFile.h CpptrajStdio.h CurveFit.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h ProgressBar.h PubFFT.h Random.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SimplexMin.h Spline.h StringRoutines.h TextFormat.h Timer.h Topology.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 AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h TextFormat.h Timer.h Topology.h Vec3.h From a0a3f0daf64de422613a9b6ce06c84f15d16e787 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 12 Apr 2022 09:13:30 -0400 Subject: [PATCH 08/38] Fix comments --- src/qcprot.c | 1 + src/qcprot.h | 29 +++++++++++++---------------- 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/src/qcprot.c b/src/qcprot.c index 077ad4ecfb..e96e542e3b 100644 --- a/src/qcprot.c +++ b/src/qcprot.c @@ -245,6 +245,7 @@ FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double 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) diff --git a/src/qcprot.h b/src/qcprot.h index 19e371f4f1..f355f39fa2 100644 --- a/src/qcprot.h +++ b/src/qcprot.h @@ -128,22 +128,19 @@ extern "C" { */ //double InnerProduct(double *A, double **coords1, double **coords2, const int len, const double *weight); -/* Calculate the RMSD, and/or the optimal rotation matrix. - - Input: - A[9] -- the inner product of two structures - E0 -- (G1 + G2) * 0.5 - len -- the size of the system - minScore-- if( minScore > 0 && rmsd < minScore) then calculate only the rmsd; - otherwise, calculate both the RMSD & the rotation matrix - Output: - rot[9] -- the rotation matrix in the order of xx, xy, xz, yx, yy, yz, zx, zy, zz - rmsd -- the RMSD value - Return: - only the rmsd was calculated if < 0 - both the RMSD & rotational matrix calculated if > 0 -*/ -int FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double wsum, double minScore); +/** 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, double minScore); /* Center the coordinates. From b8251a9622fb235553d5f271baeea1329dad13ab Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 12 Apr 2022 09:50:41 -0400 Subject: [PATCH 09/38] Add version where rotation matrix is not calculated. --- src/QuaternionRMSD.cpp | 136 +++++++++++++++++++++++++++++++++++------ 1 file changed, 119 insertions(+), 17 deletions(-) diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp index 554bc9ac0d..3516d56430 100644 --- a/src/QuaternionRMSD.cpp +++ b/src/QuaternionRMSD.cpp @@ -3,11 +3,112 @@ #include "Constants.h" #include "CpptrajStdio.h" +/** Calculate the coordinate covariance matrix and perform the + * quaternion RMSD calculation. + */ +double do_quaternion_rmsd(Frame const& Ref, Frame& Tgt, + double* U, double* Trans, + bool useMass, double minScore) +{ + 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 = + FastCalcRMSDAndRotation(U, rot.Dptr(), &rmsd, mwss, total_mass, minScore); + //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, double minScore) { - Trans.Zero(); + return do_quaternion_rmsd(Ref, Tgt, U.Dptr(), Trans.Dptr(), useMass, minScore); +/* Trans.Zero(); int ncoord = Ref.size(); double total_mass; if (useMass) { @@ -87,26 +188,27 @@ double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, } } mwss *= 0.5; // E0 = 0.5*Sum(xn^2+yn^2) -/* - // Save coordinates - int ix = 0; - for (int im = 0; im < len_; im++, ix += 3) - { - Xtgt_[0][im] = X_[ix ]; - Xtgt_[1][im] = X_[ix+1]; - Xtgt_[2][im] = X_[ix+2]; - Xref_[0][im] = R_[ix ]; - Xref_[1][im] = R_[ix+1]; - Xref_[2][im] = R_[ix+2]; - } +// +//// Save coordinates +//int ix = 0; +//for (int im = 0; im < len_; im++, ix += 3) +//{ +// Xtgt_[0][im] = X_[ix ]; +// Xtgt_[1][im] = X_[ix+1]; +// Xtgt_[2][im] = X_[ix+2]; +// Xref_[0][im] = R_[ix ]; +// Xref_[1][im] = R_[ix+1]; +// Xref_[2][im] = R_[ix+2]; +//} - return CalcRMSDRotationalMatrix( Xref_, Xtgt_, len_, U.Dptr(), M_ ); -*/ +//return CalcRMSDRotationalMatrix( Xref_, Xtgt_, len_, U.Dptr(), M_ ); +// double rmsd; //int err = FastCalcRMSDAndRotation(U.Dptr(), rot.Dptr(), &rmsd, mwss, total_mass, minScore); //mprintf("DEBUG: qcprot returned %i\n", err); return rmsd; +*/ } double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, @@ -116,9 +218,9 @@ double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, return QuaternionRMSD_CenteredRef(Ref, Tgt, U, Trans, useMass, -1); } +/** Calculate quaternion RMSD only. */ double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, bool useMass) { - Matrix_3x3 U; Vec3 Trans; - return QuaternionRMSD_CenteredRef(Ref, Tgt, U, Trans, useMass, 99999); + return do_quaternion_rmsd(Ref, Tgt, 0, Trans.Dptr(), useMass, -1); } From 8a4baca925645c1f4eac00199f7751d6f6749319 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Sun, 24 Apr 2022 07:12:58 -0400 Subject: [PATCH 10/38] Start CoordCovarMatrix --- src/CoordCovarMatrix.h | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 src/CoordCovarMatrix.h diff --git a/src/CoordCovarMatrix.h b/src/CoordCovarMatrix.h new file mode 100644 index 0000000000..da82287457 --- /dev/null +++ b/src/CoordCovarMatrix.h @@ -0,0 +1,37 @@ +#ifndef INC_COORDCOVARMATRIX_H +#define INC_COORDCOVARMATRIX_H +template class CoordCovarMatrix { + public: + CoordCovarMatrix() {} + /// Calculate covariance matrix of centered Coords and Reference (R = Xt * Ref) + T CalcCovariance_MassWt(int ncoord, T const* Ref, T const* Tgt, T const* Mass) { + T mwss = 0.0; + for (unsigned int idx = 0; idx != 9; idx++) + rot_[idx] = 0; + int im = 0; + for (int i = 0; i < ncoord; i += 3, im++) + { + T xt = Tgt[i ]; + T yt = Tgt[i+1]; + T zt = Tgt[i+2]; + T xr = Ref[i ]; + T yr = Ref[i+1]; + T zr = Ref[i+2]; + T atom_mass = 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; + } + return mwss; + } + private: + T rot_[9]; ///< Hold coordinate covariance matrix, row-major +}; +#endif From bfa4d31cdf19bf9212103514eae69052e5c01010 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Sun, 24 Apr 2022 07:19:32 -0400 Subject: [PATCH 11/38] Use atom mask --- src/CoordCovarMatrix.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/CoordCovarMatrix.h b/src/CoordCovarMatrix.h index da82287457..a40853c620 100644 --- a/src/CoordCovarMatrix.h +++ b/src/CoordCovarMatrix.h @@ -4,20 +4,21 @@ template class CoordCovarMatrix { public: CoordCovarMatrix() {} /// Calculate covariance matrix of centered Coords and Reference (R = Xt * Ref) - T CalcCovariance_MassWt(int ncoord, T const* Ref, T const* Tgt, T const* Mass) { + T CalcCovariance_MassWt(int nselected, const int* imask, T const* Ref, T const* Tgt, T const* Mass) { T mwss = 0.0; for (unsigned int idx = 0; idx != 9; idx++) rot_[idx] = 0; - int im = 0; - for (int i = 0; i < ncoord; i += 3, im++) + for (int aidx = 0; aidx != nselected; aidx++) { + int at = imask[aidx]; + int i = at * 3; T xt = Tgt[i ]; T yt = Tgt[i+1]; T zt = Tgt[i+2]; T xr = Ref[i ]; T yr = Ref[i+1]; T zr = Ref[i+2]; - T atom_mass = Mass[im]; + T atom_mass = Mass[at]; 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; From e1f3531786af225cbeb6bc2e488a15cb762c2192 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Sun, 24 Apr 2022 07:21:13 -0400 Subject: [PATCH 12/38] Update comments --- src/CoordCovarMatrix.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CoordCovarMatrix.h b/src/CoordCovarMatrix.h index a40853c620..16d8d2bcd1 100644 --- a/src/CoordCovarMatrix.h +++ b/src/CoordCovarMatrix.h @@ -3,7 +3,7 @@ template class CoordCovarMatrix { public: CoordCovarMatrix() {} - /// Calculate covariance matrix of centered Coords and Reference (R = Xt * Ref) + /// Calculate mass-weighted covariance matrix of centered Coords and Reference (R = Xt * Ref) T CalcCovariance_MassWt(int nselected, const int* imask, T const* Ref, T const* Tgt, T const* Mass) { T mwss = 0.0; for (unsigned int idx = 0; idx != 9; idx++) From 5c49bc0e9b20b5dde906ddafe41dd123067a8dcc Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Sun, 24 Apr 2022 13:07:25 -0400 Subject: [PATCH 13/38] Add in translations to origin for tgt and ref --- src/CoordCovarMatrix.h | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/CoordCovarMatrix.h b/src/CoordCovarMatrix.h index 16d8d2bcd1..b04a820f75 100644 --- a/src/CoordCovarMatrix.h +++ b/src/CoordCovarMatrix.h @@ -4,7 +4,9 @@ template class CoordCovarMatrix { public: CoordCovarMatrix() {} /// Calculate mass-weighted covariance matrix of centered Coords and Reference (R = Xt * Ref) - T CalcCovariance_MassWt(int nselected, const int* imask, T const* Ref, T const* Tgt, T const* Mass) { + T CalcCovariance_MassWt(int nselected, const int* imask, T const* Ref, T const* Tgt, T const* Mass, + T const* refTrans, T const* tgtTrans) + { T mwss = 0.0; for (unsigned int idx = 0; idx != 9; idx++) rot_[idx] = 0; @@ -12,12 +14,12 @@ template class CoordCovarMatrix { { int at = imask[aidx]; int i = at * 3; - T xt = Tgt[i ]; - T yt = Tgt[i+1]; - T zt = Tgt[i+2]; - T xr = Ref[i ]; - T yr = Ref[i+1]; - T zr = Ref[i+2]; + T xt = Tgt[i ] + tgtTrans[0]; + T yt = Tgt[i+1] + tgtTrans[1]; + T zt = Tgt[i+2] + tgtTrans[2]; + T xr = Ref[i ] + refTrans[0]; + T yr = Ref[i+1] + refTrans[1]; + T zr = Ref[i+2] + refTrans[2]; T atom_mass = Mass[at]; mwss += atom_mass * ( (xt*xt)+(yt*yt)+(zt*zt)+(xr*xr)+(yr*yr)+(zr*zr) ); rot_[0] += atom_mass*xt*xr; From 184ea4ab42b473fc6cd9a06c86b69661018c42a5 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Sat, 18 Jun 2022 12:51:16 -0400 Subject: [PATCH 14/38] Put coordinate covariance matrix calc into Matrix_3x3 --- src/Matrix_3x3.cpp | 33 +++++++++++++++++++++++++++++++++ src/Matrix_3x3.h | 2 ++ 2 files changed, 35 insertions(+) diff --git a/src/Matrix_3x3.cpp b/src/Matrix_3x3.cpp index 30e58ed386..464f3a04b5 100644 --- a/src/Matrix_3x3.cpp +++ b/src/Matrix_3x3.cpp @@ -532,6 +532,39 @@ Vec3 Matrix_3x3::AxisOfRotation(double theta) { return Vec3(0.0, 0.0, 0.0); } +/** Calculate coordinate covariance matrix between Ref and Tgt. */ +double Matrix_3x3::CalcCovariance(int nselected, const int* imask, double const* Ref, double const* Tgt, + double const* Mass, + Vec3 const& refTrans, Vec3 const& tgtTrans) +{ + double mwss = 0.0; + for (unsigned int idx = 0; idx != 9; idx++) + M_[idx] = 0; + for (int aidx = 0; aidx != nselected; aidx++) + { + int at = imask[aidx]; + int i = at * 3; + double xt = Tgt[i ] + tgtTrans[0]; + double yt = Tgt[i+1] + tgtTrans[1]; + double zt = Tgt[i+2] + tgtTrans[2]; + double xr = Ref[i ] + refTrans[0]; + double yr = Ref[i+1] + refTrans[1]; + double zr = Ref[i+2] + refTrans[2]; + double atom_mass = Mass[at]; + mwss += atom_mass * ( (xt*xt)+(yt*yt)+(zt*zt)+(xr*xr)+(yr*yr)+(zr*zr) ); + M_[0] += atom_mass*xt*xr; + M_[1] += atom_mass*xt*yr; + M_[2] += atom_mass*xt*zr; + M_[3] += atom_mass*yt*xr; + M_[4] += atom_mass*yt*yr; + M_[5] += atom_mass*yt*zr; + M_[6] += atom_mass*zt*xr; + M_[7] += atom_mass*zt*yr; + M_[8] += atom_mass*zt*zr; + } + return mwss; +} + #ifdef MPI /** Broadcast matrix from master to other ranks. */ void Matrix_3x3::BroadcastMatrix(Parallel::Comm const& commIn) { diff --git a/src/Matrix_3x3.h b/src/Matrix_3x3.h index d3fd97b460..9fec420521 100644 --- a/src/Matrix_3x3.h +++ b/src/Matrix_3x3.h @@ -50,6 +50,8 @@ class Matrix_3x3 { void RotationAroundY(double, double); void CalcRotationMatrix(Vec3 const&, double); void CalcRotationMatrix(double, double, double); + /// Calculate coordinate covariance matrix + double CalcCovariance(int, const int*, double const*, double const*, double const*, Vec3 const&, Vec3 const&); double RotationAngle(); Vec3 AxisOfRotation(double); /// Multiply 3x3 matrix times double[3] From 4e02135735ca1987fbd05c66b1f6d112194ae84a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Sat, 25 Jun 2022 10:51:02 -0400 Subject: [PATCH 15/38] Convert to C++ --- src/cpptrajdepend | 2 +- src/cpptrajfiles | 5 +++-- src/{qcprot.c => qcprot.cpp} | 7 ++++--- src/qcprot.h | 12 +++--------- 4 files changed, 11 insertions(+), 15 deletions(-) rename src/{qcprot.c => qcprot.cpp} (98%) diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 1c7928a9bb..af6a78e075 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -456,6 +456,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.c qcprot.h +qcprot.o : qcprot.cpp 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 82fd79ca6c..763872f07d 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -358,6 +358,7 @@ COMMON_SOURCES= \ ProgressBar.cpp \ ProgressTimer.cpp \ PubFFT.cpp \ + qcprot.cpp \ QuaternionRMSD.cpp \ Pucker_PuckerMask.cpp \ Pucker_PuckerSearch.cpp \ @@ -427,8 +428,8 @@ COMMON_SOURCES= \ vmdplugin/dtrplugin.cpp \ xoshiro128plusplus.cpp -CSOURCES= molsurf.c qcprot.c - +CSOURCES= molsurf.c + # The files below are used by cpptraj and libcpptraj but must be compiled # differently for libcpptraj. SOURCES=$(COMMON_SOURCES) \ diff --git a/src/qcprot.c b/src/qcprot.cpp similarity index 98% rename from src/qcprot.c rename to src/qcprot.cpp index e96e542e3b..f7d6bf1395 100644 --- a/src/qcprot.c +++ b/src/qcprot.cpp @@ -86,7 +86,8 @@ * ******************************************************************************/ #include "qcprot.h" - +#include +#include "CpptrajStdio.h" /* static double InnerProduct(double *A, double **coords1, double **coords2, const int len, const double *weight) @@ -163,7 +164,7 @@ InnerProduct(double *A, double **coords1, double **coords2, const int len, const */ int -FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double len, double minScore) +Cpptraj::QCPRot::FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double len, double minScore) { double Sxx, Sxy, Sxz, Syx, Syy, Syz, Szx, Szy, Szz; double Szz2, Syy2, Sxx2, Sxy2, Syz2, Sxz2, Syx2, Szy2, Szx2, @@ -239,7 +240,7 @@ FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double } if (i == 50) - fprintf(stderr,"\nMore than %d iterations needed!\n", i); + 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)); diff --git a/src/qcprot.h b/src/qcprot.h index f355f39fa2..880b8d0bf1 100644 --- a/src/qcprot.h +++ b/src/qcprot.h @@ -1,8 +1,7 @@ #ifndef INC_QCPROT_H #define INC_QCPROT_H -#ifdef __cplusplus -extern "C" { -#endif +namespace Cpptraj { +namespace QCPRot { /******************************************************************************* * -/_|:|_|_\- * @@ -79,10 +78,6 @@ extern "C" { * ******************************************************************************/ -#include -#include -#include - /* Calculate the RMSD & rotational matrix. Input: @@ -151,7 +146,6 @@ int FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, dou You can use CenterCoords() for this. */ //void CenterCoords(double **coords, const int len); -#ifdef __cplusplus } -#endif +} #endif From 2cfde45bd3c3e7602278520f9a8e3be3c6a03c61 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Sat, 25 Jun 2022 10:53:07 -0400 Subject: [PATCH 16/38] Add namespace to call --- src/QuaternionRMSD.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp index 3516d56430..dc8090ba1b 100644 --- a/src/QuaternionRMSD.cpp +++ b/src/QuaternionRMSD.cpp @@ -97,7 +97,7 @@ double do_quaternion_rmsd(Frame const& Ref, Frame& Tgt, double rmsd; //int err = - FastCalcRMSDAndRotation(U, rot.Dptr(), &rmsd, mwss, total_mass, minScore); + Cpptraj::QCPRot::FastCalcRMSDAndRotation(U, rot.Dptr(), &rmsd, mwss, total_mass, minScore); //mprintf("DEBUG: qcprot returned %i\n", err); return rmsd; } From 7a3291ae9b4a315af90c74f91dc896ee2e75cb88 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Sat, 25 Jun 2022 10:56:16 -0400 Subject: [PATCH 17/38] Add test for 2drms with quaternion --- test/Test_2DRMS/RunTest.sh | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/test/Test_2DRMS/RunTest.sh b/test/Test_2DRMS/RunTest.sh index b1744de3b7..b80e927ada 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 < Date: Sat, 25 Jun 2022 10:58:01 -0400 Subject: [PATCH 18/38] Add feedback for quat rmsd --- src/Analysis_Rms2d.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Analysis_Rms2d.cpp b/src/Analysis_Rms2d.cpp index 33cf6db80f..fe06f1c6c5 100644 --- a/src/Analysis_Rms2d.cpp +++ b/src/Analysis_Rms2d.cpp @@ -22,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 | quat] [mass]\n" "\t[reftraj [parm | parmindex <#>] []]\n" "\t[corr ]\n" " Calculate RMSD between all frames in , or between frames in\n" @@ -31,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() From aacb8d7f41a4d47c23688c0148135fb846f6b5c8 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 28 Jul 2022 11:22:17 -0400 Subject: [PATCH 19/38] Remove unused CalcCovariance --- src/Matrix_3x3.cpp | 33 --------------------------------- src/Matrix_3x3.h | 4 ---- 2 files changed, 37 deletions(-) diff --git a/src/Matrix_3x3.cpp b/src/Matrix_3x3.cpp index 011405b3ee..b77d6d23a9 100644 --- a/src/Matrix_3x3.cpp +++ b/src/Matrix_3x3.cpp @@ -584,39 +584,6 @@ Vec3 Matrix_3x3::AxisOfRotation(double theta) const { } // ----------------------------------------------------------------------------- -/** Calculate coordinate covariance matrix between Ref and Tgt. */ -double Matrix_3x3::CalcCovariance(int nselected, const int* imask, double const* Ref, double const* Tgt, - double const* Mass, - Vec3 const& refTrans, Vec3 const& tgtTrans) -{ - double mwss = 0.0; - for (unsigned int idx = 0; idx != 9; idx++) - M_[idx] = 0; - for (int aidx = 0; aidx != nselected; aidx++) - { - int at = imask[aidx]; - int i = at * 3; - double xt = Tgt[i ] + tgtTrans[0]; - double yt = Tgt[i+1] + tgtTrans[1]; - double zt = Tgt[i+2] + tgtTrans[2]; - double xr = Ref[i ] + refTrans[0]; - double yr = Ref[i+1] + refTrans[1]; - double zr = Ref[i+2] + refTrans[2]; - double atom_mass = Mass[at]; - mwss += atom_mass * ( (xt*xt)+(yt*yt)+(zt*zt)+(xr*xr)+(yr*yr)+(zr*zr) ); - M_[0] += atom_mass*xt*xr; - M_[1] += atom_mass*xt*yr; - M_[2] += atom_mass*xt*zr; - M_[3] += atom_mass*yt*xr; - M_[4] += atom_mass*yt*yr; - M_[5] += atom_mass*yt*zr; - M_[6] += atom_mass*zt*xr; - M_[7] += atom_mass*zt*yr; - M_[8] += atom_mass*zt*zr; - } - return mwss; -} - #ifdef MPI /** Broadcast matrix from master to other ranks. */ void Matrix_3x3::BroadcastMatrix(Parallel::Comm const& commIn) { diff --git a/src/Matrix_3x3.h b/src/Matrix_3x3.h index 9926188416..20104949ee 100644 --- a/src/Matrix_3x3.h +++ b/src/Matrix_3x3.h @@ -72,10 +72,6 @@ class Matrix_3x3 { Matrix_3x3& operator*=(double); /// \return Result of multiplying this matrix times given scalar Matrix_3x3 operator*(double) const; - - /// Calculate coordinate covariance matrix - double CalcCovariance(int, const int*, double const*, double const*, double const*, Vec3 const&, Vec3 const&); - /// Multiply 3x3 matrix times double[3] void TimesVec(double* result, const double* rhs) const { double x = rhs[0]; From a6dba01af891ccb4a6198ef04ed4ddaca3b3f5f5 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 28 Jul 2022 11:23:35 -0400 Subject: [PATCH 20/38] Remove CoordCovarMatrix class, unused --- src/CoordCovarMatrix.h | 40 ---------------------------------------- 1 file changed, 40 deletions(-) delete mode 100644 src/CoordCovarMatrix.h diff --git a/src/CoordCovarMatrix.h b/src/CoordCovarMatrix.h deleted file mode 100644 index b04a820f75..0000000000 --- a/src/CoordCovarMatrix.h +++ /dev/null @@ -1,40 +0,0 @@ -#ifndef INC_COORDCOVARMATRIX_H -#define INC_COORDCOVARMATRIX_H -template class CoordCovarMatrix { - public: - CoordCovarMatrix() {} - /// Calculate mass-weighted covariance matrix of centered Coords and Reference (R = Xt * Ref) - T CalcCovariance_MassWt(int nselected, const int* imask, T const* Ref, T const* Tgt, T const* Mass, - T const* refTrans, T const* tgtTrans) - { - T mwss = 0.0; - for (unsigned int idx = 0; idx != 9; idx++) - rot_[idx] = 0; - for (int aidx = 0; aidx != nselected; aidx++) - { - int at = imask[aidx]; - int i = at * 3; - T xt = Tgt[i ] + tgtTrans[0]; - T yt = Tgt[i+1] + tgtTrans[1]; - T zt = Tgt[i+2] + tgtTrans[2]; - T xr = Ref[i ] + refTrans[0]; - T yr = Ref[i+1] + refTrans[1]; - T zr = Ref[i+2] + refTrans[2]; - T atom_mass = Mass[at]; - 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; - } - return mwss; - } - private: - T rot_[9]; ///< Hold coordinate covariance matrix, row-major -}; -#endif From 2aa4dde3939fd3fd30e538e230f54091781e8764 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 28 Jul 2022 11:25:50 -0400 Subject: [PATCH 21/38] Update dependencies --- src/cpptrajdepend | 2 +- src/cpptrajfiles | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cpptrajdepend b/src/cpptrajdepend index a24d6886bf..1e7636e73f 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -457,6 +457,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 qcprot.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 a31c7e4eb9..8f58c5b8ae 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -430,7 +430,7 @@ COMMON_SOURCES= \ xoshiro128plusplus.cpp CSOURCES= molsurf.c - + # The files below are used by cpptraj and libcpptraj but must be compiled # differently for libcpptraj. SOURCES=$(COMMON_SOURCES) \ From 7a7644c632417805f700c74881ddae3bc3ba3f5c Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 28 Jul 2022 12:23:52 -0400 Subject: [PATCH 22/38] Initial incarnation of the quaternion rmsd cluster metric --- src/Cluster/CMakeLists.txt | 1 + src/Cluster/Metric.h | 2 +- src/Cluster/Metric_QuatRMSD.cpp | 134 ++++++++++++++++++++++++++++++++ src/Cluster/Metric_QuatRMSD.h | 40 ++++++++++ src/Cluster/clusterdepend | 1 + src/Cluster/clusterfiles | 1 + 6 files changed, 178 insertions(+), 1 deletion(-) create mode 100644 src/Cluster/Metric_QuatRMSD.cpp create mode 100644 src/Cluster/Metric_QuatRMSD.h 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/Metric_QuatRMSD.cpp b/src/Cluster/Metric_QuatRMSD.cpp new file mode 100644 index 0000000000..4a03418097 --- /dev/null +++ b/src/Cluster/Metric_QuatRMSD.cpp @@ -0,0 +1,134 @@ +#include "Metric_QuatRMSD.h" +#include "Centroid_Coord.h" +#include "Cframes.h" +#include "../CpptrajStdio.h" +#include "../QuaternionRMSD.h" + +int Cpptraj::Cluster::Metric_QuatRMSD::Init(DataSet_Coords* dIn, AtomMask const& maskIn, + bool nofit, bool useMass, int debugIn) +{ + // TODO better error handles + if (dIn == 0) { + mprinterr("Internal Error: Metric_QuatRMSD::Init() called with null data set.\n"); + return 1; + } + coords_ = dIn; + mask_ = maskIn; + useMass_ = useMass; + + return 0; +} + +int Cpptraj::Cluster::Metric_QuatRMSD::Setup() { + if (coords_->Top().SetupIntegerMask( mask_ )) return 1; + mask_.MaskInfo(); +# ifdef DEBUG_CLUSTER + mprintf("DEBUG: QuatRMSD metric topology: %s %s %i\n", coords_->legend(), + coords_->Top().c_str(), coords_->Top().Natom()); +# endif + frm1_.SetupFrameFromMask(mask_, coords_->Top().Atoms()); + frm2_ = frm1_; + return 0; +} + +double Cpptraj::Cluster::Metric_QuatRMSD::FrameDist(int f1, int f2) { + coords_->GetFrame( f1, frm1_, mask_ ); + coords_->GetFrame( f2, frm2_, mask_ ); + return QuaternionRMSD_CenteredRef(frm2_, frm1_, useMass_); +} + +double Cpptraj::Cluster::Metric_QuatRMSD::CentroidDist(Centroid* c1, Centroid* c2) { + // Centroid is already at origin. + return QuaternionRMSD_CenteredRef( ((Centroid_Coord*)c2)->Cframe(), + ((Centroid_Coord*)c1)->Cframe(), useMass_ ); +} + +double Cpptraj::Cluster::Metric_QuatRMSD::FrameCentroidDist(int f1, Centroid* c1) { + coords_->GetFrame( f1, frm1_, mask_ ); + // 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 (SRMSD_.Fit()) + cent->Cframe().CenterOnOrigin(useMass_); + } else { + Matrix_3x3 Rot; + Vec3 TgtTrans; + QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, TgtTrans, useMass_ ); + //if (SRMSD_.Fit()) { + //frm1_.Translate( TgtTrans ); + frm1_.Rotate( Rot ); + //} + cent->Cframe() += frm1_; + } + } + 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]); +} + +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; +} + +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_ ); + QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, Trans, useMass_ ); + //if (SRMSD_.Fit()) { + //frm2_.Translate( SRMSD_.TgtTrans() ); + frm1_.Rotate( 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 quaternion RMS calc. */ +std::string Cpptraj::Cluster::Metric_QuatRMSD::Description() const { + std::string description("qrmsd " + mask_.MaskExpression()); + //if (!SRMSD_.Fit()) 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 (!SRMSD_.Fit()) + // 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..39461cb381 --- /dev/null +++ b/src/Cluster/Metric_QuatRMSD.h @@ -0,0 +1,40 @@ +#ifndef INC_CLUSTER_METRIC_QUATRMSD_H +#define INC_CLUSTER_METRIC_QUATRMSD_H +#include "Metric.h" +#include "../AtomMask.h" +#include "../DataSet_Coords.h" +namespace Cpptraj { +namespace Cluster { + +/// Symmetry-corrected RMS distance calc for Coords DataSet. +class Metric_QuatRMSD : public Metric { + public: + Metric_QuatRMSD() : Metric(QRMSD), useMass_(false) {} + int Init(DataSet_Coords*,AtomMask const&,bool,bool,int); + /// \return whether RMSD is mass-weighted + bool UseMass() const { return useMass_; } + /// \return Atom mask + AtomMask const& Mask() const { return mask_; } + // ----- Metric ------------------------------ + int Setup(); + double FrameDist(int, int); + double CentroidDist( Centroid*, Centroid* ); + double FrameCentroidDist(int, Centroid*); + void CalculateCentroid(Centroid*, Cframes const&); + Centroid* NewCentroid(Cframes const&); + void FrameOpCentroid(int, Centroid*, double, CentOpType); + Metric* Copy() { return new Metric_QuatRMSD( * this ); } + std::string Description() const; + void Info() const; + unsigned int Ntotal() const { return (unsigned int)coords_->Size(); } + private: + DataSet_Coords* coords_; + AtomMask mask_; + Frame frm1_; + Frame frm2_; + bool useMass_; +}; + +} +} +#endif diff --git a/src/Cluster/clusterdepend b/src/Cluster/clusterdepend index ae2b0552ac..853b68e7ef 100644 --- a/src/Cluster/clusterdepend +++ b/src/Cluster/clusterdepend @@ -16,6 +16,7 @@ DynamicMatrix.o : DynamicMatrix.cpp ../ArrayIterator.h ../CpptrajStdio.h ../Matr 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 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 \ From b621fded6aac696ac38f72b4c384d0c266742b42 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 28 Jul 2022 13:27:12 -0400 Subject: [PATCH 23/38] Add qrmsd keyword --- src/Cluster/MetricArray.cpp | 11 ++++++++--- src/Cluster/clusterdepend | 2 +- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/Cluster/MetricArray.cpp b/src/Cluster/MetricArray.cpp index 5a2a6b3dc8..04fab3f94d 100644 --- a/src/Cluster/MetricArray.cpp +++ b/src/Cluster/MetricArray.cpp @@ -23,6 +23,7 @@ #include "Metric_Scalar.h" #include "Metric_SRMSD.h" #include "Metric_Torsion.h" +#include "Metric_QuatRMSD.h" /** CONSTRUCTOR */ Cpptraj::Cluster::MetricArray::MetricArray() : @@ -300,7 +301,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 +312,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,14 +331,15 @@ 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 if (usesrms) coordsMetricType = Metric::QRMSD; else coordsMetricType = Metric::RMS; // default // For each input set, set up the appropriate metric @@ -373,6 +376,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), nofit, useMass, debug_); break; case Metric::SCALAR : err = ((Metric_Scalar*)met)->Init((DataSet_1D*)*ds); break; case Metric::TORSION : diff --git a/src/Cluster/clusterdepend b/src/Cluster/clusterdepend index 853b68e7ef..9d60ca7856 100644 --- a/src/Cluster/clusterdepend +++ b/src/Cluster/clusterdepend @@ -14,7 +14,7 @@ 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 ../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 From 3d14652d5a527565be581ea0ec924ec1c1f6a005 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 28 Jul 2022 14:07:29 -0400 Subject: [PATCH 24/38] Try to use quat rmsd for clustering. Doesnt work yet. --- src/Cluster/MetricArray.cpp | 8 ++-- src/Cluster/Metric_QuatRMSD.cpp | 81 +++++++++++++++++++++++---------- src/Cluster/Metric_QuatRMSD.h | 24 +++++----- src/cpptrajdepend | 3 +- 4 files changed, 77 insertions(+), 39 deletions(-) diff --git a/src/Cluster/MetricArray.cpp b/src/Cluster/MetricArray.cpp index 04fab3f94d..f0e95bc4d4 100644 --- a/src/Cluster/MetricArray.cpp +++ b/src/Cluster/MetricArray.cpp @@ -278,7 +278,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; @@ -292,6 +293,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"); @@ -339,7 +341,7 @@ int Cpptraj::Cluster::MetricArray::initMetricArray(DataSetList const& setsToClus if (usedme) coordsMetricType = Metric::DME; else if (userms) coordsMetricType = Metric::RMS; else if (usesrms) coordsMetricType = Metric::SRMSD; - else if (usesrms) coordsMetricType = Metric::QRMSD; + else if (useqrms) coordsMetricType = Metric::QRMSD; else coordsMetricType = Metric::RMS; // default // For each input set, set up the appropriate metric @@ -377,7 +379,7 @@ int Cpptraj::Cluster::MetricArray::initMetricArray(DataSetList const& setsToClus 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), nofit, useMass, debug_); break; + 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 index 4a03418097..c2d27e852e 100644 --- a/src/Cluster/Metric_QuatRMSD.cpp +++ b/src/Cluster/Metric_QuatRMSD.cpp @@ -4,49 +4,81 @@ #include "../CpptrajStdio.h" #include "../QuaternionRMSD.h" +/** Initialize the metric. */ int Cpptraj::Cluster::Metric_QuatRMSD::Init(DataSet_Coords* dIn, AtomMask const& maskIn, - bool nofit, bool useMass, int debugIn) + 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: QuatRMSD metric topology: %s %s %i\n", coords_->legend(), + mprintf("DEBUG: QRMSD metric topology: %s %s %i\n", coords_->legend(), coords_->Top().c_str(), coords_->Top().Natom()); # endif - frm1_.SetupFrameFromMask(mask_, coords_->Top().Atoms()); + 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_ ); - return QuaternionRMSD_CenteredRef(frm2_, frm1_, useMass_); +# ifdef DEBUG_CLUSTER + double rms; + frm1_.printAtomCoord(0); + frm2_.printAtomCoord(0); + if (nofit_) + rms = frm1_.RMSD_NoFit( frm2_, useMass_ ); + else + rms = frm1_.RMSD( frm2_, useMass_ ); + mprintf("\tMetric_QuatRMSD::FrameDist(%i, %i)= %g\n", f1, f2, rms); + return rms; +# else + //if (nofit_) + // return frm1_.RMSD_NoFit( frm2_, useMass_ ); + //else + return QuaternionRMSD_CenteredRef( frm2_, frm1_, useMass_ ); +# endif } +/** \return RMSD between two given centroids. */ double Cpptraj::Cluster::Metric_QuatRMSD::CentroidDist(Centroid* c1, Centroid* c2) { - // Centroid is already at origin. - return QuaternionRMSD_CenteredRef( ((Centroid_Coord*)c2)->Cframe(), - ((Centroid_Coord*)c1)->Cframe(), useMass_ ); + //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_ ); - // Centroid is already at origin. - return QuaternionRMSD_CenteredRef( ((Centroid_Coord*)c1)->Cframe(), frm1_, useMass_ ); + //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 @@ -63,24 +95,23 @@ void Cpptraj::Cluster::Metric_QuatRMSD::CalculateCentroid(Centroid* centIn, Cfr coords_->GetFrame( *frm, frm1_, mask_ ); if (cent->Cframe().empty()) { cent->Cframe() = frm1_; - //if (SRMSD_.Fit()) + //if (!nofit_) cent->Cframe().CenterOnOrigin(useMass_); } else { - Matrix_3x3 Rot; - Vec3 TgtTrans; - QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, TgtTrans, useMass_ ); - //if (SRMSD_.Fit()) { - //frm1_.Translate( TgtTrans ); + //if (!nofit_) { + QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, Trans, useMass_ ); frm1_.Rotate( 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]); + //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() ); @@ -88,16 +119,18 @@ Cpptraj::Cluster::Centroid* Cpptraj::Cluster::Metric_QuatRMSD::NewCentroid( Cfra 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) + CentOpType OP) { Matrix_3x3 Rot; Vec3 Trans; Centroid_Coord* cent = (Centroid_Coord*)centIn; coords_->GetFrame( frame, frm1_, mask_ ); - QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, Trans, useMass_ ); - //if (SRMSD_.Fit()) { - //frm2_.Translate( SRMSD_.TgtTrans() ); + //if (!nofit_) { + QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, Trans, useMass_ ); frm1_.Rotate( Rot ); //} cent->Cframe().Multiply( oldSize ); @@ -110,10 +143,10 @@ void Cpptraj::Cluster::Metric_QuatRMSD::FrameOpCentroid(int frame, Centroid* cen } } -/** \return Description of quaternion RMS calc. */ +/** \return Description of RMS calc. */ std::string Cpptraj::Cluster::Metric_QuatRMSD::Description() const { std::string description("qrmsd " + mask_.MaskExpression()); - //if (!SRMSD_.Fit()) description.append(" nofit"); + //if (nofit_) description.append(" nofit"); if (useMass_) description.append(" mass"); return description; } @@ -126,7 +159,7 @@ void Cpptraj::Cluster::Metric_QuatRMSD::Info() const { mprintf(" (mask '%s')", mask_.MaskString()); if (useMass_) mprintf(", mass-weighted"); - //if (!SRMSD_.Fit()) + //if (nofit_) // mprintf(", no fitting"); //else mprintf(" best-fit"); diff --git a/src/Cluster/Metric_QuatRMSD.h b/src/Cluster/Metric_QuatRMSD.h index 39461cb381..cf959d08e4 100644 --- a/src/Cluster/Metric_QuatRMSD.h +++ b/src/Cluster/Metric_QuatRMSD.h @@ -3,37 +3,39 @@ #include "Metric.h" #include "../AtomMask.h" #include "../DataSet_Coords.h" +#include "../Frame.h" namespace Cpptraj { namespace Cluster { -/// Symmetry-corrected RMS distance calc for Coords DataSet. +/// RMS cluster distance calc for Coords DataSet class Metric_QuatRMSD : public Metric { public: - Metric_QuatRMSD() : Metric(QRMSD), useMass_(false) {} - int Init(DataSet_Coords*,AtomMask const&,bool,bool,int); - /// \return whether RMSD is mass-weighted - bool UseMass() const { return useMass_; } - /// \return Atom mask - AtomMask const& Mask() const { return mask_; } - // ----- Metric ------------------------------ + 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); - Metric* Copy() { return new Metric_QuatRMSD( * this ); } 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_; - bool useMass_; -}; +}; } } diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 1e7636e73f..03c313266b 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -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 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 From 3b272d16d92fe15c986a254d55c5a7d7dee3e1dd Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 28 Jul 2022 14:11:31 -0400 Subject: [PATCH 25/38] Fix DEBUG_CLUSTER --- src/Cluster/Metric_QuatRMSD.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Cluster/Metric_QuatRMSD.cpp b/src/Cluster/Metric_QuatRMSD.cpp index c2d27e852e..26751a12d9 100644 --- a/src/Cluster/Metric_QuatRMSD.cpp +++ b/src/Cluster/Metric_QuatRMSD.cpp @@ -49,10 +49,10 @@ double Cpptraj::Cluster::Metric_QuatRMSD::FrameDist(int f1, int f2) { double rms; frm1_.printAtomCoord(0); frm2_.printAtomCoord(0); - if (nofit_) - rms = frm1_.RMSD_NoFit( frm2_, useMass_ ); - else - rms = frm1_.RMSD( frm2_, useMass_ ); + //if (nofit_) + // rms = frm1_.RMSD_NoFit( frm2_, useMass_ ); + //else + rms = QuaternionRMSD_CenteredRef( frm2_, frm1_, useMass_ ); mprintf("\tMetric_QuatRMSD::FrameDist(%i, %i)= %g\n", f1, f2, rms); return rms; # else From 43ba9624af89eaf85a9873279ff050a460fec0f7 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 08:30:06 -0400 Subject: [PATCH 26/38] Ensure reference structures are centered --- src/Cluster/Metric_QuatRMSD.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Cluster/Metric_QuatRMSD.cpp b/src/Cluster/Metric_QuatRMSD.cpp index 26751a12d9..39b872c33b 100644 --- a/src/Cluster/Metric_QuatRMSD.cpp +++ b/src/Cluster/Metric_QuatRMSD.cpp @@ -52,6 +52,7 @@ double Cpptraj::Cluster::Metric_QuatRMSD::FrameDist(int f1, int f2) { //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; @@ -59,6 +60,7 @@ double Cpptraj::Cluster::Metric_QuatRMSD::FrameDist(int f1, int f2) { //if (nofit_) // return frm1_.RMSD_NoFit( frm2_, useMass_ ); //else + frm2_.CenterOnOrigin(useMass_); return QuaternionRMSD_CenteredRef( frm2_, frm1_, useMass_ ); # endif } From d292db48927b95eaa700ebefa766d9a2c3c59c85 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 08:48:42 -0400 Subject: [PATCH 27/38] Start quaternion rmsd clustering test. Info not matching yet. --- test/Test_Cluster_QuatRMSD/RunTest.sh | 39 +++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100755 test/Test_Cluster_QuatRMSD/RunTest.sh diff --git a/test/Test_Cluster_QuatRMSD/RunTest.sh b/test/Test_Cluster_QuatRMSD/RunTest.sh new file mode 100755 index 0000000000..95174c8ef3 --- /dev/null +++ b/test/Test_Cluster_QuatRMSD/RunTest.sh @@ -0,0 +1,39 @@ +#!/bin/bash + +. ../MasterTest.sh + +# Clean +CleanFiles cluster.in C?.cnumvtime.dat C?.cinfo.dat C?.summary.dat \ + C?.cmatrix.dat + +TESTNAME='Clustering with quaternion RMSD tests' + +INPUT="-i cluster.in" + +UNITNAME='HA clustering with regular RMSD calc.' +cat > cluster.in < cluster.in < Date: Fri, 29 Jul 2022 09:08:01 -0400 Subject: [PATCH 28/38] QuaternionRMSD rotation vectors are in columns; need to apply inverse rotation --- src/Cluster/Metric_QuatRMSD.cpp | 5 +++-- src/Cluster/Metric_RMS.cpp | 1 + src/Frame.h | 13 +++++++++++++ 3 files changed, 17 insertions(+), 2 deletions(-) diff --git a/src/Cluster/Metric_QuatRMSD.cpp b/src/Cluster/Metric_QuatRMSD.cpp index 39b872c33b..924ec81d85 100644 --- a/src/Cluster/Metric_QuatRMSD.cpp +++ b/src/Cluster/Metric_QuatRMSD.cpp @@ -102,7 +102,8 @@ void Cpptraj::Cluster::Metric_QuatRMSD::CalculateCentroid(Centroid* centIn, Cfr } else { //if (!nofit_) { QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, Trans, useMass_ ); - frm1_.Rotate( Rot ); + //Rot.Print("CalculateCentroid"); // DEBUG + frm1_.InverseRotate( Rot ); //} cent->Cframe() += frm1_; } @@ -133,7 +134,7 @@ void Cpptraj::Cluster::Metric_QuatRMSD::FrameOpCentroid(int frame, Centroid* cen coords_->GetFrame( frame, frm1_, mask_ ); //if (!nofit_) { QuaternionRMSD_CenteredRef( cent->Cframe(), frm1_, Rot, Trans, useMass_ ); - frm1_.Rotate( Rot ); + frm1_.InverseRotate( Rot ); //} cent->Cframe().Multiply( oldSize ); if (OP == ADDFRAME) { 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/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]; From 4bd90957a5ed05d9c59165fcd66237bce9851681 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 09:21:11 -0400 Subject: [PATCH 29/38] Get rid of minScore (unused and unneeded); add some more code docs. --- src/QuaternionRMSD.cpp | 22 ++++++++++------------ src/QuaternionRMSD.h | 1 - src/qcprot.cpp | 10 +++++----- src/qcprot.h | 2 +- 4 files changed, 16 insertions(+), 19 deletions(-) diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp index dc8090ba1b..e05e7c95db 100644 --- a/src/QuaternionRMSD.cpp +++ b/src/QuaternionRMSD.cpp @@ -5,10 +5,15 @@ /** 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, double minScore) + bool useMass) { Trans[0] = 0; Trans[1] = 0; @@ -97,7 +102,7 @@ double do_quaternion_rmsd(Frame const& Ref, Frame& Tgt, double rmsd; //int err = - Cpptraj::QCPRot::FastCalcRMSDAndRotation(U, rot.Dptr(), &rmsd, mwss, total_mass, minScore); + Cpptraj::QCPRot::FastCalcRMSDAndRotation(U, rot.Dptr(), &rmsd, mwss, total_mass); //mprintf("DEBUG: qcprot returned %i\n", err); return rmsd; } @@ -105,9 +110,9 @@ double do_quaternion_rmsd(Frame const& Ref, Frame& Tgt, /** Calculate quaternion RMSD with translation vector and rotation matrix. */ double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, Matrix_3x3& U, Vec3& Trans, - bool useMass, double minScore) + bool useMass) { - return do_quaternion_rmsd(Ref, Tgt, U.Dptr(), Trans.Dptr(), useMass, minScore); + return do_quaternion_rmsd(Ref, Tgt, U.Dptr(), Trans.Dptr(), useMass); /* Trans.Zero(); int ncoord = Ref.size(); double total_mass; @@ -211,16 +216,9 @@ double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, */ } -double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, - Matrix_3x3& U, Vec3& Trans, - bool useMass) -{ - return QuaternionRMSD_CenteredRef(Ref, Tgt, U, Trans, useMass, -1); -} - /** 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, -1); + return do_quaternion_rmsd(Ref, Tgt, 0, Trans.Dptr(), useMass); } diff --git a/src/QuaternionRMSD.h b/src/QuaternionRMSD.h index 02f884b126..34f6d2b442 100644 --- a/src/QuaternionRMSD.h +++ b/src/QuaternionRMSD.h @@ -1,7 +1,6 @@ #ifndef INC_QUATERNIONRMSD_H #define INC_QUATERNIONRMSD_H #include "Frame.h" -double QuaternionRMSD_CenteredRef(Frame const&, Frame&, Matrix_3x3&, Vec3&, bool, double); double QuaternionRMSD_CenteredRef(Frame const&, Frame&, Matrix_3x3&, Vec3&, bool); double QuaternionRMSD_CenteredRef(Frame const&, Frame&, bool); #endif diff --git a/src/qcprot.cpp b/src/qcprot.cpp index f7d6bf1395..72404c5349 100644 --- a/src/qcprot.cpp +++ b/src/qcprot.cpp @@ -163,8 +163,8 @@ InnerProduct(double *A, double **coords1, double **coords2, const int len, const } */ -int -Cpptraj::QCPRot::FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double len, double minScore) +/// 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, @@ -248,9 +248,9 @@ Cpptraj::QCPRot::FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, d /* 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. +// 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; diff --git a/src/qcprot.h b/src/qcprot.h index 880b8d0bf1..ed435c6be0 100644 --- a/src/qcprot.h +++ b/src/qcprot.h @@ -135,7 +135,7 @@ namespace QCPRot { * \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, double minScore); +int FastCalcRMSDAndRotation(double *rot, double *A, double *rmsd, double E0, double len); /* Center the coordinates. From 91b4f6b9855cb07e6e837f1e37c21ce0519f5392 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 09:24:37 -0400 Subject: [PATCH 30/38] Remove unused code --- src/QuaternionRMSD.cpp | 102 +---------------------------------------- 1 file changed, 1 insertion(+), 101 deletions(-) diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp index e05e7c95db..1040f0f0bc 100644 --- a/src/QuaternionRMSD.cpp +++ b/src/QuaternionRMSD.cpp @@ -15,6 +15,7 @@ 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; @@ -113,107 +114,6 @@ double QuaternionRMSD_CenteredRef(Frame const& Ref, Frame& Tgt, bool useMass) { return do_quaternion_rmsd(Ref, Tgt, U.Dptr(), Trans.Dptr(), useMass); -/* Trans.Zero(); - 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.Neg(); - 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) -// -//// Save coordinates -//int ix = 0; -//for (int im = 0; im < len_; im++, ix += 3) -//{ -// Xtgt_[0][im] = X_[ix ]; -// Xtgt_[1][im] = X_[ix+1]; -// Xtgt_[2][im] = X_[ix+2]; -// Xref_[0][im] = R_[ix ]; -// Xref_[1][im] = R_[ix+1]; -// Xref_[2][im] = R_[ix+2]; -//} - -//return CalcRMSDRotationalMatrix( Xref_, Xtgt_, len_, U.Dptr(), M_ ); -// - double rmsd; - //int err = - FastCalcRMSDAndRotation(U.Dptr(), rot.Dptr(), &rmsd, mwss, total_mass, minScore); - //mprintf("DEBUG: qcprot returned %i\n", err); - return rmsd; -*/ } /** Calculate quaternion RMSD only. */ From e8fe62ba1f1c0ae008f7bb3d5eb6d48e45ef5a6a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 09:27:48 -0400 Subject: [PATCH 31/38] Add quaternion cluster test to cluster target --- test/Makefile | 1 + 1 file changed, 1 insertion(+) 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) From 746cbedc781a1f2db69988cc591841a0bb87c93d Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 10:04:07 -0400 Subject: [PATCH 32/38] Fix typo --- src/Cluster/MetricArray.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cluster/MetricArray.cpp b/src/Cluster/MetricArray.cpp index f0e95bc4d4..4b72a3779e 100644 --- a/src/Cluster/MetricArray.cpp +++ b/src/Cluster/MetricArray.cpp @@ -333,7 +333,7 @@ int Cpptraj::Cluster::MetricArray::initMetricArray(DataSetList const& setsToClus std::string wgtArgStr = analyzeArgs.GetStringKey("wgt"); // Check args - if (usedme + userms + usesrms + useqrms> 1) { + if (usedme + userms + usesrms + useqrms > 1) { mprinterr("Error: Specify either 'dme', 'rms', 'srmsd', or 'qrmsd'.\n"); return 1; } From 2afb924b15aa8b7eb8fb812a790e5f7a68c91840 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 10:14:46 -0400 Subject: [PATCH 33/38] Add qrmsd citation info --- src/Analysis_Rms2d.cpp | 2 ++ src/Cluster/MetricArray.cpp | 7 +++++-- src/Cluster/clusterdepend | 2 +- src/QuaternionRMSD.cpp | 11 +++++++++++ src/QuaternionRMSD.h | 4 ++++ src/cpptrajdepend | 2 +- 6 files changed, 24 insertions(+), 4 deletions(-) diff --git a/src/Analysis_Rms2d.cpp b/src/Analysis_Rms2d.cpp index fe06f1c6c5..95dc16f708 100644 --- a/src/Analysis_Rms2d.cpp +++ b/src/Analysis_Rms2d.cpp @@ -141,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; } diff --git a/src/Cluster/MetricArray.cpp b/src/Cluster/MetricArray.cpp index 4b72a3779e..f667ce56c4 100644 --- a/src/Cluster/MetricArray.cpp +++ b/src/Cluster/MetricArray.cpp @@ -17,6 +17,7 @@ #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" @@ -341,8 +342,10 @@ int Cpptraj::Cluster::MetricArray::initMetricArray(DataSetList const& setsToClus if (usedme) coordsMetricType = Metric::DME; else if (userms) coordsMetricType = Metric::RMS; else if (usesrms) coordsMetricType = Metric::SRMSD; - else if (useqrms) coordsMetricType = Metric::QRMSD; - 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) diff --git a/src/Cluster/clusterdepend b/src/Cluster/clusterdepend index 9d60ca7856..d731f3ce12 100644 --- a/src/Cluster/clusterdepend +++ b/src/Cluster/clusterdepend @@ -14,7 +14,7 @@ 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_QuatRMSD.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 diff --git a/src/QuaternionRMSD.cpp b/src/QuaternionRMSD.cpp index 1040f0f0bc..2d0a28b84b 100644 --- a/src/QuaternionRMSD.cpp +++ b/src/QuaternionRMSD.cpp @@ -3,6 +3,17 @@ #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. diff --git a/src/QuaternionRMSD.h b/src/QuaternionRMSD.h index 34f6d2b442..83170b3e91 100644 --- a/src/QuaternionRMSD.h +++ b/src/QuaternionRMSD.h @@ -1,6 +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/cpptrajdepend b/src/cpptrajdepend index 03c313266b..8e704c1e10 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -164,7 +164,7 @@ 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_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 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 From c42e00eb5132cd41572d98a874c9178c580c309e Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 10:17:37 -0400 Subject: [PATCH 34/38] Add blurb about qcprot.c and authors to README --- README.md | 2 ++ 1 file changed, 2 insertions(+) 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). From 2361d06e2280082e00b27227907ca91c72c8d528 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 10:20:00 -0400 Subject: [PATCH 35/38] Use qrmsd keyword to be consistent with cluster --- src/Analysis_Rms2d.cpp | 4 ++-- test/Test_2DRMS/RunTest.sh | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Analysis_Rms2d.cpp b/src/Analysis_Rms2d.cpp index 95dc16f708..1a4ae98ff1 100644 --- a/src/Analysis_Rms2d.cpp +++ b/src/Analysis_Rms2d.cpp @@ -22,7 +22,7 @@ Analysis_Rms2d::Analysis_Rms2d() : void Analysis_Rms2d::Help() const { mprintf("\t[crdset ] [] [] [out ]\n" - "\t[dme | nofit | srmsd | quat] [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" @@ -53,7 +53,7 @@ Analysis::RetType Analysis_Rms2d::Setup(ArgList& analyzeArgs, AnalysisSetup& set mode_ = DME; else if (analyzeArgs.hasKey("srmsd")) mode_ = SRMSD; - else if (analyzeArgs.hasKey("quat")) + else if (analyzeArgs.hasKey("qrmsd")) mode_ = QUAT; else mode_ = RMS_FIT; diff --git a/test/Test_2DRMS/RunTest.sh b/test/Test_2DRMS/RunTest.sh index b80e927ada..69c8c6373d 100755 --- a/test/Test_2DRMS/RunTest.sh +++ b/test/Test_2DRMS/RunTest.sh @@ -92,7 +92,7 @@ cat > rms.in < Date: Fri, 29 Jul 2022 10:25:46 -0400 Subject: [PATCH 36/38] Add qrmsd keyword to 2drms and cluster entries --- doc/cpptraj.lyx | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) 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 From 15c1b103b558fdff85afb14ef5d12f781f3b1912 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 10:26:50 -0400 Subject: [PATCH 37/38] Remove unneeded debug print --- src/cuda_kernels/kernel_wrappers.cu | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/cuda_kernels/kernel_wrappers.cu b/src/cuda_kernels/kernel_wrappers.cu index a3c40af5bb..fb772a0b28 100644 --- a/src/cuda_kernels/kernel_wrappers.cu +++ b/src/cuda_kernels/kernel_wrappers.cu @@ -68,7 +68,6 @@ void Action_Closest_Center(const double *SolventMols_, double *D_, const double* printf("NMols = %d, NAtoms = %d\n", NMols, NAtoms); printf("active_size = %d\n", active_size); printf("NBlocks = %d\n", NBlocks); - printf("BLOCKDIM = %d\n", BLOCKDIM); printf("sizeof(double) = %d\n", sizeof(double)); printf("About to launch kernel.\n"); @@ -170,7 +169,6 @@ void Action_Closest_NoCenter(const double *SolventMols_, double *D_, const doubl printf("NMols = %d, NAtoms = %d\n", NMols, NAtoms); printf("active_size = %d\n", active_size); printf("NBlocks = %d\n", NBlocks); - printf("BLOCKDIM = %d\n", BLOCKDIM); printf("sizeof(double) = %d\n", sizeof(double)); printf("About to launch kernel.\n"); From f3bc2e530c0362dd6b252b99f459466a9af3d676 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 29 Jul 2022 10:28:25 -0400 Subject: [PATCH 38/38] 6.12.2. Revision bump for adding quaternion rmsd calc to 2drms and cluster. --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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