SDAIonicStepper.C 3.5 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17 18 19 20 21 22
// SDAIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////

#include "SDAIonicStepper.h"
using namespace std;

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
23
void SDAIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
Francois Gygi committed
24
{
25 26
  double fp0;
  bool wolfe1, wolfe2;
27

28 29 30 31
  // 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++ )
Francois Gygi committed
32
  {
33 34 35 36
    for ( int i = 0; i < r0_[is].size(); i++ )
    {
      largest_force = max(largest_force,fabs(f0[is][i]));
    }
Francois Gygi committed
37
  }
38

39
  if ( largest_force > max_force )
Francois Gygi committed
40
  {
41 42 43 44 45 46
    if ( s_.ctxt_.onpe0() )
      cout << "  SDAIonicStepper: 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++ )
Francois Gygi committed
47
    {
48 49 50 51
      for ( int i = 0; i < r0_[is].size(); i++ )
      {
        rp_[is][i] = r0_[is][i] + alpha_sd * f0[is][i];
      }
Francois Gygi committed
52
    }
53 54 55 56 57 58 59
    constraints_.enforce_r(r0_,rp_);
    rm_ = r0_;
    r0_ = rp_;
    atoms_.set_positions(r0_);
    // reset the SDA algorithm
    first_step_ = true;
    return;
Francois Gygi committed
60
  }
61

62 63 64
  // SDA algorithm

  if ( !first_step_ )
Francois Gygi committed
65
  {
66 67 68 69 70 71 72 73 74 75 76
    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_);
Francois Gygi committed
77
    if ( s_.ctxt_.onpe0() )
78 79 80 81 82 83 84 85 86
    {
      cout << "  SDAIonicStepper: fpc = " << fpc_ << endl;
      cout << "  SDAIonicStepper: fp0 = " << fp0 << endl;
      cout << "  SDAIonicStepper: ec = " << ec_ << " e0 = " << e0 <<  endl;
      cout << "  SDAIonicStepper: ec_ + fpc_ * sigma1_ * alpha_ ="
           << ec_ + fpc_ * sigma1_ * alpha_ << endl;
      cout << "  SDAIonicStepper: wolfe1/wolfe2 = "
           << wolfe1 << "/" << wolfe2 << endl;
    }
Francois Gygi committed
87
  }
88

89
  if ( first_step_ || (wolfe1 && wolfe2) || linmin_.done() )
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
  {
    // set new descent direction
    // pc = f0
    fc_ = f0;
    pc_ = fc_;
    // 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();
  }

111
  alpha_ = linmin_.next_alpha(alpha_,e0,fp0);
Francois Gygi committed
112

113 114 115 116
  if ( s_.ctxt_.onpe0() )
    cout << "  SDAIonicStepper: alpha = " << alpha_ << endl;

  // rp = rc + alpha * pc
Francois Gygi committed
117 118 119 120
  for ( int is = 0; is < r0_.size(); is++ )
  {
    for ( int i = 0; i < r0_[is].size(); i++ )
    {
121
      rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i];
Francois Gygi committed
122 123
    }
  }
124

Francois Gygi committed
125
  constraints_.enforce_r(r0_,rp_);
Francois Gygi committed
126 127 128
  rm_ = r0_;
  r0_ = rp_;
  atoms_.set_positions(r0_);
129 130

  first_step_ = false;
Francois Gygi committed
131
}
132