Change underlying data type in radial action. (#868)

* Use unsigned long instead of int to avoid overflow

* Add missing noefile and orderparamfile keywords to help

* Fix IRED manual entry

* Revision bump for changing radial int to unsigned long

* Ensure volume is summed across processes as well

* Remove commented out code

* Revision bump for changing radial hist type from int to unsigned long

Co-authored-by: Daniel R. Roe <daniel.roe@nih.gov>

Author: drroe

doc/cpptraj.lyx

@@ -41458,7 +41458,11 @@ ired
 \end_layout
 
 \begin_layout LyX-Code
-ired [relax freq <MHz> [NHdist <distnh>]] [order <order>]
+ired [relax freq <MHz> [NHdist <distnh>] [noefile <noefilename>]]
+\end_layout
+
+\begin_layout LyX-Code
+     [order <order>] [orderparamfile <orderfilename>]
 \end_layout
 
 \begin_layout LyX-Code
@@ -41485,15 +41489,7 @@ freq
 \begin_inset space ~
 \end_inset
 
-<MHz>
-\begin_inset space ~
-\end_inset
-
-[NHdist
-\begin_inset space ~
-\end_inset
-
-<distnh>]]
+<MHz>]
 \series default
  Should only be used when ired vectors represent N-H bonds; calculate correlatio
 n times 
@@ -41519,55 +41515,82 @@ n times
 \noun default
 \color inherit
 \lang american
- for each eigenmode and relaxation rates and NOEs for each N-H vector.
+ 
+\lang english
+for each eigenmode and relaxation rates and NOEs for each N-H vector.
  'freq <MHz>' (required) is the Lamor frequency of the measurement.
- 'NHdist <distnh>' specifies the length of the NH bond in Angstroms (default
- is 1.02).
 \end_layout
 
+\begin_deeper
 \begin_layout Description
-order
+[NHdist
+\begin_inset space ~
+\end_inset
+
+<distnh>] Specifies the length of the NH bond in Angstroms (default is 1.02).
+\end_layout
+
+\begin_layout Description
+[noefile
+\begin_inset space ~
+\end_inset
+
+<noefilename>] File to write the T1, T2, and NOE data to.
+\end_layout
+
+\end_deeper
+\begin_layout Description
+[order
 \begin_inset space ~
 \end_inset
 
-<order> Order of the Legendre polynomials to use when calculating spherical
+<order>] Order of the Legendre polynomials to use when calculating spherical
  harmonics (default 2).
 \end_layout
 
+\begin_layout Description
+[orderparamfile
+\begin_inset space ~
+\end_inset
+
+<orderfilename>] File to write the S2 data to.
+\end_layout
+
 \begin_layout Description
 
 \series bold
-tstep
+[tstep
 \begin_inset space ~
 \end_inset
 
-<tstep>
+<tstep>]
 \series default
  Time between snapshots in ps (default 1.0).
 \end_layout
 
 \begin_layout Description
 
 \series bold
-tcorr
+[tcorr
 \begin_inset space ~
 \end_inset
 
-<tcorr>
+<tcorr>]
 \series default
  Maximum time to calculate correlation functions for in ps (default 10000.0).
 \end_layout
 
 \begin_layout Description
 
 \series bold
-out
+[out
 \begin_inset space ~
 \end_inset
 
 <filename>
 \series default
- Name of file to write output to.
+ Name of file to write plateau and TauM data.
+ Also the prefix for the .cmt and .cjt files (see below).
 \end_layout
 
 \begin_layout Description
@@ -41698,19 +41721,15 @@ literal "true"
 \noun default
 \color inherit
 \lang american
- will be written to 
-\shape italic
-filename.cmt
-\shape default
-.
- Autocorrelation functions for each vector will be written to 
-\shape italic
-filename.cjt
-\shape default
-.
+ 
+\lang english
+will be written to <filename>.cmt.
+ Autocorrelation functions for each vector will be written to <filename>.cjt.
  Relaxation rates and NOEs for each N-H vector will be written to <filename>
  or added to the the end of the standard output.
- For the calculation of 
+ For the calculation of
+\lang american
+ 
 \family roman
 \series medium
 \shape up
@@ -41734,16 +41753,18 @@ filename.cjt
 \noun default
 \color inherit
 \lang american
- the normalized correlation functions and only the first third of the analyzed
+ 
+\lang english
+the normalized correlation functions and only the first third of the analyzed
  time steps will be used.
  For further information on the convergence of correlation functions see
  [Schneider, Brünger, Nilges, 
-\shape italic
+\emph on
 J.
  Mol.
  Biol.
 
-\shape default
+\emph default
 
 \series bold
 285

src/Action_Radial.cpp

@@ -9,8 +9,9 @@
 
 // CONSTRUCTOR
 Action_Radial::Action_Radial() :
-  RDF_(0),
-  rdf_thread_(0),
+# ifdef _OPENMP
+  threadsCombined_(false),
+# endif
   rmode_(NORMAL),
   currentParm_(0),
   intramol_distances_(0),
@@ -43,18 +44,7 @@ void Action_Radial::Help() const {
           "  of residues/molecules selected by mask1 or mask2.\n");
 }
 
-// DESTRUCTOR
-Action_Radial::~Action_Radial() {
-  //fprintf(stderr,"Radial Destructor.\n");
-  if (RDF_!=0) delete[] RDF_;
-  if (rdf_thread_!=0) {
-    for (int i=0; i < numthreads_; i++)
-      delete[] rdf_thread_[i];
-    delete[] rdf_thread_;
-  }
-}
-
-inline Action::RetType RDF_ERR(const char* msg) {
+inline Action::RetType Rdf_Err(const char* msg) {
   mprinterr("Error: %s\n", msg);
   return Action::ERR;
 }
@@ -136,7 +126,7 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
   // Set up output dataset.
   Dset_ = init.DSL().AddSet( DataSet::DOUBLE, MetaData(actionArgs.GetStringNext(), "",
                                                        MetaData::NOT_TS), "g(r)");
-  if (Dset_ == 0) return RDF_ERR("Could not allocate RDF data set.");
+  if (Dset_ == 0) return Rdf_Err("Could not allocate RDF data set.");
   DataFile* outfile = init.DFL().AddDataFile(outfilename, actionArgs);
   if (outfile != 0) outfile->AddDataSet( Dset_ );
   // Make default precision a little higher than normal
@@ -157,7 +147,7 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
   if (intrdfFile != 0) {
     intrdf_ = init.DSL().AddSet( DataSet::DOUBLE, MetaData(Dset_->Meta().Name(), "int",
                                                            MetaData::NOT_TS) );
-    if (intrdf_ == 0) return RDF_ERR("Could not allocate RDF integral data set.");
+    if (intrdf_ == 0) return Rdf_Err("Could not allocate RDF integral data set.");
     intrdf_->SetupFormat().SetFormatWidthPrecision(12,6);
     intrdf_->SetLegend("Int[" + Mask2_.MaskExpression() + "]");
     intrdf_->SetDim(Dimension::X, Rdim);
@@ -168,7 +158,7 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
   if (rawrdfFile != 0) {
     rawrdf_ = init.DSL().AddSet( DataSet::DOUBLE, MetaData(Dset_->Meta().Name(), "raw",
                                                            MetaData::NOT_TS) );
-    if (rawrdf_ == 0) return RDF_ERR("Could not allocate raw RDF data set.");
+    if (rawrdf_ == 0) return Rdf_Err("Could not allocate raw RDF data set.");
     rawrdf_->SetupFormat().SetFormatWidthPrecision(12,6);
     rawrdf_->SetLegend("Raw[" + Dset_->Meta().Legend() + "]");
     rawrdf_->SetDim(Dimension::X, Rdim);
@@ -182,9 +172,9 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
   if (rawrdf_ != 0) rawrdf_->SetNeedsSync( false );
 # endif
   // Set up histogram
-  RDF_ = new int[ numBins_ ];
-  std::fill(RDF_, RDF_ + numBins_, 0);
+  RDF_.assign( numBins_, 0 );
 # ifdef _OPENMP
+  threadsCombined_ = false;
   // Since RDF is shared by all threads and we cant guarantee that a given
   // bin in RDF wont be accessed at the same time by the same thread,
   // each thread needs its own bin space.
@@ -193,12 +183,10 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
   if (omp_get_thread_num()==0)
     numthreads_ = omp_get_num_threads();
 }
-  rdf_thread_ = new int*[ numthreads_ ];
-  for (int i=0; i < numthreads_; i++) {
-    rdf_thread_[i] = new int[ numBins_ ];
-    std::fill(rdf_thread_[i], rdf_thread_[i] + numBins_, 0);
-  }
-# endif
+  rdf_thread_.resize( numthreads_ );
+  for (int i=0; i < numthreads_; i++)
+    rdf_thread_[i].assign( numBins_, 0 );
+# endif /* _OPENMP */
 
   mprintf("    RADIAL: Calculating RDF for atoms in mask [%s]",Mask1_.MaskString());
   if (!mask2.empty()) 
@@ -590,31 +578,31 @@ int Action_Radial::SyncAction() {
 # ifdef _OPENMP
   CombineRdfThreads();
 # endif
-  int total_frames = 0;
-  trajComm_.ReduceMaster( &total_frames, &numFrames_, 1, MPI_INT, MPI_SUM );
+  double total_volume = 0;
+  trajComm_.ReduceMaster( &total_volume, &volume_, 1, MPI_DOUBLE, MPI_SUM );
+  unsigned long total_frames = 0;
+  trajComm_.ReduceMaster( &total_frames, &numFrames_, 1, MPI_UNSIGNED_LONG, MPI_SUM );
   if (trajComm_.Master()) {
+    volume_ = total_volume;
     numFrames_ = total_frames;
-    int* sum_bins = new int[ numBins_ ];
-    trajComm_.ReduceMaster( sum_bins, RDF_, numBins_, MPI_INT, MPI_SUM );
-    std::copy( sum_bins, sum_bins + numBins_, RDF_ );
-    delete[] sum_bins;
+    Iarray sum_bins( numBins_ );
+    trajComm_.ReduceMaster( &sum_bins[0], &RDF_[0], numBins_, MPI_UNSIGNED_LONG, MPI_SUM );
+    RDF_ = sum_bins;
   } else
-    trajComm_.ReduceMaster( 0,        RDF_, numBins_, MPI_INT, MPI_SUM );
+    trajComm_.ReduceMaster( 0,            &RDF_[0], numBins_, MPI_UNSIGNED_LONG, MPI_SUM );
   return 0;
 }
 #endif
 
 #ifdef _OPENMP
 /** Combine results from each rdf_thread into rdf. */
 void Action_Radial::CombineRdfThreads() {
-  if (rdf_thread_ == 0) return;
+  if (threadsCombined_) return;
   for (int thread = 0; thread < numthreads_; thread++) { 
     for (int bin = 0; bin < numBins_; bin++)
       RDF_[bin] += rdf_thread_[thread][bin];
-    delete[] rdf_thread_[thread];
   }
-  delete[] rdf_thread_;
-  rdf_thread_ = 0;
+  threadsCombined_ = true;
 }
 #endif
 
@@ -629,7 +617,7 @@ void Action_Radial::Print() {
 # ifdef _OPENMP
   CombineRdfThreads(); 
 # endif
-  mprintf("    RADIAL: %i frames,", numFrames_);
+  mprintf("    RADIAL: %lu frames,", numFrames_);
   double nmask1 = (double)Mask1_.Nselected();
   double nmask2 = (double)Mask2_.Nselected();
   int numSameAtoms = 0;
@@ -666,7 +654,7 @@ void Action_Radial::Print() {
 
   // If useVolume, calculate the density from the average volume
   if (useVolume_) {
-    double avgVol = volume_ / numFrames_;
+    double avgVol = volume_ / (double)numFrames_;
     mprintf("\tAverage volume is %f Ang^3.\n",avgVol);
     density_ = (nmask1 * nmask2 - (double)numSameAtoms) / avgVol;
     mprintf("\tAverage density is %f distances / Ang^3.\n",density_);
@@ -684,7 +672,7 @@ void Action_Radial::Print() {
   for (int bin = 0; bin < numBins_; bin++) {
     //mprintf("DBG:\tNumBins= %i\n",rdf[bin]); 
     // Number of particles in this volume slice over all frames.
-    double N = (double) RDF_[bin];
+    double N = (double)RDF_[bin];
     if (rawrdf_ != 0)
       rawrdf_->Add(bin, &N);
     // r, r + dr
@@ -696,7 +684,7 @@ void Action_Radial::Print() {
     double expectedD = dv * density_;
     if (debug_>0)
       mprintf("    \tBin %f->%f <Pop>=%f, V=%f, D=%f, norm %f distances.\n",
-              R,Rdr,N/numFrames_,dv,density_,expectedD);
+              R,Rdr,N/(double)numFrames_,dv,density_,expectedD);
     // Divide by # frames
     double norm = expectedD * (double)numFrames_;
     N /= norm;

src/Action_Radial.h

@@ -6,7 +6,6 @@
 class Action_Radial: public Action {
   public:
     Action_Radial();
-    ~Action_Radial();
     DispatchObject* Alloc() const { return (DispatchObject*)new Action_Radial(); }
     void Help() const;
   private:
@@ -21,18 +20,23 @@ class Action_Radial: public Action {
 #   ifdef _OPENMP
     void CombineRdfThreads();
 #   endif
-    ImagedAction image_;      ///< Image routines.
-    int* RDF_;                ///< Hold bin counts.
-    int** rdf_thread_;        ///< Hold bin count on each thread.
-    AtomMask Mask1_;          ///< Atoms to calculate RDF for.
-    AtomMask Mask2_;          ///< Optional mask to calc RDF to atoms in Mask1.
-    AtomMask OuterMask_;      ///< Mask with the most atoms.
-    AtomMask InnerMask_;      ///< Mask with the fewest atoms.
+    typedef std::vector<unsigned long> Iarray;
+
+#   ifdef _OPENMP
+    bool threadsCombined_;           ///< True if CombineRdfThreads() has been called.
+#   endif
+    ImagedAction image_;             ///< Image routines.
+    Iarray RDF_;                     ///< Hold bin counts.
+    std::vector<Iarray> rdf_thread_; ///< Hold bin count on each thread.
+    AtomMask Mask1_;                 ///< Atoms to calculate RDF for.
+    AtomMask Mask2_;                 ///< Optional mask to calc RDF to atoms in Mask1.
+    AtomMask OuterMask_;             ///< Mask with the most atoms.
+    AtomMask InnerMask_;             ///< Mask with the fewest atoms.
     typedef std::vector<AtomMask> Marray;
     Marray Sites1_;
     Marray Sites2_;
     enum RmodeType { NORMAL=0, NO_INTRAMOL, CENTER1, CENTER2, BYSITE };
-    RmodeType rmode_;         ///< Type of calculation to perform.
+    RmodeType rmode_;                ///< Type of calculation to perform.
     enum SmodeType { OFF = 0, BYRES, BYMOL };
     SmodeType siteMode1_;
     SmodeType siteMode2_;
@@ -45,7 +49,7 @@ class Action_Radial: public Action {
     double one_over_spacing_; ///< 1/spacing, used to avoid man division ops.
     int numBins_;             ///< The number of bins.
     int numthreads_;          ///< Number of threads.
-    int numFrames_;           ///< Number of frames for which RDF is calcd.
+    unsigned long numFrames_; ///< Number of frames for which RDF is calcd.
     double density_;          ///< Particle density (molecules/Ang^3).
     DataSet* Dset_;
     DataSet* intrdf_;