Amber-MD · drroe · Jun 17, 2022 · Jun 16, 2022 · Jun 17, 2022 · Jun 17, 2022
diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx
@@ -34750,20 +34750,24 @@ radial
 
 \end_inset
 
- [out <outfilename>] <spacing> <maximum> <solvent mask1>
+ [out <outfilename>] <spacing> <maximum> <solvent mask1> [<solute mask2>]
 \end_layout
 
 \begin_layout LyX-Code
-       [<solute mask2>] [noimage] [density <density> | volume] [<dataset
- name>]
+       [noimage]
 \end_layout
 
 \begin_layout LyX-Code
-       [intrdf <file>] [rawrdf <file>]
+   	[density <density> | volume] [<dataset name>] [intrdf <file>] [rawrdf
+ <file>]
 \end_layout
 
 \begin_layout LyX-Code
-       [{{center1|center2|nointramol} | [byres1] [byres2] [bymol1] [bymol2]}]]
+   	[{{center1|center2|nointramol|toxyz <x>,<y>,<z>} |
+\end_layout
+
+\begin_layout LyX-Code
+         [byres1] [byres2] [bymol1] [bymol2]}]
 \begin_inset Separator latexpar
 \end_inset
 
@@ -34831,7 +34835,25 @@ mask2>] (Optional) If specified calculate RDF of all atoms in <solvent mask1>
 \end_layout
 
 \begin_layout Description
-[<setname>] Name of output data sets.
+[<dataset name>] Name of output data sets.
+\end_layout
+
+\begin_layout Description
+[intrdf
+\begin_inset space ~
+\end_inset
+
+<file>] Calculate integral of RDF bin values (averaged over # of frames
+ but otherwise not normalized) and write to <file> (can be same as <output_filen
+ame>).
+\end_layout
+
+\begin_layout Description
+[rawrdf
+\begin_inset space ~
+\end_inset
+
+<file>] Write raw (non-normalized) RDF values to <file>.
 \end_layout
 
 \begin_layout Description
@@ -34849,21 +34871,12 @@ mask2>] (Optional) If specified calculate RDF of all atoms in <solvent mask1>
 \end_layout
 
 \begin_layout Description
-[intrdf
+[toxyz
 \begin_inset space ~
 \end_inset
 
-<file>] Calculate integral of RDF bin values (averaged over # of frames
- but otherwise not normalized) and write to <file> (can be same as <output_filen
-ame>).
-\end_layout
-
-\begin_layout Description
-[rawrdf
-\begin_inset space ~
-\end_inset
-
-<file>] Write raw (non-normalized) RDF values to <file>.
+<x>,<y>,<z>] Calculate RDF from center of atoms in <solvent mask1> to point
+ specified by <x> <y> and <z> (in Ang.).
 \end_layout
 
 \begin_layout Description

diff --git a/src/Action_Radial.cpp b/src/Action_Radial.cpp
@@ -33,9 +33,11 @@ Action_Radial::Action_Radial() :
 {} 
 
 void Action_Radial::Help() const {
-  mprintf("\t[out <outfilename>] <spacing> <maximum> <solvent mask1> [<solute mask2>] [noimage]\n"
+  mprintf("\t[out <outfilename>] <spacing> <maximum> <solvent mask1> [<solute mask2>]\n"
+          "\t[noimage]\n"
           "\t[density <density> | volume] [<dataset name>] [intrdf <file>] [rawrdf <file>]\n"
-          "\t[{{center1|center2|nointramol} | [byres1] [byres2] [bymol1] [bymol2]}]\n"
+          "\t[{{center1|center2|nointramol|toxyz <x>,<y>,<z>} |\n"
+          "\t  [byres1] [byres2] [bymol1] [bymol2]}]\n"
           "  Calculate the radial distribution function (RDF) of atoms in <solvent mask1>.\n"
           "  If <solute mask2> is given calculate RDF of all atoms in <solvent mask1>\n"
           "  to each atom in <solute mask2>.\n"
@@ -62,7 +64,9 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
   std::string outfilename = actionArgs.GetStringKey("out");
   // Default particle density (mols/Ang^3) for water based on 1.0 g/mL
   density_ = actionArgs.getKeyDouble("density",0.033456);
+
   // Determine mode, by site TODO better integrate with other modes
+  bool needMask2 = true;
   siteMode1_ = OFF;
   siteMode2_ = OFF;
   if (actionArgs.hasKey("byres1")) siteMode1_ = BYRES;
@@ -76,8 +80,21 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
     rmode_ = CENTER2;
   else if (actionArgs.hasKey("nointramol"))
     rmode_ = NO_INTRAMOL;
-  else
+  else if (actionArgs.Contains("toxyz")) {
+    std::string toxyz = actionArgs.GetStringKey("toxyz");
+    ArgList toxyzArg( toxyz, "," );
+    if (toxyzArg.Nargs() != 3) {
+      mprinterr("Error: Expected a comma-separated list of 3 coordinates after 'toxyz'.\n");
+      return Action::ERR;
+    }
+    specified_xyz_[0] = toxyzArg.getNextDouble(0);
+    specified_xyz_[1] = toxyzArg.getNextDouble(0);
+    specified_xyz_[2] = toxyzArg.getNextDouble(0);
+    rmode_ = SPECIFIED;
+    needMask2 = false;
+  } else
     rmode_ = NORMAL;
+
   // Check for mode incompatibility
   if (siteMode1_ != OFF || siteMode2_ != OFF) {
     if (rmode_ != NORMAL) {
@@ -86,6 +103,7 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
     }
     rmode_ = BYSITE;
   }
+
   useVolume_ = actionArgs.hasKey("volume");
   DataFile* intrdfFile = init.DFL().AddDataFile(actionArgs.GetStringKey("intrdf"));
   DataFile* rawrdfFile = init.DFL().AddDataFile(actionArgs.GetStringKey("rawrdf"));
@@ -114,11 +132,13 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
   if (Mask1_.SetMaskString(mask1)) return Action::ERR;
 
   // Check for second mask - if none specified use first mask
-  std::string mask2 = actionArgs.GetMaskNext();
-  if (!mask2.empty()) {
-    if (Mask2_.SetMaskString(mask2)) return Action::ERR;
-  } else {
-    if (Mask2_.SetMaskString(mask1)) return Action::ERR;
+  if (needMask2) {
+    std::string mask2 = actionArgs.GetMaskNext();
+    if (!mask2.empty()) {
+      if (Mask2_.SetMaskString(mask2)) return Action::ERR;
+    } else {
+      if (Mask2_.SetMaskString(mask1)) return Action::ERR;
+    }
   }
   // If filename not yet specified check for backwards compat.
   if (outfilename.empty() && actionArgs.Nargs() > 1 && !actionArgs.Marked(1))
@@ -133,7 +153,10 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
   // Make default precision a little higher than normal
   Dset_->SetupFormat().SetFormatWidthPrecision(12,6);
   // Set DataSet legend from mask strings.
-  Dset_->SetLegend(Mask1_.MaskExpression() + " => " + Mask2_.MaskExpression());
+  if (needMask2)
+    Dset_->SetLegend(Mask1_.MaskExpression() + " => " + Mask2_.MaskExpression());
+  else
+    Dset_->SetLegend(Mask1_.MaskExpression());
   // TODO: Set Yaxis label in DataFile
   // Calculate number of bins
   one_over_spacing_ = 1 / spacing_;
@@ -150,7 +173,10 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
                                                            MetaData::NOT_TS) );
     if (intrdf_ == 0) return Rdf_Err("Could not allocate RDF integral data set.");
     intrdf_->SetupFormat().SetFormatWidthPrecision(12,6);
-    intrdf_->SetLegend("Int[" + Mask2_.MaskExpression() + "]");
+    if (needMask2)
+      intrdf_->SetLegend("Int[" + Mask2_.MaskExpression() + "]");
+    else
+      intrdf_->SetLegend("Int[" + Mask1_.MaskExpression() + "]");
     intrdf_->SetDim(Dimension::X, Rdim);
     intrdfFile->AddDataSet( intrdf_ );
   } else
@@ -189,9 +215,9 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
     rdf_thread_[i].assign( numBins_, 0 );
 # endif /* _OPENMP */
 
-  mprintf("    RADIAL: Calculating RDF for atoms in mask [%s]",Mask1_.MaskString());
-  if (!mask2.empty()) 
-    mprintf(" to atoms in mask [%s]",Mask2_.MaskString());
+  mprintf("    RADIAL: Calculating RDF for atoms in mask1 [%s]",Mask1_.MaskString());
+  if (Mask2_.MaskStringSet()) 
+    mprintf(" to atoms in mask2 [%s]",Mask2_.MaskString());
   mprintf("\n");
   if (outfile != 0)
     mprintf("\tOutput to %s.\n", outfile->DataFilename().full());
@@ -217,6 +243,9 @@ Action::RetType Action_Radial::Init(ArgList& actionArgs, ActionInit& init, int d
       mprintf("\tUsing center of all atoms selected by mask2.\n");
     else if (rmode_==NO_INTRAMOL)
       mprintf("\tIgnoring intramolecular distances.\n");
+    else if (rmode_ == SPECIFIED)
+      mprintf("\tCalculating RDF of atoms selected by mask1 to point %g %g %g\n",
+              specified_xyz_[0], specified_xyz_[1], specified_xyz_[2]);
   }
   mprintf("\tHistogram max %f, spacing %f, bins %i.\n",maximum,
           spacing_,numBins_);
@@ -318,10 +347,12 @@ Action::RetType Action_Radial::Setup(ActionSetup& setup) {
     mprintf("Warning: First mask has no atoms.\n");
     return Action::SKIP;
   }
-  if (setup.Top().SetupIntegerMask( Mask2_ ) ) return Action::ERR;
-  if (Mask2_.None()) {
-    mprintf("Warning: Second mask has no atoms.\n");
-    return Action::SKIP;
+  if (Mask2_.MaskStringSet()) {
+    if (setup.Top().SetupIntegerMask( Mask2_ ) ) return Action::ERR;
+    if (Mask2_.None()) {
+      mprintf("Warning: Second mask has no atoms.\n");
+      return Action::SKIP;
+    }
   }
   imageOpt_.SetupImaging( setup.CoordInfo().TrajBox().HasBox() );
 
@@ -341,6 +372,8 @@ Action::RetType Action_Radial::Setup(ActionSetup& setup) {
   } else if (rmode_ == CENTER2) {
     OuterMask_ = Mask2_;
     InnerMask_ = Mask1_;
+  } else if (rmode_ == SPECIFIED) {
+    InnerMask_ = Mask1_;
   } else if (rmode_ == BYSITE) {
     // One or both masks will be by residue.
     int err = 0;
@@ -393,6 +426,8 @@ Action::RetType Action_Radial::Setup(ActionSetup& setup) {
   if (rmode_ == BYSITE) {
     mprintf("\t%zu sites selected by Mask1 (%i atoms), %zu sites selected by Mask2 (%i atoms)\n",
             Sites1_.size(), Mask1_.Nselected(), Sites2_.size(), Mask2_.Nselected());
+  } else if (rmode_ == SPECIFIED) {
+    mprintf("\t%i atoms in Mask1.\n", Mask1_.Nselected());
   } else {
     mprintf("\t%i atoms in Mask1, %i atoms in Mask2\n",
             Mask1_.Nselected(), Mask2_.Nselected());
@@ -532,9 +567,14 @@ Action::RetType Action_Radial::DoAction(int frameNum, ActionFrame& frm) {
     }
 #   endif
   // ---------------------------------------------
-  } else { // CENTER1 || CENTER2
+  } else { // CENTER1 || CENTER2 || SPECIFIED
     // Calculation of center of one Mask to all atoms in other Mask
-    Vec3 coord_center = frm.Frm().VGeometricCenter(OuterMask_);
+    // or specified point to all atoms in a mask (InnerMask_).
+    Vec3 coord_center;
+    if (rmode_ == SPECIFIED)
+      coord_center = specified_xyz_;
+    else
+      coord_center = frm.Frm().VGeometricCenter(OuterMask_);
     int mask2_max = InnerMask_.Nselected();
 #   ifdef _OPENMP
 #   pragma omp parallel private(nmask2,atom2,D,idx,mythread)
@@ -633,6 +673,10 @@ void Action_Radial::Print() {
     // from mask 2. Assume COM of mask 2 != atom(s) in mask1.
     nmask2 = 1.0;
     numSameAtoms = 0;
+  } else if (rmode_ == SPECIFIED) {
+    // No mask 2; calculated distances to a single point.
+    nmask2 = 1.0;
+    numSameAtoms = 0;
   } else if (rmode_ == BYSITE) {
     // Count sites in common
     nmask1 = (double)Sites1_.size();

diff --git a/src/Action_Radial.h b/src/Action_Radial.h
@@ -35,7 +35,7 @@ class Action_Radial: public Action {
     typedef std::vector<AtomMask> Marray;
     Marray Sites1_;
     Marray Sites2_;
-    enum RmodeType { NORMAL=0, NO_INTRAMOL, CENTER1, CENTER2, BYSITE };
+    enum RmodeType { NORMAL=0, NO_INTRAMOL, CENTER1, CENTER2, BYSITE, SPECIFIED };
     RmodeType rmode_;                ///< Type of calculation to perform.
     enum SmodeType { OFF = 0, BYRES, BYMOL };
     SmodeType siteMode1_;
@@ -55,6 +55,7 @@ class Action_Radial: public Action {
     DataSet* intrdf_;
     DataSet* rawrdf_;
     int debug_;
+    Vec3 specified_xyz_;      ///< XYZ coordinates for SPECIFIED
 
     int SetupSiteArrayByAtom(Marray&, AtomMask const&) const;
     int SetupSiteArrayByRes(Marray&, Topology const&, AtomMask const&) const;

diff --git a/src/Version.h b/src/Version.h
@@ -12,7 +12,7 @@
  * Whenever a number that precedes <revision> is incremented, all subsequent
  * numbers should be reset to 0.
  */
-#define CPPTRAJ_INTERNAL_VERSION "V6.8.0"
+#define CPPTRAJ_INTERNAL_VERSION "V6.8.1"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif
diff --git a/test/Test_Radial/RunTest.sh b/test/Test_Radial/RunTest.sh
@@ -4,7 +4,8 @@
 
 CleanFiles radial.in Radial.agr cRadial.agr WatO-Trp4.agr WatO-Trp4.raw.agr \
            WatO-Trp4.byres.agr WatO-Trp.agr WatO-Trp.volume.agr \
-           WatO-Glu5CD.agr noimage.WatO-Glu5CD.agr
+           WatO-Glu5CD.agr noimage.WatO-Glu5CD.agr point.dat \
+           point?.agr wat.origin.agr
 
 TESTNAME='Radial tests'
 Requires netcdf maxthreads 10
@@ -44,6 +45,32 @@ RunCpptraj "$UNITNAME"
 DoTest WatO-Glu5CD.agr.save WatO-Glu5CD.agr
 DoTest noimage.WatO-Glu5CD.agr.save noimage.WatO-Glu5CD.agr
 
+CheckFor maxthreads 1
+if [ $? -eq 0 ] ; then
+  UNITNAME='Radial test, specified point'
+  cat > radial.in <<EOF
+parm ../tz2.truncoct.parm7
+trajin ../tz2.truncoct.nc 1 1
+
+#radial point0.agr    0.5 10.0 :WAT@O :5@CD
+radial point1.agr    0.5 10.0 :WAT@O toxyz 11.5742,4.3807,-13.7675
+#vector center :5@CD out point.dat
+EOF
+  RunCpptraj "$UNITNAME"
+  DoTest point1.agr.save point1.agr
+fi
+
+UNITNAME='Radial test, specified point after centering'
+cat > radial.in <<EOF
+parm ../tz2.truncoct.parm7
+trajin ../tz2.truncoct.nc
+
+autoimage origin
+radial out wat.origin.agr 0.25 15.0 :WAT@O toxyz 0,0,0
+EOF
+RunCpptraj "$UNITNAME"
+DoTest wat.origin.agr.save wat.origin.agr
+
 EndTest
 
 exit 0
diff --git a/test/Test_Radial/point1.agr.save b/test/Test_Radial/point1.agr.save
@@ -0,0 +1,28 @@
+@with g0
+@  xaxis label "Distance (Ang)"
+@  yaxis label ""
+@  legend 0.2, 0.995
+@  legend char size 0.60
+@  s0 legend ":WAT@O"
+@target G0.S0
+@type xy
+   0.250     0.000000
+   0.750     0.000000
+   1.250     0.000000
+   1.750     0.000000
+   2.250     0.000000
+   2.750     0.000000
+   3.250     1.348481
+   3.750     1.013356
+   4.250     0.263068
+   4.750     0.421297
+   5.250     1.034786
+   5.750     0.575171
+   6.250     0.608590
+   6.750     0.730530
+   7.250     0.452343
+   7.750     0.870933
+   8.250     0.838468
+   8.750     0.683289
+   9.250     0.611434
+   9.750     0.900563