Commit df62fad8 by Francois Gygi

conjugate gradient ionic stepper


git-svn-id: http://qboxcode.org/svn/qb/trunk@592 cba15fb0-1239-40c8-b417-11db7ca47a34
parent dfb790bd
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// CGIonicStepper.C // CGIonicStepper.C
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: CGIonicStepper.C,v 1.1 2008-03-26 04:57:54 fgygi Exp $ // $Id: CGIonicStepper.C,v 1.2 2008-04-06 17:46:50 fgygi Exp $
#include "CGIonicStepper.h" #include "CGIonicStepper.h"
#include <cassert> #include <cassert>
...@@ -12,41 +12,74 @@ using namespace std; ...@@ -12,41 +12,74 @@ using namespace std;
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0) void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
{ {
bool wolfe1; double fp0;
if ( first_step_ ) bool wolfe1, wolfe2;
// check if the largest component of the force f0 is larger than max_force
const double max_force = 0.1;
double largest_force = 0.0;
for ( int is = 0; is < r0_.size(); is++ )
{ {
wolfe1 = true; for ( int i = 0; i < r0_[is].size(); i++ )
alpha_ = 1.0; {
largest_force = max(largest_force,fabs(f0[is][i]));
}
} }
else
if ( largest_force > max_force )
{ {
// Check if the first Wolfe condition is satisfied if ( s_.ctxt_.onpe0() )
// compute predicted decrease for previous step cout << " CGIonicStepper: force exceeds limit, taking SD step " << endl;
// pred = pc_ * ( r0 - rc_ ) // take a steepest descent step with limited displacement and exit
double pred = 0.0; const double alpha_sd = max_force/largest_force;
// SD step
for ( int is = 0; is < r0_.size(); is++ ) for ( int is = 0; is < r0_.size(); is++ )
{ {
for ( int i = 0; i < r0_[is].size(); i++ ) for ( int i = 0; i < r0_[is].size(); i++ )
{ {
pred += pc_[is][i] * ( r0_[is][i] - rc_[is][i] ); rp_[is][i] = r0_[is][i] + alpha_sd * f0[is][i];
} }
} }
const double sigma_wolfe = 0.1; constraints_.enforce_r(r0_,rp_);
wolfe1 = ( e0 < ec_ - sigma_wolfe * pred ); rm_ = r0_;
//cout << "CGIonicStepper: pred = " << pred << endl; r0_ = rp_;
//cout << "CGIonicStepper: ec: " << ec_ << endl; atoms_.set_positions(r0_);
//cout << "CGIonicStepper: required energy: " << ec_-sigma_wolfe * pred // reset the CG algorithm
// << " actual: " << e0 << endl; first_step_ = true;
//cout << "CGIonicStepper: wolfe1 = " << wolfe1 << endl; return;
} }
if ( wolfe1 ) // CG algorithm
if ( !first_step_ )
{ {
// actual decrease from last step is sufficient wolfe1 = e0 < ec_ + fpc_ * sigma1_ * alpha_;
// accept r0 as new current point // fp0 = -proj(f0,pc)
rc_ = r0_; fp0 = 0.0;
ec_ = e0; for ( int is = 0; is < r0_.size(); is++ )
{
for ( int i = 0; i < r0_[is].size(); i++ )
{
fp0 -= f0[is][i] * pc_[is][i];
}
}
wolfe2 = fabs(fp0) < sigma2_ * fabs(fpc_);
if ( s_.ctxt_.onpe0() )
{
cout << " CGIonicStepper: fpc = " << fpc_ << endl;
cout << " CGIonicStepper: fp0 = " << fp0 << endl;
cout << " CGIonicStepper: ec = " << ec_ << " e0 = " << e0 << endl;
cout << " CGIonicStepper: ec_ + fpc_ * sigma1_ * alpha_ ="
<< ec_ + fpc_ * sigma1_ * alpha_ << endl;
cout << " CGIonicStepper: wolfe1/wolfe2 = "
<< wolfe1 << "/" << wolfe2 << endl;
}
}
if ( first_step_ || (wolfe1 && wolfe2) )
{
// set new descent direction
// pc = f0
// define new descent direction // define new descent direction
if ( first_step_ ) if ( first_step_ )
{ {
...@@ -66,8 +99,10 @@ void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0) ...@@ -66,8 +99,10 @@ void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
den += fctmp * fctmp; den += fctmp * fctmp;
} }
} }
const double beta = den > 0.0 ? num/den : 0.0; double beta = den > 0.0 ? num/den : 0.0;
//cout << "CGIonicStepper: beta = " << beta << endl; beta = max(beta,0.0);
if ( s_.ctxt_.onpe0() )
cout << " CGIonicStepper: beta = " << beta << endl;
for ( int is = 0; is < r0_.size(); is++ ) for ( int is = 0; is < r0_.size(); is++ )
{ {
for ( int i = 0; i < r0_[is].size(); i++ ) for ( int i = 0; i < r0_[is].size(); i++ )
...@@ -77,14 +112,28 @@ void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0) ...@@ -77,14 +112,28 @@ void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
} }
} }
fc_ = f0; fc_ = f0;
alpha_ *= 1.05; // fpc = d_e / d_alpha in direction pc
} fpc_ = 0.0;
else for ( int is = 0; is < r0_.size(); is++ )
{ {
// backtrack for ( int i = 0; i < r0_[is].size(); i++ )
alpha_ *= 0.8; {
fpc_ -= fc_[is][i] * pc_[is][i];
}
}
ec_ = e0;
rc_ = r0_;
fp0 = fpc_;
// reset line minimizer
linmin_.reset();
} }
//cout << "CGIonicStepper: alpha = " << alpha_ << endl;
alpha_ = linmin_.newalpha(alpha_,e0,fp0);
if ( s_.ctxt_.onpe0() )
cout << " CGIonicStepper: alpha = " << alpha_ << endl;
// rp = rc + alpha * pc
for ( int is = 0; is < r0_.size(); is++ ) for ( int is = 0; is < r0_.size(); is++ )
{ {
for ( int i = 0; i < r0_[is].size(); i++ ) for ( int i = 0; i < r0_[is].size(); i++ )
...@@ -92,6 +141,7 @@ void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0) ...@@ -92,6 +141,7 @@ void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i]; rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i];
} }
} }
constraints_.enforce_r(r0_,rp_); constraints_.enforce_r(r0_,rp_);
rm_ = r0_; rm_ = r0_;
r0_ = rp_; r0_ = rp_;
......
...@@ -3,12 +3,13 @@ ...@@ -3,12 +3,13 @@
// CGIonicStepper.h: // CGIonicStepper.h:
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: CGIonicStepper.h,v 1.1 2008-03-26 04:57:54 fgygi Exp $ // $Id: CGIonicStepper.h,v 1.2 2008-04-06 17:46:50 fgygi Exp $
#ifndef CGIONICSTEPPER_H #ifndef CGIONICSTEPPER_H
#define CGIONICSTEPPER_H #define CGIONICSTEPPER_H
#include "IonicStepper.h" #include "IonicStepper.h"
#include "LineMinimizer.h"
#include <vector> #include <vector>
class CGIonicStepper : public IonicStepper class CGIonicStepper : public IonicStepper
...@@ -19,12 +20,14 @@ class CGIonicStepper : public IonicStepper ...@@ -19,12 +20,14 @@ class CGIonicStepper : public IonicStepper
std::vector<std::vector< double> > rc_; std::vector<std::vector< double> > rc_;
std::vector<std::vector< double> > pc_; std::vector<std::vector< double> > pc_;
std::vector<std::vector< double> > fc_; std::vector<std::vector< double> > fc_;
double ec_; double ec_, fpc_;
double alpha_; double alpha_, sigma1_, sigma2_;
LineMinimizer linmin_;
public: public:
CGIonicStepper(Sample& s) : IonicStepper(s), first_step_(true) {} CGIonicStepper(Sample& s) : IonicStepper(s), first_step_(true),
sigma1_(0.1), sigma2_(0.5) { linmin_.set_sigma1(sigma1_); }
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