Logo Search packages:      
Sourcecode: ball version File versions

conjugateGradient.C

// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
// $Id: conjugateGradient.C,v 1.38.8.7 2007/08/07 09:13:05 aleru Exp $
//
// Minimize the potential energy of a system using a nonlinear conjugate 
// gradient method with  line search

#include <BALL/MOLMEC/MINIMIZATION/conjugateGradient.h>
#include <BALL/MOLMEC/MINIMIZATION/lineSearch.h>
#include <BALL/COMMON/limits.h>
#include <BALL/MOLMEC/COMMON/forceField.h>

//#define BALL_DEBUG
#undef BALL_DEBUG

// The default method to use for the CG direction update
// (FLETCHER_REEVES | POLAK_RIBIERE | SHANNO)
#define DEFAULT_METHOD SHANNO

namespace BALL 
{
      using namespace std; 
      
      // Set the default values for options of this class 
      const char* ConjugateGradientMinimizer::Option::UPDATE_METHOD = "update_method"; 
      
      const Size ConjugateGradientMinimizer::Default::UPDATE_METHOD = DEFAULT_METHOD; 
      
      // Default constructor 
      // It does nothing but calling its base class constructor 
00032   ConjugateGradientMinimizer::ConjugateGradientMinimizer()
            : EnergyMinimizer(),
                  line_search_(),
                  unscaled_direction_(),
                  number_of_atoms_(0),
                  updt_method_(DEFAULT_METHOD),
                  first_iter_(true),
                  old_gtg_(0.),
                  a_i_(),
                  b_i_(),
                  p_t_(),
                  y_t_(),
                  p_i_(),
                  y_i_(),
                  D_1_(0.),
                  D_4_(0.),
                  restart_frequency_(1),
                  last_restart_iter_(0)
      {
            line_search_.setMinimizer(*this);

            options.setDefaultInteger(Option::UPDATE_METHOD, Default::UPDATE_METHOD);
      }

      // Constructor initialized with a force field
00057       ConjugateGradientMinimizer::ConjugateGradientMinimizer(ForceField& force_field)
            : EnergyMinimizer(),
                  line_search_(),
                  unscaled_direction_(),
                  number_of_atoms_(0),
                  updt_method_(DEFAULT_METHOD),
                  first_iter_(true),
                  old_gtg_(0.),
                  a_i_(),
                  b_i_(),
                  p_t_(),
                  y_t_(),
                  p_i_(),
                  y_i_(),
                  D_1_(0.),
                  D_4_(0.),
                  restart_frequency_(1),
                  last_restart_iter_(0)
      {
            line_search_.setMinimizer(*this);

            options.setDefaultInteger(Option::UPDATE_METHOD, Default::UPDATE_METHOD);
            valid_ = setup(force_field);
            
            if (!valid_)
            {
                  Log.error() << "ConjugateGradientMinimizer: setup failed! " << endl;
            }
      }

      // Constructor initialized with a force field and a snapshot manager 
      ConjugateGradientMinimizer::ConjugateGradientMinimizer
00089             (ForceField& force_field, SnapShotManager* ssm)
            : EnergyMinimizer(),
                  line_search_(),
                  unscaled_direction_(),
                  number_of_atoms_(0),
                  updt_method_(DEFAULT_METHOD),
                  first_iter_(true),
                  old_gtg_(0.),
                  a_i_(),
                  b_i_(),
                  p_t_(),
                  y_t_(),
                  p_i_(),
                  y_i_(),
                  D_1_(0.),
                  D_4_(0.),
                  restart_frequency_(1),
                  last_restart_iter_(0)
      {
            line_search_.setMinimizer(*this);

            options.setDefaultInteger(Option::UPDATE_METHOD, Default::UPDATE_METHOD);
            valid_ = setup(force_field, ssm);
            
            if (!valid_)
            {
                  Log.error() << "ConjugateGradientMinimizer: setup failed! " << endl;
            }
      }

      
      // Constructor initialized with a force field and a set of options
      ConjugateGradientMinimizer::ConjugateGradientMinimizer
00122             (ForceField& force_field, const Options& new_options) 
            : EnergyMinimizer(),
                  line_search_(),
                  unscaled_direction_(),
                  number_of_atoms_(0),
                  updt_method_(DEFAULT_METHOD),
                  first_iter_(true),
                  old_gtg_(0.),
                  a_i_(),
                  b_i_(),
                  p_t_(),
                  y_t_(),
                  p_i_(),
                  y_i_(),
                  D_1_(0.),
                  D_4_(0.),
                  restart_frequency_(1),
                  last_restart_iter_(0)
      {
            line_search_.setMinimizer(*this);

            options.setDefaultInteger(Option::UPDATE_METHOD, Default::UPDATE_METHOD);
            
            // The actual work is done in setup 
            valid_ = setup(force_field, new_options);
            
            if (!valid_)
            {
                  Log.error() << " Setup of conjugate gradient minimizer has failed! " << endl; 
            }
      }
      
      // Constructor initialized with a force field, a snapshot manager, and a set of options
      ConjugateGradientMinimizer::ConjugateGradientMinimizer
00156             (ForceField& force_field, SnapShotManager* ssm, const Options& new_options)
            : EnergyMinimizer(),
                  line_search_(),
                  unscaled_direction_(),
                  number_of_atoms_(0),
                  updt_method_(DEFAULT_METHOD),
                  first_iter_(true),
                  old_gtg_(0.),
                  a_i_(),
                  b_i_(),
                  p_t_(),
                  y_t_(),
                  p_i_(),
                  y_i_(),
                  D_1_(0.),
                  D_4_(0.),
                  restart_frequency_(1),
                  last_restart_iter_(0)
      {
            line_search_.setMinimizer(*this);

            options.setDefaultInteger(Option::UPDATE_METHOD, Default::UPDATE_METHOD);
            
            // The actual work is done in setup 
            valid_ = setup(force_field, ssm, new_options);
            
            if (!valid_)
            {
                  Log.error() << " Setup of conjugate gradient minimizer has failed! " << endl; 
            }
      }
      
      // Destructor
00189       ConjugateGradientMinimizer::~ConjugateGradientMinimizer()
      {
      }
      
      // The copy constructor 
      ConjugateGradientMinimizer::ConjugateGradientMinimizer 
00195             (const ConjugateGradientMinimizer& rhs)
            : EnergyMinimizer(rhs),
                  line_search_(rhs.line_search_),
                  unscaled_direction_(rhs.unscaled_direction_),
                  number_of_atoms_(rhs.number_of_atoms_),
                  updt_method_(rhs.updt_method_),
                  first_iter_(rhs.first_iter_),
                  old_gtg_(rhs.old_gtg_),
                  a_i_(rhs.a_i_),
                  b_i_(rhs.b_i_),
                  p_t_(rhs.p_t_),
                  y_t_(rhs.y_t_),
                  p_i_(rhs.p_i_),
                  y_i_(rhs.y_i_),
                  D_1_(rhs.D_1_),
                  D_4_(rhs.D_4_),
                  restart_frequency_(rhs.restart_frequency_),
                  last_restart_iter_(rhs.last_restart_iter_)
      {
            line_search_.setMinimizer(*this);
      }
      
      // The assignment operator
      const ConjugateGradientMinimizer& ConjugateGradientMinimizer::operator =
00219             (const ConjugateGradientMinimizer& rhs)
      {
            EnergyMinimizer::operator = (rhs);
            line_search_        = rhs.line_search_;
            updt_method_        = rhs.updt_method_;
            unscaled_direction_ = rhs.unscaled_direction_;
            number_of_atoms_    = rhs.number_of_atoms_;
            first_iter_         = rhs.first_iter_;
            old_gtg_            = rhs.old_gtg_;
            a_i_                = rhs.a_i_;
            b_i_                = rhs.b_i_;
            p_t_                = rhs.p_t_;
            y_t_                = rhs.y_t_;
            p_i_                = rhs.p_i_;
            y_i_                = rhs.y_i_;
            D_1_                = rhs.D_1_;
            D_4_                = rhs.D_4_;
            restart_frequency_  = rhs.restart_frequency_;
            last_restart_iter_  = rhs.last_restart_iter_;
            line_search_.setMinimizer(*this);
            return *this;
      }
      
      // This method is responsible for doing the specific setup of this class       
00243       bool ConjugateGradientMinimizer::specificSetup()
      {
            // Make sure the force field is assigned and valid!
            if (force_field_ == 0 || !force_field_->isValid())
            {
                  return false;
            }
            
            // Invalidate the initial gradient in order to ensure
            // its re-evaluation at the start of minimize().
            initial_grad_.invalidate();
            old_grad_.invalidate();
            
            // Get the options
            updt_method_ = (Size) options.getInteger(Option::UPDATE_METHOD);
            
            // determine the number of atoms
            number_of_atoms_ = (Size)force_field_->getNumberOfMovableAtoms();
            
            // Reset the restart counter
            last_restart_iter_ = 0;
            
            return true;
      }
      
      // Set explicitly the update method
00269       void ConjugateGradientMinimizer::setUpdateMethod(UpdateMethod updt)
      {
            updt_method_ = updt;
            options.setInteger(Option::UPDATE_METHOD, (Size)updt);
      }
      
      // Return the update method
00276       ConjugateGradientMinimizer::UpdateMethod ConjugateGradientMinimizer::getUpdateMethod() const
      {
            return (UpdateMethod)updt_method_;
      }
      
      // This method determines the new search direction. Along this
      // vector we will try to find the next solution.
00283       void ConjugateGradientMinimizer::updateDirection()
      {
            if (!initial_grad_.isValid())
            {
                  // Compute the initial energy and the initial forces
                  updateEnergy();
                  updateForces();
                  
                  // Store them for later use
                  storeGradientEnergy();
            }
            
            if (first_iter_)
            {
                  // This is the first iteration. We use the
                  // normalized negative gradient as first direction
                  direction_ = initial_grad_;
                  direction_.negate();
                  if (updt_method_ != SHANNO)
                  {
                        old_gtg_ = 0.;
                        for(Size i = 0; i < number_of_atoms_; ++i)
                        {
                              // We do the calculation in 'double' accuracy
                              old_gtg_ += (double)initial_grad_[i].x*(double)initial_grad_[i].x;
                              old_gtg_ += (double)initial_grad_[i].y*(double)initial_grad_[i].y;
                              old_gtg_ += (double)initial_grad_[i].z*(double)initial_grad_[i].z;
                        }
                        unscaled_direction_ = direction_;
                  }
                  direction_.normalize();
                  first_iter_        = false;
                  last_restart_iter_ = 0;
                  return;
            }
            
            // Just in case: check whether
            //  (1) (enough) memory is allocated for our search direction
            //  (2) the old gradient are invalid.
            if ((direction_.size() == 0) || !old_grad_.isValid())
            {
                  // This mustn't happen at all. If this routine is called all data(-structures) should be fine.
                  // But we are here! Something went wrong and we don't know what. Most kindly the data is
                  // modified by an external program/thread, so we do the best we can do: first we update the forces
                  // and set the direction to the negative gradient. Additionally, we force a restart.
                  Log.error() << "dir: " << direction_.isValid()
                                                      << "  initial_grad: " << initial_grad_.isValid()
                                                      << "  old_grad: " << old_grad_.isValid() << endl;
                  Log.error() << "ConjugateGradient::updateDirection: invalid gradient or direction  -  cannot use CG." << endl;
                  
                  // Calculate the current gradient and determine the current
                  // direction as the direction of steepest descent.
                  updateForces();
                  direction_ = current_grad_;
                  direction_.negate();
                  unscaled_direction_ = direction_;
                  direction_.normalize();
                  last_restart_iter_ = 0;
                  return;
            }
            
            /* There are three different conjugate gradient methods implemented:
                  
                        - Fletcher-Reeves (FR)
                        - Polak-Ribiere (PR)
                        - Shanno
                        
                   For FR and PR methods, the new direction vector is calculated as follows:
                        
                              d_k = - g_k + \beta d_{k-1}
                        
                   Where d_k and d_{k-1} are the new and the old search direction 
                   and g_k is the current gradient. \beta is determined by one of 
                   methods stated above.
            */
            switch (updt_method_)
            {
                  case FLETCHER_REEVES:
                  {
                        // We compute the new direction by the Fletcher-Reeves update.
                        
                        // \beta = \frac{<g_{k}, g_{k}>}{<g_{k-1}, g_{k-1}>}, see [1].
                        
                        double gtg = 0.;
                        
                        // In order to use 'double' accuracy the calculation is done
                        // 'by hand'.
                        for(Size i = 0; i < number_of_atoms_; ++i)
                        {
                              gtg += (double)initial_grad_[i].x*(double)initial_grad_[i].x;
                              gtg += (double)initial_grad_[i].y*(double)initial_grad_[i].y;
                              gtg += (double)initial_grad_[i].z*(double)initial_grad_[i].z;
                        }
                        
                        // Check whether we should proceed by a restart
                        if ((last_restart_iter_ == restart_frequency_) || (same_energy_counter_ > 0))
                        {
                              direction_ = initial_grad_;
                              direction_.negate();
                              unscaled_direction_ = direction_;
                              direction_.normalize();
                              old_gtg_ = gtg;
                              last_restart_iter_ = 0;
                              return;
                        }
                        
                        double beta = gtg;
                        if (old_gtg_ > cutlo_)
                        {
                              beta    /= old_gtg_;
                              old_gtg_ = gtg;
                        }
                        else
                        {
                              // Something went wrong. The recent gradient is about 0 or
                              // is something wrong with the stored value???
                              // We force a restart.
                              direction_ = initial_grad_;
                              direction_.negate();
                              unscaled_direction_ = direction_;
                              direction_.normalize();
                              old_gtg_ = gtg;
                              last_restart_iter_ = 0;
                              return;
                        }
                        
                        // Calculate the new conjugate gradient search direction:
                        // direction_ = - initial_gradient_ + beta * direction_;
                        // and a simple angle test for restart, see [5]
                        direction_.norm = 0.0;
                        double stg = 0.;
                        for (Size i = 0; i < number_of_atoms_; ++i)
                        {
                              // Same case as above, use 'double' accuracy
                              direction_[i].x  = -initial_grad_[i].x + beta*unscaled_direction_[i].x;
                              direction_[i].y  = -initial_grad_[i].y + beta*unscaled_direction_[i].y;
                              direction_[i].z  = -initial_grad_[i].z + beta*unscaled_direction_[i].z;
                              
                              // We don't need the following values computed in 'double' accuracy
                              direction_.norm += direction_[i] * direction_[i];
                              stg             += direction_[i]*initial_grad_[i];
                        }
                        direction_.norm = sqrt(direction_.norm);
                        if (-stg/direction_.norm/initial_grad_.norm < 1.e-3)
                        {
                              // Case: beta = 0, i.e. a restart
                              direction_ = initial_grad_;
                              direction_.negate();
                              unscaled_direction_ = direction_;
                              direction_.normalize();
                              last_restart_iter_ = 0;
                              return;
                        }
                        else
                        {
                              // Assign the norm of the new direction
                              direction_.rms = direction_.norm / (3.0 * (double)number_of_atoms_);
                        
                              // Don't risc a "NaN"
                              if (direction_.norm >= cutlo_)
                              {
                                    direction_.inv_norm = 1.0 / direction_.norm;
                              }
                              else
                              {
                                    direction_.inv_norm = sqrt(Limits<float>::max());
                              }
                              unscaled_direction_ = direction_;
                        }
                        break;
                  }
                  case POLAK_RIBIERE:
                  {
                        // We compute the new direction by the Polak-Ribiere update.
                        
                        // \beta = \frac{<g_k, (g_k - g_{k-1})>}{<g_{k-1}, g_{k-1}>}
                        double gtg = 0.;
                        double gg  = 0.;
                        for (Size i = 0; i < number_of_atoms_; i++)
                        {
                              // We want 'double' accuracy
                              gtg += (double)initial_grad_[i].x*(double)initial_grad_[i].x;
                              gg  += (double)initial_grad_[i].x*(double)old_grad_[i].x;
                              gtg += (double)initial_grad_[i].y*(double)initial_grad_[i].y;
                              gg  += (double)initial_grad_[i].y*(double)old_grad_[i].y;
                              gtg += (double)initial_grad_[i].z*(double)initial_grad_[i].z;
                              gg  += (double)initial_grad_[i].z*(double)old_grad_[i].z;
                        }
                        
                        // Check whether we should proceed by a restart.
                        // We also force a restart if the energy difference is too small
                        if ((last_restart_iter_ == restart_frequency_) || (same_energy_counter_ > 0))
                        {
                              direction_ = initial_grad_;
                              direction_.negate();
                              unscaled_direction_ = direction_;
                              direction_.normalize();
                              old_gtg_ = gtg;
                              last_restart_iter_ = 0;
                              return;
                        }
                        
                        double beta = gtg - gg;
                        if (old_gtg_ > cutlo_)
                        {
                              beta    /= old_gtg_;
                              old_gtg_ = gtg;
                        }
                        else
                        {
                              // Something went wrong. The recent gradient is about 0 or
                              // is something wrong with the stored value???
                              // We force a restart.
                              direction_ = initial_grad_;
                              direction_.negate();
                              unscaled_direction_ = direction_;
                              direction_.normalize();
                              old_gtg_ = gtg;
                              last_restart_iter_ = 0;
                              return;
                        }
                        
                        if (beta < 0.)
                        {
                              // This is the proposed criterion in [4]
                              // beta = 0.;
                              direction_ = initial_grad_;
                              direction_.negate();
                              unscaled_direction_ = direction_;
                              direction_.normalize();
                              last_restart_iter_ = 0;
                        }
                        else
                        {
                              // Calculate the new conjugate gradient search direction:
                              // direction_ = - initial_gradient_ + beta * direction_;
                              // and a smiple angle test for restart
                              direction_.norm = 0.0;
                              float stg = 0.;
                              for (Size i = 0; i < number_of_atoms_; i++)
                              {
                                    // Same case as above, use 'double' accuracy
                                    direction_[i].x  = -initial_grad_[i].x + beta*unscaled_direction_[i].x;
                                    direction_[i].y  = -initial_grad_[i].y + beta*unscaled_direction_[i].y;
                                    direction_[i].z  = -initial_grad_[i].z + beta*unscaled_direction_[i].z;
                                    
                                    // We don't need the following values computed in 'double' accuracy
                                    direction_.norm += direction_[i] * direction_[i];
                                    stg             += direction_[i]*initial_grad_[i];
                              }
                              // Assign the norm of the new direction
                              direction_.norm = sqrt(direction_.norm);
                              if (-stg/direction_.norm/initial_grad_.norm < 1.e-3)
                              {
                                    // beta = 0.;
                                    direction_ = initial_grad_;
                                    direction_.negate();
                                    unscaled_direction_ = direction_;
                                    direction_.normalize();
                                    last_restart_iter_ = 0;
                              }
                              else
                              {
                                    direction_.rms = direction_.norm / (3.0 * (double)number_of_atoms_);
                                    
                                    // Don't risc a "NaN"
                                    if (direction_.norm >= cutlo_)
                                    {
                                          direction_.inv_norm = 1.0 / direction_.norm;
                                    }
                                    else
                                    {
                                          direction_.inv_norm = sqrt(Limits<float>::max());
                                    }
                                    unscaled_direction_ = direction_;
                              }
                        }
                        break;
                  }
                  case SHANNO:
                  default:
                  {
                        /* Literature: see [2], [3]
                               The new search direction d_k is calculated as follows:
                                 d_i+1 = - b_i + D_8 / D_5 * a_i - (1 + D_6/D_5) * D_8/D_5 - D_7/D_5) * p_i
                               with
                               Index i = old iteration 
                               a_i = D_1 / D_4 * y_i - D_2 / D_4 * y_t + (2 D_2/D_1 - D_3/D_4) * p_t
                               b_i = D_1 / D_4 * g_i - C_1/D_4 * y_t + (2 C_1/D_1 - C_2/D_4) * p_t
                               p_i = step_length_ * d_old;
                               p_t = step_length_t * d_t
                               D_1 = p_t * y_t         D_2 = p_t * y_i
                               D_3 = y_t * y_i         C_1 = p_t * grad_new
                               C_2 = y_t * grad_new    D_4 = y_t * y_t
                               D_5 = p_i * y_i         D_6 = a_i * y_i 
                               D_7 = a_i * grad_new    D_8 = p_i * grad_new 
                               y_i = grad_new - grad_old 
                               y_t = grad_new - grad_old  in iteration t
                               d_t = search direction in iteration t
                        */
                        double sum1 = 0;
                        double sum2 = 0;
                        for (Size i = 0; i < number_of_atoms_; i++)
                        {
                              y_i_[i] = initial_grad_[i] - old_grad_[i];
                              p_i_[i] = direction_[i] * step_;
                              
                              sum1 += y_i_[i] * y_i_[i];
                              sum2 += p_i_[i] * p_i_[i];
                        }
                        
                        if (sum1 < cutlo_ || sum2 < cutlo_)
                        {
                              // Take the current gradient as the new search direction
                              direction_ = initial_grad_;
                              direction_.negate();
                              // Not necessary in Shanno case
                              //unscaled_direction_ = direction_;
                              direction_.normalize();
                              last_restart_iter_ = 0;
                              return;
                        }
                        
                        if ((last_restart_iter_ == restart_frequency_) || (number_of_iterations_ == 1))
                        {
                              double condition = 0;
                              for (Size i = 0; i < number_of_atoms_; i++)
                              {
                                    condition += initial_grad_[i] * old_grad_[i];
                              }
                              
                              // Take the absolute value
                              condition = fabs(condition);
                              
                              if ((number_of_iterations_ == 1) || (condition >= 0.2 * initial_grad_.norm * initial_grad_.norm))
                              {
                                    D_1_ = 0.;
                                    D_4_ = 0.;
                                    for (Size i = 0; i < number_of_atoms_; i++)
                                    {
                                          p_t_[i] = direction_[i]*step_;
                                          y_t_[i] = initial_grad_[i] - old_grad_[i];
                                          D_1_ += p_t_[i]*y_t_[i];
                                          D_4_ += y_t_[i]*y_t_[i];
                                    }
                                    
                                    if (fabs(D_4_) <= cutlo_)
                                    {
                                          throw Exception::DivisionByZero(__FILE__, __LINE__);
                                    }
                              }
                              last_restart_iter_ = 0;
                        }
                        
                        // Calculate all the auxiliary values
                        
                        double C_1 = 0.;
                        double C_2 = 0.;
                        double D_2 = 0.;
                        double D_3 = 0.;
                        for (Size i = 0; i < number_of_atoms_; i++)
                        {
                              C_1 += p_t_[i] * initial_grad_[i];
                              C_2 += y_t_[i] * initial_grad_[i];
                              D_2 += p_t_[i] * y_i_[i];
                              D_3 += y_t_[i] * y_i_[i];
                        }
                        
                        double factor1 = D_1_ / D_4_;
                        double factor2 = D_2 / D_4_;
                        double factor3 = (2. * D_2 / D_1_ - D_3 / D_4_);
                        
                        for (Size i = 0; i < number_of_atoms_; i++)
                        {
                              a_i_[i] = y_i_[i] * factor1 - y_t_[i] * factor2 + p_t_[i] * factor3;
                        }
                        
                        factor2 = C_1 / D_4_;
                        factor3 = (2. * C_1 / D_1_ - C_2 / D_4_);
                        
                        for (Size i = 0; i < number_of_atoms_; i++)
                        {
                              b_i_[i] = initial_grad_[i] * factor1 - y_t_[i] * factor2 + p_t_[i] * factor3;
                        }
                        
                        double D_5 = 0.;
                        double D_6 = 0.;
                        double D_7 = 0.;
                        double D_8 = 0.;
                        for (Size i = 0; i < number_of_atoms_; i++)
                        {
                              D_5 += p_i_[i] * y_i_[i];
                              D_6 += a_i_[i] * y_i_[i];
                              D_7 += a_i_[i] * initial_grad_[i];
                              D_8 += p_i_[i] * initial_grad_[i];
                        }
                        
                        if (fabs(D_5) <= cutlo_)
                        {
                              // Take the current gradient as the new search direction
                              direction_ = initial_grad_;
                              direction_.negate();
                              // Not necessary in Shanno case
                              //unscaled_direction_ = direction_;
                              direction_.normalize();
                              last_restart_iter_ = 0;
                              return;
                        }
                        
                        // Set the new search direction
                        direction_.norm = 0.0;
                        
                        factor1 = D_8 / D_5;
                        factor2 = (1. + D_6 / D_5) * factor1 - D_7 / D_5;
                        
                        for (Size i = 0; i < number_of_atoms_; i++)
                        {
                              direction_[i] = -b_i_[i] + a_i_[i] * factor1 - p_i_[i] * factor2;
                              direction_.norm += direction_[i] * direction_[i];
                        }
                        
                        // Assign the norm of the new direction
                        direction_.norm = sqrt(direction_.norm);
                        direction_.rms = direction_.norm / (3.0 * (double)number_of_atoms_);
                        
                        // Don't risc a "NaN"
                        if (direction_.norm >= cutlo_)
                        {
                              direction_.inv_norm = 1.0 / direction_.norm;
                        }
                        else
                        {
                              direction_.inv_norm = sqrt(Limits<float>::max());
                        }
                        // No necessary in Shanno case
                        //unscaled_direction_ = direction_;
                  }
            }
      } // end of method 'updateDirection'

      // The minimizer optimizes the energy of the system 
      // by using a conjugate gradient method. 
      // Return value is true when no further steps can be taken!
00726       bool ConjugateGradientMinimizer::minimize(Size iterations, bool resume)
      {
            aborted_ = false;
            
            // Check for validity of minimizer and force field
            if (!isValid() || getForceField() == 0 || !getForceField()->isValid())
            {
                  Log.error() << "ConjugateGradientMinimizer: is not initialized correctly!" << std::endl;
                  aborted_ = true;
                  return false;
            }
            
            Size noatoms = force_field_->getNumberOfMovableAtoms();
            // (1) Make sure we have something worth moving.
            if (noatoms == 0)
            {
                  return true;
            }
            // (2) Check that outside of this minimizer the number of movable atoms didn't change
            if (number_of_atoms_ != noatoms)
            {
                  // The number of movable atoms has changed, so we must start from scratch
                  number_of_atoms_ = noatoms;
                  resume = false;
            }
            // (3) Check that outside of this minimizer the options didn't change
            if ((Size)options.getInteger(Option::UPDATE_METHOD) != updt_method_)
            {
                  updt_method_ = (Size)options.getInteger(Option::UPDATE_METHOD);
                  resume = false;
            }
            
            // Define an alias for the atom vector
            AtomVector& atoms(const_cast<AtomVector&>(getForceField()->getAtoms()));
            
            // If we start from scratch (i.e. no restart) we have to make sure
            // to calculate all the quantities we need before we start
            if (!resume)
            {
                  // Reset the number of iterations if the job is not resumed.
                  setNumberOfIterations(0);
                  same_energy_counter_ = 0;
                  initial_grad_.invalidate();
                  current_grad_.invalidate();
                  first_iter_ = true;
                  last_restart_iter_ = 0;
                  
                  restart_frequency_ = 3 * number_of_atoms_;
                  if (updt_method_ == SHANNO)
                  {
                        p_t_.resize(number_of_atoms_);
                        y_t_.resize(number_of_atoms_);
                        y_i_.resize(number_of_atoms_);
                        p_i_.resize(number_of_atoms_);
                        a_i_.resize(number_of_atoms_);
                        b_i_.resize(number_of_atoms_);
                  }
                  
                  // Obviously, we don't have "old" energies yet, so we initialize it a with 
                  // sensible value. We don't need "old" gradients here.
                  old_energy_ = Limits<float>::max();
            }
            
            Size max_iterations = std::min(getNumberOfIterations() + iterations, getMaxNumberOfIterations());
            
            #ifdef BALL_DEBUG
                  Log.info() << "CGM: minimize(" << iterations << ", " << resume << ")" << endl;
            #endif
            
            // save the current atom positions
            atoms.savePositions();
            bool converged = false; 
            
            // Iterate: while not converged and not enough iterations
            while (!converged && (getNumberOfIterations() < max_iterations))  
            {
                  // Try to take a new step
                  double stp = findStep();
                  
                  // use this step as new reference step if findStep was successful
                  if (stp > 0.)
                  {
                        atoms.savePositions();
                  }
                  
                  // Store the gradient and the energy
                  old_grad_   = initial_grad_;
                  old_energy_ = initial_energy_;
                  
                  // Store the current gradient and energy
                  storeGradientEnergy();
                  
                  #ifdef BALL_DEBUG
                        Log << "CGM: end of main: current grad RMS = " << current_grad_.rms << std::endl;
                  #endif
                  
                  // Check for convergence.
                  converged = isConverged() || (stp == 0.);
                  
                  // Increment iteration counter, take snapshots, print energy,
                  // update pair lists, and check the same-energy counter
                  finishIteration();
                  ++last_restart_iter_;
                  
                  if ((!converged) && (stp < 0.))
                  {
                        // Nasty case: No convergence and the step computation failed.
                        // We must give up:-(
                        aborted_ = true;
                        return false;
                  }
                  if (Maths::isNan(force_field_->getEnergy())) 
                  {
                        aborted_ = true;
                        return false;
                  }
                  if (Maths::isNan(getGradient().rbegin()->x) ||
                      Maths::isNan(getGradient().rbegin()->y) ||
                      Maths::isNan(getGradient().rbegin()->z)) 
                  {
                        aborted_ = true;
                        return false;
                  }
                  if (abort_by_energy_enabled_)
                  {
                        if (force_field_->getEnergy() > abort_energy_) 
                        {
                              aborted_ = true;
                              return false;
                        }
                  }
            }
            
            return converged;
      } // end of method 'minimize' 
      
      
00863       double ConjugateGradientMinimizer::findStep()
      {
            // Compute the new direction
            updateDirection();
            
            // Define an alias for the atom vector
            AtomVector& atoms(const_cast<AtomVector&>(getForceField()->getAtoms()));
            
            // We perform a line search
            
            // No need to assure the maximum displacement here since our
            // line search pays attention to this constraint.
            
            bool result = line_search_.minimize(step_);
            if (!result)
            {
                  // Something went wrong.
                  // No success. We proceed with a restart.
                  
                  // Set the search direction to the normalized negative gradient. Since we proceed
                  // with a restart, we mustn't update the stored vectors by 'updateDirection' and there is
                  // no need to compute anything by 'updateDirection'.
                  
                  // Just in case: force field update (to update the pair list)
                  atoms.resetPositions();
                  force_field_->update();
                  
                  // Compute the initial energy and the initial forces
                  initial_energy_ = force_field_->updateEnergy();
                  force_field_->updateForces();
                  initial_grad_.set(force_field_->getAtoms());
                  
                  direction_ = initial_grad_;
                  direction_.negate();
                  if (updt_method_ != SHANNO)
                  {
                        unscaled_direction_ = direction_;
                  }
                  direction_.normalize();
                  last_restart_iter_ = 0;
                  
                  Size iter = 0;
                  while ((!result) && (iter < 12))
                  {
                        result = line_search_.minimize(step_);
                        if (!result)
                        {
                              for(Size i = 0; i <  number_of_atoms_; ++i)
                              {
                                    direction_[i] *= 0.5;
                              }
                              direction_.norm     *= 0.5;
                              direction_.rms      *= 0.5;
                              direction_.inv_norm *= 2.;
                              atoms.resetPositions();
                        }
                        else
                        {
                              return step_;
                        }
                        ++iter;
                  }
                  // If we are here something went wrong
                  // Not even such scaled steepest descent steps can manage
                  // the line search to exit successfully?
                  // We must be at a local minimizer...
                  step_ = 0.;
            }
            
            #ifdef BALL_DEBUG
                  Log.info() << "LineSearch: step = " << step << " result = " << result << endl;
            #endif
            
            return step_;
      } // end of method 'findStep'
      
} // end of namespace BALL

Generated by  Doxygen 1.6.0   Back to index