MDIonicStepper.h 2.48 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
// $Id: MDIonicStepper.h,v 1.13 2010-04-16 22:41:55 fgygi Exp $
19 20 21 22 23 24

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

#ifndef MDIONICSTEPPER_H
#define MDIONICSTEPPER_H

#include "IonicStepper.h"

class MDIonicStepper : public IonicStepper
{
  private:
57

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

  public:
68

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

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

#endif