diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf index 9b3e7cd2c2..77fc309b0d 100644 Binary files a/doc/CpptrajManual.pdf and b/doc/CpptrajManual.pdf differ diff --git a/doc/DocumentChecksums.txt b/doc/DocumentChecksums.txt index 82d9c91a21..7491bdc6e0 100644 --- a/doc/DocumentChecksums.txt +++ b/doc/DocumentChecksums.txt @@ -1,3 +1,3 @@ f6f8cb1a79951d80a9d2656fd9c30f55 CpptrajDevelopmentGuide.lyx -3b89e9556662b1c6de2d0e38ab894588 cpptraj.lyx +369e9afee9918ff0797c4a4bf08473fd cpptraj.lyx 5d9b5b5ed47a3ded57b6464df99b3585 CpptrajManual.lyx diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index e5cae04ac0..9169966825 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -7501,7 +7501,7 @@ Solvation/box: \end_layout \begin_layout Description -solvatebox|solvateoct Solvate the system in either an orthorhombic or truncated octadedral box. +solvatebox|solvateoct Solvate the system in either an orthorhombic or truncated octahedral box. \end_layout \begin_deeper @@ -7551,8 +7551,7 @@ nsolvent <#> EXPERIMENTAL. Instead of a buffer, specify a target number of waters to solvate the system with. - Currently only works for orthorhombic cells (i.e. - not truncated octahedron). + \end_layout \end_deeper @@ -7769,14 +7768,23 @@ fixatomorder The build procedure is similar to that of LEaP, but CPPTRAJ's build tends to be faster and more robust, particularly for large systems and/or system containing CMAP terms. - Note that since LEaP uses an imprecise value of PI compared to CPPTRAJ/sander/pmemd, + +\end_layout + +\begin_layout Standard +Note that since LEaP uses an imprecise value of PI compared to CPPTRAJ/sander/pmemd, the topologies produced by CPPTRAJ will have slight differences in any terms involving PI (e.g. angles, dihedrals, etc). For 1:1 direct comparison to LEaP topologies, CPPTRAJ can be compiled using LEaP's value of PI by configuring with the '-leappi' configure flag, - which activates the CPPTRAJ_USE_LEAP_PI compiler define. + which activates the CPPTRAJ_USE_LEAP_PI compiler define; + +\series bold +this should only be used to compare directly to LEaP +\series default +. Note also that CPPTRAJ retains PDB information (chain ID, original residue number, Bfactor/Occupancy etc) which is printed to the output Topology. @@ -7898,6 +7906,16 @@ addionsrand ion2, or both. Note that currently CPPTRAJ does not support adding ions if no solvent has been added. + Either a +\series bold +buffer +\series default + size (in Angstroms) must be given, + or a target number of waters given with +\series bold +nsolvent +\series default +. \end_layout \begin_layout Enumerate diff --git a/src/Action_MinImage.cpp b/src/Action_MinImage.cpp index 045adb70a6..70b5185b25 100644 --- a/src/Action_MinImage.cpp +++ b/src/Action_MinImage.cpp @@ -11,7 +11,8 @@ Action_MinImage::Action_MinImage() : dist_(0), useMass_(true), - calcUsingMask_(false) + calcUsingMask_(false), + toFace_(false) // HIDDEN, TODO make not hidden after testing {} void Action_MinImage::Help() const { @@ -25,28 +26,43 @@ Action::RetType Action_MinImage::Init(ArgList& actionArgs, ActionInit& init, int // Get Keywords useMass_ = !(actionArgs.hasKey("geom")); calcUsingMask_ = actionArgs.hasKey("maskcenter"); + toFace_ = actionArgs.hasKey("toface"); + if (calcUsingMask_ && toFace_) { + mprinterr("Error: Specify either 'maskcenter' or 'toface', not both.\n"); + return Action::ERR; + } DataFile* outfile = init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs ); // Get Masks std::string mask1 = actionArgs.GetMaskNext(); - std::string mask2 = actionArgs.GetMaskNext(); - if (mask1.empty() || mask2.empty()) { - mprinterr("Error: Requires 2 masks\n"); - return Action::ERR; - } if (Mask1_.SetMaskString(mask1)) return Action::ERR; - if (Mask2_.SetMaskString(mask2)) return Action::ERR; + std::string mask2; + if (!toFace_) { + mask2 = actionArgs.GetMaskNext(); + if (mask1.empty() || mask2.empty()) { + mprinterr("Error: Requires 2 masks\n"); + return Action::ERR; + } + if (Mask2_.SetMaskString(mask2)) return Action::ERR; + } // Dataset to store distances MetaData md(actionArgs.GetStringNext()); md.SetScalarMode( MetaData::M_DISTANCE ); dist_ = init.DSL().AddSet(DataSet::DOUBLE, md, "MID"); if (dist_==0) return Action::ERR; + // Update metadata in case a default name was generated + md.SetName( dist_->Meta().Name() ); if (outfile != 0) outfile->AddDataSet( dist_ ); - if (!calcUsingMask_) { + if (toFace_ || !calcUsingMask_) { md.SetAspect("A1"); atom1_ = init.DSL().AddSet(DataSet::INTEGER, md); - md.SetAspect("A2"); - atom2_ = init.DSL().AddSet(DataSet::INTEGER, md); + if (toFace_) { + md.SetAspect("FACE"); + atom2_ = init.DSL().AddSet(DataSet::STRING, md); + } else { + md.SetAspect("A2"); + atom2_ = init.DSL().AddSet(DataSet::INTEGER, md); + } if (atom1_ == 0 || atom2_ == 0) return Action::ERR; if (outfile != 0) { outfile->AddDataSet( atom1_ ); @@ -66,7 +82,9 @@ Action::RetType Action_MinImage::Init(ArgList& actionArgs, ActionInit& init, int minAtom2_.resize( numthreads ); mprintf(" MINIMAGE: Looking for closest approach of"); - if (calcUsingMask_) { + if (toFace_) { + mprintf(" atoms in mask '%s' to unit cell faces.\n", Mask1_.MaskString()); + } else if (calcUsingMask_) { mprintf(" center of mask %s\n\tto images of center of mask %s\n", Mask1_.MaskString(), Mask2_.MaskString()); if (useMass_) @@ -92,12 +110,21 @@ Action::RetType Action_MinImage::Setup(ActionSetup& setup) { return Action::SKIP; } if (setup.Top().SetupIntegerMask( Mask1_ )) return Action::ERR; - if (setup.Top().SetupIntegerMask( Mask2_ )) return Action::ERR; - mprintf("\t%s (%i atoms) to %s (%i atoms)\n",Mask1_.MaskString(), Mask1_.Nselected(), - Mask2_.MaskString(),Mask2_.Nselected()); - if (Mask1_.None() || Mask2_.None()) { - mprintf("Warning: One or both masks have no atoms.\n"); - return Action::SKIP; + if (!toFace_) { + if (setup.Top().SetupIntegerMask( Mask2_ )) return Action::ERR; + mprintf("\t%s (%i atoms) to %s (%i atoms)\n", + Mask1_.MaskString(), Mask1_.Nselected(), + Mask2_.MaskString(), Mask2_.Nselected()); + if (Mask1_.None() || Mask2_.None()) { + mprintf("Warning: One or both masks have no atoms.\n"); + return Action::SKIP; + } + } else { + mprintf("\t%s (%i atoms)\n", Mask1_.MaskString(), Mask1_.Nselected()); + if (Mask1_.None()) { + mprintf("Warning: Mask has no atoms.\n"); + return Action::SKIP; + } } return Action::OK; @@ -116,6 +143,75 @@ static void WriteMatrix(Matrix_3x3 const& ucell, PDBfile& pdbout, const char* na } */ +/** Distance of each atom in mask to each unit cell face. */ +double Action_MinImage::minDistToCellFace(AtomMask const& maskIn, Frame const& frameIn, int frameNum) +{ + Box const& boxIn = frameIn.BoxCrd(); + // Store face coords in cartesian space + static const char* faceStr[] = { "+X", "-X", "+Y", "-Y", "+Z", "-Z" }; + Vec3 Faces[6]; + Faces[0] = boxIn.UnitCell().TransposeMult( Vec3(1.0, 0.5, 0.5) ); // +X + Faces[1] = boxIn.UnitCell().TransposeMult( Vec3(0.0, 0.5, 0.5) ); // -X + Faces[2] = boxIn.UnitCell().TransposeMult( Vec3(0.5, 1.0, 0.5) ); // +Y + Faces[3] = boxIn.UnitCell().TransposeMult( Vec3(0.5, 0.0, 0.5) ); // -Y + Faces[4] = boxIn.UnitCell().TransposeMult( Vec3(0.5, 0.5, 1.0) ); // +Z + Faces[5] = boxIn.UnitCell().TransposeMult( Vec3(0.5, 0.5, 0.0) ); // -Z + + minDist_.assign(minDist_.size(), DBL_MAX); + int mythread = 0; + for (AtomMask::const_iterator at = maskIn.begin(); at != maskIn.end(); ++at) + { + // To fractional space + Vec3 frac = boxIn.FracCell() * Vec3( frameIn.XYZ(*at) ); + // Bring back into main unit cell + frac[0] = frac[0] - floor(frac[0]); + frac[1] = frac[1] - floor(frac[1]); + frac[2] = frac[2] - floor(frac[2]); + if (frac[0] < 0.0) frac[0] += 1.0; + if (frac[1] < 0.0) frac[1] += 1.0; + if (frac[2] < 0.0) frac[2] += 1.0; + // Back to Cartesian + Vec3 cart = boxIn.UnitCell().TransposeMult( frac ); + // Distance to each face + for (unsigned int idx = 0; idx < 6; idx++) { + Vec3 dxyz = Faces[idx] - cart; + double dist2 = dxyz.Magnitude2(); +// mprintf("DEBUG: %3i Dist from at %i to face %s is %g face={%g %g %g} cart={%g %g %g}\n", frameNum+1, *at+1, faceStr[idx], sqrt(dist2), +// Faces[idx][0], Faces[idx][1], Faces[idx][2], cart[0], cart[1], cart[2]); + //minDist2 = std::min(minDist2, dist2); + if (dist2 < minDist_[mythread]) { + minDist_[mythread] = dist2; + minAtom1_[mythread] = *at; + minAtom2_[mythread] = idx; + //rprintf("DEBUG: New Min Dist: Atom %i to %i (%g)\n", + // Mask1_[m1]+1,Mask2_[m2]+1,sqrt(Dist2)); + //pdbout_.WriteHET(1, minxyz_[0], minxyz_[1], minxyz_[2]); + } + } + } + // Find lowest minDist + double globalMin = minDist_[0]; + int min1 = minAtom1_[0]; + int min2 = minAtom2_[0]; + for (unsigned int n = 1; n != minDist_.size(); n++) + if (minDist_[n] < globalMin) { + globalMin = minDist_[n]; + min1 = minAtom1_[n]; + min2 = minAtom2_[n]; + } + ++min1; + //++min2; + // By convention make atom 1 the lowest number + //int lowest_num = std::min(min1, min2); + //int highest_num = std::max(min1, min2); + //atom1_->Add(frameNum, &lowest_num); + //atom2_->Add(frameNum, &highest_num); + atom1_->Add(frameNum, &min1); + atom2_->Add(frameNum, faceStr[min2]); + + return globalMin; +} + double Action_MinImage::MinNonSelfDist2(Vec3 const& a1, Vec3 const& a2, Box const& boxIn) { // a1.Print("A1"); // a2.Print("A2"); @@ -165,7 +261,10 @@ Action::RetType Action_MinImage::DoAction(int frameNum, ActionFrame& frm) { // pdbout_.OpenWrite("minimage.pdb"); double Dist2; - if (calcUsingMask_) { + if (toFace_) { + Dist2 = minDistToCellFace(Mask1_, frm.Frm(), frameNum); + Dist2 = sqrt( Dist2 ); + } else if (calcUsingMask_) { // Use center of mask1 and mask2 Vec3 a1, a2; if (useMass_) { diff --git a/src/Action_MinImage.h b/src/Action_MinImage.h index 4dbc606c8b..1b45fdbd13 100644 --- a/src/Action_MinImage.h +++ b/src/Action_MinImage.h @@ -14,6 +14,8 @@ class Action_MinImage: public Action { Action::RetType DoAction(int, ActionFrame&); void Print() {} + double minDistToCellFace(AtomMask const&, Frame const&, int); + double MinNonSelfDist2(Vec3 const&, Vec3 const&, Box const&); DataSet* dist_; ///< Will hold DataSet of calculated distances. @@ -21,6 +23,7 @@ class Action_MinImage: public Action { DataSet* atom2_; bool useMass_; ///< If true, mass-weight distances. bool calcUsingMask_; ///< If true use center of masks + bool toFace_; AtomMask Mask1_; AtomMask Mask2_; std::vector minDist_; diff --git a/src/Structure/Solvate.cpp b/src/Structure/Solvate.cpp index ae5029ba53..702fc525fa 100644 --- a/src/Structure/Solvate.cpp +++ b/src/Structure/Solvate.cpp @@ -6,14 +6,15 @@ #include "../DataSet_Coords.h" #include "../DataSetList.h" #include "../Frame.h" +# ifdef CPPTRAJ_DEBUG_SOLVATE +#include "../StringRoutines.h" // DEBUG, integerToString +#include "../Trajout_Single.h" // DEBUG +#endif #include "../Topology.h" #include "../Parm/ParameterSet.h" #include //std::max #include // cos, sin, sqrt -// DEBUG -//#inc lude "../Trajout_Single.h" - using namespace Cpptraj::Structure; /** CONSTRUCTOR */ @@ -71,10 +72,6 @@ int Solvate::InitSolvate(ArgList& argIn, bool octIn, int debugIn) { doTruncatedOct_ = octIn; nsolvent_ = (unsigned int)argIn.getKeyInt("nsolvent", 0); if (nsolvent_ > 0) { - if (doTruncatedOct_) { - mprinterr("Error: 'nsolvent' currently does not work for truncated octahedral cells.\n"); - return 1; - } clip_ = false; mprintf("Warning: The 'nsolvent' functionality is currently EXPERIMENTAL.\n"); } @@ -338,8 +335,9 @@ int Solvate::SetVdwBoundingBox(Topology& topOut, Frame& frameOut, Cpptraj::Parm: } /** Scale buffer if needed to meet diagonal clearance. */ -void Solvate::octBoxCheck(Frame const& frameOut) { - if (frameOut.Natom() < 1) return; +double Solvate::octBoxCheck(Frame const& frameOut, double dXWidth, double dYWidth, double dZWidth, std::vector const& soluteRadii) +{ + if (frameOut.Natom() < 1) return 0; double dPBuf[4]; dPBuf[0] = bufferX_; @@ -348,7 +346,7 @@ void Solvate::octBoxCheck(Frame const& frameOut) { dPBuf[3] = bufferD_; //mprintf("dPBuf %f %f %f %f\n", dPBuf[0], dPBuf[1], dPBuf[2], dPBuf[3]); - const double* XYZ = frameOut.XYZ(0); +/* const double* XYZ = frameOut.XYZ(0); double dXmax = XYZ[0]; double dYmax = XYZ[1]; double dZmax = XYZ[2]; @@ -371,12 +369,20 @@ void Solvate::octBoxCheck(Frame const& frameOut) { dZmax = XYZ[2]; else if ( XYZ[2] < dZmin ) dZmin = XYZ[2]; - } + }*/ + double dTemp = std::max(dXWidth, dYWidth); + dTemp = std::max(dTemp, dZWidth); + dXWidth = dYWidth = dZWidth = dTemp; // calc halfbox on centers +/* double dXhalf = 0.5 * (dXmax - dXmin) + dPBuf[0]; double dYhalf = 0.5 * (dYmax - dYmin) + dPBuf[1]; double dZhalf = 0.5 * (dZmax - dZmin) + dPBuf[2]; +*/ + double dXhalf = 0.5 * dXWidth + dPBuf[0]; + double dYhalf = 0.5 * dYWidth + dPBuf[1]; + double dZhalf = 0.5 * dZWidth + dPBuf[2]; // find unit vector of diagonal double dX = dYhalf * dZhalf; @@ -393,32 +399,50 @@ void Solvate::octBoxCheck(Frame const& frameOut) { double dMax = 0.0; for (int at = 0; at < frameOut.Natom(); at++) { - XYZ = frameOut.XYZ(at); + const double* XYZ = frameOut.XYZ(at); + dX = XYZ[0]; + dY = XYZ[1]; + dZ = XYZ[2]; +/* dTmp = fabs(XYZ[0]) * dXunit + fabs(XYZ[1]) * dYunit + fabs(XYZ[2]) * dZunit; if ( dTmp > dMax ) dMax = dTmp; +*/ + double dR = soluteRadii[at]; + dTmp = fabs(dX) * dXunit + + fabs(dY) * dYunit + + fabs(dZ) * dZunit + + +dR; + + if (dTmp > dMax) + dMax = dTmp; } // calc distance of diagonal face from origin - double dBmax = 0.5 * sqrt( dXhalf*dXhalf + dYhalf*dYhalf + dZhalf*dZhalf ); + //double dBmax = 0.5 * sqrt( dXhalf*dXhalf + dYhalf*dYhalf + dZhalf*dZhalf ); + double dBmax = 1.5 / sqrt(1 / (dXhalf * dXhalf) + 1 / (dYhalf * dYhalf) + 1 / (dZhalf * dZhalf)); // see if diagonal clearance is satisfied dTmp = dMax + dPBuf[3]; if ( dTmp <= dBmax ) { if ( dPBuf[3] == 0.0 ) mprintf("\t (Diagonal clearance is %f)\n", dMax ); - return; + return 1.0; } // not satisfied: scale up box dTmp /= dBmax; mprintf("\t Scaling up box by a factor of %f to meet diagonal cut criterion\n", dTmp ); - bufferX_ *= dTmp; - bufferY_ *= dTmp; - bufferZ_ *= dTmp; + bufferX_ = dXhalf * dTmp - (0.5 * dXWidth); + bufferY_ = dYhalf * dTmp - (0.5 * dYWidth); + bufferZ_ = dZhalf * dTmp - (0.5 * dZWidth); + //bufferX_ += dXhalf * (dTmp - 1.0); //bufferX_ *= dTmp; + //bufferY_ += dYhalf * (dTmp - 1.0); //bufferY_ *= dTmp; + //bufferZ_ += dZhalf * (dTmp - 1.0); //bufferZ_ *= dTmp; + return dTmp; } /** Rotate 45 deg. around z axis, (90-tetra/2) around y axis, 90 around x axis. @@ -541,10 +565,60 @@ unsigned int Solvate::nBoxesInCube(int n) { return cube; } -static inline void checkBuffer(double& buf, double min, double max) { - if (buf < min) buf = min; - if (buf > max) buf = max; +static inline void checkBuffer(double& buf, double min, double max, double lastbuffer, double& change, int diff) +{ + //if (buf < min) buf = min; + //if (buf > max) buf = max; + bool updateChangeVal = false; + if (buf > max) { + // Pick a target halfway between last buffer and max instead + double halftgt = ((max - lastbuffer) / 2.0) + lastbuffer; + mprintf("DEBUG: New buffer %g would have been > max %g; choosing %g instead.\n", + buf, max, halftgt); + buf = halftgt; + updateChangeVal = true; + } else if (buf < min) { + // pick a target halfway between min and last buffer instead + double halftgt = lastbuffer - ((lastbuffer - min) / 2.0); + mprintf("DEBUG: New buffer %g would have been < min %g; choosing %g instead.\n", + buf, min, halftgt); + buf = halftgt; + updateChangeVal = true; + } + if (updateChangeVal) { + // Update the change value TODO is this ok? + double newchange = fabs((buf - lastbuffer) / (double)diff); + mprintf("DEBUG: Old change %g new change %g\n", change, newchange); + change = newchange; + } +} + +#ifdef CPPTRAJ_DEBUG_SOLVATE +static inline void printIntermediateStructure(int idx, CharMask const& cmask, Topology const& topOut, Frame const& frameOut) { +// cmask.MaskInfo(); + AtomMask imask( cmask.ConvertToIntMask(), topOut.Natom() ); +// imask.MaskInfo(); + Topology* newParm = topOut.modifyStateByMask( imask ); + if (newParm == 0) { + mprinterr("Error: printIntermediateStructure: Could not create topology with the desired # of solvent.\n"); + return; + } + //newParm->Brief("Topology with target # of solvent:"); + + Frame newFrame; + newFrame.SetupFrameV(newParm->Atoms(), frameOut.CoordsInfo()); + newFrame.SetFrame( frameOut, imask ); + + // DEBUG + std::string fnameOut("nsolvent." + integerToString(idx) + ".mol2"); + Trajout_Single trajout; + trajout.InitTrajWrite(fnameOut, ArgList(), DataSetList(), TrajectoryFile::MOL2FILE); + trajout.SetupTrajWrite(newParm, newFrame.CoordsInfo(), 1); + trajout.WriteSingle(0, newFrame); + trajout.EndTraj(); + // END DEBUG } +#endif /** Create box, fill with target number of solvent molecules. */ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptraj::Parm::ParameterSet const& set0, @@ -600,13 +674,13 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra return 1; } // Check if buffers need to be increased for trunc oct. - if (doTruncatedOct_) - octBoxCheck( frameOut ); + //if (doTruncatedOct_) + // octBoxCheck( frameOut, boxX, boxY, boxZ, soluteRadii ); mprintf("\t Solute vdw bounding box: %-5.3f %-5.3f %-5.3f\n", boxX, boxY, boxZ); // Ratio of solute box size to solvent box size double soluteFac = (boxX*boxY*boxZ) / solventBoxVol; - mprintf("\t Solute is %g times the solvent box.\n", soluteFac); + mprintf("\t Solute box is %g times the solvent box.\n", soluteFac); // Estimate how many solvent molecules will be removed by solvent unsigned int nMolsRemovedBySolute = (unsigned int)ceil( soluteFac * (double)SOLVENTBOX.Top().Nmol() ); mprintf("\t Estimating approximately %u solvent molecules will be removed by solute.\n", nMolsRemovedBySolute); @@ -617,7 +691,8 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra while (nSolventAdded < nsolvent_) { nSolventAdded = (long int)(nBoxesInCube(layer) * (unsigned int)SOLVENTBOX.Top().Nmol()); nSolventAdded = nSolventAdded - (long int)nMolsRemovedBySolute; - mprintf("DEBUG: Layer %i, # solvent added = %li\n", layer, nSolventAdded); + if (debug_ > 0) + mprintf("DEBUG: Layer %i, # solvent added = %li\n", layer, nSolventAdded); layer++; } // NOTE: layer is actually layer+1 right now @@ -628,7 +703,8 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra double dXWidth = ((double)layer * solventX); double dYWidth = ((double)layer * solventY); double dZWidth = ((double)layer * solventZ); - mprintf("DEBUG: Estimated size with buffer: %f %f %f\n", dXWidth, dYWidth, dZWidth); + if (debug_ > 0) + mprintf("DEBUG: Estimated size with buffer: %f %f %f\n", dXWidth, dYWidth, dZWidth); // See how many solvent boxes required in each dimension int iX = (int)( dXWidth / solventX ) + 1; @@ -649,14 +725,15 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra addSolventUnits(iX, iY, iZ, soluteMaxR, dXStart, dYStart, dZStart, solventX, solventY, solventZ, solventFrame, SOLVENTBOX.Top(), frameOut, topOut, soluteRadii, solventRadii); +# ifdef CPPTRAJ_DEBUG_SOLVATE // DEBUG -// Trajout_Single trajout; -// trajout.InitTrajWrite("BigCube.mol2", ArgList(), DataSetList(), TrajectoryFile::MOL2FILE); -// trajout.SetupTrajWrite(&topOut, frameOut.CoordsInfo(), 1); -// trajout.WriteSingle(0, frameOut); -// trajout.EndTraj(); + Trajout_Single trajout; + trajout.InitTrajWrite("BigCube.mol2", ArgList(), DataSetList(), TrajectoryFile::MOL2FILE); + trajout.SetupTrajWrite(&topOut, frameOut.CoordsInfo(), 1); + trajout.WriteSingle(0, frameOut); + trajout.EndTraj(); // END DEBUG - +# endif // Define the size of the new solvent/solute system double maxX, maxY, maxZ; soluteRadii = getAtomRadii(soluteMaxR, topOut, set0) ; @@ -664,11 +741,19 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra mprinterr("Error: Setting vdw bounding box for solute/solvent system failed.\n", topOut.c_str()); return 1; } + double scaleFac = 1.0; + if (doTruncatedOct_) + scaleFac = octBoxCheck( frameOut, maxX, maxY, maxZ, soluteRadii ); + maxX *= scaleFac; + maxY *= scaleFac; + maxZ *= scaleFac; + mprintf("\t Initial trial cell size is %g %g %g\n", maxX, maxY, maxZ); if (firstSolventRes >= topOut.Nres()) { mprinterr("Internal Error: Solvate::SolvateBoxWithExactNumber(): No solvent added.\n"); return 1; } + mprintf("\t %i solvent residues in initial cell.\n", topOut.Nres() - firstSolventRes); //unsigned int nBoxesForSolvent = (unsigned int)ceil( (double)nsolvent / (double)SOLVENTBOX.Top().Nmol() ); // // Calculate the solvent box density @@ -692,9 +777,8 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra //unsigned int nSolventBoxesNeeded = (unsigned int)ceil( (double)neededNumberOfSolvent / (double)SOLVENTBOX.Top().Nmol() ); //mprintf("\t Estimating %u solvent boxes will be needed.\n", nSolventBoxesNeeded); - frameOut.CenterOnOrigin(false); +/* frameOut.CenterOnOrigin(false);*/ - // FIXME make user-specifiable double alpha, beta, gamma; if (doTruncatedOct_) { alpha = Box::TruncatedOctAngle(); @@ -705,6 +789,7 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra beta = 90.0; gamma = 90.0; } + // FIXME make user-specifiable int negtol = -5; double bufX = boxX + ((maxX - boxX) / 2.0); @@ -712,9 +797,9 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra double bufZ = boxZ + ((maxZ - boxZ) / 2.0); if (doTruncatedOct_) { - if (bufX >= bufY && bufX >= bufZ) { bufY = bufX; bufZ = bufX; } - if (bufY >= bufX && bufY >= bufZ) { bufX = bufY; bufZ = bufY; } - if (bufZ >= bufX && bufZ >= bufY) { bufX = bufZ; bufY = bufZ; } + double dTemp = std::max(bufX, bufY); + dTemp = std::max(dTemp, bufZ); + bufX = bufY = bufZ = dTemp; } // Mask selecting everything that should be kept. Always keep solute @@ -739,12 +824,56 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra int lastdiff = 0; Box newBox; while (loop) { - mprintf("DEBUG: %i) Buffer %f %f %f\n", ntries, bufX, bufY, bufZ); + // Set clipping + clipX_ = 0.5 * bufX; + clipY_ = 0.5 * bufY; + clipZ_ = 0.5 * bufZ; + //mprintf("cCriteria: %f %f %f\n", clipX_, clipY_, clipZ_); // DEBUG + mprintf("\t %i) Trial cell: %f %f %f\n", ntries, bufX, bufY, bufZ); // Define the unit cell newBox.SetupFromXyzAbg( bufX, bufY, bufZ, alpha, beta, gamma ); - newBox.PrintInfo(); - // Convert to fractional + if (debug_ > 0) + newBox.PrintInfo(); int nSolventInCell = 0; + // ------------------------------------------- + for (int ires = firstSolventRes; ires < topOut.Nres(); ires++) { + Residue const& solventRes = topOut.Res(ires); + bool collision = false; + for (int vat = solventRes.FirstAtom(); vat != solventRes.LastAtom(); vat++) + { + // Check for clipping + const double* VXYZ = frameOut.XYZ(vat); + double dXabs = fabs(VXYZ[0]); + double dYabs = fabs(VXYZ[1]); + double dZabs = fabs(VXYZ[2]); + if ( dXabs >= clipX_ ) { collision = true; break; } + if ( dYabs >= clipY_ ) { collision = true; break; } + if ( dZabs >= clipZ_ ) { collision = true; break; } + // Check if atom falls outside of oct/diagonal clip + if (doTruncatedOct_) { + double dX = 0.5 * (clipX_ - dXabs) / clipX_; + dX = fabs( dX - 0.5 ); + double dY = 0.5 * (clipY_ - dYabs) / clipY_; + dY = fabs( dY - 0.5 ); + double dZ = 0.5 * (clipZ_ - dZabs) / clipZ_; + dZ = fabs( dZ - 0.5 ); + if ( ( dX + dY + dZ ) > 0.75 ) { collision = true; break; } + } + } // END loop over solvent residue atoms + if (collision) { + // Remove solvent + for (int at = solventRes.FirstAtom(); at != solventRes.LastAtom(); at++) + cmask.SelectAtom(at, false); + } else { + // Keep solvent + for (int at = solventRes.FirstAtom(); at != solventRes.LastAtom(); at++) + cmask.SelectAtom(at, true); + nSolventInCell++; + } + } // END loop over solvent residues + // ------------------------------------------- +/* + // Convert to fractional Matrix_3x3 const& recip = newBox.FracCell(); for (int ires = firstSolventRes; ires < topOut.Nres(); ires++) { bool inCell = true; @@ -773,12 +902,14 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra cmask.SelectAtom(at, false); } } // END loop over solvent residues +*/ + // ------------------------------------------- //cmask.MaskInfo(); // DEBUG //mprintf("DEBUG: Should select %i atoms.\n", firstSolventAtom + (nSolventInCell*topOut.Res(firstSolventRes).NumAtoms())); // How far off is it? int diff = nsolvent_ - nSolventInCell; - mprintf("DEBUG:\t%i solvent residues in cell. Diff: %i\n", nSolventInCell, diff); + mprintf("\t\t%i solvent residues in trial cell. Diff: %i\n", nSolventInCell, diff); if (diff < 0) { if (smallestNegDiff == 0 || diff > smallestNegDiff) { @@ -790,8 +921,9 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra // Update smallest buffer; presumably is closer smallestBuf = Vec3(bufX, bufY, bufZ); } - } - mprintf("DEBUG: SmallestNegDiff= %i nTimesSmallestSeen= %i\n", smallestNegDiff, nTimesSmallestSeen); + } + if (debug_ > 0) + mprintf("DEBUG: SmallestNegDiff= %i nTimesSmallestSeen= %i\n", smallestNegDiff, nTimesSmallestSeen); // If this is the first time through choose an appropriate change val if (ntries == 0) { change = 0.001; @@ -800,11 +932,11 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra if ((diff < 0 && nTimesSmallestSeen > 2) || (diff < 0 && diff >= negtol)) { // Close enough, just remove the last water residue(s) int nToRemove = -diff; - mprintf("\tOnly %i solvent molecules off, removing them.\n", nToRemove); + mprintf("\t Only %i solvent molecules off, removing them.\n", nToRemove); for (int ires = topOut.Nres() - 1; ires > firstSolventRes-1; ires--) { Residue const& currentRes = topOut.Res(ires); if (cmask.AtomInCharMask( currentRes.FirstAtom() )) { - mprintf("Removing %s\n", topOut.TruncResNameNum(ires).c_str()); + mprintf("\t\tRemoving %s\n", topOut.TruncResNameNum(ires).c_str()); for (int at = currentRes.FirstAtom(); at != currentRes.LastAtom(); at++) cmask.SelectAtom(at, false); nToRemove--; @@ -824,7 +956,7 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra loop = false; } else if (diff == 0) { loop = false; - mprintf("DEBUG: Found.\n"); + if (debug_ > 0) mprintf("DEBUG: Found.\n"); } else { double CHANGE_DIR = 0; if ( lastdiff > 0 && diff < 0) { @@ -853,26 +985,39 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra lastdiff = diff; // Choose a new buffer value double CD = change * (double)diff; // TODO should be minus? + //mprintf("DEBUG: CD= %16.8E\n", CD); double newbufX = bufX + CD; // TODO should be minus double newbufY = bufY + CD; // TODO should be minus double newbufZ = bufZ + CD; // TODO should be minus + //mprintf("DEBUG: New proposed buffer sizes: %f %f %f\n", newbufX, newbufY, newbufZ); //CheckBuffer $newbuffer; //lastbuffer=$buffer - checkBuffer(newbufX, boxX, maxX); - checkBuffer(newbufY, boxY, maxY); - checkBuffer(newbufZ, boxZ, maxZ); + checkBuffer(newbufX, boxX, maxX, bufX, change, diff); + if (doTruncatedOct_) { + newbufY = newbufZ = newbufX; + } else { + checkBuffer(newbufY, boxY, maxY, bufY, change, diff); + checkBuffer(newbufZ, boxZ, maxZ, bufZ, change, diff); + } + bufX = newbufX; bufY = newbufY; bufZ = newbufZ; - mprintf("DEBUG: Change: %G (%g) newBuf= %f %f %f\n", change, CHANGE_DIR, bufX, bufY, bufZ); + if (debug_ > 0) + mprintf("DEBUG: Change: %G (%g) newBuf= %f %f %f\n", change, CHANGE_DIR, bufX, bufY, bufZ); } // END change calc - +# ifdef CPPTRAJ_DEBUG_SOLVATE + // DEBUG + printIntermediateStructure(ntries, cmask, topOut, frameOut); + // END DEBUG +# endif // Safety valve ntries++; if (ntries > 100) { mprintf("\tTaking too long - moving on...\n"); loop = false; - mprintf("DEBUG: SmallestNegDiff %i (%i) : %f %f %f\n", smallestNegDiff, nTimesSmallestSeen, smallestBuf[0], smallestBuf[1], smallestBuf[2]); + if (debug_ > 0) + mprintf("DEBUG: SmallestNegDiff %i (%i) : %f %f %f\n", smallestNegDiff, nTimesSmallestSeen, smallestBuf[0], smallestBuf[1], smallestBuf[2]); } } // END cell loop @@ -887,6 +1032,15 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra topOut = *newParm; delete newParm; topOut.Brief("Topology with target # of solvent:"); + // Reset the solvent residue numbers + int solventResNum = topOut.Res(firstSolventRes - 1).OriginalResNum() + 1; + for (int ires = firstSolventRes; ires < topOut.Nres(); ires++) + topOut.SetRes(ires).SetOriginalNum( solventResNum++ ); + // Reset the solvent atom numbers TODO zero these? + if (!topOut.PdbSerialNum().empty()) { + for (int at = topOut.Res(firstSolventRes).FirstAtom(); at < topOut.Natom(); at++) + topOut.ModifyPdbSerialNum()[at] = at+1; + } Frame newFrame; newFrame.SetupFrameV(topOut.Atoms(), frameOut.CoordsInfo()); @@ -904,7 +1058,7 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra ewald_rotate(frameOut, dAngle); //mprintf("EwaldRotate: %f\n", dAngle); // Add an angstrom to the desired box size rather than using the bounding box size - dAngle = clipX_ + .5; // NOTE: If truncoct is enabled for nsolvent, clipX_ will be 0 + dAngle = clipX_ + .5; boxX = boxY = boxZ = dAngle * sqrt(3.0) * 0.5; } @@ -924,6 +1078,7 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra mprintf("\t Volume: %5.3lf A^3\n", frameOut.BoxCrd().CellVolume()); // Reimage FIXME create separate init routine for autoimage +/* Action_AutoImage autoimage; ArgList actionargs; DataSetList tempdsl; @@ -945,8 +1100,21 @@ int Solvate::SolvateBoxWithExactNumber(Topology& topOut, Frame& frameOut, Cpptra mprinterr("Error: Could not autoimage post-solvate.\n"); return 1; } - +*/ calculateDensity(topOut, frameOut, set0); + // Center if needed + if (center_) { + double dX2 = boxX * 0.5; + double dY2 = boxY * 0.5; + double dZ2 = boxZ * 0.5; + double* xptr = frameOut.xAddress(); + for (int at = 0; at < frameOut.Natom(); at++, xptr += 3) + { + xptr[0] += dX2; + xptr[1] += dY2; + xptr[2] += dZ2; + } + } return 0; } @@ -990,7 +1158,7 @@ int Solvate::SolvateBox(Topology& topOut, Frame& frameOut, Cpptraj::Parm::Parame } // Check if buffers need to be increased for trunc oct. if (doTruncatedOct_) - octBoxCheck( frameOut ); + octBoxCheck( frameOut, boxX, boxY, boxZ, soluteRadii ); mprintf("\t Solute vdw bounding box: %-5.3f %-5.3f %-5.3f\n", boxX, boxY, boxZ); double dXWidth = boxX + bufferX_ * 2; diff --git a/src/Structure/Solvate.h b/src/Structure/Solvate.h index a3a2c5df3f..e21effc70c 100644 --- a/src/Structure/Solvate.h +++ b/src/Structure/Solvate.h @@ -66,7 +66,7 @@ class Solvate { Frame const&, Topology const&, Frame const&, std::vector const&, std::vector const&) const; /// Scale up buffer sizes for oct box if needed - void octBoxCheck(Frame const&); + double octBoxCheck(Frame const&,double,double,double,std::vector const&); /// Ewald rotate for trun. oct system static void ewald_rotate(Frame&, double&); /// Add solvent unit boxes diff --git a/src/Structure/structuredepend b/src/Structure/structuredepend index 219597caff..df56e6dee6 100644 --- a/src/Structure/structuredepend +++ b/src/Structure/structuredepend @@ -14,7 +14,7 @@ Model.o : Model.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../Coordinat OldBuilder.o : OldBuilder.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h BuildAtom.h Chirality.h InternalCoords.h OldBuilder.h StructureEnum.h Zmatrix.h PdbCleaner.o : PdbCleaner.cpp ../ArgList.h ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h PdbCleaner.h RingFinder.o : RingFinder.cpp ../Atom.h ../AtomMask.h ../Box.h ../CharMask.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h RingFinder.h -Solvate.o : Solvate.cpp ../Action.h ../ActionState.h ../Action_AutoImage.h ../ArgList.h ../AssociatedData.h ../Atom.h ../AtomMask.h ../BaseIOtype.h ../Box.h ../CharMask.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DataFile.h ../DataFileList.h ../DataSet.h ../DataSetList.h ../DataSet_Coords.h ../DataSet_Coords_REF.h ../Dimension.h ../DispatchObject.h ../FileIO.h ../FileName.h ../FileTypes.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../MetaData.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Parm/../AtomType.h ../Parm/../CmapParmHolder.h ../Parm/../Constants.h ../Parm/../FileName.h ../Parm/../NameType.h ../Parm/../ParameterTypes.h ../Parm/../Parm/ParmEnum.h ../Parm/../TypeNameHolder.h ../Parm/DihedralParmHolder.h ../Parm/ImproperParmHolder.h ../Parm/ParameterSet.h ../Parm/ParmEnum.h ../Parm/ParmHolder.h ../Range.h ../ReferenceFrame.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../TextFormat.h ../Timer.h ../Topology.h ../Unit.h ../Vec3.h Solvate.h +Solvate.o : Solvate.cpp ../Action.h ../ActionFrameCounter.h ../ActionState.h ../Action_AutoImage.h ../ArgList.h ../AssociatedData.h ../Atom.h ../AtomMask.h ../BaseIOtype.h ../Box.h ../CharMask.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DataFile.h ../DataFileList.h ../DataSet.h ../DataSetList.h ../DataSet_Coords.h ../DataSet_Coords_REF.h ../Dimension.h ../DispatchObject.h ../FileIO.h ../FileName.h ../FileTypes.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../MetaData.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../OutputTrajCommon.h ../Parallel.h ../ParameterTypes.h ../Parm/../AtomType.h ../Parm/../CmapParmHolder.h ../Parm/../Constants.h ../Parm/../FileName.h ../Parm/../NameType.h ../Parm/../ParameterTypes.h ../Parm/../Parm/ParmEnum.h ../Parm/../TypeNameHolder.h ../Parm/DihedralParmHolder.h ../Parm/ImproperParmHolder.h ../Parm/ParameterSet.h ../Parm/ParmEnum.h ../Parm/ParmHolder.h ../Range.h ../ReferenceFrame.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../StringRoutines.h ../SymbolExporting.h ../TextFormat.h ../Timer.h ../Topology.h ../TrajectoryFile.h ../Trajout_Single.h ../Unit.h ../Vec3.h Solvate.h StructureEnum.o : StructureEnum.cpp StructureEnum.h StructureRoutines.o : StructureRoutines.cpp ../Atom.h ../CpptrajStdio.h ../NameType.h ../Residue.h ../SymbolExporting.h StructureRoutines.h Sugar.o : Sugar.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h Sugar.h SugarToken.h diff --git a/src/Version.h b/src/Version.h index d138368205..ddbee0927b 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 "V7.4.0" +#define CPPTRAJ_INTERNAL_VERSION "V7.5.0" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index b540b9c32e..bf47da9c56 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -505,7 +505,7 @@ Structure/Model.o : Structure/Model.cpp Atom.h AtomMask.h Box.h Constants.h Coor Structure/OldBuilder.o : Structure/OldBuilder.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/BuildAtom.h Structure/Chirality.h Structure/InternalCoords.h Structure/OldBuilder.h Structure/StructureEnum.h Structure/Zmatrix.h SymbolExporting.h Topology.h Unit.h Vec3.h Structure/PdbCleaner.o : Structure/PdbCleaner.cpp ArgList.h Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/PdbCleaner.h SymbolExporting.h Topology.h Unit.h Vec3.h Structure/RingFinder.o : Structure/RingFinder.cpp Atom.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/RingFinder.h SymbolExporting.h Topology.h Unit.h Vec3.h -Structure/Solvate.o : Structure/Solvate.cpp Action.h ActionState.h Action_AutoImage.h ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h Structure/../Parm/../AtomType.h Structure/../Parm/../CmapParmHolder.h Structure/../Parm/../Constants.h Structure/../Parm/../FileName.h Structure/../Parm/../NameType.h Structure/../Parm/../ParameterTypes.h Structure/../Parm/../Parm/ParmEnum.h Structure/../Parm/../TypeNameHolder.h Structure/../Parm/DihedralParmHolder.h Structure/../Parm/ImproperParmHolder.h Structure/../Parm/ParameterSet.h Structure/../Parm/ParmEnum.h Structure/../Parm/ParmHolder.h Structure/Solvate.h SymbolExporting.h TextFormat.h Timer.h Topology.h Unit.h Vec3.h +Structure/Solvate.o : Structure/Solvate.cpp Action.h ActionFrameCounter.h ActionState.h Action_AutoImage.h ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h Structure/../Parm/../AtomType.h Structure/../Parm/../CmapParmHolder.h Structure/../Parm/../Constants.h Structure/../Parm/../FileName.h Structure/../Parm/../NameType.h Structure/../Parm/../ParameterTypes.h Structure/../Parm/../Parm/ParmEnum.h Structure/../Parm/../TypeNameHolder.h Structure/../Parm/DihedralParmHolder.h Structure/../Parm/ImproperParmHolder.h Structure/../Parm/ParameterSet.h Structure/../Parm/ParmEnum.h Structure/../Parm/ParmHolder.h Structure/Solvate.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajectoryFile.h Trajout_Single.h Unit.h Vec3.h Structure/StructureEnum.o : Structure/StructureEnum.cpp Structure/StructureEnum.h Structure/StructureRoutines.o : Structure/StructureRoutines.cpp Atom.h CpptrajStdio.h NameType.h Residue.h Structure/StructureRoutines.h SymbolExporting.h Structure/Sugar.o : Structure/Sugar.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/Sugar.h Structure/SugarToken.h SymbolExporting.h Topology.h Unit.h Vec3.h diff --git a/test/Test_MinImage/RunTest.sh b/test/Test_MinImage/RunTest.sh index 265713be4b..5a40cfc479 100755 --- a/test/Test_MinImage/RunTest.sh +++ b/test/Test_MinImage/RunTest.sh @@ -2,7 +2,7 @@ . ../MasterTest.sh -CleanFiles min.in min.dat +CleanFiles min.in min.dat toface.dat test.toface.dat test.pdb TESTNAME='Minimum non-self imaged distance test' Requires netcdf maxthreads 10 @@ -14,8 +14,46 @@ parm ../tz2.truncoct.parm7 trajin ../tz2.truncoct.nc minimage m1 :1-13 :1-13 out min.dat minimage m2 :1-13 :1-13 maskcenter out min.dat +minimage m3 :1-13 toface out toface.dat EOF RunCpptraj "Minimum Image test." DoTest min.dat.save min.dat +DoTest toface.dat.save toface.dat + +cat > test.pdb < min.in <