Merge pull request #657 from drroe/scalevel

Add ability to scale velocities to 'setvelocity'

Author: drroe

src/Action_SetVelocity.cpp

@@ -6,12 +6,17 @@
 Action_SetVelocity::Action_SetVelocity() : tempi_(0.0), zeroMomentum_(false) {}
 
 void Action_SetVelocity::Help() const {
-  mprintf("\t[<mask>] [{tempi <temperature> | modify}] [ig <random seed>]\n"
+  mprintf("\t[<mask>] [{ tempi <temperature> |\n"
+          "\t            scale [factor <fac>] [sx <xfac>] [sy <yfac>] [sz <zfac>] |\n"
+          "\t            none | \n"
+          "\t            modify}] [ig <random seed>]\n"
           "\t[%s] [%s]\n", Constraints::constraintArgs, Constraints::rattleArgs);
-  mprintf("\t[zeromomentum]\n");
+  mprintf("\t[zeromomentum] [ig <random seed>]\n");
   mprintf("  Set velocities in frame for atoms in <mask> using Maxwellian distribution\n" 
           "  based on given temperature; default 300.0 K. If tempi is 0.0 set\n"
           "  velocities of atoms in mask to 0.0.\n"
+          "  If 'scale' is specified scale the velocities by the given factor(s).\n"
+          "  If 'none' or a temperature of 0.0 is specified, velocities will be zeroed.\n"
           "  If 'modify' is specified do not set; only modify existing velocities.\n"
           "  This is useful e.g. with 'ntc' or 'zeromomentum'.\n"
           "  If 'ntc' is specified attempt to correct velocities for constraints.\n"
@@ -24,11 +29,17 @@ Action::RetType Action_SetVelocity::Init(ArgList& actionArgs, ActionInit& init,
 {
   // Keywords
   tempi_ = actionArgs.getKeyDouble("tempi", 300.0);
-  if (tempi_ < Constants::SMALL)
+  if (actionArgs.hasKey("none") || tempi_ < Constants::SMALL)
     mode_ = ZERO;
   else if (actionArgs.hasKey("modify"))
     mode_ = MODIFY;
-  else
+  else if (actionArgs.hasKey("scale")) {
+    double sf = actionArgs.getKeyDouble("factor", 1.0);
+    scaleFac_[0] = actionArgs.getKeyDouble("sx", sf);
+    scaleFac_[1] = actionArgs.getKeyDouble("sy", sf);
+    scaleFac_[2] = actionArgs.getKeyDouble("sz", sf);
+    mode_ = SCALE;
+  } else
     mode_ = SET;
   int ig_ = actionArgs.getKeyInt("ig", -1);
   RN_.rn_set( ig_ );
@@ -56,6 +67,9 @@ Action::RetType Action_SetVelocity::Init(ArgList& actionArgs, ActionInit& init,
       mprintf("\tRandom seed is %i\n", ig_);
   } else if (mode_ == MODIFY)
     mprintf(" Modifying any existing velocities for atoms in mask '%s'\n", Mask_.MaskString());
+  else if (mode_ == SCALE)
+    mprintf(" Scaling velocities by X= %g, Y= %g, Z= %g\n",
+            scaleFac_[0], scaleFac_[1], scaleFac_[2]);
   else if (mode_ == ZERO)
     mprintf(" Zeroing velocities for atoms in mask '%s'\n", Mask_.MaskString());
   if (cons_.Type() != Constraints::OFF) {
@@ -101,6 +115,11 @@ Action::RetType Action_SetVelocity::Setup(ActionSetup& setup) {
     mprintf("Warning: 'modify' specified but no velocity info, skipping.\n");
     return Action::SKIP;
   }
+  // If scale need to have existing velocity info
+  if (mode_ == SCALE && !setup.CoordInfo().HasVel()) {
+    mprintf("Warning: 'scale' specified but no velocity info, skipping.\n");
+    return Action::SKIP;
+  }
   // Always add velocity info even if not strictly necessary
   cInfo_ = setup.CoordInfo();
   cInfo_.SetVelocity( true );
@@ -132,6 +151,14 @@ Action::RetType Action_SetVelocity::DoAction(int frameNum, ActionFrame& frm) {
       V[1] = RN_.rn_gauss(0.0, *sd);
       V[2] = RN_.rn_gauss(0.0, *sd);
     }
+  } else if (mode_ == SCALE) {
+    for (AtomMask::const_iterator atom = Mask_.begin(); atom != Mask_.end(); ++atom)
+    {
+      double* V = newFrame_.vAddress() + (*atom * 3);
+      V[0] *= scaleFac_[0];
+      V[1] *= scaleFac_[1];
+      V[2] *= scaleFac_[2];
+    }
   }
   // Correct velocities for constraints
   if (cons_.Type() != Constraints::OFF)

src/Action_SetVelocity.h

@@ -15,13 +15,14 @@ class Action_SetVelocity : public Action {
     Action::RetType DoAction(int, ActionFrame&);
     void Print() {}
 
-    enum ModeType { SET = 0, ZERO, MODIFY };
+    enum ModeType { SET = 0, ZERO, MODIFY, SCALE };
 
     typedef std::vector<double> Darray;
 
     AtomMask Mask_;    ///< Atoms to set/modify velocities of
     Darray SD_;        ///< Hold sqrt(kB*(1/mass)) for each atom
     double tempi_;     ///< Temperature to generate velocity distribution at.
+    Vec3 scaleFac_;    ///< Velocity scaling factors.
     ModeType mode_;    ///< Set velocity, zero velocity, modify existing velocity.
     Constraints cons_; ///< Hold constraint info 
     Random_Number RN_;

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.7"
+#define CPPTRAJ_INTERNAL_VERSION "V4.10.8"
 #endif

test/Test_SetVelocity/RunTest.sh

@@ -2,12 +2,14 @@
 
 . ../MasterTest.sh
 
-CleanFiles setvel.in tz2.vel.rst7 V1.dat
+CleanFiles setvel.in tz2.vel.rst7 V1.dat tz2.scale.rst7
 
 INPUT='-i setvel.in'
 
-TESTNAME='Set Velocity test'
+TESTNAME='Set Velocity tests'
 Requires maxthreads 1
+
+UNITNAME='Set Velocity test'
 cat > setvel.in <<EOF
 parm ../tz2.parm7
 trajin ../tz2.rst7
@@ -25,9 +27,20 @@ setvelocity modify zeromomentum
 vector V1 momentum out V1.dat
 go
 EOF
-RunCpptraj "$TESTNAME"
+RunCpptraj "$UNITNAME"
 DoTest tz2.vel.rst7.save tz2.vel.rst7
 DoTest V1.dat.save V1.dat
 
+UNITNAME='Scale velocity test'
+cat > setvel.in <<EOF
+parm ../tz2.parm7
+trajin tz2.vel.rst7.save
+strip !:13
+setvelocity scale factor 2.0
+trajout tz2.scale.rst7
+EOF
+RunCpptraj "$UNITNAME"
+DoTest tz2.scale.rst7.save tz2.scale.rst7
+
 EndTest
 exit 0

test/Test_SetVelocity/tz2.scale.rst7.save

@@ -0,0 +1,8 @@
+Cpptraj Generated Restart                                                       
+    6
+  -4.8320000 -10.2020000  -0.3440000  -4.6470000  -9.2390000  -0.5870000
+  -5.7200000 -10.8660000  -1.2770000  -5.5910000 -11.9460000  -1.3400000
+  -6.7570000 -10.8860000  -0.9420000  -5.7770000 -10.3580000  -2.2400000
+   0.3643370   0.0717144   0.1748470   1.7579932   1.0938414   2.7569718
+  -0.7964882   0.7930536  -0.1924086   1.0052930  -0.5603402   0.6513296
+  -0.3944468   2.7520096   3.7998472  -0.7440500  -1.4602778  -0.5855208