BMDIonicStepper.C 2.67 KB
Newer Older
Francois Gygi committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// BMDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
#include "BMDIonicStepper.h"
using namespace std;

////////////////////////////////////////////////////////////////////////////////
void BMDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
{
  // f0 contains forces at r0
  // e0 contains energy at r0
  // Compute new positions rp using the velocity Verlet algorithm
  // enforce constraints for rp
  // update rm <- r0, r0 <- rp, and update atomset

30 31 32 33 34 35 36 37 38 39
  // 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_.set_positions(r0_);
Francois Gygi committed
40 41 42 43 44 45 46 47 48 49 50 51
}

////////////////////////////////////////////////////////////////////////////////
void BMDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
{
  // compute velocities v0_ using r0_, rm_ and f0(r0)
  // enforce constraints for vp
  // adjust velocities with the thermostat

  e0_ = e0;
  for ( int is = 0; is < v0_.size(); is++ )
    for ( int i = 0; i < v0_[is].size(); i++ )
52
      v0_[is][i] = r0_[is][i] - rm_[is][i] + bmd_fac_ * f0[is][i];
Francois Gygi committed
53 54 55 56

  // check if energy increased
  if ( e0_ > em_ )
  {
57
    // rescale velocity or reset to zero
Francois Gygi committed
58 59 60
    for ( int is = 0; is < r0_.size(); is++ )
      for ( int i = 0; i < r0_[is].size(); i++ )
      {
61 62 63
        // stop if velocity component is opposed to force
        if ( v0_[is][i] * f0[is][i] < 0.0 )
          v0_[is][i] = 0.0;
Francois Gygi committed
64 65 66 67 68 69 70
      }
  }
  else
  {
    // accelerate
    for ( int is = 0; is < r0_.size(); is++ )
      for ( int i = 0; i < r0_[is].size(); i++ )
71
        v0_[is][i] *= 1.05;
Francois Gygi committed
72
  }
73
  constraints_.enforce_v(r0_,v0_);
Francois Gygi committed
74 75 76 77 78 79 80 81 82 83 84 85
  compute_ekin();
}

////////////////////////////////////////////////////////////////////////////////
void BMDIonicStepper::compute_ekin(void)
{
  ekin_ = 0.0;
  for ( int is = 0; is < v0_.size(); is++ )
  {
    for ( int i = 0; i < v0_[is].size(); i++ )
    {
      const double v = v0_[is][i];
86
      ekin_ += v * v / ( 4.0 * bmd_fac_ );
Francois Gygi committed
87 88 89
    }
  }
}