diff --git a/src/CGIonicStepper.C b/src/CGIonicStepper.C index 29e5b73..f3c4f52 100644 --- a/src/CGIonicStepper.C +++ b/src/CGIonicStepper.C @@ -3,7 +3,7 @@ // 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 @@ -12,41 +12,74 @@ using namespace std; //////////////////////////////////////////////////////////////////////////////// void CGIonicStepper::compute_r(double e0, const vector >& f0) { - bool wolfe1; - if ( first_step_ ) + double fp0; + 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; - alpha_ = 1.0; + for ( int i = 0; i < r0_[is].size(); i++ ) + { + largest_force = max(largest_force,fabs(f0[is][i])); + } } - else + + if ( largest_force > max_force ) { - // Check if the first Wolfe condition is satisfied - // compute predicted decrease for previous step - // pred = pc_ * ( r0 - rc_ ) - double pred = 0.0; + if ( s_.ctxt_.onpe0() ) + cout << " CGIonicStepper: force exceeds limit, taking SD step " << endl; + // take a steepest descent step with limited displacement and exit + const double alpha_sd = max_force/largest_force; + // SD step 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] ); + rp_[is][i] = r0_[is][i] + alpha_sd * f0[is][i]; } } - const double sigma_wolfe = 0.1; - wolfe1 = ( e0 < ec_ - sigma_wolfe * pred ); - //cout << "CGIonicStepper: pred = " << pred << endl; - //cout << "CGIonicStepper: ec: " << ec_ << endl; - //cout << "CGIonicStepper: required energy: " << ec_-sigma_wolfe * pred - // << " actual: " << e0 << endl; - //cout << "CGIonicStepper: wolfe1 = " << wolfe1 << endl; + constraints_.enforce_r(r0_,rp_); + rm_ = r0_; + r0_ = rp_; + atoms_.set_positions(r0_); + // reset the CG algorithm + first_step_ = true; + return; } - if ( wolfe1 ) + // CG algorithm + + if ( !first_step_ ) { - // actual decrease from last step is sufficient - // accept r0 as new current point - rc_ = r0_; - ec_ = e0; + wolfe1 = e0 < ec_ + fpc_ * sigma1_ * alpha_; + // fp0 = -proj(f0,pc) + fp0 = 0.0; + 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 if ( first_step_ ) { @@ -66,8 +99,10 @@ void CGIonicStepper::compute_r(double e0, const vector >& f0) den += fctmp * fctmp; } } - const double beta = den > 0.0 ? num/den : 0.0; - //cout << "CGIonicStepper: beta = " << beta << endl; + double beta = den > 0.0 ? num/den : 0.0; + beta = max(beta,0.0); + if ( s_.ctxt_.onpe0() ) + cout << " CGIonicStepper: beta = " << beta << endl; for ( int is = 0; is < r0_.size(); is++ ) { for ( int i = 0; i < r0_[is].size(); i++ ) @@ -77,14 +112,28 @@ void CGIonicStepper::compute_r(double e0, const vector >& f0) } } fc_ = f0; - alpha_ *= 1.05; - } - else - { - // backtrack - alpha_ *= 0.8; + // fpc = d_e / d_alpha in direction pc + fpc_ = 0.0; + for ( int is = 0; is < r0_.size(); is++ ) + { + for ( int i = 0; i < r0_[is].size(); i++ ) + { + 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 i = 0; i < r0_[is].size(); i++ ) @@ -92,6 +141,7 @@ void CGIonicStepper::compute_r(double e0, const vector >& f0) rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i]; } } + constraints_.enforce_r(r0_,rp_); rm_ = r0_; r0_ = rp_; diff --git a/src/CGIonicStepper.h b/src/CGIonicStepper.h index c05ac8f..12a4e31 100644 --- a/src/CGIonicStepper.h +++ b/src/CGIonicStepper.h @@ -3,12 +3,13 @@ // 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 #define CGIONICSTEPPER_H #include "IonicStepper.h" +#include "LineMinimizer.h" #include class CGIonicStepper : public IonicStepper @@ -19,12 +20,14 @@ class CGIonicStepper : public IonicStepper std::vector > rc_; std::vector > pc_; std::vector > fc_; - double ec_; - double alpha_; + double ec_, fpc_; + double alpha_, sigma1_, sigma2_; + LineMinimizer linmin_; 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 >& f0); void compute_v(double e0, const std::vector >& f0) {}