BMDIonicStepper.C 2.67 KB
 Francois Gygi committed Apr 30, 2009 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 . // //////////////////////////////////////////////////////////////////////////////// // // BMDIonicStepper.C // //////////////////////////////////////////////////////////////////////////////// #include "BMDIonicStepper.h" using namespace std; //////////////////////////////////////////////////////////////////////////////// void BMDIonicStepper::compute_r(double e0, const vector >& 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  Francois Gygi committed Sep 18, 2011 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 Apr 30, 2009 40 41 42 43 44 45 46 47 48 49 50 51 } //////////////////////////////////////////////////////////////////////////////// void BMDIonicStepper::compute_v(double e0, const vector >& 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++ )  Francois Gygi committed Sep 18, 2011 52  v0_[is][i] = r0_[is][i] - rm_[is][i] + bmd_fac_ * f0[is][i];  Francois Gygi committed Apr 30, 2009 53 54 55 56  // check if energy increased if ( e0_ > em_ ) {  Francois Gygi committed Sep 18, 2011 57  // rescale velocity or reset to zero  Francois Gygi committed Apr 30, 2009 58 59 60  for ( int is = 0; is < r0_.size(); is++ ) for ( int i = 0; i < r0_[is].size(); i++ ) {  Francois Gygi committed Sep 18, 2011 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 Apr 30, 2009 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++ )  Francois Gygi committed Sep 18, 2011 71  v0_[is][i] *= 1.05;  Francois Gygi committed Apr 30, 2009 72  }  Francois Gygi committed Sep 18, 2011 73  constraints_.enforce_v(r0_,v0_);  Francois Gygi committed Apr 30, 2009 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];  Francois Gygi committed Sep 18, 2011 86  ekin_ += v * v / ( 4.0 * bmd_fac_ );  Francois Gygi committed Apr 30, 2009 87 88 89  } } }