diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index b719c6d5f2..b676fe3513 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 @@ -33336,129 +33401,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 +33564,229 @@ 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 + +\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 +\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. + 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 +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 +<# 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 @@ -35068,7 +35233,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 +35249,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 +35272,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 +35284,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 +35396,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 +35423,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 +35438,7 @@ Specify the grid center explicitly. Default is the origin. \end_layout +\end_deeper \begin_layout Description \family roman @@ -35204,6 +35453,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 +35470,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 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/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/DataFile.cpp b/src/DataFile.cpp index 3fefcf445a..d7befd8462 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() : @@ -55,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 }, @@ -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..33de0ad57f 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,34 @@ 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"); + 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 +72,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 +180,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 +195,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 +255,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 +272,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 +328,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 +346,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 +360,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 +385,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..5786159a30 100644 --- a/src/DataIO_OpenDx.h +++ b/src/DataIO_OpenDx.h @@ -4,10 +4,11 @@ /// 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 ReadHelp(); static void WriteHelp(); int processWriteArgs(ArgList&); int WriteData(FileName const&,DataSetList const&); @@ -18,9 +19,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) { 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 diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 3ce10297e1..e661447540 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 @@ -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 @@ -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 \