Commit ac685f95 by Francois Gygi

reverted to simple SD algorithm usins dt^2/mass


git-svn-id: http://qboxcode.org/svn/qb/trunk@594 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 0797b6b3
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// SDIonicStepper.C // SDIonicStepper.C
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: SDIonicStepper.C,v 1.5 2008-03-26 04:57:54 fgygi Exp $ // $Id: SDIonicStepper.C,v 1.6 2008-04-06 17:47:36 fgygi Exp $
#include "SDIonicStepper.h" #include "SDIonicStepper.h"
using namespace std; using namespace std;
...@@ -11,61 +11,17 @@ using namespace std; ...@@ -11,61 +11,17 @@ using namespace std;
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void SDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0) void SDIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
{ {
bool wolfe1; // Steepest descent step
if ( first_step_ )
{
wolfe1 = true;
alpha_ = 1.0;
}
else
{
// Check if the first Wolfe condition is satisfied
// compute predicted decrease for previous step
// pred = pc_ * ( r0 - rc_ )
double pred = 0.0;
for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
pred += pc_[is][i] * ( r0_[is][i] - rc_[is][i]);
}
}
const double sigma_wolfe = 0.1;
wolfe1 = ( e0 < ec_ - sigma_wolfe * pred );
//cout << "SDIonicStepper: pred = " << pred << endl;
//cout << "SDIonicStepper: ec: " << ec_ << endl;
//cout << "SDIonicStepper: required energy: " << ec_-sigma_wolfe * pred
// << " actual: " << e0 << endl;
//cout << "SDIonicStepper: wolfe1 = " << wolfe1 << endl;
}
if ( wolfe1 )
{
// actual decrease from last step is sufficient
// accept r0 as new current point
rc_ = r0_;
ec_ = e0;
// define new descent direction
pc_ = f0;
alpha_ *= 1.05;
}
else
{
// backtrack
alpha_ *= 0.8;
}
//cout << "SDIonicStepper: alpha = " << alpha_ << endl;
for ( int is = 0; is < r0_.size(); is++ ) for ( int is = 0; is < r0_.size(); is++ )
{ {
const double dt2bym = dt_ * dt_ / pmass_[is];
for ( int i = 0; i < r0_[is].size(); i++ ) for ( int i = 0; i < r0_[is].size(); i++ )
{ {
rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i]; rp_[is][i] = r0_[is][i] + dt2bym * f0[is][i];
} }
} }
constraints_.enforce_r(r0_,rp_); constraints_.enforce_r(r0_,rp_);
rm_ = r0_; rm_ = r0_;
r0_ = rp_; r0_ = rp_;
atoms_.set_positions(r0_); atoms_.set_positions(r0_);
first_step_ = false;
} }
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// SDIonicStepper.h: // SDIonicStepper.h:
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: SDIonicStepper.h,v 1.7 2008-03-26 04:57:54 fgygi Exp $ // $Id: SDIonicStepper.h,v 1.8 2008-04-06 17:47:36 fgygi Exp $
#ifndef SDIONICSTEPPER_H #ifndef SDIONICSTEPPER_H
#define SDIONICSTEPPER_H #define SDIONICSTEPPER_H
...@@ -15,15 +15,9 @@ class SDIonicStepper : public IonicStepper ...@@ -15,15 +15,9 @@ class SDIonicStepper : public IonicStepper
{ {
private: private:
bool first_step_;
std::vector<std::vector< double> > rc_;
std::vector<std::vector< double> > pc_;
double ec_;
double alpha_;
public: public:
SDIonicStepper(Sample& s) : IonicStepper(s), first_step_(true) {} SDIonicStepper(Sample& s) : IonicStepper(s) {}
void compute_r(double e0, const std::vector<std::vector< double> >& f0); void compute_r(double e0, const std::vector<std::vector< double> >& f0);
void compute_v(double e0, const std::vector<std::vector< double> >& f0) {} void compute_v(double e0, const std::vector<std::vector< double> >& f0) {}
......
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