Commit 25c5a3f5 by Francois Gygi

New thermostat options: SCALING ANDERSON LOWE OFF


git-svn-id: http://qboxcode.org/svn/qb/trunk@468 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 60604bea
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// MDIonicStepper.C // MDIonicStepper.C
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.C,v 1.9 2005-06-27 22:28:33 fgygi Exp $ // $Id: MDIonicStepper.C,v 1.10 2006-08-22 15:16:10 fgygi Exp $
#include "MDIonicStepper.h" #include "MDIonicStepper.h"
...@@ -50,7 +50,8 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0) ...@@ -50,7 +50,8 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
} }
} }
compute_ekin(); compute_ekin();
if ( thermostat_ )
if ( thermostat_ == "SCALING" )
{ {
eta_ = tanh ( ( temp() - th_temp_ ) / th_width_ ) / th_time_; eta_ = tanh ( ( temp() - th_temp_ ) / th_width_ ) / th_time_;
if ( s_.ctxt_.onpe0() ) if ( s_.ctxt_.onpe0() )
...@@ -70,6 +71,97 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0) ...@@ -70,6 +71,97 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
} }
ekin_ *= fac * fac; ekin_ *= fac * fac;
} }
else if ( thermostat_ == "ANDERSON" )
{
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
for ( int is = 0; is < nsp_; is++ )
{
const double m = pmass_[is];
const double width = sqrt( boltz * th_temp_ / m );
for ( int ia = 0; ia < na_[is]; ia++ )
{
// if th_time is zero, set probability to one
if ( th_time_ == 0.0 || drand48() < dt_ / th_time_ )
{
// draw gaussian random variables
double xi0 = drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() - 6.0;
double xi1 = drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() - 6.0;
double xi2 = drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() - 6.0;
v0_[is][3*ia+0] = width * xi0;
v0_[is][3*ia+1] = width * xi1;
v0_[is][3*ia+2] = width * xi2;
}
}
}
}
else if ( thermostat_ == "LOWE" )
{
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
// scan all atom pairs in the space (is,ia)
//int npairs = 0;
for ( int is1 = 0; is1 < nsp_; is1++ )
{
for ( int is2 = is1; is2 < nsp_; is2++ )
{
const double m1 = pmass_[is1];
const double m2 = pmass_[is2];
const double mu = m1*m2/(m1+m2);
const double width = sqrt(boltz * th_temp_ / mu);
for ( int ia1 = 0; ia1 < na_[is1]; ia1++ )
{
int ia2min = 0;
if ( is1 == is2 ) ia2min = ia1 + 1;
for ( int ia2 = ia2min; ia2 < na_[is2]; ia2++ )
{
// if th_time is zero, set probability to one
if ( th_time_ == 0.0 || drand48() < dt_ / th_time_ )
{
//cout << " pair " << is1 << " " << ia1 << " "
// << is2 << " " << ia2 << endl;
D3vector r1(&r0_[is1][3*ia1]);
D3vector r2(&r0_[is2][3*ia2]);
D3vector r12 = r1 - r2;
s_.wf.cell().fold_in_ws(r12);
D3vector e12 = normalized(r12);
D3vector v1(&v0_[is1][3*ia1]);
D3vector v2(&v0_[is2][3*ia2]);
D3vector v12 = v1 - v2;
// draw a gaussian random variable
double xi = drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() +
drand48() + drand48() + drand48() - 6.0;
double lambda = xi * width;
D3vector dv12 = mu * ( lambda - v12*e12 ) * e12;
D3vector v1p = v1 + (1.0/m1) * dv12;
D3vector v2p = v2 - (1.0/m2) * dv12;
v0_[is1][3*ia1+0] = v1p.x;
v0_[is1][3*ia1+1] = v1p.y;
v0_[is1][3*ia1+2] = v1p.z;
v0_[is2][3*ia2+0] = v2p.x;
v0_[is2][3*ia2+1] = v2p.y;
v0_[is2][3*ia2+2] = v2p.z;
}
//npairs++;
}
}
}
}
//cout << " npairs: " << npairs << endl;
compute_ekin();
}
constraints_.enforce_v(r0_,v0_); constraints_.enforce_v(r0_,v0_);
atoms_.set_velocities(v0_); atoms_.set_velocities(v0_);
} }
......
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// MDIonicStepper.h: // MDIonicStepper.h:
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.h,v 1.7 2006-05-13 05:38:17 fgygi Exp $ // $Id: MDIonicStepper.h,v 1.8 2006-08-22 15:16:10 fgygi Exp $
// //
// IonicStepper is used in the following way // IonicStepper is used in the following way
...@@ -48,14 +48,14 @@ class MDIonicStepper : public IonicStepper ...@@ -48,14 +48,14 @@ class MDIonicStepper : public IonicStepper
double th_width_; double th_width_;
double ekin_; double ekin_;
double eta_; double eta_;
bool thermostat_; string thermostat_;
void compute_ekin(void); void compute_ekin(void);
public: public:
MDIonicStepper(Sample& s) : IonicStepper(s) MDIonicStepper(Sample& s) : IonicStepper(s)
{ {
thermostat_ = ( s.ctrl.thermostat == "ON" ); thermostat_ = s.ctrl.thermostat;
th_temp_ = s.ctrl.th_temp; th_temp_ = s.ctrl.th_temp;
th_time_ = s.ctrl.th_time; th_time_ = s.ctrl.th_time;
th_width_ = s.ctrl.th_width; th_width_ = s.ctrl.th_width;
......
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// Thermostat.h // Thermostat.h
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: Thermostat.h,v 1.1 2003-08-28 00:28:51 fgygi Exp $ // $Id: Thermostat.h,v 1.2 2006-08-22 15:16:10 fgygi Exp $
#ifndef THERMOSTAT_H #ifndef THERMOSTAT_H
#define THERMOSTAT_H #define THERMOSTAT_H
...@@ -33,10 +33,11 @@ class Thermostat : public Var ...@@ -33,10 +33,11 @@ class Thermostat : public Var
} }
string v = argv[1]; string v = argv[1];
if ( !( v == "ON" || v == "OFF" ) ) if ( !( v == "SCALING" || v == "ANDERSON" || v == "LOWE" || v == "OFF" ) )
{ {
if ( ui->onpe0() ) if ( ui->onpe0() )
cout << " thermostat must be ON or OFF" << endl; cout << " thermostat must be SCALING or ANDERSON or LOWE or OFF"
<< endl;
return 1; return 1;
} }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment