Fix solvating with truncated octahedral cells (Amber-MD#1174)

* Truncated octahedron diagonal cut scaling fix suggested by Nikolai Skrynnikov and Stephan Schott.

* Fix minimage action setup when no data set name specified

* Try to implement toface keyword to check closest approach to unit cell face

* Add a toface test

* Add another toface test

* Store face a string. Ensure negative frac coords are wrapped back into 0.0 to 1.0 range.

* Fix up solvateoct diagonal check to be more like the original fixed code

* Original test save calculated when LEAPPI active, update for regular PI

* Add debug print of intermediate structures

* Try to handle the trunc oct box better in nsolvent

* Hide some debug info. Ensure all box lengths are properly enforced to be equal for truncoct

* Do clipping just like the original solvate command - results in nicer looking trunc oct cells

* Add the centering step

* Hide some debug info. Clean up output

* Remove unused code

* Add ifdef, update depends

* nsolvent now works with truncoct

* Remove old include

* Improve build manual entry

* V7.5.0. Minor version bump for fixed truncated octahedral cell calculation (thanks to N. Skrynnikov for the fix). Also enables truncated octahedron for 'nsolvent' keyword.

* Update manual PDF

---------

Co-authored-by: Daniel R. Roe <daniel.roe@nih.gov>

Author: drroe

doc/CpptrajManual.pdf

@@ -1,3 +1,3 @@
 f6f8cb1a79951d80a9d2656fd9c30f55  CpptrajDevelopmentGuide.lyx
-3b89e9556662b1c6de2d0e38ab894588  cpptraj.lyx
+369e9afee9918ff0797c4a4bf08473fd  cpptraj.lyx
 5d9b5b5ed47a3ded57b6464df99b3585  CpptrajManual.lyx

doc/DocumentChecksums.txt

@@ -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

doc/cpptraj.lyx

@@ -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_) {

src/Action_MinImage.cpp

@@ -14,13 +14,16 @@ 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.
     DataSet* atom1_;
     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<double> minDist_;