diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index d37a598123..3985e5aa19 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -34750,20 +34750,24 @@ radial \end_inset - [out ] + [out ] [] \end_layout \begin_layout LyX-Code - [] [noimage] [density | volume] [] + [noimage] \end_layout \begin_layout LyX-Code - [intrdf ] [rawrdf ] + [density | volume] [] [intrdf ] [rawrdf + ] \end_layout \begin_layout LyX-Code - [{{center1|center2|nointramol} | [byres1] [byres2] [bymol1] [bymol2]}]] + [{{center1|center2|nointramol|toxyz ,,} | +\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 \end_layout \begin_layout Description -[] Name of output data sets. +[] Name of output data sets. +\end_layout + +\begin_layout Description +[intrdf +\begin_inset space ~ +\end_inset + +] Calculate integral of RDF bin values (averaged over # of frames + but otherwise not normalized) and write to (can be same as ). +\end_layout + +\begin_layout Description +[rawrdf +\begin_inset space ~ +\end_inset + +] Write raw (non-normalized) RDF values to . \end_layout \begin_layout Description @@ -34849,21 +34871,12 @@ mask2>] (Optional) If specified calculate RDF of all atoms in \end_layout \begin_layout Description -[intrdf +[toxyz \begin_inset space ~ \end_inset -] Calculate integral of RDF bin values (averaged over # of frames - but otherwise not normalized) and write to (can be same as ). -\end_layout - -\begin_layout Description -[rawrdf -\begin_inset space ~ -\end_inset - -] Write raw (non-normalized) RDF values to . +,,] Calculate RDF from center of atoms in to point + specified by and (in Ang.). \end_layout \begin_layout Description diff --git a/src/Action_Radial.cpp b/src/Action_Radial.cpp index c040ac3458..30cb36cf91 100644 --- 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 ] [] [noimage]\n" + mprintf("\t[out ] []\n" + "\t[noimage]\n" "\t[density | volume] [] [intrdf ] [rawrdf ]\n" - "\t[{{center1|center2|nointramol} | [byres1] [byres2] [bymol1] [bymol2]}]\n" + "\t[{{center1|center2|nointramol|toxyz ,,} |\n" + "\t [byres1] [byres2] [bymol1] [bymol2]}]\n" " Calculate the radial distribution function (RDF) of atoms in .\n" " If is given calculate RDF of all atoms in \n" " to each atom in .\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 index 24e83df6b3..7982d7c636 100644 --- a/src/Action_Radial.h +++ b/src/Action_Radial.h @@ -35,7 +35,7 @@ class Action_Radial: public Action { typedef std::vector 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 index d856148aca..bcd2c18817 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V6.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 index 2a08456133..9d2b00ef9b 100755 --- 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 < radial.in <