Major changes to unit cell handling. Fix bugs in GIST, LIE, PairDist actions. Miscellaneous other fixes.

LICENSE

@@ -1,6 +1,6 @@
   _______________________________________________________________________
 
-                 CPPTRAJ: 2019 - LICENSE AGREEMENT
+                 CPPTRAJ: 2021 - LICENSE AGREEMENT
   _______________________________________________________________________
 
   cpptraj is free software; you can redistribute it and/or modify it

README.md

@@ -46,7 +46,7 @@ CPPTRAJ see the following publication:
 Disclaimer and Copyright
 ========================
 
-CPPTRAJ is Copyright (c) 2010-2020 Daniel R. Roe.
+CPPTRAJ is Copyright (c) 2010-2021 Daniel R. Roe.
 The terms for using, copying, modifying, and distributing CPPTRAJ are
 specified in the file LICENSE.
 

doc/CpptrajDevelopmentGuide.lyx

@@ -1912,6 +1912,73 @@ selectedTgt.SetCoordinates( Tgt2, mask );
 rmsd = selectedTgt.RMSD_CenteredRef( selectedRef, ...
 \end_layout
 
+\begin_layout Subsection
+Box
+\end_layout
+
+\begin_layout Standard
+The Box class holds all unit cell information.
+ It has the unit cell vectors (stored in a Matrix_3x3 class in row-major
+ order), as well as the fractional cell matrix (for converting from Cartesian
+ to fractional coordinates), the box lengths and angles, and the unit cell
+ volume.
+ There are in general two ways to set up Box with unit cell information:
+ the SetupX routines (e.g.
+ SetupFromXyzAbg()) and AssignX routines (e.g.
+ AssignFromXyzAbg()).
+ The SetupX routines perform some extra checks and will return 1 if something
+ is wrong with the box.
+ The AssignX routines do not perform these checks and so should be used
+ where performance is needed (or the state of the box doesn't matter).
+\end_layout
+
+\begin_layout Standard
+Both the SetupX and AssignX will set the unit and fractional cell matrices,
+ the box lengths and angles, and the cell volume each time they are called.
+ The ImageOption class can then be used to determine what kind of imaging
+ can be performed on the unit cell using the distance calculation routines
+ (DIST2(), DIST()) found in DistRoutines.h.
+ For example, this is how distances are calculated using ImageOption imageOpt_
+ in the distance action.
+\end_layout
+
+\begin_layout Standard
+1.
+ In Init(), determine if imaging should be used if possible:
+\end_layout
+
+\begin_layout LyX-Code
+imageOpt_.InitImaging( !(actionArgs.hasKey("noimage")) );
+\end_layout
+
+\begin_layout Standard
+2.
+ In Setup(), determine if imaging will be possible based on if box info
+ is present:
+\end_layout
+
+\begin_layout LyX-Code
+imageOpt_.SetupImaging( setup.CoordInfo().TrajBox().HasBox() );
+\end_layout
+
+\begin_layout Standard
+3.
+ In DoAction(), determine what imaging to use and then calculate distance
+ using DIST() from DistRoutines.h:
+\end_layout
+
+\begin_layout LyX-Code
+if (imageOpt_.ImagingEnabled())
+\end_layout
+
+\begin_layout LyX-Code
+  imageOpt_.SetImageType( frm.Frm().BoxCrd().Is_X_Aligned_Ortho() );
+\end_layout
+
+\begin_layout LyX-Code
+double Dist = DIST(imageOpt_.ImagingType(), a1, a2_, frm.Frm().BoxCrd());
+\end_layout
+
 \begin_layout Section
 Console and File Input/Output
 \end_layout

doc/cpptraj.lyx

@@ -15595,6 +15595,100 @@ trajout output.crd onlyframes 2,5-7
 only frames 2, 5, 6, and 7 from input.crd will be written to output.crd.
 \end_layout
 
+\begin_layout Subsubsection*
+Cell not X-aligned Warning
+\end_layout
+
+\begin_layout Standard
+Certain Actions (e.g.
+
+\series bold
+\emph on
+align
+\series default
+\emph default
+, 
+\series bold
+\emph on
+rms
+\series default
+\emph default
+, 
+\series bold
+\emph on
+principal
+\series default
+\emph default
+, etc.) can rotate the unit cell vectors (i.e.
+ the box) if they are present.
+ Some trajectory formats do not support writing out box coordinates if the
+ unit cell is not 
+\begin_inset Quotes eld
+\end_inset
+
+X-aligned
+\begin_inset Quotes erd
+\end_inset
+
+; in other words, if the unit cell 
+\begin_inset Quotes eld
+\end_inset
+
+A
+\begin_inset Quotes erd
+\end_inset
+
+ vector is not aligned with the coordinate X-axis and the 
+\begin_inset Quotes eld
+\end_inset
+
+B
+\begin_inset Quotes erd
+\end_inset
+
+ vector is not in the X-Y plane.
+ If this is the case, the following warnings may appear:
+\end_layout
+
+\begin_layout LyX-Code
+Warning: Unit cell is not X-aligned.
+ Box cannot be properly stored as <format>.
+\end_layout
+
+\begin_layout LyX-Code
+Warning: Set <#>; unit cell is not X-aligned.
+ Box cannot be properly stored as <format>.
+\end_layout
+
+\begin_layout Standard
+This means that the frame will be written with the X-aligned unit cell instead
+ of the actual unit cell.
+ Imaging will not be possible with a trajectory written this way.
+ Currently the only trajectory formats that support writing non-X-aligned
+ cells are the Gromacs TRR and XTC formats.
+\end_layout
+
+\begin_layout Standard
+If unit cell information is no longer needed, it can be removed (via e.g.
+ the 
+\series bold
+\emph on
+box
+\series default
+\emph default
+ action, the 
+\series bold
+\emph on
+strip
+\series default
+\emph default
+ action with the 
+\series bold
+'nobox'
+\series default
+ keyword, etc.) to prevent these warnings from triggering.
+\end_layout
+
 \begin_layout Subsubsection
 
 \shape italic
@@ -16020,6 +16114,13 @@ veltraj Write velocity trajectory instead of coordinates.
 \end_layout
 
 \end_deeper
+\begin_layout Standard
+Note that by default CPPTRAJ will try to write the symmetric shape matrix
+ if box information is present.
+ If this is not possible, CPPTRAJ will fall back to writing unit cell parameters
+ (lengths and angles) as long as the cell is X-aligned.
+\end_layout
+
 \begin_layout Subsubsection
 
 \emph on
@@ -16037,6 +16138,11 @@ dt Time step tp multiply set numbers by (default 1.0).
 \end_layout
 
 \end_deeper
+\begin_layout Standard
+Note: these formats can write rotated (i.e.
+ non-X-aligned) unit cells.
+\end_layout
+
 \begin_layout Subsubsection
 
 \shape italic

src/Action.h

@@ -60,8 +60,5 @@ class Action : public DispatchObject {
     /// Process array of frames before start required by ranks > 0 (not called by master)
     virtual int ParallelPreloadFrames(FArray const&) { return 0; } // TODO: pure virtual
 #   endif
-  protected:
-    /// Print a warning that imaging will not be possible if coordinates are being rotated.
-    static void CheckImageRotationWarning(ActionSetup const&, const char*);
 };
 #endif

src/Action_Align.cpp

@@ -84,9 +84,6 @@ Action::RetType Action_Align::Setup(ActionSetup& setup) {
   if (REF_.SetupRef(setup.Top(), tgtMask_.Nselected()))
     return Action::SKIP;
 
-  // Warn if PBC and rotating
-  Action::CheckImageRotationWarning(setup, "the alignment");
-
   return Action::OK;
 }
 
@@ -98,6 +95,7 @@ Action::RetType Action_Align::DoAction(int frameNum, ActionFrame& frm) {
   tgtFrame_.SetCoordinates(frm.Frm(), tgtMask_);
   tgtFrame_.RMSD_CenteredRef(REF_.SelectedRef(), rot_, tgtTrans_, useMass_);
   frm.ModifyFrm().Trans_Rot_Trans(movMask_, tgtTrans_, rot_, REF_.RefTrans());
+  frm.ModifyFrm().ModifyBox().RotateUcell( rot_ );
   REF_.PreviousRef( frm.Frm() );
   return Action::MODIFY_COORDS;
 }

src/Action_AreaPerMol.cpp

@@ -58,14 +58,14 @@ Action::RetType Action_AreaPerMol::Init(ArgList& actionArgs, ActionInit& init, i
   */
 Action::RetType Action_AreaPerMol::Setup(ActionSetup& setup) {
   // Needs box info
-  if (setup.CoordInfo().TrajBox().Type() == Box::NOBOX) {
+  if (!setup.CoordInfo().TrajBox().HasBox()) {
     mprintf("Warning: No box information for '%s', cannot calculate area.\n",
             setup.Top().c_str());
     return Action::SKIP;
   }
   // Probably will not work for non-orthorhombic cells
-  if (setup.CoordInfo().TrajBox().Type() != Box::ORTHO)
-    mprintf("Warning: Box is not orthorhombic, calculated area may not be correct.\n");
+  if (!setup.CoordInfo().TrajBox().Is_X_Aligned_Ortho())
+    mprintf("Warning: Box is not X-aligned orthorhombic, calculated area may not be correct.\n");
   // Determine how many molecules are selected
   if (Mask1_.MaskStringSet()) {
     if (setup.Top().SetupCharMask(Mask1_)) return Action::ERR;
@@ -93,11 +93,11 @@ Action::RetType Action_AreaPerMol::Setup(ActionSetup& setup) {
 Action::RetType Action_AreaPerMol::DoAction(int frameNum, ActionFrame& frm) {
   double area;
   if (areaType_ == XY)
-    area = frm.Frm().BoxCrd().BoxX() * frm.Frm().BoxCrd().BoxY();
+    area = frm.Frm().BoxCrd().Param(Box::X) * frm.Frm().BoxCrd().Param(Box::Y);
   else if (areaType_ == XZ) 
-    area = frm.Frm().BoxCrd().BoxX() * frm.Frm().BoxCrd().BoxZ();
+    area = frm.Frm().BoxCrd().Param(Box::X) * frm.Frm().BoxCrd().Param(Box::Z);
   else // if areaType_ == YZ
-    area = frm.Frm().BoxCrd().BoxY() * frm.Frm().BoxCrd().BoxZ();
+    area = frm.Frm().BoxCrd().Param(Box::Y) * frm.Frm().BoxCrd().Param(Box::Z);
 
   area = area / Nmols_;
 

src/Action_AtomMap.cpp

@@ -232,7 +232,6 @@ Action::RetType Action_AtomMap::Setup(ActionSetup& setup) {
   }
   if (rmsfit_) {
     mprintf("\trmsfit specified, %i atoms.\n",rmsRefFrame_.Natom());
-    Action::CheckImageRotationWarning(setup, "the RMS-fit");
     return Action::OK;
   }
   mprintf("\tMap for parm %s -> %s (%i atom).\n",TgtFrame_->Top().c_str(),
@@ -257,6 +256,7 @@ Action::RetType Action_AtomMap::DoAction(int frameNum, ActionFrame& frm) {
     Vec3 Trans, refTrans;
     double R = rmsTgtFrame_.RMSD(rmsRefFrame_, Rot, Trans, refTrans, false);
     frm.ModifyFrm().Trans_Rot_Trans(Trans, Rot, refTrans);
+    frm.ModifyFrm().ModifyBox().RotateUcell( Rot );
     if (rmsdata_!=0)
       rmsdata_->Add(frameNum, &R);
   } else {