modified velocity rescaling alg, enforce constraints on velocity,

use a common ionic mass


git-svn-id: http://qboxcode.org/svn/qb/trunk@988 cba15fb0-1239-40c8-b417-11db7ca47a34

src/BMDIonicStepper.C

@@ -15,8 +15,6 @@
 // BMDIonicStepper.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: BMDIonicStepper.C,v 1.1 2009-04-30 23:13:14 fgygi Exp $
-
 #include "BMDIonicStepper.h"
 using namespace std;
 
@@ -29,41 +27,17 @@ void BMDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
   // enforce constraints for rp
   // update rm <- r0, r0 <- rp, and update atomset
 
-  if ( e0 > em_ )
-  {
-    // energy has increased, reset r0 to rm and make a steepest descent step
-    r0_ = rm_;
-    for ( int is = 0; is < r0_.size(); is++ )
-    {
-      const double dt2by2m = dt_ * dt_ / ( 2.0 * pmass_[is] );
-      for ( int i = 0; i < r0_[is].size(); i++ )
-      {
-        rp_[is][i] = rm_[is][i] + dt2by2m * fm_[is][i];
-      }
-    }
-    constraints_.enforce_r(r0_,rp_);
-    r0_ = rp_;
-    e0_ = em_;
-    atoms_.set_positions(r0_);
-  }
-  else
-  {
-    // compute rp
-    for ( int is = 0; is < r0_.size(); is++ )
-    {
-      const double dt2by2m = dt_ * dt_ / ( 2.0 * pmass_[is] );
-      for ( int i = 0; i < r0_[is].size(); i++ )
-      {
-        rp_[is][i] = r0_[is][i] + v0_[is][i] * dt_ + dt2by2m * f0[is][i];
-      }
-    }
-    constraints_.enforce_r(r0_,rp_);
-    rm_ = r0_;
-    fm_ = f0;
-    em_ = e0_;
-    r0_ = rp_;
-    atoms_.set_positions(r0_);
-  }
+  // compute rp
+  for ( int is = 0; is < r0_.size(); is++ )
+    for ( int i = 0; i < r0_[is].size(); i++ )
+      rp_[is][i] = r0_[is][i] + v0_[is][i] + bmd_fac_ * f0[is][i];
+  constraints_.enforce_r(r0_,rp_);
+  rm_ = r0_;
+  fm_ = f0;
+  em_ = e0_;
+  r0_ = rp_;
+  atoms_.sync_positions(r0_);
+  atoms_.set_positions(r0_);
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -74,41 +48,30 @@ void BMDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
   // adjust velocities with the thermostat
 
   e0_ = e0;
-  assert(dt_ != 0.0);
   for ( int is = 0; is < v0_.size(); is++ )
-  {
-    assert(pmass_[is] > 0.0);
-    const double dtby2m = dt_ / ( 2.0 * pmass_[is] );
     for ( int i = 0; i < v0_[is].size(); i++ )
-    {
-      const double vhalf = ( r0_[is][i] - rm_[is][i] ) / dt_;
-      v0_[is][i] = vhalf + dtby2m * f0[is][i];
-    }
-  }
+      v0_[is][i] = r0_[is][i] - rm_[is][i] + bmd_fac_ * f0[is][i];
 
   // check if energy increased
   if ( e0_ > em_ )
   {
-    // reset the velocity to zero
+    // rescale velocity or reset to zero
     for ( int is = 0; is < r0_.size(); is++ )
-    {
       for ( int i = 0; i < r0_[is].size(); i++ )
       {
-        v0_[is][i] = 0.0;
+        // stop if velocity component is opposed to force
+        if ( v0_[is][i] * f0[is][i] < 0.0 )
+          v0_[is][i] = 0.0;
       }
-    }
   }
   else
   {
     // accelerate
     for ( int is = 0; is < r0_.size(); is++ )
-    {
       for ( int i = 0; i < r0_[is].size(); i++ )
-      {
-        v0_[is][i] *= 1.1;
-      }
-    }
+        v0_[is][i] *= 1.05;
   }
+  constraints_.enforce_v(r0_,v0_);
   compute_ekin();
 }
 
@@ -121,7 +84,7 @@ void BMDIonicStepper::compute_ekin(void)
     for ( int i = 0; i < v0_[is].size(); i++ )
     {
       const double v = v0_[is][i];
-      ekin_ += 0.5 * pmass_[is] * v * v;
+      ekin_ += v * v / ( 4.0 * bmd_fac_ );
     }
   }
 }

src/BMDIonicStepper.h

@@ -15,8 +15,6 @@
 // BMDIonicStepper.h:
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: BMDIonicStepper.h,v 1.1 2009-04-30 23:13:14 fgygi Exp $
-
 //
 // IonicStepper is used in the following way
 //
@@ -55,6 +53,7 @@ class BMDIonicStepper : public IonicStepper
 {
   private:
 
+  const double bmd_fac_;
   double e0_,em_;
   double ekin_;
   std::vector<std::vector< double> > fm_;
@@ -62,7 +61,7 @@ class BMDIonicStepper : public IonicStepper
 
   public:
 
-  BMDIonicStepper(Sample& s) : IonicStepper(s)
+  BMDIonicStepper(Sample& s) : IonicStepper(s), bmd_fac_(0.025)
   {
     e0_ = 0.0;
     em_ = 0.0;