From 52a945570f92f65539bf841f4b77bce451614c5a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 21 Aug 2020 13:29:36 -0400 Subject: [PATCH 1/7] Start importing changes from the fixspam branch. Improved formatting for OpenDX files, option to read grids in as double precision, enable DIV and MULT for grids, warn if extra args remain after 'purewater' specified for SPAM, use DataSet/DataFile framework for peaks in SPAM, add DataIO for read and write of peaks data, fix Tinker file recognition, make sure SPAM is not hidden, skip peaks where bin calc could overflow, avoid using zero bandwidth for KDE hist when only 1 value available. --- src/Action_Spam.cpp | 192 ++++++++++++++++++++-------------- src/Action_Spam.h | 6 +- src/DataFile.cpp | 5 + src/DataFile.h | 2 +- src/DataIO_OpenDx.cpp | 97 ++++++++++++++--- src/DataIO_OpenDx.h | 8 +- src/DataIO_Peaks.cpp | 166 +++++++++++++++++++++++++++++ src/DataIO_Peaks.h | 17 +++ src/DataSet.cpp | 3 +- src/DataSet.h | 2 +- src/DataSetList.cpp | 2 + src/DataSet_3D.h | 9 +- src/DataSet_GridDbl.h | 32 +++--- src/DataSet_GridFlt.h | 34 +++--- src/DataSet_Vector.cpp | 1 + src/DataSet_Vector_Scalar.cpp | 80 ++++++++++++++ src/DataSet_Vector_Scalar.h | 38 +++++++ src/GridBin.h | 1 + src/KDE.cpp | 4 +- src/RPNcalc.cpp | 6 +- src/TextFormat.h | 6 ++ src/TinkerFile.cpp | 1 + 22 files changed, 583 insertions(+), 129 deletions(-) create mode 100644 src/DataIO_Peaks.cpp create mode 100644 src/DataIO_Peaks.h create mode 100644 src/DataSet_Vector_Scalar.cpp create mode 100644 src/DataSet_Vector_Scalar.h diff --git a/src/Action_Spam.cpp b/src/Action_Spam.cpp index 2368a3e141..47e6163972 100644 --- a/src/Action_Spam.cpp +++ b/src/Action_Spam.cpp @@ -12,6 +12,8 @@ #include "KDE.h" #include "OnlineVarT.h" // Stats #include "DataSet_Mesh.h" +#include "DataSet_Vector_Scalar.h" +#include "DataIO_Peaks.h" // CONSTRUCTOR Action_Spam::Action_Spam() : @@ -32,33 +34,67 @@ Action_Spam::Action_Spam() : ds_dh_(0), ds_ds_(0), Nframes_(0), - overflow_(false) + overflow_(false), + peaksData_(0) { } +/** Search for DataSet with peaks data. If that fails, try to load peaks + * from a file. + */ +int Action_Spam::GetPeaks(std::string const& name, DataSetList const& dsl) +{ + // Check for peaks DataSet. + DataSet* ds = dsl.FindSetOfType(name, DataSet::VECTOR_SCALAR); + if (ds == 0) { + // No set found. See if file exists. + FileName fname(name); + if (!File::Exists(fname)) { + File::ErrorMsg( fname.full() ); + mprinterr("Error: No peak data or file with name '%s'.\n", name.c_str()); + return 1; + } + // Try to load peaks from file. + DataIO_Peaks infile; + if (infile.ReadData(fname, peaksdsl_, fname.Base())) { + mprinterr("Error: Could not load peaks data from %s\n", fname.full()); + return 1; + } + // Sanity check + if (peaksdsl_.size() < 1 || peaksdsl_[0]->Type() != DataSet::VECTOR_SCALAR) { + mprinterr("Error: Could not allocate peaks data set for file.\n"); + return 1; + } + peaksData_ = (DataSet_Vector_Scalar*)peaksdsl_[0]; + } else { + peaksData_ = (DataSet_Vector_Scalar*)ds; + } + return 0; +} + void Action_Spam::Help() const { mprintf("\t[name ] [out ] [cut ] [solv ]\n" "\t{ purewater |\n" - "\t [reorder] [info ] [summary ]\n" + "\t [reorder] [info ] [summary ]\n" "\t [site_size ] [sphere] [temperature ]\n" "\t [dgbulk ] [dhbulk ] }\n" " Perform SPAM water analysis. If 'purewater' is specified calculate\n" " bulk energy values for a pure water system. Otherwise determine SPAM\n" " energies from peaks previously identified from the 'volmap' action.\n" - " : Output data set name.\n" - " : Data file with all SPAM energies for each snapshot.\n" - " : Non-bonded cutoff for energy evaluation\n" - " : Name of the solvent residues\n" - " [purewater]: The system is pure water---used to parametrize the bulk values.\n" - " : File with the peak locations present (XYZ- format)\n" - " [reorder] : The solvent should be re-ordered so the same solvent molecule\n" - " is always in the same site.\n" - " : File with stats about which sites are occupied when.\n" - " : File with the summary of all SPAM results.\n" - " : Size of the water site around each density peak.\n" - " [sphere] : Treat each site like a sphere.\n" - " : Temperature at which SPAM calculation was run.\n" - " : SPAM free energy of the bulk solvent in kcal/mol\n" - " : SPAM enthalpy of the bulk solvent in kcal/mol\n"); + " : Output data set name.\n" + " : Data file with all SPAM energies for each snapshot.\n" + " : Non-bonded cutoff for energy evaluation\n" + " : Name of the solvent residues\n" + " [purewater] : The system is pure water---used to parametrize the bulk values.\n" + " : Dataset/File (XYZ format) with the peak locations present.\n" + " [reorder] : The solvent should be re-ordered so the same solvent molecule\n" + " is always in the same site.\n" + " : File with stats about which sites are occupied when.\n" + " : File with the summary of all SPAM results.\n" + " : Size of the water site around each density peak.\n" + " [sphere] : Treat each site like a sphere.\n" + " : Temperature at which SPAM calculation was run.\n" + " : SPAM free energy of the bulk solvent in kcal/mol\n" + " : SPAM enthalpy of the bulk solvent in kcal/mol\n"); } // Action_Spam::Init() @@ -71,7 +107,7 @@ Action::RetType Action_Spam::Init(ArgList& actionArgs, ActionInit& init, int deb // Always use imaged distances image_.InitImaging(true); // This is needed everywhere in this function scope - FileName filename; + std::string peaksname; // See if we're doing pure water. If so, we don't need a peak file purewater_ = actionArgs.hasKey("purewater"); @@ -108,14 +144,14 @@ Action::RetType Action_Spam::Init(ArgList& actionArgs, ActionInit& init, int deb myDSL_.push_back( ds ); DG_BULK_ = 0.0; DH_BULK_ = 0.0; + // Shouldn't need any more arguments. + if (actionArgs.NremainingArgs() > 0) + mprintf("Warning: 'purewater' specified but more arguments remain.\n"); } else { - // Get the file name with the peaks defined in it - filename.SetFileName( actionArgs.GetStringNext() ); - if (filename.empty()) { - mprinterr("Error: No Peak file specified.\n"); - return Action::ERR; - } else if (!File::Exists(filename)) { - File::ErrorMsg( filename.full() ); + // Get the file/dataset name with the peaks defined in it + peaksname = actionArgs.GetStringNext(); + if (peaksname.empty()) { + mprinterr("Error: No Peak dataset/file specified.\n"); return Action::ERR; } // Get the remaining optional arguments @@ -139,39 +175,13 @@ Action::RetType Action_Spam::Init(ArgList& actionArgs, ActionInit& init, int deb // If it's a sphere, square the radius to compare with if (sphere_) site_size_ *= site_size_; - // Parse through the peaks file and extract the peaks - CpptrajFile peakfile; - if (peakfile.OpenRead(filename)) { - mprinterr("SPAM: Error: Could not open %s for reading!\n", filename.full()); + // Get or load the peaks data + if (GetPeaks(peaksname, init.DSL())) { + mprinterr("Error: Could not get peaks.\n"); return Action::ERR; } - std::string line = peakfile.GetLine(); - int npeaks = 0; - while (!line.empty()) { - if (sscanf(line.c_str(), "%d", &npeaks) != 1) { - line = peakfile.GetLine(); - continue; - } - line = peakfile.GetLine(); - break; - } - while (!line.empty()) { - double x, y, z, dens; - if (sscanf(line.c_str(), "C %lg %lg %lg %lg", &x, &y, &z, &dens) != 4) { - line = peakfile.GetLine(); - continue; - } - line = peakfile.GetLine(); - peaks_.push_back(Vec3(x, y, z)); - } - peakfile.CloseFile(); - // Check that our initial number of peaks matches our parsed peaks. Warn - // otherwise - if (npeaks != (int)peaks_.size()) - mprinterr("SPAM: Warning: %s claims to have %d peaks, but really has %zu!\n", - filename.full(), npeaks, peaks_.size()); // Now add all of the individual peak energy data sets - for (int i = 0; i < (int)peaks_.size(); i++) { + for (unsigned int i = 0; i < peaksData_->Size(); i++) { DataSet* ds = init.DSL().AddSet(DataSet::DOUBLE, MetaData(ds_name,i+1)); if (ds == 0) return Action::ERR; myDSL_.push_back( ds ); @@ -198,7 +208,7 @@ Action::RetType Action_Spam::Init(ArgList& actionArgs, ActionInit& init, int deb # endif // peakFrameData will keep track of omitted frames for each peak. peakFrameData_.clear(); - peakFrameData_.resize( peaks_.size() ); + peakFrameData_.resize( peaksData_->Size() ); } // Determine if energy calculation needs to happen calcEnergy_ = (summaryfile != 0 || datafile != 0); @@ -223,10 +233,8 @@ Action::RetType Action_Spam::Init(ArgList& actionArgs, ActionInit& init, int deb mprintf("\tPrinting solvent SPAM summary to %s\n", summaryfile->DataFilename().full()); } else { - mprintf("\tSolvent [%s] density peaks taken from %s.\n", - solvname_.c_str(), filename.base()); - mprintf("\t%zu density peaks will be analyzed from %s.\n", - peaks_.size(), filename.base()); + mprintf("\tSolvent [%s], %zu density peaks taken from %s.\n", + solvname_.c_str(), peaksData_->Size(), peaksData_->legend()); mprintf("\tOccupation information printed to %s.\n", infofile_->Filename().full()); mprintf("\tSites are "); if (sphere_) @@ -536,17 +544,17 @@ Action::RetType Action_Spam::DoSPAM(int frameNum, Frame& frameIn) { t_assign_.Start(); // Loop through each peak and then scan through every residue, and assign a // solvent residue to each peak - int pknum = 0; - for (Varray::const_iterator pk = peaks_.begin(); pk != peaks_.end(); ++pk, ++pknum) + for (unsigned int pknum = 0; pknum < peaksData_->Size(); pknum++) { + Vec3 const& pk = peaksData_->Vec( pknum ); for (unsigned int resnum = 0; resnum != comlist_.size(); resnum++) { // If we're inside, make sure this residue is not already `claimed'. If it // is, assign it to the closer peak center - if ((this->*Inside_)(*pk, comlist_[resnum], site_size_)) { + if ((this->*Inside_)(pk, comlist_[resnum], site_size_)) { if (resPeakNum_[resnum] > 0) { - Vec3 diff1 = comlist_[resnum] - *pk; - Vec3 diff2 = comlist_[resnum] - peaks_[ resPeakNum_[resnum] ]; + Vec3 diff1 = comlist_[resnum] - pk; + Vec3 diff2 = comlist_[resnum] - peaksData_->Vec( resPeakNum_[resnum] ); // If we are closer, update. Otherwise do nothing if (diff1.Magnitude2() < diff2.Magnitude2()) resPeakNum_[resnum] = pknum; @@ -563,8 +571,8 @@ Action::RetType Action_Spam::DoSPAM(int frameNum, Frame& frameIn) { * -frameNum to this peak's data set in peakFrameData_. */ typedef std::vector Barray; - Barray occupied(peaks_.size(), false); - Barray doubled(peaks_.size(), false); // to avoid double-additions + Barray occupied(peaksData_->Size(), false); + Barray doubled(peaksData_->Size(), false); // to avoid double-additions for (Iarray::const_iterator it = resPeakNum_.begin(); it != resPeakNum_.end(); it++) { @@ -578,12 +586,12 @@ Action::RetType Action_Spam::DoSPAM(int frameNum, Frame& frameIn) { } } // Now loop through and add all non-occupied sites - for (unsigned int i = 0; i < peaks_.size(); i++) + for (unsigned int i = 0; i < peaksData_->Size(); i++) if (!occupied[i]) peakFrameData_[i].push_back(frameNum); // Now adjust the occupied vectors to only contain 'true' for sites we need to // analyze (i.e., make all doubled points 'unoccupied') - for (unsigned int i = 0; i < peaks_.size(); i++) + for (unsigned int i = 0; i < peaksData_->Size(); i++) if (doubled[i]) occupied[i] = false; t_occupy_.Stop(); @@ -591,7 +599,7 @@ Action::RetType Action_Spam::DoSPAM(int frameNum, Frame& frameIn) { // If we have to calculate energies, do that here if (calcEnergy_) { int peak; - int npeaks = (int)peaks_.size(); + int npeaks = (int)peaksData_->Size(); const double ZERO = 0.0; // Loop through every peak, then loop through the water molecules to find // which one is in that site, and calculate the LJ and EEL energies for that @@ -611,6 +619,8 @@ Action::RetType Action_Spam::DoSPAM(int frameNum, Frame& frameIn) { // provide some time savings. double ene = Calculate_Energy(frameIn, solvent_residues_[i]); myDSL_[peak]->Add(frameNum, &ene); + //mprintf("DEBUG: Frm %6i peak %6i residx %6u solvres %6i ene %g\n", + // frameNum+1, peak+1, i, solvent_residues_[i].OriginalResNum(), ene); break; } } else @@ -628,7 +638,7 @@ Action::RetType Action_Spam::DoSPAM(int frameNum, Frame& frameIn) { /* Loop over every occupied site and swap the atoms so the same solvent * residue is always in the same site */ - for (int i = 0; i < (int)peaks_.size(); i++) { + for (int i = 0; i < (int)peaksData_->Size(); i++) { // Skip unoccupied sites if (!occupied[i]) continue; for (unsigned int j = 0; j < solvent_residues_.size(); j++) { @@ -700,16 +710,33 @@ int Action_Spam::Calc_G_Wat(DataSet* dsIn, unsigned int peaknum) // Calculate distribution of energy values using KDE. Get the bandwidth // factor here since we already know the SD. double BWfac = KDE::BandwidthFactor( enevec.Size() ); - //mprintf("DEBUG:\tNvals=%zu min=%g max=%g BWfac=%g\n", enevec.Size(), min, max, BWfac); + if (debug_ > 0) + mprintf("DEBUG:\tNvals=%zu min=%g max=%g BWfac=%g\n", enevec.Size(), min, max, BWfac); // Estimate number of bins the same way spamstats.py does. int nbins = (int)(((max - min) / BWfac) + 0.5) + 100; + if (nbins < 0) { + // Probably an overflow due to extremely large energy. + mprintf("Warning: Large magnitude energy observed for peak %u (min=%g max=%g)\n", + peaknum+1, min, max); + mprintf("Warning: Skipping peak.\n"); + return -1; + } HistBin Xdim(nbins, min - (50*BWfac), BWfac, "P(Ewat)"); //Xdim.CalcBinsOrStep(min - Havg.variance(), max + Havg.variance(), 0.0, nbins, "P(Ewat)"); - if (debug_ > 0) Xdim.PrintHistBin(); + if (debug_ > 0) { + mprintf("DEBUG:"); + Xdim.PrintHistBin(); + } DataSet_double kde1; KDE gkde; - if (gkde.CalcKDE( kde1, enevec, Xdim, 1.06 * sqrt(Havg.variance()) * BWfac )) { + double bandwidth; + if (enevec.Size() == 1) { + // Special case. Juse use BWfac to avoid a zero bandwidth. + bandwidth = BWfac; + } else + bandwidth = 1.06 * sqrt(Havg.variance()) * BWfac; + if (gkde.CalcKDE( kde1, enevec, Xdim, bandwidth )) { mprinterr("Error: Could not calculate E KDE histogram.\n"); return -1; } @@ -722,23 +749,26 @@ int Action_Spam::Calc_G_Wat(DataSet* dsIn, unsigned int peaknum) double Ewat = kde1.Xcrd(i); double PEwat = kde1.Dval(i); sumQ += (PEwat * exp( -Ewat * KB )); + //mprintf("DEBUG:\t\tEwat %20.10E PEwat %20.10E sumQ %20.10E\n", Ewat, PEwat, sumQ); } - //mprintf("DEBUG: sumQ= %20.10E\n", sumQ); + if (debug_ > 0) + mprintf("DEBUG: peak %6u sumQ= %20.10E\n", peaknum+1, sumQ); double DG = -RT * log(BWfac * sumQ); double adjustedDG = DG - DG_BULK_; double adjustedDH = Havg.mean() - DH_BULK_; double ntds = adjustedDG - adjustedDH; + if (ds_dg_ == 0) { - mprintf("\tSPAM bulk energy values:\n"); - mprintf("\t = %g, = %g +/- %g, -TdS= %g\n", adjustedDG, adjustedDH, + mprintf("\tSPAM bulk energy values:\n" + "\t = %g, = %g +/- %g, -TdS= %g\n", adjustedDG, adjustedDH, sqrt(Havg.variance()), ntds); } else { ((DataSet_Mesh*)ds_dg_)->AddXY(peaknum+1, adjustedDG); ((DataSet_Mesh*)ds_dh_)->AddXY(peaknum+1, adjustedDH); ((DataSet_Mesh*)ds_ds_)->AddXY(peaknum+1, ntds); } - + // DEBUG if (debug_ > 1) { FileName rawname("dbgraw." + integerToString(peaknum+1) + ".dat"); @@ -819,8 +849,8 @@ void Action_Spam::Print() { mprinterr("Warning: SPAM: Some frames had a box too small for the cutoff.\n"); // Print information about each missing peak - infofile_->Printf("# There are %d density peaks and %d frames\n\n", - (int)peaks_.size(), Nframes_); + infofile_->Printf("# There are %zu density peaks and %d frames\n\n", + peaksData_->Size(), Nframes_); // Loop over every Data set for (unsigned int i = 0; i < peakFrameData_.size(); i++) { // Skip peaks with 0 unoccupied sites @@ -850,7 +880,7 @@ void Action_Spam::Print() { if (err == 1) n_peaks_no_energy++; else if (err == -1) - mprintf("Warning: Error calculating SPAM energies for peak %zu\n", ds - myDSL_.begin()); + mprintf("Warning: Error calculating SPAM energies for peak %u\n", p + 1); } if (n_peaks_no_energy > 0) mprintf("Warning: No energies for %i peaks.\n", n_peaks_no_energy); diff --git a/src/Action_Spam.h b/src/Action_Spam.h index 491e9574f1..ffbc5a7758 100644 --- a/src/Action_Spam.h +++ b/src/Action_Spam.h @@ -5,6 +5,8 @@ #include "Vec3.h" #include "Timer.h" #include "PairList.h" +// Forward declares +class DataSet_Vector_Scalar; /** SPAM is a water profiling technique developed by Guanglei Cui at GlaxoSmithKline (GSK). The original implementation involved a set of specialized @@ -52,6 +54,7 @@ class Action_Spam: public Action { Action::RetType DoPureWater(int, Frame const&); Action::RetType DoSPAM(int, Frame&); + int GetPeaks(std::string const&, DataSetList const&); typedef bool (Action_Spam::*FxnType)(Vec3, Vec3, double) const; bool inside_box(Vec3, Vec3, double) const; bool inside_sphere(Vec3, Vec3, double) const; @@ -86,11 +89,12 @@ class Action_Spam: public Action { DataSet* ds_ds_; ///< Hold final -T*S values for each peak Parray peakFrameData_; ///< A list of all omitted frames for each peak DSarray myDSL_; ///< Hold energy data sets - Varray peaks_; ///< List of each peak location Varray comlist_; ///< For given frame, each residue C.O.M. coords. Rarray solvent_residues_; ///< List of each solvent residue int Nframes_; ///< Total number of frames bool overflow_; ///< True if cutoff overflowed our box coordinates + DataSetList peaksdsl_; ///< Will allocate DataSet for peaks data if loading from a file. + DataSet_Vector_Scalar* peaksData_; ///< Hold peaks DataSet // Timers Timer t_action_; Timer t_resCom_; diff --git a/src/DataFile.cpp b/src/DataFile.cpp index 3fefcf445a..6074f00472 100644 --- a/src/DataFile.cpp +++ b/src/DataFile.cpp @@ -28,6 +28,7 @@ #include "DataIO_CharmmOutput.h" #include "DataIO_Cpout.h" #include "DataIO_CharmmRtfPrm.h" +#include "DataIO_Peaks.h" // CONSTRUCTOR DataFile::DataFile() : @@ -73,6 +74,7 @@ const FileTypes::AllocToken DataFile::DF_AllocArray[] = { { "CHARMM Output", 0, 0, DataIO_CharmmOutput::Alloc}, { "Amber CPOUT", DataIO_Cpout::ReadHelp, DataIO_Cpout::WriteHelp, DataIO_Cpout::Alloc}, { "CHARMM RTF/PRM", 0, 0, DataIO_CharmmRtfPrm::Alloc }, + { "Peaks", 0, 0, DataIO_Peaks::Alloc }, { "Unknown Data file", 0, 0, 0 } }; @@ -94,6 +96,7 @@ const FileTypes::KeyToken DataFile::DF_KeyArray[] = { { CHARMMREPD, "charmmrepd",".exch" }, { CHARMMOUT, "charmmout", ".charmmout"}, { CHARMMRTFPRM, "charmmrtfprm", ".rtfprm"}, + { PEAKS, "peaks", ".peaks" }, { UNKNOWN_DATA, 0, 0 } }; @@ -113,6 +116,7 @@ const FileTypes::KeyToken DataFile::DF_WriteKeyArray[] = { { NCCMATRIX, "nccmatrix", ".nccmatrix" }, { CPOUT, "cpout", ".cpout" }, { CHARMMRTFPRM, "charmmrtfprm", ".prm" }, + { PEAKS, "peaks", ".peaks" }, { UNKNOWN_DATA, 0, 0 } }; @@ -268,6 +272,7 @@ int DataFile::SetupDatafile(FileName const& fnameIn, ArgList& argIn, // Set up DataIO based on format. dataio_ = (DataIO*)FileTypes::AllocIO( DF_AllocArray, dfType_, false ); if (dataio_ == 0) return Error("Error: Data file allocation failed.\n"); + dataio_->SetDebug( debug_ ); # ifdef MPI // Default to TrajComm master can write. threadCanWrite_ = Parallel::TrajComm().Master(); diff --git a/src/DataFile.h b/src/DataFile.h index 4da2653e06..0517c00774 100644 --- a/src/DataFile.h +++ b/src/DataFile.h @@ -17,7 +17,7 @@ class DataFile { enum DataFormatType { DATAFILE=0, XMGRACE, GNUPLOT, XPLOR, OPENDX, REMLOG, MDOUT, EVECS, VECTRAJ, XVG, CCP4, CMATRIX, NCCMATRIX, CHARMMREPD, CHARMMFASTREP, - CHARMMOUT, CPOUT, CHARMMRTFPRM, UNKNOWN_DATA + CHARMMOUT, CPOUT, CHARMMRTFPRM, PEAKS, UNKNOWN_DATA }; DataFile(); ~DataFile(); diff --git a/src/DataIO_OpenDx.cpp b/src/DataIO_OpenDx.cpp index bbbd63662b..74d4b0d973 100644 --- a/src/DataIO_OpenDx.cpp +++ b/src/DataIO_OpenDx.cpp @@ -2,10 +2,25 @@ #include // atof #include "DataIO_OpenDx.h" #include "CpptrajStdio.h" -#include "DataSet_GridFlt.h" +#include "DataSet_GridFlt.h" // Included to get default format +#include "DataSet_GridDbl.h" // Included to get default format #include "BufferedLine.h" #include "ProgressBar.h" +DataIO_OpenDx::DataIO_OpenDx() : + DataIO(false, false, true), // Valid for 3D only + gridWriteMode_(BIN_CORNER), + gridReadType_(DataSet::GRID_FLT) +{ + // Get the default formats for float/double grid datasets + DataSet* ds = new DataSet_GridFlt(); + fmt_gridFlt_ = ds->Format(); + delete ds; + ds = new DataSet_GridDbl(); + fmt_gridDbl_ = ds->Format(); + delete ds; +} + bool DataIO_OpenDx::ID_DataFormat( CpptrajFile& infile ) { bool isDX = false; if (!infile.OpenFile()) { @@ -17,13 +32,29 @@ bool DataIO_OpenDx::ID_DataFormat( CpptrajFile& infile ) { return isDX; } +/** Process read options. */ +int DataIO_OpenDx::processReadArgs(ArgList& argIn) { + std::string typestr = argIn.GetStringKey("type"); + if (!typestr.empty()) { + if (typestr == "float") + gridReadType_ = DataSet::GRID_FLT; + else if (typestr == "double") + gridReadType_ = DataSet::GRID_DBL; + else { + mprinterr("Error: Unrecognized grid type: %s\n", typestr.c_str()); + return 1; + } + } + return 0; +} + // DataIO_OpenDx::ReadData() int DataIO_OpenDx::ReadData(FileName const& fname, DataSetList& datasetlist, std::string const& dsname) { // TODO append? // Add grid data set. Default to float for now. - DataSet* ds = datasetlist.AddSet( DataSet::GRID_FLT, dsname, "GRID" ); + DataSet* ds = datasetlist.AddSet( gridReadType_, dsname, "GRID" ); if (ds==0) return 1; if (LoadGrid(fname.full(), *ds)) { // Load failed. Erase grid data set. @@ -36,8 +67,7 @@ int DataIO_OpenDx::ReadData(FileName const& fname, // DataIO_OpenDx::LoadGrid() int DataIO_OpenDx::LoadGrid(const char* filename, DataSet& ds) { - // TODO: This may need to be changed if new 3D types introduced. - DataSet_GridFlt& grid = static_cast( ds ); + DataSet_3D& grid = static_cast( ds ); // Open file BufferedLine infile; if (infile.OpenFileRead(filename)) return 1; @@ -145,7 +175,7 @@ int DataIO_OpenDx::LoadGrid(const char* filename, DataSet& ds) mprintf("Warning: Check that data region ends with a newline.\n"); break; } - grid[ndata++] = (float)atof(infile.NextToken()); + grid.SetGrid(ndata++, atof(infile.NextToken())); } progress.Update( ndata ); } @@ -160,6 +190,7 @@ void DataIO_OpenDx::WriteHelp() { "\tgridext: Like 'bincenter', but also print extra layer of empty bins.\n"); } +/** Process opendx write options. */ int DataIO_OpenDx::processWriteArgs(ArgList& argIn) { if (argIn.hasKey("bincenter")) gridWriteMode_ = BIN_CENTER; else if (argIn.hasKey("gridwrap")) gridWriteMode_ = WRAP; @@ -219,6 +250,7 @@ int DataIO_OpenDx::WriteSet3D(DataSet const& setIn, CpptrajFile& outfile) const return err; } +/** Write the header for an OpenDX file. */ void DataIO_OpenDx::WriteDxHeader(CpptrajFile& outfile, size_t NX, size_t NY, size_t NZ, double LX, double LY, double LZ, @@ -235,8 +267,36 @@ void DataIO_OpenDx::WriteDxHeader(CpptrajFile& outfile, NX, NY, NZ, NX*NY*NZ); } +/** Create a text format string using the given TextFormat. If the default + * %12.4f is in use, upgrade it to a general format. + */ +std::string DataIO_OpenDx::CreateFmtString(DataSet::DataType typeIn, + TextFormat const& fmtIn) const +{ + TextFormat fmtOut; + bool usingDefaultFmt = false; + if ( typeIn == DataSet::GRID_FLT ) { + usingDefaultFmt = (fmtIn == fmt_gridFlt_); + } else if ( typeIn == DataSet::GRID_DBL ) { + usingDefaultFmt = (fmtIn == fmt_gridDbl_); + } else { + mprinterr("Internal Error: Unhandled type in CreateFmtString(): %s\n", + DataSet::description(typeIn)); + return std::string(); + } + if (usingDefaultFmt) + fmtOut = TextFormat(TextFormat::GDOUBLE, -1, -1); + else + fmtOut = TextFormat(fmtIn.FormatType(), fmtIn.Width(), fmtIn.Precision()); + return fmtOut.Fmt(); +} + +/** Write the grid wrapped. */ int DataIO_OpenDx::WriteGridWrap(DataSet const& setIn, CpptrajFile& outfile) const { DataSet_3D const& set = static_cast( setIn ); + std::string myFmt = CreateFmtString( setIn.Type(), setIn.Format() ); + if (myFmt.empty()) return 1; + myFmt = " " + myFmt; // Need to construct a grid mesh around bins, with points centered on the bins. int mesh_x = set.NX(); int mesh_y = set.NY(); @@ -263,7 +323,7 @@ int DataIO_OpenDx::WriteGridWrap(DataSet const& setIn, CpptrajFile& outfile) con if (ik < 0 ) bk = mesh_z - 1; else if (ik == mesh_z) bk = 0; else bk = ik; - outfile.Printf(" %g", set.GetElement(bi, bj, bk)); + outfile.Printf(myFmt.c_str(), set.GetElement(bi, bj, bk)); ++nvals; if (nvals == 5) { outfile.Printf("\n"); @@ -281,7 +341,7 @@ int DataIO_OpenDx::WriteGridWrap(DataSet const& setIn, CpptrajFile& outfile) con if (zero_x || zero_y || ik < 0 || ik == mesh_z) outfile.Printf(" 0"); else - outfile.Printf(" %g", set.GetElement(ii, ij, ik)); + outfile.Printf(myFmt.c_str(), set.GetElement(ii, ij, ik)); ++nvals; if (nvals == 5) { outfile.Printf("\n"); @@ -295,8 +355,21 @@ int DataIO_OpenDx::WriteGridWrap(DataSet const& setIn, CpptrajFile& outfile) con return 0; } +/** Write the grid normally. */ int DataIO_OpenDx::WriteGrid(DataSet const& setIn, CpptrajFile& outfile) const { DataSet_3D const& set = static_cast( setIn ); + std::string myFmt = CreateFmtString( setIn.Type(), setIn.Format() ); + if (myFmt.empty()) return 1; + std::string myFmt2 = myFmt + " " + myFmt; + std::string myFmt3 = myFmt + " " + myFmt + " " + myFmt; + if (debug_ > 0) { + mprintf("DEBUG: Grid data %s format is \"%s\"\n", set.legend(), set.Format().fmt()); + mprintf("DEBUG: New formats are \"%s\" \"%s\" \"%s\"\n", + myFmt.c_str(), myFmt2.c_str(), myFmt3.c_str()); + } + myFmt.append("\n"); + myFmt2.append("\n"); + myFmt3.append("\n"); Vec3 oxyz = set.Bin().GridOrigin(); if (gridWriteMode_ == BIN_CENTER) // Origin needs to be shifted to center of bin located at 0,0,0 @@ -307,18 +380,18 @@ int DataIO_OpenDx::WriteGrid(DataSet const& setIn, CpptrajFile& outfile) const { // Now print out the data. size_t gridsize = set.Size(); if (gridsize == 1) - outfile.Printf("%g\n", set[0]); + outfile.Printf(myFmt.c_str(), set[0]); else if (gridsize == 2) - outfile.Printf("%g %g\n", set[0], set[1]); + outfile.Printf(myFmt2.c_str(), set[0], set[1]); else if (gridsize > 2) { // Data is already in row-major form (z-axis changes // fastest), so no need to do any kind of data adjustment for (size_t i = 0UL; i < gridsize - 2UL; i += 3UL) - outfile.Printf("%g %g %g\n", set[i], set[i+1], set[i+2]); + outfile.Printf(myFmt3.c_str(), set[i], set[i+1], set[i+2]); // Print out any points we may have missed switch (gridsize % 3) { - case 2: outfile.Printf("%g %g\n", set[gridsize-2], set[gridsize-1]); break; - case 1: outfile.Printf("%g\n", set[gridsize-1]); break; + case 2: outfile.Printf(myFmt2.c_str(), set[gridsize-2], set[gridsize-1]); break; + case 1: outfile.Printf(myFmt.c_str(), set[gridsize-1]); break; } } return 0; diff --git a/src/DataIO_OpenDx.h b/src/DataIO_OpenDx.h index 306b50dd6a..cae0cfcf00 100644 --- a/src/DataIO_OpenDx.h +++ b/src/DataIO_OpenDx.h @@ -4,9 +4,9 @@ /// Write OpenDx format data files. class DataIO_OpenDx : public DataIO { public: - DataIO_OpenDx() : DataIO(false, false, true), gridWriteMode_(BIN_CORNER) {} // Valid for 3D only + DataIO_OpenDx(); static BaseIOtype* Alloc() { return (BaseIOtype*)new DataIO_OpenDx(); } - int processReadArgs(ArgList&) { return 0; } + int processReadArgs(ArgList&); int ReadData(FileName const&, DataSetList&, std::string const&); static void WriteHelp(); int processWriteArgs(ArgList&); @@ -18,9 +18,13 @@ class DataIO_OpenDx : public DataIO { int WriteSet3D(DataSet const&, CpptrajFile&) const; void WriteDxHeader(CpptrajFile&, size_t, size_t, size_t, double, double, double, Matrix_3x3 const&, Vec3 const&) const; + std::string CreateFmtString(DataSet::DataType, TextFormat const&) const; int WriteGridWrap(DataSet const&, CpptrajFile&) const; int WriteGrid(DataSet const&, CpptrajFile&) const; GridWriteType gridWriteMode_; + DataSet::DataType gridReadType_; + TextFormat fmt_gridFlt_; ///< Format for DataSet_GridFlt + TextFormat fmt_gridDbl_; ///< Format for DataSet_GridDbl }; #endif diff --git a/src/DataIO_Peaks.cpp b/src/DataIO_Peaks.cpp new file mode 100644 index 0000000000..3bf7c5eabc --- /dev/null +++ b/src/DataIO_Peaks.cpp @@ -0,0 +1,166 @@ +#include "DataIO_Peaks.h" +#include "CpptrajStdio.h" +#include "DataSet_Vector_Scalar.h" +#include "StringRoutines.h" +#include // sscanf + +/// CONSTRUCTOR +DataIO_Peaks::DataIO_Peaks() +{ + SetValid( DataSet::VECTOR_SCALAR ); +} + +// DataIO_Peaks::ID_DataFormat() +bool DataIO_Peaks::ID_DataFormat(CpptrajFile& infile) +{ + if (infile.OpenFile()) return false; + std::string firstLine = NoTrailingWhitespace(infile.GetLine()); + std::string secondLine = NoTrailingWhitespace(infile.GetLine()); + std::string thirdLine = infile.GetLine(); + //mprintf("DEBUG: ID peaks.\n"); + //mprintf("'%s'\n", firstLine.c_str()); + //mprintf("'%s'\n", secondLine.c_str()); + //mprintf("'%s'\n", thirdLine.c_str()); + infile.CloseFile(); + if (validInteger(firstLine)) { + //mprintf("DEBUG: Line 1 int.\n"); + if (secondLine.empty()) { + //mprintf("DEBUG: Line 2 empty.\n"); + ArgList line3(thirdLine); + if (line3.Nargs() == 5) { + if (line3[0] == "C") { + for (int col = 1; col < 5; col++) + if (!validDouble(line3[col])) return false; + //mprintf("DEBUG: Line 3 C and 4 doubles.\n"); + return true; + } + } + } + } + return false; +} + +// DataIO_Peaks::ReadHelp() +void DataIO_Peaks::ReadHelp() +{ + +} + +// DataIO_Peaks::processReadArgs() +int DataIO_Peaks::processReadArgs(ArgList& argIn) +{ + + return 0; +} + +// DataIO_Peaks::ReadData() +int DataIO_Peaks::ReadData(FileName const& fname, DataSetList& dsl, std::string const& dsname) +{ + // Figure out the data set + unsigned int peakidx = 0; + DataSet* ds = dsl.CheckForSet( dsname ); + if (ds == 0) { + // New set + ds = dsl.AddSet( DataSet::VECTOR_SCALAR, dsname, "peaks" ); + if (ds == 0) return 1; + } else { + // Appending. + if (ds->Type() != DataSet::VECTOR_SCALAR) { + mprinterr("Error: Set '%s' is not vector with scalar, cannot append peaks to it.\n", + ds->legend()); + return 1; + } + peakidx = ds->Size(); + mprintf("\tAppending peaks to set '%s', offset %u\n", ds->legend(), peakidx); + } + // Parse through the peaks file and extract the peaks + CpptrajFile peakfile; + if (peakfile.OpenRead(fname)) { + mprinterr("Error: Could not open %s for reading!\n", fname.full()); + return 1; + } + // FORMAT: + // Line 1 : # peaks + int npeaks_in_file = 0; + const char* ptr = peakfile.NextLine(); + if (ptr == 0) { + mprinterr("Error: Unexpected EOF when reading # peaks..\n"); + return 1; + } + if (sscanf(ptr, "%i", &npeaks_in_file) != 1) { + mprinterr("Error: Could not read number of peaks in file.\n"); + return 1; + } + mprintf("\tAttempting to read %i peaks.\n", npeaks_in_file); + // Line 2 : Blank + ptr = peakfile.NextLine(); + // Line 3-X : C + int npeaks_read = 0; + ptr = peakfile.NextLine(); + double pbuf[4]; + while (ptr != 0) { + if (sscanf(ptr, "C %lg %lg %lg %lg", pbuf, pbuf+1, pbuf+2, pbuf+3) != 4) { + mprintf("Warning: Unexpected format for line, skipping: '%s'\n", + NoTrailingWhitespace(std::string(ptr)).c_str()); + } else { + npeaks_read++; + ds->Add(peakidx++, pbuf); + } + ptr = peakfile.NextLine(); + } + peakfile.CloseFile(); + mprintf("\tRead %i peaks.\n", npeaks_read); + if (npeaks_read != npeaks_in_file) + mprintf("Warning: Expected %i peaks, actually got %i peaks.\n", npeaks_in_file, npeaks_read); + + return 0; +} + +// DataIO_Peaks::WriteHelp() +void DataIO_Peaks::WriteHelp() +{ + +} + +// DataIO_Peaks::processWriteArgs() +int DataIO_Peaks::processWriteArgs(ArgList& argIn) +{ + + return 0; +} + +// DataIO_Peaks::WriteData() +int DataIO_Peaks::WriteData(FileName const& fname, DataSetList const& dsl) +{ + if (dsl.size() > 1) + mprintf("Warning: Writing multiple sets to peak file may result in invalid format.\n"); + CpptrajFile outfile; + // Loop over sets. Only write VECTOR_SCALAR sets. + for (DataSetList::const_iterator it = dsl.begin(); it != dsl.end(); ++it) + { + if ((*it)->Type() != DataSet::VECTOR_SCALAR) { + mprintf("Warning: Set '%s' is not vector with scalar, cannot be used for peaks file.\n", + (*it)->legend()); + } else { + DataSet_Vector_Scalar const& ds = static_cast( *(*it) ); + if (ds.Size() > 0) { + // Only open the file when a valid set is found to match old Action_Volmap behavior. + if (!outfile.IsOpen()) { + if (outfile.OpenWrite( fname )) { + mprinterr("Error: Could not open %s for write.\n", fname.full()); + return 1; + } + } + outfile.Printf("%zu\n\n", ds.Size()); + for (unsigned int i = 0; i < ds.Size(); i++) + { + Vec3 const& vxyz = ds.Vec(i); + outfile.Printf("C %16.8f %16.8f %16.8f %16.8f\n", + vxyz[0], vxyz[1], vxyz[2], ds.Val(i)); + } + } + } + } // END loop over sets. + + return 0; +} diff --git a/src/DataIO_Peaks.h b/src/DataIO_Peaks.h new file mode 100644 index 0000000000..07ecd3ca3a --- /dev/null +++ b/src/DataIO_Peaks.h @@ -0,0 +1,17 @@ +#ifndef INC_DATAIO_PEAKS_H +#define INC_DATAIO_PEAKS_H +#include "DataIO.h" +/// +class DataIO_Peaks : public DataIO { + public: + DataIO_Peaks(); + static void ReadHelp(); + static void WriteHelp(); + static BaseIOtype* Alloc() { return (BaseIOtype*)new DataIO_Peaks(); } + int processReadArgs(ArgList&); + int ReadData(FileName const&, DataSetList&, std::string const&); + int processWriteArgs(ArgList&); + int WriteData(FileName const&, DataSetList const&); + bool ID_DataFormat(CpptrajFile&); +}; +#endif diff --git a/src/DataSet.cpp b/src/DataSet.cpp index c462a52007..f054fcafd8 100644 --- a/src/DataSet.cpp +++ b/src/DataSet.cpp @@ -30,7 +30,8 @@ const char* DataSet::Descriptions_[] = { "pH REMD (implicit)", // PH_IMPL "parameters", // PARAMETERS "tensor", // TENSOR - "string variable" // STRINGVAR + "string variable", // STRINGVAR + "vector with scalar" // VECTOR_SCALAR }; // CONSTRUCTOR diff --git a/src/DataSet.h b/src/DataSet.h index 152141823d..76a8ec7c0f 100644 --- a/src/DataSet.h +++ b/src/DataSet.h @@ -28,7 +28,7 @@ class DataSet { UNKNOWN_DATA=0, DOUBLE, FLOAT, INTEGER, STRING, MATRIX_DBL, MATRIX_FLT, COORDS, VECTOR, MODES, GRID_FLT, GRID_DBL, REMLOG, XYMESH, TRAJ, REF_FRAME, MAT3X3, TOPOLOGY, CMATRIX, CMATRIX_NOMEM, CMATRIX_DISK, PH, PH_EXPL, PH_IMPL, - PARAMETERS, TENSOR, STRINGVAR + PARAMETERS, TENSOR, STRINGVAR, VECTOR_SCALAR }; /// Group DataSet belongs to. enum DataGroup { diff --git a/src/DataSetList.cpp b/src/DataSetList.cpp index c96e3b545c..413a411c7e 100644 --- a/src/DataSetList.cpp +++ b/src/DataSetList.cpp @@ -33,6 +33,7 @@ #include "DataSet_Parameters.h" #include "DataSet_Tensor.h" #include "DataSet_StringVar.h" +#include "DataSet_Vector_Scalar.h" bool DataSetList::useDiskCache_ = false; @@ -81,6 +82,7 @@ DataSet* DataSetList::NewSet(DataSet::DataType typeIn) { case DataSet::PARAMETERS : ds = DataSet_Parameters::Alloc(); break; case DataSet::TENSOR : ds = DataSet_Tensor::Alloc(); break; case DataSet::STRINGVAR : ds = DataSet_StringVar::Alloc(); break; + case DataSet::VECTOR_SCALAR : ds = DataSet_Vector_Scalar::Alloc(); break; // Sanity check default: mprinterr("Internal Error: No allocator for DataSet type '%s'\n", diff --git a/src/DataSet_3D.h b/src/DataSet_3D.h index 0070a4840a..7967f48772 100644 --- a/src/DataSet_3D.h +++ b/src/DataSet_3D.h @@ -1,9 +1,8 @@ #ifndef INC_DATASET_3D_H #define INC_DATASET_3D_H #include "DataSet.h" -#include "CpptrajFile.h" -#include "Box.h" #include "GridBin.h" +class Box; /// Interface for 3D DataSets. // FIXME: Use DataSet Dims? class DataSet_3D : public DataSet { @@ -19,6 +18,10 @@ class DataSet_3D : public DataSet { int Allocate(SizeArray const&) { return 1; } // TODO enable? /// \return Data from grid at x/y/z point. virtual double GetElement(size_t, size_t, size_t) const = 0; + /// Set grid to value + virtual void SetGrid(size_t, double) = 0; + /// Increment the specified grid point by value + virtual void IncrementElement(size_t, size_t, size_t, double) = 0; /// \return Data from grid. virtual double operator[](size_t) const = 0; /// \return size of X dimension. @@ -33,6 +36,8 @@ class DataSet_3D : public DataSet { virtual void ReverseIndex(long int, size_t&, size_t&, size_t&) const = 0; /// Increment specified voxel by given amount. virtual void UpdateVoxel(long int, double) = 0; + /// Divide all elements by the given scalar + virtual void operator/=(double) = 0; // ------------------------------------------- // TODO: Remove this. Only needed by DataSet_1D.h void Add(size_t,const void*) { } diff --git a/src/DataSet_GridDbl.h b/src/DataSet_GridDbl.h index 3fb10cd26a..4ebe9614d3 100644 --- a/src/DataSet_GridDbl.h +++ b/src/DataSet_GridDbl.h @@ -11,6 +11,18 @@ class DataSet_GridDbl : public DataSet_3D { double& operator[](size_t idx) { return grid_[idx]; } static DataSet* Alloc() { return (DataSet*)new DataSet_GridDbl(); } Grid const& InternalGrid() const { return grid_; } + void SetElement(size_t x,size_t y,size_t z,double v) { grid_.setGrid(x,y,z,v); } + /// Type definition of iterator over grid elements. + typedef Grid::iterator iterator; + iterator begin() { return grid_.begin(); } + iterator end() { return grid_.end(); } + /// Increment grid bin corresponding to point by given value. + inline long int Increment(Vec3 const&, double); + inline long int Increment(const double*, double); + /// Increment grid bin by given value. + inline long int Increment(size_t,size_t,size_t,double); + /// \return grid value at specified bin. + double GridVal(size_t x,size_t y,size_t z) const { return grid_.element(x,y,z); } // ----- DataSet functions ------------------- size_t Size() const { return grid_.size(); } # ifdef MPI @@ -23,6 +35,8 @@ class DataSet_GridDbl : public DataSet_3D { // ----- DataSet_3D functions ---------------- int Allocate3D(size_t x,size_t y,size_t z) { return grid_.resize(x,y,z); } double GetElement(size_t x,size_t y,size_t z) const { return grid_.element(x,y,z); } + void SetGrid(size_t idx, double d) { grid_[idx] = d; } + void IncrementElement(size_t x, size_t y, size_t z, double val) { Increment(x,y,z,val); } double operator[](size_t idx) const { return grid_[idx]; } size_t NX() const { return grid_.NX(); } size_t NY() const { return grid_.NY(); } @@ -32,19 +46,11 @@ class DataSet_GridDbl : public DataSet_3D { void ReverseIndex(long int n, size_t& i, size_t& j, size_t& k) const { return grid_.ReverseIndex(n,i,j,k); } void UpdateVoxel(long int i, double val) { grid_[i] += val; } - // ------------------------------------------- - void SetElement(size_t x,size_t y,size_t z,double v) { grid_.setGrid(x,y,z,v); } - /// Type definition of iterator over grid elements. - typedef Grid::iterator iterator; - iterator begin() { return grid_.begin(); } - iterator end() { return grid_.end(); } - /// Increment grid bin corresponding to point by given value. - inline long int Increment(Vec3 const&, double); - inline long int Increment(const double*, double); - /// Increment grid bin by given value. - inline long int Increment(size_t,size_t,size_t,double); - /// \return grid value at specified bin. - double GridVal(size_t x,size_t y,size_t z) const { return grid_.element(x,y,z); } + /// Divide all elements by val + void operator/=(double val) { + for (Grid::iterator it = grid_.begin(); it != grid_.end(); ++it) + (*it) /= val; + } private: Grid grid_; }; diff --git a/src/DataSet_GridFlt.h b/src/DataSet_GridFlt.h index 4f8ca9aedf..cb7169b364 100644 --- a/src/DataSet_GridFlt.h +++ b/src/DataSet_GridFlt.h @@ -11,6 +11,18 @@ class DataSet_GridFlt : public DataSet_3D { float& operator[](size_t idx) { return grid_[idx]; } static DataSet* Alloc() { return (DataSet*)new DataSet_GridFlt(); } Grid const& InternalGrid() const { return grid_; } + void SetElement(size_t x,size_t y,size_t z,float v) { grid_.setGrid(x,y,z,v); } + /// Type definition of iterator over grid elements. + typedef Grid::iterator iterator; + iterator begin() { return grid_.begin(); } + iterator end() { return grid_.end(); } + /// Increment grid bin corresponding to point by given value. + inline long int Increment(Vec3 const&, float); + inline long int Increment(const double*, float); + /// Increment grid bin by given value. + inline long int Increment(size_t,size_t,size_t,float); + /// \return grid value at specified bin. + float GridVal(size_t x,size_t y,size_t z) const { return grid_.element(x,y,z); } // ----- DataSet functions ------------------- size_t Size() const { return grid_.size(); } # ifdef MPI @@ -23,6 +35,8 @@ class DataSet_GridFlt : public DataSet_3D { // ----- DataSet_3D functions ---------------- int Allocate3D(size_t x, size_t y, size_t z) { return grid_.resize(x,y,z); } double GetElement(size_t x, size_t y, size_t z) const { return (double)grid_.element(x,y,z); } + void SetGrid(size_t idx, double d) { grid_[idx] = d; } + void IncrementElement(size_t x, size_t y, size_t z, double val) { Increment(x, y, z, (float)val); } double operator[](size_t idx) const { return (double)grid_[idx]; } size_t NX() const { return grid_.NX(); } size_t NY() const { return grid_.NY(); } @@ -32,19 +46,13 @@ class DataSet_GridFlt : public DataSet_3D { void ReverseIndex(long int n, size_t& i, size_t& j, size_t& k) const { return grid_.ReverseIndex(n,i,j,k); } void UpdateVoxel(long int i, double val) { grid_[i] += (float)val; } - // ------------------------------------------- - void SetElement(size_t x,size_t y,size_t z,float v) { grid_.setGrid(x,y,z,v); } - /// Type definition of iterator over grid elements. - typedef Grid::iterator iterator; - iterator begin() { return grid_.begin(); } - iterator end() { return grid_.end(); } - /// Increment grid bin corresponding to point by given value. - inline long int Increment(Vec3 const&, float); - inline long int Increment(const double*, float); - /// Increment grid bin by given value. - inline long int Increment(size_t,size_t,size_t,float); - /// \return grid value at specified bin. - float GridVal(size_t x,size_t y,size_t z) const { return grid_.element(x,y,z); } + /// Divide all elements by val + void operator/=(double val) { + float fval = (float)val; + for (Grid::iterator it = grid_.begin(); it != grid_.end(); ++it) + (*it) /= fval; + } + private: Grid grid_; }; diff --git a/src/DataSet_Vector.cpp b/src/DataSet_Vector.cpp index 694649981a..0041b31127 100644 --- a/src/DataSet_Vector.cpp +++ b/src/DataSet_Vector.cpp @@ -68,6 +68,7 @@ void DataSet_Vector::WriteBuffer(CpptrajFile &cbuffer, SizeArray const& pIn) con // DataSet_Vector::Append() int DataSet_Vector::Append(DataSet* dsIn) { + if (dsIn == 0) return 1; if (dsIn->Empty()) return 0; if (dsIn->Type() != VECTOR) return 1; Varray const& vIn = ((DataSet_Vector*)dsIn)->vectors_; diff --git a/src/DataSet_Vector_Scalar.cpp b/src/DataSet_Vector_Scalar.cpp new file mode 100644 index 0000000000..351eaf2cec --- /dev/null +++ b/src/DataSet_Vector_Scalar.cpp @@ -0,0 +1,80 @@ +#include "DataSet_Vector_Scalar.h" +#include "CpptrajStdio.h" +#include // std::copy + +/// CONSTRUCTOR TODO should this be in VECTOR_1D group? +DataSet_Vector_Scalar::DataSet_Vector_Scalar() : + // TODO marking this as 0-dimension so that it is forced to use + // the DataIO_Peaks format, but may want to change this down the + // road. + DataSet(VECTOR_SCALAR, GENERIC, TextFormat(TextFormat::DOUBLE, 8, 4, 4), 0) +{ } + +/** Reserve space in vector array. */ +int DataSet_Vector_Scalar::Allocate(SizeArray const& Nin) { + if (!Nin.empty()) { + vecs_.reserve( Nin[0] ); + vals_.reserve( Nin[0] ); + } + return 0; +} + +/** Write vector + scalar to output file. */ +void DataSet_Vector_Scalar::WriteBuffer(CpptrajFile& cbuffer, SizeArray const& pIn) const { + if (pIn[0] >= Size()) { + cbuffer.Printf(format_.fmt(), 0.0, 0.0, 0.0, 0.0); /// VXYZ SCALAR + } else { + Vec3 const& vxyz = vecs_[pIn[0]]; + cbuffer.Printf(format_.fmt(), vxyz[0], vxyz[1], vxyz[2], vals_[pIn[0]]); + } +} + +/** Append vector scalar set to this set. */ +int DataSet_Vector_Scalar::Append(DataSet* dsIn) { + if (dsIn == 0) return 1; + if (dsIn->Empty()) return 0; + if (dsIn->Type() != VECTOR_SCALAR) return 1; + Varray const& vIn = ((DataSet_Vector_Scalar*)dsIn)->vecs_; + Darray const& dIn = ((DataSet_Vector_Scalar*)dsIn)->vals_; + size_t oldsize = vecs_.size(); + vecs_.resize( oldsize + vIn.size() ); + std::copy( vIn.begin(), vIn.end(), vecs_.begin() + oldsize ); + vals_.resize( oldsize + dIn.size() ); + std::copy( dIn.begin(), dIn.end(), vals_.begin() + oldsize ); + return 0; +} + +/** \return memory usage of the set in bytes. */ +size_t DataSet_Vector_Scalar::MemUsageInBytes() const { + return ((vecs_.size()*Vec3::DataSize()) + vals_.size()); +} + +#ifdef MPI +int DataSet_Vector_Scalar::Sync(size_t total, std::vector const& rank_frames, + Parallel::Comm const& commIn) +{ + if (commIn.Size()==1) return 0; + double buf[4]; + if (commIn.Master()) { + // Resize to accept data from other ranks. + vecs_.resize( total ); + vals_.resize( total ); + int vidx = rank_frames[0]; + for (int rank = 1; rank < commIn.Size(); rank++) { + for (int ridx = 0; ridx != rank_frames[rank]; ridx++, vidx++) { + commIn.SendMaster( buf, 4, rank, MPI_DOUBLE ); + std::copy( buf, buf+3, vecs_[vidx].Dptr() ); + vals_[vidx] = buf[3]; + } + } + } else { + // Send data to master + for (unsigned int ridx = 0; ridx != vecs_.size(); ++ridx) { + std::copy( vecs_[ridx].Dptr(), vecs_[ridx].Dptr()+3, buf ); + buf[3] = vals_[ridx]; + commIn.SendMaster( buf, 4, commIn.Rank(), MPI_DOUBLE ); + } + } + return 0; +} +#endif diff --git a/src/DataSet_Vector_Scalar.h b/src/DataSet_Vector_Scalar.h new file mode 100644 index 0000000000..1f56fa9018 --- /dev/null +++ b/src/DataSet_Vector_Scalar.h @@ -0,0 +1,38 @@ +#ifndef INC_DATASET_VECTOR_SCALAR_H +#define INC_DATASET_VECTOR_SCALAR_H +#include "DataSet.h" +#include "Vec3.h" +/// Hold XYZ value and associated scalar. +class DataSet_Vector_Scalar : public DataSet { + public: + DataSet_Vector_Scalar(); + static DataSet* Alloc() { return (DataSet*)new DataSet_Vector_Scalar(); } + // ----- DataSet functions ------------------- + size_t Size() const { return vecs_.size(); } + void Info() const { return; } + int Allocate(SizeArray const&); + inline void Add(size_t, const void*); + void WriteBuffer(CpptrajFile&, SizeArray const&) const; + int Append(DataSet*); + size_t MemUsageInBytes() const; +# ifdef MPI + int Sync(size_t, std::vector const&, Parallel::Comm const&); +# endif + // ------------------------------------------- + Vec3 const& Vec(unsigned int i) const { return vecs_[i]; } + double Val(unsigned int i) const { return vals_[i]; } + private: + typedef std::vector Varray; + typedef std::vector Darray; + + Varray vecs_; ///< Hold XYZ values. + Darray vals_; ///< Hold scalar values. +}; + +/** Add XYZ scalar to the DataSet */ +void DataSet_Vector_Scalar::Add(size_t idx, const void* ptrIn) { + const double* ptr = (const double*)ptrIn; + vecs_.push_back( Vec3(ptr[0], ptr[1], ptr[2]) ); + vals_.push_back( ptr[3] ); +} +#endif diff --git a/src/GridBin.h b/src/GridBin.h index 93ae7b66ff..5e4307ef43 100644 --- a/src/GridBin.h +++ b/src/GridBin.h @@ -1,6 +1,7 @@ #ifndef INC_GRIDBIN_H #define INC_GRIDBIN_H #include "Matrix_3x3.h" +#include "Box.h" /// Class used to perform binning on/get voxel coords of 3D grids. class GridBin { public: diff --git a/src/KDE.cpp b/src/KDE.cpp index cddf7fa21a..33f544ef3f 100644 --- a/src/KDE.cpp +++ b/src/KDE.cpp @@ -179,7 +179,9 @@ int KDE::CalcKDE(DataSet_double& Out, DataSet_1D const& Pdata, omp_set_num_threads( original_num_threads ); # endif // Normalize - for (unsigned int j = 0; j < Out.Size(); j++) + for (unsigned int j = 0; j < Out.Size(); j++) { + //mprintf("DEBUG:\t\t\tkde %g %g %g\n", Out[j], total, bandwidth); Out[j] /= (total * bandwidth); + } return 0; } diff --git a/src/RPNcalc.cpp b/src/RPNcalc.cpp index 9fe5d3b061..d591a6022a 100644 --- a/src/RPNcalc.cpp +++ b/src/RPNcalc.cpp @@ -763,7 +763,11 @@ int RPNcalc::TokenLoop(DataSetList& DSL) const { { DataSet_3D const& G1 = static_cast( *ds1 ); DataSet_3D const& G2 = static_cast( *ds2 ); - if (T->Type() == OP_MINUS || T->Type() == OP_PLUS) { + if (T->Type() == OP_MINUS || + T->Type() == OP_PLUS || + T->Type() == OP_DIV || + T->Type() == OP_MULT) + { if (G1.NX() != G2.NX() || G1.NY() != G2.NY() || G1.NZ() != G2.NZ()) { mprinterr("Error: Grid operation '%s' requires both grids have" " same dimensions.\n", T->Description()); diff --git a/src/TextFormat.h b/src/TextFormat.h index 6fa91eb521..7998b3624a 100644 --- a/src/TextFormat.h +++ b/src/TextFormat.h @@ -51,9 +51,15 @@ class TextFormat { /// \return pointer to format string. const char* fmt() const { return fmt_.c_str(); } std::string const& Fmt() const { return fmt_; } + FmtType FormatType() const { return type_; } int Width() const { return width_; } int Precision() const { return precision_; } int ColumnWidth() const { return colwidth_; } + /// \return True if given TextFormat matches this format type, width, and precision. + bool operator==(TextFormat const& rhs) const { + return ( (type_ == rhs.type_) && (width_ == rhs.width_) && (precision_ == rhs.precision_) ); + } + /// \return Description of the given format type. static const char* typeDescription(FmtType f) { return TypeDesc_[f]; } private: static inline bool IsDoubleType(FmtType); diff --git a/src/TinkerFile.cpp b/src/TinkerFile.cpp index ba492cf717..a83ad40cab 100644 --- a/src/TinkerFile.cpp +++ b/src/TinkerFile.cpp @@ -31,6 +31,7 @@ static inline int SetNatomAndTitle(ArgList& lineIn, int& natom, std::string& tit } static inline bool IsAtomLine(ArgList& lineIn) { + if (lineIn.Nargs() < 1) return false; for (int i = 0; i < lineIn.Nargs(); i++) { std::string item = lineIn.GetStringNext(); if (i == 0 || i >= 5) { From 4e1b96564f8ac4749623f4ea8066e26540ddd543 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 21 Aug 2020 15:32:26 -0400 Subject: [PATCH 2/7] Update dependencies. --- src/cpptrajdepend | 10 ++++++---- src/cpptrajfiles | 2 ++ 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 3ce10297e1..14438aee66 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -73,7 +73,7 @@ Action_RunningAvg.o : Action_RunningAvg.cpp Action.h ActionState.h Action_Runnin 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 AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h 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 ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_Scale.o : Action_Scale.cpp Action.h ActionState.h Action_Scale.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_SetVelocity.o : Action_SetVelocity.cpp Action.h ActionState.h Action_SetVelocity.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h Constraints.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Random.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h -Action_Spam.o : Action_Spam.cpp Action.h ActionState.h Action_Spam.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mesh.h DataSet_double.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h HistBin.h ImagedAction.h KDE.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OnlineVarT.h PairList.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h +Action_Spam.o : Action_Spam.cpp Action.h ActionState.h Action_Spam.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataIO_Peaks.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mesh.h DataSet_Vector_Scalar.h DataSet_double.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h HistBin.h ImagedAction.h KDE.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OnlineVarT.h PairList.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_Strip.o : Action_Strip.cpp Action.h ActionState.h ActionTopWriter.h Action_Strip.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_Surf.o : Action_Surf.cpp Action.h ActionState.h Action_Surf.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_SymmetricRmsd.o : Action_SymmetricRmsd.cpp Action.h ActionState.h Action_SymmetricRmsd.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.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 Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h @@ -171,7 +171,7 @@ CpptrajFile.o : CpptrajFile.cpp CpptrajFile.h CpptrajStdio.h FileIO.h FileIO_Bzi CpptrajState.o : CpptrajState.cpp Action.h ActionList.h ActionState.h Action_CreateCrd.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_Topology.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleNavigator.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 Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Vec3.h CpptrajStdio.o : CpptrajStdio.cpp Parallel.h CurveFit.o : CurveFit.cpp CurveFit.h -DataFile.o : DataFile.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h ClusterDist.h ClusterSieve.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataFile.h DataIO.h DataIO_CCP4.h DataIO_CharmmFastRep.h DataIO_CharmmOutput.h DataIO_CharmmRepLog.h DataIO_CharmmRtfPrm.h DataIO_Cmatrix.h DataIO_Cpout.h DataIO_Evecs.h DataIO_Gnuplot.h DataIO_Grace.h DataIO_Mdout.h DataIO_NC_Cmatrix.h DataIO_OpenDx.h DataIO_RemLog.h DataIO_Std.h DataIO_VecTraj.h DataIO_XVG.h DataIO_Xplor.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Cmatrix.h DataSet_Cmatrix_MEM.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_RemLog.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 NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajectoryFile.h TypeNameHolder.h Vec3.h +DataFile.o : DataFile.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h ClusterDist.h ClusterSieve.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataFile.h DataIO.h DataIO_CCP4.h DataIO_CharmmFastRep.h DataIO_CharmmOutput.h DataIO_CharmmRepLog.h DataIO_CharmmRtfPrm.h DataIO_Cmatrix.h DataIO_Cpout.h DataIO_Evecs.h DataIO_Gnuplot.h DataIO_Grace.h DataIO_Mdout.h DataIO_NC_Cmatrix.h DataIO_OpenDx.h DataIO_Peaks.h DataIO_RemLog.h DataIO_Std.h DataIO_VecTraj.h DataIO_XVG.h DataIO_Xplor.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Cmatrix.h DataSet_Cmatrix_MEM.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_RemLog.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 NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajectoryFile.h TypeNameHolder.h Vec3.h DataFileList.o : DataFileList.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h PDBfile.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h DataIO.o : DataIO.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataSet.h DataSetList.h DataSet_2D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_MatrixDbl.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h DataIO_CCP4.o : DataIO_CCP4.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h ByteRoutines.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_CCP4.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h @@ -186,14 +186,15 @@ DataIO_Gnuplot.o : DataIO_Gnuplot.cpp ArgList.h Array1D.h AssociatedData.h Atom. DataIO_Grace.o : DataIO_Grace.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Grace.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h DataSet_string.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 ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h DataIO_Mdout.o : DataIO_Mdout.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Mdout.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.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 ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h DataIO_NC_Cmatrix.o : DataIO_NC_Cmatrix.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h ClusterDist.h ClusterSieve.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_NC_Cmatrix.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Cmatrix.h DataSet_Cmatrix_MEM.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h FileIO.h FileName.h Frame.h Hungarian.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h -DataIO_OpenDx.o : DataIO_OpenDx.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_OpenDx.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h +DataIO_OpenDx.o : DataIO_OpenDx.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_OpenDx.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridDbl.h DataSet_GridFlt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h +DataIO_Peaks.o : DataIO_Peaks.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Peaks.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Vector_Scalar.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 ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h DataIO_RemLog.o : DataIO_RemLog.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_RemLog.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_RemLog.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 ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h DataIO_Std.o : DataIO_Std.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h ClusterDist.h ClusterSieve.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Std.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Cmatrix_MEM.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mat3x3.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_string.h Dimension.h FileIO.h FileName.h Frame.h GridBin.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 ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h DataIO_VecTraj.o : DataIO_VecTraj.cpp ActionFrameCounter.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_VecTraj.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Vector.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 SymbolExporting.h TextFormat.h Timer.h Topology.h TrajectoryFile.h Trajout_Single.h TypeNameHolder.h Vec3.h DataIO_XVG.o : DataIO_XVG.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_XVG.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.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 ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h DataIO_Xplor.o : DataIO_Xplor.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Xplor.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h DataSet.o : DataSet.cpp AssociatedData.h CpptrajFile.h CpptrajStdio.h DataSet.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h -DataSetList.o : DataSetList.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h Box.h ClusterDist.h ClusterSieve.h ComplexArray.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Cmatrix_DISK.h DataSet_Cmatrix_MEM.h DataSet_Cmatrix_NOMEM.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_GridDbl.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_PHREMD.h DataSet_PHREMD_Explicit.h DataSet_PHREMD_Implicit.h DataSet_Parameters.h DataSet_RemLog.h DataSet_StringVar.h DataSet_Tensor.h DataSet_Topology.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_disk.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h Hungarian.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h SymmetricTensor.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TypeNameHolder.h Vec3.h +DataSetList.o : DataSetList.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h Box.h ClusterDist.h ClusterSieve.h ComplexArray.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Cmatrix_DISK.h DataSet_Cmatrix_MEM.h DataSet_Cmatrix_NOMEM.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_GridDbl.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_PHREMD.h DataSet_PHREMD_Explicit.h DataSet_PHREMD_Implicit.h DataSet_Parameters.h DataSet_RemLog.h DataSet_StringVar.h DataSet_Tensor.h DataSet_Topology.h DataSet_Vector.h DataSet_Vector_Scalar.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_disk.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h Hungarian.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Spline.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h SymmetricTensor.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TypeNameHolder.h Vec3.h DataSet_1D.o : DataSet_1D.cpp ArrayIterator.h AssociatedData.h ComplexArray.h Constants.h Corr.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h PubFFT.h Range.h TextFormat.h DataSet_3D.o : DataSet_3D.cpp AssociatedData.h Box.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_3D.h Dimension.h FileIO.h FileName.h GridBin.h Matrix_3x3.h MetaData.h Parallel.h Range.h TextFormat.h Vec3.h DataSet_Cmatrix.o : DataSet_Cmatrix.cpp ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AtomType.h Box.h ClusterDist.h ClusterSieve.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h DataSet_Cmatrix.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 SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Topology.h TypeNameHolder.h Vec3.h @@ -219,6 +220,7 @@ DataSet_StringVar.o : DataSet_StringVar.cpp AssociatedData.h CpptrajFile.h Cpptr DataSet_Tensor.o : DataSet_Tensor.cpp AssociatedData.h CpptrajFile.h DataSet.h DataSet_Tensor.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h SymmetricTensor.h TextFormat.h DataSet_Topology.o : DataSet_Topology.cpp AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Topology.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ParmFile.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Topology.h TypeNameHolder.h Vec3.h DataSet_Vector.o : DataSet_Vector.cpp ArrayIterator.h AssociatedData.h ComplexArray.h Constants.h Corr.h CpptrajFile.h DataSet.h DataSet_1D.h DataSet_Vector.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h PubFFT.h Range.h SymbolExporting.h TextFormat.h Vec3.h +DataSet_Vector_Scalar.o : DataSet_Vector_Scalar.cpp AssociatedData.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Vector_Scalar.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h Vec3.h DataSet_double.o : DataSet_double.cpp AssociatedData.h CpptrajFile.h DataSet.h DataSet_1D.h DataSet_double.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h DataSet_float.o : DataSet_float.cpp AssociatedData.h CpptrajFile.h DataSet.h DataSet_1D.h DataSet_float.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h DataSet_integer_disk.o : DataSet_integer_disk.cpp AssociatedData.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h DataSet_integer.h DataSet_integer_disk.h Dimension.h FileIO.h FileName.h File_TempName.h MetaData.h NC_Routines.h Parallel.h Range.h TextFormat.h diff --git a/src/cpptrajfiles b/src/cpptrajfiles index 717427b8f4..9128409fde 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -186,6 +186,7 @@ COMMON_SOURCES= \ DataIO_Gnuplot.cpp \ DataIO_Grace.cpp \ DataIO_OpenDx.cpp \ + DataIO_Peaks.cpp \ DataIO_Mdout.cpp \ DataIO_NC_Cmatrix.cpp \ DataIO_RemLog.cpp \ @@ -219,6 +220,7 @@ COMMON_SOURCES= \ DataSet_RemLog.cpp \ DataSet_Topology.cpp \ DataSet_Vector.cpp \ + DataSet_Vector_Scalar.cpp \ DataSet_double.cpp \ DataSet_float.cpp \ DataSet_integer_disk.cpp \ From 5bf3abf4de8648dffc1e988ff85146b817d7d004 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 21 Aug 2020 15:38:59 -0400 Subject: [PATCH 3/7] Changes from fixspam. Peaks are stored in a data set now. Allow volmap to be compiled double precision if desired. Fix some uninitialized vars. --- src/Action_Volmap.cpp | 218 +++++++++++++++++++++++++++++------------- src/Action_Volmap.h | 25 +++-- src/cpptrajdepend | 2 +- 3 files changed, 170 insertions(+), 75 deletions(-) diff --git a/src/Action_Volmap.cpp b/src/Action_Volmap.cpp index 1aaff46c4b..76d5e1ece1 100644 --- a/src/Action_Volmap.cpp +++ b/src/Action_Volmap.cpp @@ -3,11 +3,21 @@ #include "Action_Volmap.h" #include "Constants.h" // PI #include "CpptrajStdio.h" +#ifdef VOLMAP_DOUBLE +# include "DataSet_GridDbl.h" +# define VOLMAP_DS_T DataSet_GridDbl +# define VOLMAP_T double +#else +# include "DataSet_GridFlt.h" +# define VOLMAP_DS_T DataSet_GridFlt +# define VOLMAP_T float +#endif #ifdef _OPENMP # include #endif -const double Action_Volmap::sqrt_8_pi_cubed = sqrt(8.0*Constants::PI*Constants::PI*Constants::PI); +const double Action_Volmap::sqrt_8_pi_cubed_ = sqrt(8.0*Constants::PI*Constants::PI*Constants::PI); + // CONSTRUCTOR Action_Volmap::Action_Volmap() : radiiType_(UNSPECIFIED), @@ -17,11 +27,15 @@ Action_Volmap::Action_Volmap() : xmin_(0.0), ymin_(0.0), zmin_(0.0), + debug_(0), Nframes_(0), setupGridOnMask_(false), spheremode_(false), grid_(0), + total_volume_(0), + calcpeaks_(false), peakfile_(0), + peakdata_(0), peakcut_(0.05), buffer_(3.0), radscale_(1.0), @@ -40,25 +54,27 @@ void Action_Volmap::Help() const { } void Action_Volmap::RawHelp() const { - mprintf("\t[out ] [radscale ] [stepfac ] [sphere]\n" + mprintf("\t[out ] [radscale ] [stepfac ]\n" + "\t[sphere] [radii {vdw | element}]\n" + "\t[calcpeaks] [peakcut ] [peakfile ]\n" "\t{ data |\n" "\t name [ ]\n" "\t { size [center ] |\n" "\t centermask [buffer ] |\n" "\t boxref }\n" "\t}\n" - "\t[radii {vdw | element}] [peakcut ] [peakfile ]\n"); + ); } // Action_Volmap::Init() Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int debugIn) { + debug_ = debugIn; # ifdef MPI trajComm_ = init.TrajComm(); # endif // Get specific keywords peakcut_ = actionArgs.getKeyDouble("peakcut", 0.05); - peakfile_ = init.DFL().AddCpptrajFile(actionArgs.GetStringKey("peakfile"), "Volmap Peaks"); spheremode_ = actionArgs.hasKey("sphere"); radscale_ = 1.0; stepfac_ = 4.1; @@ -101,6 +117,12 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d return Action::ERR; } +# ifdef VOLMAP_DOUBLE + static const DataSet::DataType GridDataType = DataSet::GRID_DBL; +# else + static const DataSet::DataType GridDataType = DataSet::GRID_FLT; +# endif + if (mode != DATASET) { Vec3 Sizes(0.0); Vec3 Center(0.0); @@ -154,7 +176,7 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d return Action::ERR; } // Allocate grid dataset - grid_ = (DataSet_GridFlt*)init.DSL().AddSet(DataSet::GRID_FLT, setname, "VOLMAP"); + grid_ = (VOLMAP_DS_T*)init.DSL().AddSet(GridDataType, setname, "VOLMAP"); if (grid_ == 0) { mprinterr("Error: Could not create grid dataset '%s'\n", setname.c_str()); return Action::ERR; @@ -174,7 +196,7 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d } } else { // Get existing grid dataset - grid_ = (DataSet_GridFlt*)init.DSL().FindSetOfType( setup_arg, DataSet::GRID_FLT ); + grid_ = (VOLMAP_DS_T*)init.DSL().FindSetOfType( setup_arg, GridDataType ); if (grid_ == 0) { mprinterr("Error: Could not find grid data set with name '%s'\n", setup_arg.c_str()); @@ -201,12 +223,36 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d return Action::ERR; } if (densitymask_.SetMaskString(reqmask)) return Action::ERR; - // Get output filename + // See if peaks are being determined + calcpeaks_ = actionArgs.hasKey("calcpeaks"); + peakfile_ = 0; + std::string pfilename = actionArgs.GetStringKey("peakfile"); + if (!pfilename.empty()) { + peakfile_ = init.DFL().AddDataFile(pfilename, actionArgs); + if (peakfile_ == 0) { + mprinterr("Error: Unable to allocate peak file.\n"); + return Action::ERR; + } + calcpeaks_ = true; + } + if (calcpeaks_) { + peakdata_ = init.DSL().AddSet(DataSet::VECTOR_SCALAR, MetaData(grid_->Meta().Name(), "peaks")); + if (peakdata_ == 0) { + mprinterr("Error: Unable to allocate peak data set.\n"); + return Action::ERR; + } +# ifdef MPI + peakdata_->SetNeedsSync( false ); +# endif + if (peakfile_ != 0) peakfile_->AddDataSet( peakdata_ ); + } + // Get output filename (newer syntax) std::string outfilename = actionArgs.GetStringKey("out"); if (outfilename.empty()) outfilename = actionArgs.GetStringNext(); // Backwards compat. DataFile* outfile = init.DFL().AddDataFile( outfilename, actionArgs ); if (outfile != 0) outfile->AddDataSet( grid_ ); + // Create total volume set total_volume_ = init.DSL().AddSet(DataSet::DOUBLE, MetaData(grid_->Meta().Name(), "totalvol")); if (total_volume_ == 0) return Action::ERR; @@ -228,6 +274,11 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d // Info mprintf(" VOLMAP: Grid spacing will be %.2fx%.2fx%.2f Angstroms\n", dx_, dy_, dz_); +# ifdef VOLMAP_DOUBLE + mprintf("\tUsing double precision grid.\n"); +# else + mprintf("\tUsing single precision grid.\n"); +# endif if (setupGridOnMask_) mprintf("\tGrid centered around mask %s with %.2f Ang. clearance\n", centermask_.MaskExpression().c_str(), buffer_); @@ -244,13 +295,17 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d mprintf("\tWhen smearing Gaussian, voxels farther than radii/2 will be skipped.\n"); mprintf("\tDividing radii by %f\n", 1.0/radscale_); mprintf("\tFactor for determining number of bins to smear Gaussian is %f\n", stepfac_); + if (outfile != 0) mprintf("\tDensity will wrtten to '%s'\n", outfile->DataFilename().full()); mprintf("\tGrid dataset name is '%s'\n", grid_->legend()); mprintf("\tTotal grid volume dataset name is '%s'\n", total_volume_->legend()); - if (peakfile_ != 0) - mprintf("\tDensity peaks above %.3f will be printed to %s in XYZ-format\n", - peakcut_, peakfile_->Filename().full()); + if (calcpeaks_) { + mprintf("\tDensity peaks above %.3f will be saved to %s\n", peakcut_, peakdata_->legend()); + if (peakfile_ != 0) + mprintf("\tDensity peaks will be printed to %s in XYZ format.\n", + peakfile_->DataFilename().full()); + } # ifdef _OPENMP if (GRID_THREAD_.size() > 1) mprintf("\tParallelizing calculation with %zu threads.\n", GRID_THREAD_.size()); @@ -299,6 +354,7 @@ Action::RetType Action_Volmap::Setup(ActionSetup& setup) { else radiiToUse = ELEMENT; } + double maxRad = 0; for (AtomMask::const_iterator atom = densitymask_.begin(); atom != densitymask_.end(); ++atom) { double rad = 0.0; @@ -307,10 +363,27 @@ Action::RetType Action_Volmap::Setup(ActionSetup& setup) { else if (radiiToUse == ELEMENT) rad = setup.Top()[*atom].ElementRadius(); if (rad > 0.0) { - halfradii_.push_back( (float)(rad * radscale_ / 2) ); + halfradii_.push_back( (double)(rad * radscale_ / 2) ); Atoms_.push_back( *atom ); + // For determining function range. + if (halfradii_.back() > maxRad) maxRad = halfradii_.back(); } } + // Try to determine the max value we may encounter for the exponential. + //mprintf("\tMax observed radius: %g Ang\n", maxRad); + int nxstep = (int) ceil(stepfac_ * maxRad / dx_)+1; + int nystep = (int) ceil(stepfac_ * maxRad / dy_)+1; + int nzstep = (int) ceil(stepfac_ * maxRad / dz_)+1; + double maxx = (double)nxstep * dx_; + double maxy = (double)nystep * dy_; + double maxz = (double)nzstep * dz_; + double maxDist = maxx*maxx + maxy*maxy + maxz*maxz; + //mprintf("DEBUG: nx= %i ny= %i nz= %i\n", nxstep, nystep, nzstep); + //mprintf("DEBUG: %g %g %g %g\n", maxx, maxy, maxz, maxDist); + maxDist *= (-1.0 / (2.0 * maxRad * maxRad)); + if (debug_ > 1) + mprintf("DEBUG: maxDist= %g\n", maxDist); + if ((int)Atoms_.size() < densitymask_.Nselected()) mprintf("Warning: %i atoms have 0.0 radii and will be skipped.\n", densitymask_.Nselected() - (int)Atoms_.size()); @@ -400,48 +473,53 @@ Action::RetType Action_Volmap::DoAction(int frameNum, ActionFrame& frm) { # endif for (midx = 0; midx < maxidx; midx++) { double rhalf = (double)halfradii_[midx]; - double rcut2; - if (spheremode_) - rcut2 = rhalf*rhalf; - else - rcut2 = 99999999.0; - atom = Atoms_[midx]; - Vec3 pt = Vec3(frm.Frm().XYZ(atom)); - int ix = (int) ( floor( (pt[0]-xmin_) / dx_ + 0.5 ) ); + double rcut2; + if (spheremode_) + rcut2 = rhalf*rhalf; + else + rcut2 = 99999999.0; + atom = Atoms_[midx]; + Vec3 pt = Vec3(frm.Frm().XYZ(atom)); + int ix = (int) ( floor( (pt[0]-xmin_) / dx_ + 0.5 ) ); + /* See how many steps in X dimension we smear out our Gaussian. This + * formula is taken to be consistent with VMD's volmap tool. + */ + int nxstep = (int) ceil(stepfac_ * rhalf / dx_); + if (ix >= -nxstep && ix <= nX + nxstep) { int iy = (int) ( floor( (pt[1]-ymin_) / dy_ + 0.5 ) ); - int iz = (int) ( floor( (pt[2]-zmin_) / dz_ + 0.5 ) ); - /* See how many steps in each dimension we smear out our Gaussian. This - * formula is taken to be consistent with VMD's volmap tool - */ - int nxstep = (int) ceil(stepfac_ * rhalf / dx_); + // See how many steps in Y dim we smear out Gaussian. int nystep = (int) ceil(stepfac_ * rhalf / dy_); - int nzstep = (int) ceil(stepfac_ * rhalf / dz_); - if (ix < -nxstep || ix > nX + nxstep || - iy < -nystep || iy > nY + nystep || - iz < -nzstep || iz > nZ + nzstep) - continue; - // Calculate the gaussian normalization factor (in 3 dimensions with the - // given half-radius) - double norm = 1 / (sqrt_8_pi_cubed * rhalf*rhalf*rhalf); - double exfac = -1.0 / (2.0 * rhalf * rhalf); - //mprintf("DBG: Atom %i norm %g exfac %g\n", atom+1, norm, exfac); + if (iy >= -nystep && iy <= nY + nystep) { + int iz = (int) ( floor( (pt[2]-zmin_) / dz_ + 0.5 ) ); + // See how many steps in Z dim we smear out Gaussian. + int nzstep = (int) ceil(stepfac_ * rhalf / dz_); + if (iz >= -nzstep && iz <= nZ + nzstep) { + // Calculate the gaussian normalization factor (in 3 dimensions with the + // given half-radius) + double norm = 1 / (sqrt_8_pi_cubed_ * rhalf*rhalf*rhalf); + double exfac = -1.0 / (2.0 * rhalf * rhalf); + //mprintf("DBG: Atom %i norm %g exfac %g\n", atom+1, norm, exfac); - int xend = std::min(ix+nxstep, nX); - int yend = std::min(iy+nystep, nY); - int zend = std::min(iz+nzstep, nZ); - for (int xval = std::max(ix-nxstep, 0); xval < xend; xval++) - for (int yval = std::max(iy-nystep, 0); yval < yend; yval++) - for (int zval = std::max(iz-nzstep, 0); zval < zend; zval++) { - Vec3 gridpt = Vec3(xmin_+xval*dx_, ymin_+yval*dy_, zmin_+zval*dz_) - pt; - double dist2 = gridpt.Magnitude2(); - if (dist2 < rcut2) { -# ifdef _OPENMP - GRID_THREAD_[mythread].incrementBy(xval, yval, zval, norm * exp(exfac * dist2)); -# else - grid_->Increment(xval, yval, zval, norm * exp(exfac * dist2)); -# endif - } - } + int xend = std::min(ix+nxstep, nX); + int yend = std::min(iy+nystep, nY); + int zend = std::min(iz+nzstep, nZ); + for (int xval = std::max(ix-nxstep, 0); xval < xend; xval++) + for (int yval = std::max(iy-nystep, 0); yval < yend; yval++) + for (int zval = std::max(iz-nzstep, 0); zval < zend; zval++) { + Vec3 gridpt = Vec3(xmin_+xval*dx_, ymin_+yval*dy_, zmin_+zval*dz_) - pt; + double dist2 = gridpt.Magnitude2(); + if (dist2 < rcut2) { + //mprintf("DEBUG: rhalf= %g dist2= %g exfac= %g exp= %g\n", rhalf, dist2, exfac, exfac*dist2); +# ifdef _OPENMP + GRID_THREAD_[mythread].incrementBy(xval, yval, zval, norm * exp(exfac * dist2)); +# else /* OPENMP */ + grid_->Increment(xval, yval, zval, norm * exp(exfac * dist2)); +# endif /* OPENMP */ + } + } // END loop over zval + } // END z dim is valid + } // END y dim is valid + } // END x dim is valid } // END loop over atoms in densitymask_ # ifdef _OPENMP } // END pragma omp parallel @@ -492,12 +570,12 @@ void Action_Volmap::Print() { CombineGridThreads(); # endif // Divide our grid by the number of frames - float nf = (float)Nframes_; - for (DataSet_GridFlt::iterator gval = grid_->begin(); gval != grid_->end(); ++gval) + VOLMAP_T nf = (VOLMAP_T)Nframes_; + for (VOLMAP_DS_T::iterator gval = grid_->begin(); gval != grid_->end(); ++gval) *gval /= nf; // Print volume estimate unsigned int nOccupiedVoxels = 0; - for (DataSet_GridFlt::iterator gval = grid_->begin(); gval != grid_->end(); ++gval) { + for (VOLMAP_DS_T::iterator gval = grid_->begin(); gval != grid_->end(); ++gval) { if (*gval > 0.0) { ++nOccupiedVoxels; //mprintf("DBG: %16.8e\n", *gval); @@ -512,7 +590,7 @@ void Action_Volmap::Print() { // "rdparm generated grid density" ); // See if we need to write the peaks out somewhere - if (peakfile_ != 0) { + if (calcpeaks_) { // Extract peaks from the current grid, setup another Grid instance. This // works by taking every grid point and analyzing all grid points adjacent // to it (including diagonals). If any of those grid points have a higher @@ -520,11 +598,11 @@ void Action_Volmap::Print() { // that value is _not_ a maximum. Any density peaks less than the minimum // filter are discarded. The result is a Grid instance with all non-peak // grid points zeroed-out. - Grid peakgrid = grid_->InternalGrid(); + Grid peakgrid = grid_->InternalGrid(); for (size_t i = 0; i < grid_->NX(); i++) for (size_t j = 0; j < grid_->NY(); j++) for (size_t k = 0; k < grid_->NZ(); k++) { - float val = grid_->GridVal(i, j, k); + VOLMAP_T val = grid_->GridVal(i, j, k); if (val < peakcut_) { peakgrid.setGrid(i, j, k, 0.0f); continue; @@ -541,29 +619,33 @@ void Action_Volmap::Print() { } } int npeaks = 0; - std::vector peakdata; + //std::vector peakdata; + double pbuf[4]; for (size_t i = 0; i < peakgrid.NX(); i++) for (size_t j = 0; j < peakgrid.NY(); j++) for (size_t k = 0; k < peakgrid.NZ(); k++) { double gval = peakgrid.element(i, j, k); if (gval > 0) { - npeaks++; - peakdata.push_back(xmin_+dx_*i); - peakdata.push_back(ymin_+dy_*j); - peakdata.push_back(zmin_+dz_*k); - peakdata.push_back(gval); + //npeaks++; + //peakdata.push_back(xmin_+dx_*i); + //peakdata.push_back(ymin_+dy_*j); + //peakdata.push_back(zmin_+dz_*k); + //peakdata.push_back(gval); + pbuf[0] = xmin_+dx_*(double)i; + pbuf[1] = ymin_+dy_*(double)j; + pbuf[2] = zmin_+dz_*(double)k; + pbuf[3] = gval; + peakdata_->Add(npeaks++, pbuf); } } - // If we have peaks, open up our peak data and print it + // Report on how many peaks found. if (npeaks > 0) { - peakfile_->Printf("%d\n\n", npeaks); - for (int i = 0; i < npeaks; i++) - peakfile_->Printf("C %16.8f %16.8f %16.8f %16.8f\n", peakdata[4*i], - peakdata[4*i+1], peakdata[4*i+2], peakdata[4*i+3]); mprintf("Volmap: %d density peaks found with higher density than %.4lf\n", npeaks, peakcut_); - }else{ + } else { mprintf("No peaks found with a density greater than %.3lf\n", peakcut_); } } } +#undef VOLMAP_DS_T +#undef VOLMAP_T diff --git a/src/Action_Volmap.h b/src/Action_Volmap.h index 5103c3c076..77041a0e26 100644 --- a/src/Action_Volmap.h +++ b/src/Action_Volmap.h @@ -1,7 +1,15 @@ #ifndef INC_ACTION_VOLMAP_H #define INC_ACTION_VOLMAP_H #include "Action.h" -#include "DataSet_GridFlt.h" +#include "Grid.h" +#ifdef VOLMAP_DOUBLE +# define VOLMAP_DS_T DataSet_GridDbl +# define VOLMAP_T double +#else +# define VOLMAP_DS_T DataSet_GridFlt +# define VOLMAP_T float +#endif +class VOLMAP_DS_T; class Action_Volmap : public Action { public: Action_Volmap(); @@ -25,25 +33,30 @@ class Action_Volmap : public Action { double dx_, dy_, dz_; /// minimum values in the x-, y-, and z-dimensions double xmin_, ymin_, zmin_; + int debug_; int Nframes_; ///< Number of frames we analyzed so we can average at the end bool setupGridOnMask_; ///< If true, set up the grid on first frame based on centermask. bool spheremode_; ///< If true, grid points farther than rhalf^2 will be skipped AtomMask centermask_; ///< Mask to center the grid on AtomMask densitymask_; ///< Max of atoms to grid. - DataSet_GridFlt* grid_; ///< Hold the grid. + VOLMAP_DS_T* grid_; ///< Hold the grid. DataSet* total_volume_; ///< Hold total grid volume. - CpptrajFile* peakfile_; ///< file name with the peak locations as Carbons in XYZ file format + bool calcpeaks_; ///< If true, calculate peaks + DataFile* peakfile_; ///< Optional file to write peak data to as Carbons in XYZ file format + DataSet* peakdata_; ///< Data set holding peak locations along with densities double peakcut_; ///< The value below which to ignore all peaks std::vector Atoms_; ///< Atoms with radii > 0.0 - std::vector halfradii_; ///< 1/2 the atomic radii of each atom in the gridded selection + std::vector halfradii_; ///< 1/2 the atomic radii of each atom in the gridded selection double buffer_; ///< Clearance between the edges of our grid and centermask_ double radscale_; ///< The scaling factor to divide all radii by double stepfac_; ///< Factor for determining how many steps to smear Gaussian - static const double sqrt_8_pi_cubed; + static const double sqrt_8_pi_cubed_; # ifdef _OPENMP - typedef std::vector< Grid > Garray; + typedef std::vector< Grid > Garray; Garray GRID_THREAD_; void CombineGridThreads(); # endif }; +#undef VOLMAP_DS_T +#undef VOLMAP_T #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 14438aee66..e661447540 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -84,7 +84,7 @@ Action_Unstrip.o : Action_Unstrip.cpp Action.h ActionState.h Action_Unstrip.h Ar Action_Unwrap.o : Action_Unwrap.cpp Action.h ActionState.h Action_Unwrap.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h ImageRoutines.h ImageTypes.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_Vector.o : Action_Vector.cpp Action.h ActionState.h Action_Vector.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h ComplexArray.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_Vector.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_VelocityAutoCorr.o : Action_VelocityAutoCorr.cpp Action.h ActionState.h Action_VelocityAutoCorr.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h ComplexArray.h Constants.h CoordinateInfo.h Corr.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Vector.h DataSet_double.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h PubFFT.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h -Action_Volmap.o : Action_Volmap.cpp Action.h ActionState.h Action_Volmap.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h +Action_Volmap.o : Action_Volmap.cpp Action.h ActionState.h Action_Volmap.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridDbl.h DataSet_GridFlt.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_Volume.o : Action_Volume.cpp Action.h ActionState.h Action_Volume.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h 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 ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_Watershell.o : Action_Watershell.cpp Action.h ActionState.h Action_Watershell.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageRoutines.h ImageTypes.h ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h Action_XtalSymm.o : Action_XtalSymm.cpp Action.h ActionState.h Action_XtalSymm.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SpaceGroup.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h From b89c3bd29b7cefc7a9f172284c29240cfdc33f7d Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 21 Aug 2020 16:01:24 -0400 Subject: [PATCH 4/7] Clean up the volmap and spam manual entries. --- doc/cpptraj.lyx | 448 +++++++++++++++++++++++++++++++----------------- 1 file changed, 290 insertions(+), 158 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index b719c6d5f2..762c92d96c 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -33336,129 +33336,6 @@ secstruct :1-22 out dssp.gnu sumout dssp.agr xmgrace dssp.agr \end_layout -\begin_layout Subsection -spam -\end_layout - -\begin_layout LyX-Code -spam [solv ] [reorder] [name ] -\end_layout - -\begin_layout LyX-Code - [purewater] [cut ] [info ] [summary ] -\end_layout - -\begin_layout LyX-Code - [site_size ] [sphere] [out ] -\end_layout - -\begin_layout LyX-Code - [dgbulk ] [dhbulk ] [temperature ] -\end_layout - -\begin_deeper -\begin_layout Description - File with the peak locations present (XYZ- format) -\end_layout - -\begin_layout Description - Name of the solvent residues -\end_layout - -\begin_layout Description - Non-bonded cutoff for energy evaluation -\end_layout - -\begin_layout Description - SPAM free energy of the bulk solvent in kcal/mol; default is -30.3 - kcal/mol (SPC/E water). -\end_layout - -\begin_layout Description - SPAM enthalpy of the bulk solvent in kcal/mol; default is -22.2 - kcal/mol (SPC/E water). -\end_layout - -\begin_layout Description - Temperature at which SPAM calculation was run. -\end_layout - -\begin_layout Description - File with stats about which sites are occupied when. -\end_layout - -\begin_layout Description - Size of the water site around each density peak. -\end_layout - -\begin_layout Description -[sphere] Treat each site like a sphere. -\end_layout - -\begin_layout Description -[purewater] The system is pure water---used to parametrize the bulk values. -\end_layout - -\begin_layout Description -[reorder] The solvent should be re-ordered so the same solvent molecule - is always in the same site. -\end_layout - -\begin_layout Description - File with the summary of all SPAM results. - If not specified, no SPAM energies will be calculated. -\end_layout - -\begin_layout Description - Data file with all SPAM energies for each snapshot. -\end_layout - -\end_deeper -\begin_layout Standard -Perform profiling of bound water molecules via SPAM analysis -\begin_inset CommandInset citation -LatexCommand citep -key "Cui2013" -literal "true" - -\end_inset - -. - Briefly, this method identifies and estimates the free energy profiles - of bound waters via calculation of the distribution of interaction energies - between the water and it's environment from explicit solvent MD trajectories. - The interaction energies are calculated using a force- and energy-shifted - electrostatic term with a hard cutoff. -\end_layout - -\begin_layout Standard -Prior to this command, the -\series bold -\emph on -volmap -\series default -\emph default - command should be run with the -\series bold -peakfile -\series default - keyword (see -\begin_inset CommandInset ref -LatexCommand vref -reference "subsec:cpptraj_volmap" - -\end_inset - -) to generate the peaks file. - If not using peaks from the -\series bold -\emph on -volmap -\series default -\emph default - command, the peaks file should have one line per peak with format: -\end_layout - \begin_layout LyX-Code C \end_layout @@ -33622,6 +33499,178 @@ Set velocities in frame for atoms in using Maxwellian distribution based on given temperature. \end_layout +\begin_layout Subsection +spam +\end_layout + +\begin_layout LyX-Code +spam [name ] [out ] [cut ] [solv ] +\end_layout + +\begin_layout LyX-Code + { purewater | +\end_layout + +\begin_layout LyX-Code + [reorder] [info ] [summary ] +\end_layout + +\begin_layout LyX-Code + [site_size ] [sphere] [temperature ] +\end_layout + +\begin_layout LyX-Code + [dgbulk ] [dhbulk ] } +\end_layout + +\begin_deeper +\begin_layout Description +name +\begin_inset space ~ +\end_inset + + Output data sets name. +\end_layout + +\begin_layout Description +out +\begin_inset space ~ +\end_inset + + Data file with all SPAM energies for each snapshot. +\end_layout + +\begin_layout Description +cut +\begin_inset space ~ +\end_inset + + Non-bonded cutoff for energy evaluation +\end_layout + +\begin_layout Description +solv +\begin_inset space ~ +\end_inset + + Name of the solvent residues. +\end_layout + +\begin_layout Description +[purewater] The system is pure water. + Used to parametrize the bulk values. + If this is specified, none of the below options are relevant. +\end_layout + +\begin_layout Description + Data set or file (XYZ- format: see below) with the peak locations + present . + +\end_layout + +\begin_layout Description +[reorder] The solvent should be re-ordered so the same solvent molecule + is always in the same site. +\end_layout + +\begin_layout Description +info +\begin_inset space ~ +\end_inset + + File with stats about which sites are occupied when. +\end_layout + +\begin_layout Description +summary +\begin_inset space ~ +\end_inset + + File with the summary of all SPAM results. + If not specified, no SPAM energies will be calculated. +\end_layout + +\begin_layout Description +site_size Size of the water site around each density peak (sphere + diameter/box edge length) in Ang. +\end_layout + +\begin_layout Description +[sphere] Treat each site like a sphere. +\end_layout + +\begin_layout Description +temperature +\begin_inset space ~ +\end_inset + + Temperature at which SPAM calculation was run. +\end_layout + +\begin_layout Description +dgbulk +\begin_inset space ~ +\end_inset + + SPAM free energy of the bulk solvent in kcal/mol; default is -30.3 + kcal/mol (SPC/E water). +\end_layout + +\begin_layout Description +dhbulk +\begin_inset space ~ +\end_inset + + SPAM enthalpy of the bulk solvent in kcal/mol; default is -22.2 + kcal/mol (SPC/E water). +\end_layout + +\end_deeper +\begin_layout Standard +Perform profiling of bound water molecules via SPAM analysis +\begin_inset CommandInset citation +LatexCommand citep +key "Cui2013" +literal "true" + +\end_inset + +. + Briefly, this method identifies and estimates the free energy profiles + of bound waters via calculation of the distribution of interaction energies + between the water and it's environment from explicit solvent MD trajectories. + The interaction energies are calculated using a force- and energy-shifted + electrostatic term with a hard cutoff. +\end_layout + +\begin_layout Standard +Prior to this command, the +\series bold +\emph on +volmap +\series default +\emph default + command should be run with the +\series bold +peakfile +\series default + keyword (see +\begin_inset CommandInset ref +LatexCommand vref +reference "subsec:cpptraj_volmap" + +\end_inset + +) to generate the peaks file. + If not using peaks from the +\series bold +\emph on +volmap +\series default +\emph default + command, the peaks file should have one line per peak with format: +\end_layout + \begin_layout Subsection stfcdiffusion \end_layout @@ -35068,7 +35117,15 @@ name "subsec:cpptraj_volmap" \end_layout \begin_layout LyX-Code -volmap filename dx dy dz [radscale ] +volmap [out ] [radscale ] [stepfac ] +\end_layout + +\begin_layout LyX-Code + [sphere] [radii {vdw | element}] +\end_layout + +\begin_layout LyX-Code + [calcpeaks] [peakcut ] [peakfile ] \end_layout \begin_layout LyX-Code @@ -35076,15 +35133,19 @@ volmap filename dx dy dz [radscale ] \end_layout \begin_layout LyX-Code - name { size [center ] | + name [ ] +\end_layout + +\begin_layout LyX-Code + { size [center ] | \end_layout \begin_layout LyX-Code - centermask [buffer ] } } + centermask [buffer ] | \end_layout \begin_layout LyX-Code - [peakcut ] [peakfile ] + boxref } } \begin_inset Separator latexpar \end_inset @@ -35095,7 +35156,11 @@ volmap filename dx dy dz [radscale ] \begin_layout Description \family roman -filename +out +\begin_inset space ~ +\end_inset + + \family default The name of the output file with the grid density. \end_layout @@ -35103,40 +35168,90 @@ The name of the output file with the grid density. \begin_layout Description \family roman -dx, + +\family default +The atom selection from which to calculate the number density. +\end_layout + +\begin_layout Description + +\family roman +radscale \begin_inset space ~ \end_inset -dy, + +\family default +Factor by which to scale radii (by division). + To match the atomic radius of Oxygen used by the VMD volmap tool, a scaling + factor of 1.36 should be used. + Default 1.0. +\end_layout + +\begin_layout Description +stepfac \begin_inset space ~ \end_inset -dz -\family default -The grid spacing (Angstroms) in the X-, Y-, and Z-dimensions, respectively + Factor for determining how many voxels to smear Gaussian (default + 4.1, 1.0 for +\series bold +sphere +\series default +). +\end_layout + +\begin_layout Description +sphere When smearing Gaussian, skip voxels farther than radii/2. +\end_layout + +\begin_layout Description +radii +\begin_inset space ~ +\end_inset + +{vdw|element} Specify either van der Waals radii (default) or elemental + radii. +\end_layout + +\begin_layout Description +calcpeaks If specified, peaks in the grid density will be calculated and + saved to set with aspect +\begin_inset Quotes eld +\end_inset + +peaks +\begin_inset Quotes erd +\end_inset + +. \end_layout \begin_layout Description \family roman - +peakcut +\begin_inset space ~ +\end_inset + + \family default -The atom selection from which to calculate the number density. +The minimum density required to consider a local maximum a 'density peak' + in the outputted peak file (default 0.05). \end_layout \begin_layout Description \family roman -radscale +peakfile \begin_inset space ~ \end_inset - + \family default -Factor by which to scale radii (by division). - To match the atomic radius of Oxygen used by the VMD volmap tool, a scaling - factor of 1.36 should be used. - Default 1.0. +A file in XYZ-format that contains a carbon atom centered at the grid point + of every local density maximum. + This file is necessary input to the spam action command. \end_layout \begin_layout Description @@ -35165,6 +35280,22 @@ centermask/buffer \begin_layout Description +\family roman +dx, +\begin_inset space ~ +\end_inset + +dy, +\begin_inset space ~ +\end_inset + +dz +\family default +The grid spacing (Angstroms) in the X-, Y-, and Z-dimensions, respectively. +\end_layout + +\begin_layout Description + \family roman size \begin_inset space ~ @@ -35176,6 +35307,7 @@ Specify the size of the grid in the X-, Y-, and Z-dimensions. Must be used alongside the center argument. \end_layout +\begin_deeper \begin_layout Description \family roman @@ -35190,6 +35322,7 @@ Specify the grid center explicitly. Default is the origin. \end_layout +\end_deeper \begin_layout Description \family roman @@ -35204,6 +35337,7 @@ The mask around which the grid should be centered (via geometric center). entered (see above) is used in its place. \end_layout +\begin_deeper \begin_layout Description \family roman @@ -35220,31 +35354,29 @@ A buffer distance, in Angstroms, by which the edges of the grid should clear The buffer is ignored if the center and size are specified (see below). \end_layout +\end_deeper \begin_layout Description - -\family roman -peakcut +boxref \begin_inset space ~ \end_inset - -\family default -The minimum density required to consider a local maximum a 'density peak' - in the outputted peak file (default 0.05). + Set up the grid using the unit cell info in the specified reference. \end_layout -\begin_layout Description +\begin_layout Standard +Data Sets Created: +\end_layout -\family roman -peakfile -\begin_inset space ~ -\end_inset +\begin_layout Description + The 3D grid. +\end_layout - -\family default -A file in XYZ-format that contains a carbon atom centered at the grid point - of every local density maximum. - This file is necessary input to the spam action command. +\begin_layout Description +[peaks] The density peaks if +\series bold +calcpeaks +\series default + specified. \end_layout \end_deeper From 7274169771b685965190e548b0695e8325f076b1 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 21 Aug 2020 16:11:43 -0400 Subject: [PATCH 5/7] Add type keyword to OpenDX read help and manual entry --- doc/cpptraj.lyx | 87 ++++++++++++++++++++++++++++++++++++++++++- src/DataFile.cpp | 2 +- src/DataIO_OpenDx.cpp | 5 +++ src/DataIO_OpenDx.h | 1 + 4 files changed, 93 insertions(+), 2 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 762c92d96c..f9bba6510e 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -2759,7 +2759,7 @@ status open \begin_layout Plain Layout \align center \begin_inset Tabular - + @@ -3250,6 +3250,53 @@ pH data only \begin_inset Text +\begin_layout Plain Layout +Density Peaks +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +.peaks +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +peaks +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +3D density peaks (spam/volmap) +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout + +\end_layout + +\end_inset + + + + +\begin_inset Text + \begin_layout Plain Layout Vector pseudo-traj \end_layout @@ -5409,6 +5456,24 @@ parmout OpenDX file options \end_layout +\begin_layout Subsubsection* +Read +\end_layout + +\begin_layout LyX-Code +[type {float|double}] +\end_layout + +\begin_deeper +\begin_layout Description +type +\begin_inset space ~ +\end_inset + +{float|double} Precision to read in 3D grid (default float). +\end_layout + +\end_deeper \begin_layout Subsubsection* Write \end_layout @@ -33671,6 +33736,26 @@ volmap command, the peaks file should have one line per peak with format: \end_layout +\begin_layout LyX-Code +<# of peaks> +\end_layout + +\begin_layout LyX-Code + +\end_layout + +\begin_layout LyX-Code +C +\end_layout + +\begin_layout LyX-Code +... +\end_layout + +\begin_layout Standard +With a 'C' line for each peak. +\end_layout + \begin_layout Subsection stfcdiffusion \end_layout diff --git a/src/DataFile.cpp b/src/DataFile.cpp index 6074f00472..d7befd8462 100644 --- a/src/DataFile.cpp +++ b/src/DataFile.cpp @@ -56,7 +56,7 @@ const FileTypes::AllocToken DataFile::DF_AllocArray[] = { { "Grace File", 0, DataIO_Grace::WriteHelp, DataIO_Grace::Alloc }, { "Gnuplot File", 0, DataIO_Gnuplot::WriteHelp,DataIO_Gnuplot::Alloc}, { "Xplor File", 0, 0, DataIO_Xplor::Alloc }, - { "OpenDX File", 0, DataIO_OpenDx::WriteHelp, DataIO_OpenDx::Alloc }, + { "OpenDX File", DataIO_OpenDx::ReadHelp, DataIO_OpenDx::WriteHelp, DataIO_OpenDx::Alloc }, { "Amber REM log", DataIO_RemLog::ReadHelp, 0, DataIO_RemLog::Alloc }, { "Amber MDOUT file", 0, 0, DataIO_Mdout::Alloc }, { "Evecs file", DataIO_Evecs::ReadHelp, 0, DataIO_Evecs::Alloc }, diff --git a/src/DataIO_OpenDx.cpp b/src/DataIO_OpenDx.cpp index 74d4b0d973..33de0ad57f 100644 --- a/src/DataIO_OpenDx.cpp +++ b/src/DataIO_OpenDx.cpp @@ -32,6 +32,11 @@ bool DataIO_OpenDx::ID_DataFormat( CpptrajFile& infile ) { return isDX; } +/** Read help options. */ +void DataIO_OpenDx::ReadHelp() { + mprintf("\ttype {float|double} : Precision to read grid (default float).\n"); +} + /** Process read options. */ int DataIO_OpenDx::processReadArgs(ArgList& argIn) { std::string typestr = argIn.GetStringKey("type"); diff --git a/src/DataIO_OpenDx.h b/src/DataIO_OpenDx.h index cae0cfcf00..5786159a30 100644 --- a/src/DataIO_OpenDx.h +++ b/src/DataIO_OpenDx.h @@ -8,6 +8,7 @@ class DataIO_OpenDx : public DataIO { static BaseIOtype* Alloc() { return (BaseIOtype*)new DataIO_OpenDx(); } int processReadArgs(ArgList&); int ReadData(FileName const&, DataSetList&, std::string const&); + static void ReadHelp(); static void WriteHelp(); int processWriteArgs(ArgList&); int WriteData(FileName const&,DataSetList const&); From d10bb6be867f5773a4d988f73e66b7743e979a07 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 21 Aug 2020 16:21:46 -0400 Subject: [PATCH 6/7] Add description of the sets created by SPAM. --- doc/cpptraj.lyx | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index f9bba6510e..b676fe3513 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -33690,6 +33690,34 @@ dhbulk kcal/mol (SPC/E water). \end_layout +\begin_layout Standard +Data Sets Created for 'purewater': +\end_layout + +\begin_layout Description + Energies for each water at each frame. +\end_layout + +\begin_layout Standard +Data Sets Created otherwise: +\end_layout + +\begin_layout Description +:<#> SPAM energies for peak <#> starting from 1. +\end_layout + +\begin_layout Description +[DG] SPAM delta G values for valid peaks. +\end_layout + +\begin_layout Description +[DH] SPAM delta H values for valid peaks. +\end_layout + +\begin_layout Description +[-TDS] SPAM -T * delta S values for valid peaks. +\end_layout + \end_deeper \begin_layout Standard Perform profiling of bound water molecules via SPAM analysis @@ -33706,6 +33734,9 @@ literal "true" between the water and it's environment from explicit solvent MD trajectories. The interaction energies are calculated using a force- and energy-shifted electrostatic term with a hard cutoff. + For a given peak, SPAM energies will only be calculated for peaks where + the peak is singly-occupied (i.e. + a multiple-occupied peak is not considered valid). \end_layout \begin_layout Standard From 0b1bb0167cd87e9f0c45456681380228d92502bb Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 21 Aug 2020 16:23:13 -0400 Subject: [PATCH 7/7] Increment revision for SPAM/volmap improvements. --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index ab3b47ce35..c0835e7b52 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 "V4.29.4" +#define CPPTRAJ_INTERNAL_VERSION "V4.29.5" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif