Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow atommap rmsfit to proceed with partial atom maps #1114

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Allow atommap rmsfit to proceed with partial atom maps #1114

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a0042c9
Add some checks for atommap with rmsfit. Fix atommap rmsfit when not all
drroe Nov 7, 2024
ab888ee
Version 6.29.8. Revision bump for improved handling of atommap rmsfit
drroe Nov 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 21 additions & 6 deletions src/Action_AtomMap.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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() );
}
}
}
}
}
Expand All @@ -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())
{
Expand All @@ -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());

Expand Down
44 changes: 28 additions & 16 deletions src/Frame.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -774,16 +774,16 @@ void Frame::SetCoordinatesByMap(Frame const& tgtIn, std::vector<int> const& mapI
* mapped: Map[newatom] = refatom (refatom -> this frame).
*/
void Frame::StripUnmappedAtoms(Frame const& refIn, std::vector<int> 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_;
Expand All @@ -796,12 +796,18 @@ void Frame::StripUnmappedAtoms(Frame const& refIn, std::vector<int> const& mapIn

double* newXptr = X_;
double* refptr = refIn.X_;
natom_ = 0;
for (std::vector<int>::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;
}
Expand All @@ -814,11 +820,11 @@ void Frame::StripUnmappedAtoms(Frame const& refIn, std::vector<int> const& mapIn
* according to the given atom map: Map[newatom] = oldatom (oldatom -> frameIn)
*/
void Frame::ModifyByMap(Frame const& frameIn, std::vector<int> 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_;
Expand All @@ -830,16 +836,22 @@ void Frame::ModifyByMap(Frame const& frameIn, std::vector<int> const& mapIn) {
remd_indices_ = frameIn.remd_indices_;

double* Xptr = X_;
natom_ = 0;
for (std::vector<int>::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 -------------------------------------------------
Expand Down
2 changes: 1 addition & 1 deletion src/Version.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.29.7"
#define CPPTRAJ_INTERNAL_VERSION "V6.29.8"
/// PYTRAJ relies on this
#define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
#endif
Loading