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