Merge pull request #658 from drroe/vectorvel

Add 'velocity' and 'force' keywords to 'vector' command.

Author: drroe

src/Action_Vector.cpp

@@ -20,20 +20,22 @@ Action_Vector::Action_Vector() :
 void Action_Vector::Help() const {
   mprintf("\t[<name>] <Type> [out <filename> [ptrajoutput]] [<mask1>] [<mask2>]\n"
           "\t[magnitude] [ired]\n"
-          "\t<Type> = { mask     | minimage  | dipole | center | corrplane | \n"
-          "\t           box      | boxcenter | ucellx | ucelly | ucellz    | \n"
-          "\t           momentum | principal [x|y|z] }\n" 
+          "\t<Type> = { mask     | minimage  | dipole | center   | corrplane | \n"
+          "\t           box      | boxcenter | ucellx | ucelly   | ucellz    | \n"
+          "\t           momentum | principal [x|y|z]  | velocity | force       }\n" 
           "  Calculate the specified coordinate vector.\n"
           "    mask: (Default) Vector from <mask1> to <mask2>.\n"
-          "    minimage: Store the minimum image vector between atoms in <mask1> and <mask2>.\n"
-          "    dipole: Dipole and center of mass of the atoms specified in <mask1>\n"
-          "    center: Store the center of mass of atoms in <mask1>.\n"
-          "    corrplane: Vector perpendicular to plane through the atoms in <mask1>.\n"
-          "    box: (No mask needed) Store the box lengths of the trajectory.\n"
-          "    boxcenter: (No mask needed) Store box center as vector.\n"
-          "    ucell{x|y|z}: (No mask needed) Store specified unit cell vector.\n"
-          "    momentum : Store total momentum vector of atoms in <mask1> (requires velocities).\n"
-          "    principal [x|y|z]: X, Y, or Z principal axis vector for atoms in <mask1>.\n");
+          "    minimage         : Store the minimum image vector between atoms in <mask1> and <mask2>.\n"
+          "    dipole           : Dipole and center of mass of the atoms specified in <mask1>\n"
+          "    center           : Store the center of mass of atoms in <mask1>.\n"
+          "    corrplane        : Vector perpendicular to plane through the atoms in <mask1>.\n"
+          "    box              : (No mask needed) Store the box lengths of the trajectory.\n"
+          "    boxcenter        : (No mask needed) Store box center as vector.\n"
+          "    ucell{x|y|z}     : (No mask needed) Store specified unit cell vector.\n"
+          "    momentum         : Store total momentum vector of atoms in <mask1> (requires velocities).\n"
+          "    principal [x|y|z]: X, Y, or Z principal axis vector for atoms in <mask1>.\n"
+          "    velocity         : Store velocity of atoms in <mask1> (requires velocities).\n"
+          "    force            : Store force of atoms in <mask1> (requires forces).\n");
 }
 
 // DESTRUCTOR
@@ -45,7 +47,7 @@ const char* Action_Vector::ModeString[] = {
   "NO_OP", "Principal X", "Principal Y", "Principal Z",
   "Dipole", "Box", "Mask", "Ired",
   "CorrPlane", "Center", "Unit cell X", "Unit cell Y", "Unit cell Z",
-  "Box Center", "MinImage", "Momentum"
+  "Box Center", "MinImage", "Momentum", "Velocity", "Force"
 };
 
 static Action::RetType WarnDeprecated() {
@@ -96,6 +98,10 @@ Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int d
     mode_ = CENTER;
   else if (actionArgs.hasKey("momentum"))
     mode_ = MOMENTUM;
+  else if (actionArgs.hasKey("velocity"))
+    mode_ = VELOCITY;
+  else if (actionArgs.hasKey("force"))
+    mode_ = FORCE;
   else if (actionArgs.hasKey("dipole"))
     mode_ = DIPOLE;
   else if (actionArgs.hasKey("box"))
@@ -184,6 +190,15 @@ Action::RetType Action_Vector::Setup(ActionSetup& setup) {
       return Action::ERR;
     }
   }
+  // Check for velocity/force
+  if ((mode_ == MOMENTUM || mode_ == VELOCITY) && !setup.CoordInfo().HasVel()) {
+    mprintf("Warning: vector %s requires velocity information. Skipping.\n", ModeString[mode_]);
+    return Action::SKIP;
+  }
+  if (mode_ == FORCE && !setup.CoordInfo().HasForce()) {
+    mprintf("Warning: vector %s requires force information. Skipping.\n", ModeString[mode_]);
+    return Action::SKIP;
+  }
   if (mask_.MaskStringSet()) {
     // Setup mask 1
     if (setup.Top().SetupIntegerMask(mask_)) return Action::ERR;
@@ -407,13 +422,29 @@ void Action_Vector::MinImage(Frame const& frm) {
   Vec3 com1 = frm.VCenterOfMass(mask_);
   Vec_->AddVxyz( MinImagedVec(com1, frm.VCenterOfMass(mask2_), ucell, recip), com1 );
 }
-                 
+
+/// \return The center of selected elements in given array.
+static inline Vec3 CalcCenter(const double* xyz, AtomMask const& maskIn) {
+  Vec3 Coord(0.0);
+  for (AtomMask::const_iterator at = maskIn.begin(); at != maskIn.end(); ++at)
+  {
+    int idx = *at * 3;
+    Coord[0] += xyz[idx  ];
+    Coord[1] += xyz[idx+1];
+    Coord[2] += xyz[idx+2];
+  }
+  Coord /= (double)maskIn.Nselected();
+  return Coord;
+}
+
 // Action_Vector::DoAction()
 Action::RetType Action_Vector::DoAction(int frameNum, ActionFrame& frm) {
   switch ( mode_ ) {
     case MASK        : Mask(frm.Frm()); break;
     case CENTER      : Vec_->AddVxyz( frm.Frm().VCenterOfMass(mask_) ); break;
-    case MOMENTUM    : Vec_->AddVxyz( frm.Frm().VMomentum(mask_) ); break; 
+    case MOMENTUM    : Vec_->AddVxyz( frm.Frm().VMomentum(mask_) ); break;
+    case VELOCITY    : Vec_->AddVxyz( CalcCenter(frm.Frm().vAddress(), mask_) ); break;
+    case FORCE       : Vec_->AddVxyz( CalcCenter(frm.Frm().fAddress(), mask_) ); break; 
     case DIPOLE      : Dipole(frm.Frm()); break;
     case PRINCIPAL_X :
     case PRINCIPAL_Y :

src/Action_Vector.h

@@ -13,7 +13,7 @@ class Action_Vector : public Action {
       NO_OP=0,   PRINCIPAL_X, PRINCIPAL_Y, PRINCIPAL_Z,
       DIPOLE,    BOX,         MASK,        IRED,
       CORRPLANE, CENTER,      BOX_X,       BOX_Y,       BOX_Z,
-      BOX_CTR,   MINIMAGE,    MOMENTUM
+      BOX_CTR,   MINIMAGE,    MOMENTUM,    VELOCITY,    FORCE
     };
     static const char* ModeString[];
 

src/Version.h

@@ -20,5 +20,5 @@
  * Whenever a number that precedes <revision> is incremented, all subsequent
  * numbers should be reset to 0.
  */
-#define CPPTRAJ_INTERNAL_VERSION "V4.10.8"
+#define CPPTRAJ_INTERNAL_VERSION "V4.10.9"
 #endif

test/Test_Vector/RunTest.sh

@@ -2,7 +2,8 @@
 
 . ../MasterTest.sh
 
-CleanFiles vector.in vtest.dat.? v8.mol2 corr.v0.v8.dat avgcoord.out res5.out
+CleanFiles vector.in vtest.dat.? vtest.dat.?? v8.mol2 corr.v0.v8.dat \
+           avgcoord.out res5.out
 
 INPUT="-i vector.in"
 # Test Vector mask, principle xyz, dipole, box 
@@ -53,18 +54,34 @@ EOF
 fi
 
 # Test momentum vector.
-UNITNAME="Momentum vector test."
+UNITNAME="Momentum/velocity vector test."
 CheckFor maxthreads 1
 if [ $? -eq 0 ] ; then
   cat > vector.in <<EOF
 parm ../tz2.parm7
 trajin ../Test_SetVelocity/tz2.vel.rst7.save
 vector v9 momentum out vtest.dat.9
+vector v10 velocity @1-3 out vtest.dat.10
 EOF
   RunCpptraj "$UNITNAME"
   DoTest vtest.dat.9.save vtest.dat.9
+  DoTest vtest.dat.10.save vtest.dat.10
 fi
 
+# Test force vector.
+UNITNAME="Force vector test."
+CheckFor maxthreads 1
+if [ $? -eq 0 ] ; then
+  cat > vector.in <<EOF
+parm ../Test_systemVF/systemVF.parm7
+trajin ../Test_systemVF/systemVF.nc 1 1
+vector v11 force @1-3 out vtest.dat.11
+EOF
+  RunCpptraj "$UNITNAME"
+  DoTest vtest.dat.11.save vtest.dat.11
+fi
+
+
 EndTest
 
 exit 0

test/Test_Vector/vtest.dat.10.save

@@ -0,0 +1,2 @@
+#Frame                                                     v10
+       1   0.8980  -0.1969   0.1939   0.0000   0.0000   0.0000

test/Test_Vector/vtest.dat.11.save

@@ -0,0 +1,2 @@
+#Frame                                                     v11
+       1   1.4218   1.0555  -1.3660   0.0000   0.0000   0.0000