MDIonicStepper.h 2.42 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17
// MDIonicStepper.h:
//
////////////////////////////////////////////////////////////////////////////////
18 19 20 21 22 23

//
// IonicStepper is used in the following way
//
// input: r0,v0
//
Francois Gygi committed
24
// compute energy e0(r0) and forces f0(r0)
25 26
// for ( k=0; k<n; k++ )
// {
Francois Gygi committed
27 28
//   // r0,v0,e0,f0 known
//   stepper->compute_r(e0,f0)
Francois Gygi committed
29 30 31 32 33
//   {
//     computes rp using r0, v0 and f0
//     restores constraints on rp using rp and r0
//     updates rm<-r0, r0<-rp and update atomset positions
//   }
34
//
Francois Gygi committed
35 36
//   compute f0(r0)
//
Francois Gygi committed
37
//   stepper->compute_v(e0,f0)
Francois Gygi committed
38 39 40 41 42 43 44 45
//   {
//     computes v0 using r0,rm,f0
//     restores constraints on v0 using r0, v0
//     modifies velocities using the thermostat (rescaling)
//     updates atomset velocities
//   }
// }
// r0,v0,f0 consistent at this point
46
//
47 48 49 50 51 52 53 54 55

#ifndef MDIONICSTEPPER_H
#define MDIONICSTEPPER_H

#include "IonicStepper.h"

class MDIonicStepper : public IonicStepper
{
  private:
56

57 58
  double th_temp_;
  double th_time_;
59
  double th_width_;
60
  double ekin_;
61
  double ekin_stepper_;
62
  double eta_;
63
  std::string thermostat_;
Francois Gygi committed
64
  void compute_ekin(void);
65 66

  public:
67

68 69
  MDIonicStepper(Sample& s) : IonicStepper(s)
  {
70
    thermostat_ = s.ctrl.thermostat;
71 72
    th_temp_ = s.ctrl.th_temp;
    th_time_ = s.ctrl.th_time;
73
    th_width_ = s.ctrl.th_width;
74
    eta_ = 0.0;
75
    ekin_ = 0.0;
76
    ekin_stepper_ = 0.0;
Francois Gygi committed
77
    atoms_.get_positions(r0_);
78
    atoms_.get_velocities(v0_);
Francois Gygi committed
79
    compute_ekin();
80 81
  }

82 83
  void compute_r(double e0, const std::vector<std::vector< double> >& f0);
  void compute_v(double e0, const std::vector<std::vector< double> >& f0);
84
  double eta(void) const { return eta_; }
85
  double ekin(void) const { return ekin_; }
86
  double ekin_stepper(void) const { return ekin_stepper_; }
87 88 89
  double temp(void) const
  {
    const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
Francois Gygi committed
90
    if ( ndofs_ > 0 )
91 92 93 94
      return 2.0 * ( ekin_ / boltz ) / ndofs_;
    else
      return 0.0;
  }
95 96 97
};

#endif