1 //////////////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (c) 2008 The Regents of the University of California 4 // 5 // This file is part of Qbox 6 // 7 // Qbox is distributed under the terms of the GNU General Public License 8 // as published by the Free Software Foundation, either version 2 of 9 // the License, or (at your option) any later version. 10 // See the file COPYING in the root directory of this distribution 11 // or <http://www.gnu.org/licenses/>. 12 // 13 //////////////////////////////////////////////////////////////////////////////// 14 // 15 // BMDIonicStepper.h: 16 // 17 //////////////////////////////////////////////////////////////////////////////// 18 // 19 // IonicStepper is used in the following way 20 // 21 // input: r0,v0 22 // 23 // compute energy e0(r0) and forces f0(r0) 24 // for ( k=0; k<n; k++ ) 25 // { 26 // // r0,v0,e0,f0 known 27 // stepper->compute_r(e0,f0) 28 // { 29 // computes rp using r0, v0 and f0 30 // restores constraints on rp using rp and r0 31 // updates rm<-r0, r0<-rp and update atomset positions 32 // } 33 // 34 // compute f0(r0) 35 // 36 // stepper->compute_v(e0,f0) 37 // { 38 // computes v0 using r0,rm,f0 39 // restores constraints on v0 using r0, v0 40 // modifies velocities using the thermostat (rescaling) 41 // updates atomset velocities 42 // } 43 // } 44 // r0,v0,f0 consistent at this point 45 // 46 47 #ifndef BMDIONICSTEPPER_H 48 #define BMDIONICSTEPPER_H 49 50 #include "IonicStepper.h" 51 52 class BMDIonicStepper : public IonicStepper 53 { 54 private: 55 56 const double bmd_fac_; 57 double e0_,em_; 58 double ekin_; 59 std::vector<std::vector< double> > fm_; 60 void compute_ekin(void); 61 62 public: 63 BMDIonicStepper(Sample & s)64 BMDIonicStepper(Sample& s) : IonicStepper(s), bmd_fac_(0.025) 65 { 66 e0_ = 0.0; 67 em_ = 0.0; 68 ekin_ = 0.0; 69 atoms_.get_positions(r0_); 70 atoms_.get_velocities(v0_); 71 fm_.resize(r0_.size()); 72 for ( int is = 0; is < fm_.size(); is++ ) 73 fm_[is].resize(r0_[is].size()); 74 compute_ekin(); 75 } 76 77 void compute_r(double e0, const std::vector<std::vector< double> >& f0); 78 void compute_v(double e0, const std::vector<std::vector< double> >& f0); ekin(void)79 double ekin(void) const { return ekin_; } temp(void)80 double temp(void) const 81 { 82 const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 ); 83 if ( ndofs_ > 0 ) 84 return 2.0 * ( ekin_ / boltz ) / ndofs_; 85 else 86 return 0.0; 87 } 88 }; 89 90 #endif 91