From a0042c924f3500038a5643b9dcc2d39070be9bf4 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 7 Nov 2024 13:11:56 -0500 Subject: [PATCH 1/2] Add some checks for atommap with rmsfit. Fix atommap rmsfit when not all atoms are mapped (allow it to proceed with whatever is mapped). --- src/Action_AtomMap.cpp | 27 ++++++++++++++++++++------ src/Frame.cpp | 44 +++++++++++++++++++++++++++--------------- 2 files changed, 49 insertions(+), 22 deletions(-) diff --git a/src/Action_AtomMap.cpp b/src/Action_AtomMap.cpp index ea747d9756..eb7fbffad2 100644 --- a/src/Action_AtomMap.cpp +++ b/src/Action_AtomMap.cpp @@ -56,6 +56,10 @@ Action::RetType Action_AtomMap::Init(ArgList& actionArgs, ActionInit& init, int maponly_ = actionArgs.hasKey("maponly"); rmsfit_ = actionArgs.hasKey("rmsfit"); renameAtoms_ = actionArgs.hasKey("changenames"); + if (rmsfit_ && renameAtoms_) { + mprinterr("Error: Specify either 'rmsfit' or 'changenames', not both.\n"); + return Action::ERR; + } std::string modestring = actionArgs.GetStringKey("mode"); if (!modestring.empty()) { if (modestring == "all") mode_ = ALL; @@ -177,6 +181,11 @@ Action::RetType Action_AtomMap::Init(ArgList& actionArgs, ActionInit& init, int // performed using all atoms that were successfully mapped from // target to reference. if (rmsfit_) { + if (Mapper.Nmapped() < 3) { + mprinterr("Error: Only %i atoms mapped, require at least 3 for RMS-fit.\n", + Mapper.Nmapped()); + return Action::ERR; + } rmsRefFrame_.SetupFrame( Mapper.Nmapped() ); rmsTgtFrame_ = rmsRefFrame_; // Set up a reference frame containing only mapped reference atoms @@ -221,13 +230,16 @@ Action::RetType Action_AtomMap::Init(ArgList& actionArgs, ActionInit& init, int newFrame_->SetupFrameM( TgtFrame_->Top().Atoms() ); // Set up new Parm newParm_ = TgtFrame_->Top().ModifyByMap(AMap_); + //newParm_->Summary(); // DEBUG if (renameAtoms_) { // newParm_ is already in the correct order for (int refatom = 0; refatom != (int)AMap_.size(); ++refatom) { if (AMap_[refatom] != -1) { - //mprintf("DEBUG: ref %i name %s becomes tgt %i\n", refatom+1, *(RefFrame_->Top()[refatom].Name()), targetatom+1); + //mprintf("DEBUG: %i name %s becomes ref name %s\n", refatom+1, + // *((*newParm_)[refatom].Name()), + // *(RefFrame_->Top()[refatom].Name())); (*newParm_).SetAtom(refatom).SetName( RefFrame_->Top()[refatom].Name() ); - } + } } } } @@ -240,10 +252,16 @@ Action::RetType Action_AtomMap::Init(ArgList& actionArgs, ActionInit& init, int * replace current parm with mapped parm. */ Action::RetType Action_AtomMap::Setup(ActionSetup& setup) { + if (rmsfit_) { + mprintf("\trmsfit specified, %i atoms.\n",rmsRefFrame_.Natom()); + return Action::OK; + } + if (maponly_) { mprintf("\tmaponly was specified, not using atom map during traj read.\n"); return Action::OK; } + if (setup.Top().Pindex() != TgtFrame_->Top().Pindex() || setup.Top().Natom() != TgtFrame_->Top().Natom()) { @@ -254,10 +272,7 @@ Action::RetType Action_AtomMap::Setup(ActionSetup& setup) { mprintf("Warning: Not using map for this topology.\n"); return Action::SKIP; } - if (rmsfit_) { - mprintf("\trmsfit specified, %i atoms.\n",rmsRefFrame_.Natom()); - return Action::OK; - } + mprintf("\tMap for parm %s -> %s (%i atom).\n",TgtFrame_->Top().c_str(), RefFrame_->Top().c_str(), TgtFrame_->Top().Natom()); diff --git a/src/Frame.cpp b/src/Frame.cpp index 7e3e69363c..475c4174c1 100644 --- a/src/Frame.cpp +++ b/src/Frame.cpp @@ -774,16 +774,16 @@ void Frame::SetCoordinatesByMap(Frame const& tgtIn, std::vector const& mapI * mapped: Map[newatom] = refatom (refatom -> this frame). */ void Frame::StripUnmappedAtoms(Frame const& refIn, std::vector const& mapIn) { - if (refIn.natom_ > maxnatom_) { - mprinterr("Error: StripUnmappedAtoms: # Input map frame atoms (%i) > max atoms (%i)\n", - refIn.natom_, maxnatom_); - return; - } - if ((int)mapIn.size() != refIn.natom_) { - mprinterr("Error: StripUnmappedAtoms: Input map size (%zu) != input frame natom (%i)\n", - mapIn.size(), refIn.natom_); - return; - } + //if (refIn.natom_ > maxnatom_) { + // mprinterr("Error: StripUnmappedAtoms: # Input map frame atoms (%i) > max atoms (%i)\n", + // refIn.natom_, maxnatom_); + // return; + //} + //if ((int)mapIn.size() != refIn.natom_) { + // mprinterr("Error: StripUnmappedAtoms: Input map size (%zu) != input frame natom (%i)\n", + // mapIn.size(), refIn.natom_); + // return; + //} step_ = refIn.step_; box_ = refIn.box_; T_ = refIn.T_; @@ -796,12 +796,18 @@ void Frame::StripUnmappedAtoms(Frame const& refIn, std::vector const& mapIn double* newXptr = X_; double* refptr = refIn.X_; + natom_ = 0; for (std::vector::const_iterator refatom = mapIn.begin(); refatom != mapIn.end(); ++refatom) { if (*refatom != -1) { + if (natom_ == maxnatom_) { + mprinterr("Error: StripUnmappedAtoms: More mapped atoms than this frame max natom (%i), map size %zu.\n", maxnatom_, mapIn.size()); + return; + } memcpy( newXptr, refptr, COORDSIZE_ ); newXptr += 3; + natom_++; } refptr += 3; } @@ -814,11 +820,11 @@ void Frame::StripUnmappedAtoms(Frame const& refIn, std::vector const& mapIn * according to the given atom map: Map[newatom] = oldatom (oldatom -> frameIn) */ void Frame::ModifyByMap(Frame const& frameIn, std::vector const& mapIn) { - if ((int)mapIn.size() > maxnatom_) { - mprinterr("Error: SetTargetByMap: Input map size (%zu) > this frame max natom (%i)\n", - mapIn.size(), maxnatom_); - return; - } + //if ((int)mapIn.size() > maxnatom_) { + // mprinterr("Error: SetTargetByMap: Input map size (%zu) > this frame max natom (%i)\n", + // mapIn.size(), maxnatom_); + // return; + //} step_ = frameIn.step_; box_ = frameIn.box_; T_ = frameIn.T_; @@ -830,16 +836,22 @@ void Frame::ModifyByMap(Frame const& frameIn, std::vector const& mapIn) { remd_indices_ = frameIn.remd_indices_; double* Xptr = X_; + natom_ = 0; for (std::vector::const_iterator oldatom = mapIn.begin(); oldatom != mapIn.end(); ++oldatom) { if (*oldatom != -1) { + if (natom_ == maxnatom_) { + mprinterr("Error: ModifyByMap: More mapped atoms than this frame max natom (%i), map size %zu.\n", maxnatom_, mapIn.size()); + return; + } memcpy( Xptr, frameIn.X_ + ((*oldatom) * 3), COORDSIZE_ ); Xptr += 3; + natom_++; } } ncoord_ = (int)(Xptr - X_); - natom_ = ncoord_ / 3; + //natom_ = ncoord_ / 3; } // ---------- BASIC ARITHMETIC ------------------------------------------------- From ab888ee3f6263869ec3ebeadfa62aece45fad8f7 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 7 Nov 2024 13:40:38 -0500 Subject: [PATCH 2/2] Version 6.29.8. Revision bump for improved handling of atommap rmsfit when not all atoms mapped. --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index dbf3238cbc..4817cfae43 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.29.7" +#define CPPTRAJ_INTERNAL_VERSION "V6.29.8" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif