diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf index 46ba898f3a..70719f82c0 100644 Binary files a/doc/CpptrajManual.pdf and b/doc/CpptrajManual.pdf differ diff --git a/doc/DocumentChecksums.txt b/doc/DocumentChecksums.txt index b8633da4e1..60e6cefb03 100644 --- a/doc/DocumentChecksums.txt +++ b/doc/DocumentChecksums.txt @@ -1,3 +1,3 @@ b37726e7a841f6fc695ecd7fb040ffbf CpptrajDevelopmentGuide.lyx -c9b19c2bade122b63363f364f86a0a69 cpptraj.lyx +736ee725ca9cc67d349f67d714f41c20 cpptraj.lyx 07c4039e732fc2eb1df0c1e0863cb949 CpptrajManual.lyx diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 54b064c15d..a97dd7af40 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -5640,7 +5640,7 @@ Read \end_layout \begin_layout LyX-Code - [vector] | + [vector [magnitude]] | \end_layout \begin_layout LyX-Code @@ -5737,7 +5737,7 @@ nosquare2d Read data as XYZ matrix (i.e. \begin_layout Description vector Read data as vector. If indices are present they will be skipped. - Assume first 3 columns after the index volumn are vector X, + Assume first 3 columns after the index column are vector X, Y, and Z, and (if present) the next 3 columns contain vector origin X, @@ -5745,6 +5745,14 @@ vector Read data as vector. and Z. \end_layout +\begin_deeper +\begin_layout Description +magnitude If specified, + assume vector file last column contains vector magnitudes; + read these into a set with Aspect 'Mag'. +\end_layout + +\end_deeper \begin_layout Description mat3x3 Read data as 3x3 matrix. If indices are present they will be skipped. @@ -55687,7 +55695,7 @@ vectormath vec1 vec2 [out ] [norm] [name // std::min, std::max +#include // sqrt #include "Analysis_VectorMath.h" #include "CpptrajStdio.h" #include "Constants.h" @@ -7,7 +8,7 @@ /// Strings corresponding to modes, used in output. const char* Analysis_VectorMath::ModeString[] = { - "Dot product", "Angle from dot product", "Cross product" }; + "Dot product", "Angle from dot product", "Cross product", "Magnitude" }; // CONSTRUCTOR Analysis_VectorMath::Analysis_VectorMath() : @@ -17,46 +18,24 @@ Analysis_VectorMath::Analysis_VectorMath() : void Analysis_VectorMath::Help() const { mprintf("\tvec1 vec2 [out ] [norm] [name ]\n" - "\t[ dotproduct | dotangle | crossproduct ]\n" + "\t[ dotproduct | dotangle | crossproduct | magnitude ]\n" " Calculate dot product, angle from dot product (degrees), or cross product\n" " for specified vectors. Either vec1 or vec2 can be size 1, otherwise they\n" - " must both be the same size.\n"); + " must both be the same size. If 'magnitude' is specified, just calculate\n" + " the magnitudes of the vectors selected by 'vec1' (no need to specify\n" + " 'vec2').\n"); } // Analysis_VectorMath::Setup() Analysis::RetType Analysis_VectorMath::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn) { - // Get Vectors - DataSetList vsets1 = setup.DSL().SelectGroupSets( analyzeArgs.GetStringKey("vec1"), - DataSet::VECTOR_1D ); - if (vsets1.empty()) { - mprinterr("Error: 'vec1' not found.\n"); - return Analysis::ERR; - } - DataSetList vsets2 = setup.DSL().SelectGroupSets( analyzeArgs.GetStringKey("vec2"), - DataSet::VECTOR_1D ); - if (vsets2.empty()) { - mprinterr("Error: 'vec2' not found.\n"); - return Analysis::ERR; - } - - if (vsets1.size() != vsets2.size()) { - mprinterr("Error: 'vec1' (%zu) and 'vec2' (%zu) do not select the same number of sets.\n", - vsets1.size(), vsets2.size()); - return Analysis::ERR; - } - - for (DataSetList::const_iterator it = vsets1.begin(); it != vsets1.end(); ++it) - vinfo1_.push_back( static_cast( *it ) ); - for (DataSetList::const_iterator it = vsets2.begin(); it != vsets2.end(); ++it) - vinfo2_.push_back( static_cast( *it ) ); - std::string setname = analyzeArgs.GetStringKey("name"); norm_ = analyzeArgs.hasKey("norm"); // Check for dotproduct/crossproduct keywords. Default is dotproduct. mode_ = DOTPRODUCT; DataSet::DataType dtype = DataSet::DOUBLE; const char* dname = "Dot"; + bool requires_two_vecs = true; if (analyzeArgs.hasKey("dotproduct")) { mode_ = DOTPRODUCT; } else if (analyzeArgs.hasKey("dotangle")) { @@ -67,7 +46,47 @@ Analysis::RetType Analysis_VectorMath::Setup(ArgList& analyzeArgs, AnalysisSetup mode_ = CROSSPRODUCT; dtype = DataSet::VECTOR; dname = "Cross"; + } else if (analyzeArgs.hasKey("magnitude")) { + mode_ = MAGNITUDE; + dname = "Mag"; + requires_two_vecs = false; + if (norm_) { + mprintf("Warning: 'norm' does not make sense with 'magnitude', ignoring.\n"); + norm_ = false; + } } + + // Get Vectors + DataSetList vsets1 = setup.DSL().SelectGroupSets( analyzeArgs.GetStringKey("vec1"), + DataSet::VECTOR_1D ); + if (vsets1.empty()) { + mprinterr("Error: 'vec1' not found.\n"); + return Analysis::ERR; + } + + DataSetList vsets2; + if (requires_two_vecs) { + vsets2 = setup.DSL().SelectGroupSets( analyzeArgs.GetStringKey("vec2"), + DataSet::VECTOR_1D ); + if (vsets2.empty()) { + mprinterr("Error: 'vec2' not found.\n"); + return Analysis::ERR; + } + + if (vsets1.size() != vsets2.size()) { + mprinterr("Error: 'vec1' (%zu) and 'vec2' (%zu) do not select the same number of sets.\n", + vsets1.size(), vsets2.size()); + return Analysis::ERR; + } + } + + for (DataSetList::const_iterator it = vsets1.begin(); it != vsets1.end(); ++it) + vinfo1_.push_back( static_cast( *it ) ); + if (requires_two_vecs) { + for (DataSetList::const_iterator it = vsets2.begin(); it != vsets2.end(); ++it) + vinfo2_.push_back( static_cast( *it ) ); + } + // Set up output file in DataFileList if necessary DataFile* outfile = setup.DFL().AddDataFile( analyzeArgs.GetStringKey("out"), analyzeArgs ); // Set up output data sets based on mode @@ -89,12 +108,22 @@ Analysis::RetType Analysis_VectorMath::Setup(ArgList& analyzeArgs, AnalysisSetup // Print Status mprintf(" VECTORMATH:"); - if (vinfo1_.size() == 1) - mprintf(" Calculating %s of vectors %s and %s\n", ModeString[mode_], vinfo1_[0]->legend(), vinfo2_[0]->legend()); - else { - mprintf(" Calculating %s of:\n", ModeString[mode_]); - for (unsigned int ii = 0; ii < vinfo1_.size(); ii++) - mprintf("\t %s and %s\n", vinfo1_[ii]->legend(), vinfo2_[ii]->legend()); + if (requires_two_vecs) { + if (vinfo1_.size() == 1) + mprintf(" Calculating %s of vectors %s and %s\n", ModeString[mode_], vinfo1_[0]->legend(), vinfo2_[0]->legend()); + else { + mprintf(" Calculating %s of:\n", ModeString[mode_]); + for (unsigned int ii = 0; ii < vinfo1_.size(); ii++) + mprintf("\t %s and %s\n", vinfo1_[ii]->legend(), vinfo2_[ii]->legend()); + } + } else { + if (vinfo1_.size() == 1) + mprintf(" Calculating %s of vector %s\n", ModeString[mode_], vinfo1_[0]->legend()); + else { + mprintf(" Calculating %s of:\n", ModeString[mode_]); + for (DVarray::const_iterator it = vinfo1_.begin(); it != vinfo1_.end(); ++it) + mprintf("\t %s\n", (*it)->legend()); + } } if (norm_) mprintf("\tVectors will be normalized.\n"); if (outfile != 0) @@ -103,6 +132,19 @@ Analysis::RetType Analysis_VectorMath::Setup(ArgList& analyzeArgs, AnalysisSetup return Analysis::OK; } +/** Calculate the magnitude of each vector. */ +int Analysis_VectorMath::Magnitude(DataSet* Dout, DataSet_Vector const& V1) +const +{ + DataSet_double& Out = static_cast( *Dout ); + Out.Resize( V1.Size() ); + for (unsigned int idx = 0; idx < V1.Size(); idx++) + { + Out[idx] = sqrt( V1[idx].Magnitude2() ); + } + return 0; +} + // Analysis_VectorMath::DotProduct() int Analysis_VectorMath::DotProduct(DataSet* Dout, DataSet_Vector& V1, DataSet_Vector& V2, unsigned int vmax, unsigned int v1inc, unsigned int v2inc) @@ -182,10 +224,30 @@ const // Analysis_VectorMath::Analyze() Analysis::RetType Analysis_VectorMath::Analyze() { - for (unsigned int ii = 0; ii < vinfo1_.size(); ii++) - { - int err = DoMath( DataOut_[ii], *(vinfo1_[ii]), *(vinfo2_[ii]) ); - if (err != 0) return Analysis::ERR; + if (vinfo2_.empty()) { + // Single set ops + for (unsigned int ii = 0; ii < vinfo1_.size(); ii++) + { + if (vinfo1_[ii]->Size() < 1) { + mprintf("Warning: Vector set '%s' is empty.\n", vinfo1_[ii]->legend()); + continue; + } + int err = 1; + if (mode_ == MAGNITUDE) + err = Magnitude( DataOut_[ii], *(vinfo1_[ii]) ); + if (err != 0) { + mprinterr("Error: A problem occurred when performing vector math for set '%s'\n", + vinfo1_[ii]->legend()); + return Analysis::ERR; + } + } + } else { + // Two set ops + for (unsigned int ii = 0; ii < vinfo1_.size(); ii++) + { + int err = DoMath( DataOut_[ii], *(vinfo1_[ii]), *(vinfo2_[ii]) ); + if (err != 0) return Analysis::ERR; + } } return Analysis::OK; } diff --git a/src/Analysis_VectorMath.h b/src/Analysis_VectorMath.h index 27365d568f..2db995b74b 100644 --- a/src/Analysis_VectorMath.h +++ b/src/Analysis_VectorMath.h @@ -11,9 +11,10 @@ class Analysis_VectorMath : public Analysis { Analysis::RetType Setup(ArgList&, AnalysisSetup&, int); Analysis::RetType Analyze(); private: - enum ModeType { DOTPRODUCT = 0, DOTANGLE, CROSSPRODUCT }; + enum ModeType { DOTPRODUCT = 0, DOTANGLE, CROSSPRODUCT, MAGNITUDE }; static const char* ModeString[]; + int Magnitude(DataSet*, DataSet_Vector const&) const; int DotProduct(DataSet*, DataSet_Vector&, DataSet_Vector&, unsigned int, unsigned int, unsigned int) const; int CrossProduct(DataSet*, DataSet_Vector&, DataSet_Vector&, diff --git a/src/DataIO_Std.cpp b/src/DataIO_Std.cpp index eb4d73be79..8f9fe8695d 100644 --- a/src/DataIO_Std.cpp +++ b/src/DataIO_Std.cpp @@ -32,7 +32,8 @@ DataIO_Std::DataIO_Std() : indexcol_(-1), isInverted_(false), hasXcolumn_(true), - writeHeader_(true), + writeHeader_(true), + read_vector_magnitude_(false), square2d_(true), sparse_(false), originSpecified_(false), @@ -92,6 +93,7 @@ void DataIO_Std::ReadHelp() { "\t\tprec {dbl|flt*} : Grid precision; double or float (default float).\n" "\t\tbin {center|corner*} : Coords specify bin centers or corners (default corners).\n" "\tvector : Read data as vector: VX VY VZ [OX OY OZ]\n" + "\t\tmagnitude : Assume vector data final column contains vector magnitude.\n" "\tmat3x3 : Read data as 3x3 matrices: M(1,1) M(1,2) ... M(3,2) M(3,3)\n"); } @@ -158,6 +160,11 @@ int DataIO_Std::processReadArgs(ArgList& argIn) { } } } + // Options for vector + if (mode_ == READVEC) { + if (argIn.hasKey("magnitude")) + read_vector_magnitude_ = true; + } // Options for 2d if (mode_ == READ2D) { if (argIn.hasKey("square2d")) square2d_ = true; @@ -853,6 +860,17 @@ int DataIO_Std::Read_Vector(std::string const& fname, return 1; } } + DataSet* magset = 0; + if (read_vector_magnitude_) { + magset = datasetlist.CheckForSet( MetaData(dsname, "Mag") ); + if (magset != 0) { + mprintf("\tAppening vector magnitude data to set '%s'\n", magset->legend()); + if (magset->Type() != DataSet::DOUBLE) { + mprinterr("Error: Cannot append magnitude data to set '%s', wrong type.\n", magset->legend()); + return 1; + } + } + } // Buffer file BufferedLine buffer; if (buffer.OpenFileRead( fname )) return 1; @@ -865,24 +883,40 @@ int DataIO_Std::Read_Vector(std::string const& fname, // 9 (VXYZ OXYZ VXYZ+OXYZ) values, optionally with indices. int ntokens = buffer.TokenizeLine( SEPARATORS ); int ncols = ntokens; // Number of columns of vector data. - int nv = 0; // Number of columns to actually read from (3 or 6). + int nexpected = 0; // Number of columns expected to actually read from (3 or 6). bool hasIndex; if (ntokens < 1) { mprinterr("Error: Could not tokenize line.\n"); return 1; } - if (ncols == 3 || ncols == 6 || ncols == 9) - hasIndex = false; - else if (ncols == 4 || ncols == 7 || ncols == 10) { - hasIndex = true; - mprintf("Warning: Not reading vector data indices.\n"); + // Try to determine what data is present based on the number of columns + if (read_vector_magnitude_) { + mprintf("\tAssuming final column contains vector magnitude data.\n"); + if (ncols == 4 || ncols == 7 || ncols == 10) + hasIndex = false; + else if (ncols == 5 || ncols == 8 || ncols == 11) { + hasIndex = true; + mprintf("Warning: Not reading vector data indices.\n"); + } else { + mprinterr("Error: Expected 3, 6, or 9 columns of vector data plus 1 column\n" + "Error: of magnitude data, got %i.\n", ncols); + return 1; + } } else { - mprinterr("Error: Expected 3, 6, or 9 columns of vector data, got %i.\n", ncols); - return 1; + if (ncols == 3 || ncols == 6 || ncols == 9) + hasIndex = false; + else if (ncols == 4 || ncols == 7 || ncols == 10) { + hasIndex = true; + mprintf("Warning: Not reading vector data indices.\n"); + } else { + mprinterr("Error: Expected 3, 6, or 9 columns of vector data, got %i.\n", ncols); + mprinterr("Error: If vector data contains a magnitude column, specify 'magnitude'.\n"); + return 1; + } } bool hasOrigins; if (ncols >= 6) { - nv = 6; + nexpected = 6; mprintf("\tReading vector X Y Z and origin X Y Z values.\n"); hasOrigins = true; // If set already exists, see if it doesnt have origins. @@ -903,7 +937,7 @@ int DataIO_Std::Read_Vector(std::string const& fname, delete oldSet; } } else { - nv = 3; + nexpected = 3; mprintf("\tReading vector X Y Z values.\n"); hasOrigins = false; // If set already exists, see if it has origins. @@ -923,26 +957,37 @@ int DataIO_Std::Read_Vector(std::string const& fname, ds = datasetlist.AddSet(dtype, dsname); if (ds == 0) return 1; } + if (read_vector_magnitude_ && magset == 0) { + magset = datasetlist.AddSet(DataSet::DOUBLE, MetaData(dsname, "Mag")); + if (magset == 0) return 1; + } // Read vector data double vec[6]; std::fill(vec, vec+6, 0.0); + double vmag = 0; + // +1 expected column if reading magnitudes + if (magset != 0) + nexpected++; size_t ndata = ds->Size(); while (linebuffer != 0) { if (hasIndex) - ntokens = sscanf(linebuffer, "%*f %lf %lf %lf %lf %lf %lf", - vec, vec+1, vec+2, vec+3, vec+4, vec+5); + ntokens = sscanf(linebuffer, "%*f %lf %lf %lf %lf %lf %lf %lf", + vec, vec+1, vec+2, vec+3, vec+4, vec+5, &vmag); else - ntokens = sscanf(linebuffer, "%lf %lf %lf %lf %lf %lf", - vec, vec+1, vec+2, vec+3, vec+4, vec+5); - if (ntokens != nv) { + ntokens = sscanf(linebuffer, "%lf %lf %lf %lf %lf %lf %lf", + vec, vec+1, vec+2, vec+3, vec+4, vec+5, &vmag); + if (ntokens != nexpected) { mprinterr("Error: In vector file, line %i: expected %i values, got %i\n", - buffer.LineNumber(), nv, ntokens); + buffer.LineNumber(), nexpected, ntokens); break; } if (hasOrigins) - ds->Add( ndata++, vec ); + ds->Add( ndata, vec ); else ((DataSet_Vector*)ds)->AddVxyz( Vec3(vec) ); // TODO: Remove if DataSet_Vector is split to with/without origins + if (magset != 0) + magset->Add(ndata, &vmag); + ndata++; linebuffer = buffer.Line(); } return 0; diff --git a/src/DataIO_Std.h b/src/DataIO_Std.h index fdd80552db..bfb6f5d5d9 100644 --- a/src/DataIO_Std.h +++ b/src/DataIO_Std.h @@ -51,6 +51,7 @@ class DataIO_Std : public DataIO { bool isInverted_; ///< For 1D writes invert X/Y. bool hasXcolumn_; bool writeHeader_; + bool read_vector_magnitude_; ///< For vector reads, assume final column contains magnitude bool square2d_; bool sparse_; ///< 3d writes, only write voxels with value > cut_ bool originSpecified_; ///< 3d reads, true if origin specified diff --git a/src/Version.h b/src/Version.h index df86dd421c..2de71ea96e 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.29.2" +#define CPPTRAJ_INTERNAL_VERSION "V6.29.3" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/test/Test_VectorMath/Magnitude.dat.save b/test/Test_VectorMath/Magnitude.dat.save new file mode 100644 index 0000000000..ea7784bea4 --- /dev/null +++ b/test/Test_VectorMath/Magnitude.dat.save @@ -0,0 +1,11 @@ +#Frame V1Mag + 1 12.0310 + 2 6.1779 + 3 5.1423 + 4 6.9181 + 5 9.0597 + 6 10.2285 + 7 7.1507 + 8 3.8188 + 9 4.4247 + 10 4.6969 diff --git a/test/Test_VectorMath/RunTest.sh b/test/Test_VectorMath/RunTest.sh index 7c2f05de38..0b538247ae 100755 --- a/test/Test_VectorMath/RunTest.sh +++ b/test/Test_VectorMath/RunTest.sh @@ -2,15 +2,18 @@ . ../MasterTest.sh -CleanFiles vectors.dat dotproduct.dat corr.in v1init_dot_v1.dat +CleanFiles vectors.dat dotproduct.dat corr.in v1init_dot_v1.dat Magnitude.dat INPUT="corr.in" TOP=../tz2.parm7 -TESTNAME='Vector dot/cross product test' -Requires netcdf +TESTNAME='Vector math tests' + # Dot/cross products -cat > corr.in < corr.in < corr.in <