Commit a65c6441 by Francois Gygi

Modif to allow negative time steps


git-svn-id: http://qboxcode.org/svn/qb/trunk@491 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 642f6c00
......@@ -3,7 +3,7 @@
// Dt.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Dt.h,v 1.1 2003-06-11 22:10:11 fgygi Exp $
// $Id: Dt.h,v 1.2 2007-03-17 01:22:50 fgygi Exp $
#ifndef DT_H
#define DT_H
......@@ -33,10 +33,10 @@ class Dt : public Var
}
double v = atof(argv[1]);
if ( v < 0.0 )
if ( v == 0.0 )
{
if ( ui->onpe0() )
cout << " dt must be non-negative" << endl;
cout << " dt must be non-zero" << endl;
return 1;
}
......
......@@ -3,7 +3,7 @@
// MDIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDIonicStepper.C,v 1.12 2007-03-17 01:14:00 fgygi Exp $
// $Id: MDIonicStepper.C,v 1.13 2007-03-17 01:22:50 fgygi Exp $
#include "MDIonicStepper.h"
using namespace std;
......@@ -39,7 +39,7 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
// enforce constraints for vp
// adjust velocities with the thermostat
assert(dt_ > 0.0);
assert(dt_ != 0.0);
for ( int is = 0; is < v0_.size(); is++ )
{
assert(pmass_[is] > 0.0);
......@@ -62,7 +62,7 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
cout << " <!-- thermostat: eta=" << eta_ << " -->" << endl;
}
const double fac = (1.0 - eta_ * dt_);
const double fac = (1.0 - eta_ * fabs(dt_));
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
......@@ -82,7 +82,7 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
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_ )
if ( th_time_ == 0.0 || drand48() < fabs(dt_) / th_time_ )
{
// draw gaussian random variables
double xi0 = drand48() + drand48() + drand48() +
......@@ -125,7 +125,7 @@ void MDIonicStepper::compute_v(double e0, const vector<vector< double> >& f0)
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_ )
if ( th_time_ == 0.0 || drand48() < fabs(dt_) / th_time_ )
{
//cout << " pair " << is1 << " " << ia1 << " "
// << is2 << " " << ia2 << endl;
......
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