CGIonicStepper.C 2.62 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
// CGIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////

#include "CGIonicStepper.h"
20
#include "CGOptimizer.h"
Francois Gygi committed
21 22 23
using namespace std;

////////////////////////////////////////////////////////////////////////////////
24
CGIonicStepper::CGIonicStepper(Sample& s) : IonicStepper(s),
25
  cgopt_(CGOptimizer(3*natoms_))
Francois Gygi committed
26
{
27
  cgopt_.set_beta_max(5.0);
28 29 30
#ifdef DEBUG
  cgopt_.set_debug_print();
#endif
31 32 33 34 35
}
////////////////////////////////////////////////////////////////////////////////
void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
{
  // CG algorithm
36

37
  valarray<double> x(3*natoms_),xp(3*natoms_),g(3*natoms_);
38 39
  vector<vector<double> > gvec;
  gvec.resize(r0_.size());
40
  for ( int is = 0, i = 0; is < r0_.size(); is++ )
41 42
  {
    gvec[is].resize(r0_[is].size());
43
    for ( int j = 0; j < r0_[is].size(); j++ )
Francois Gygi committed
44
    {
45
      x[i] = r0_[is][j];
46
      gvec[is][j] = -f0[is][j];
47
      i++;
Francois Gygi committed
48
    }
49 50 51 52 53 54 55 56 57
  }

  // enforce compatibility of the gradient with constraints
  constraints_.enforce_v(r0_,gvec);

  // copy projected gradient to g
  for ( int is = 0, i = 0; is < r0_.size(); is++ )
    for ( int j = 0; j < r0_[is].size(); j++ )
      g[i++] = gvec[is][j];
Francois Gygi committed
58

59
  cgopt_.compute_xp(x,e0,g,xp);
60

61 62 63 64 65 66 67
  // check largest displacement
  // max_disp: largest acceptable displacement
  const double max_disp = 0.05;
  double largest_disp = 0.0;
  for ( int i = 0; i < xp.size(); i++ )
    largest_disp = max(largest_disp,fabs(xp[i]-x[i]));
  if ( largest_disp > max_disp )
Francois Gygi committed
68
  {
69
    if ( s_.ctxt_.onpe0() )
70 71 72
      cout << "  CGIonicStepper: displacement exceeds limit, rescaling" << endl;
    // rescale displacement and reset the CG optimizer
    xp = x + (max_disp/largest_disp) * (xp - x);
73

74
    cgopt_.reset();
Francois Gygi committed
75
  }
76 77

  if ( s_.ctxt_.onpe0() )
78 79 80
  {
    cout << "  CGIonicStepper: alpha = " << cgopt_.alpha() << endl;
  }
81

82 83 84
  for ( int is = 0, i = 0; is < r0_.size(); is++ )
    for ( int j = 0; j < r0_[is].size(); j++ )
      rp_[is][j] = xp[i++];
85

Francois Gygi committed
86 87 88
  constraints_.enforce_r(r0_,rp_);
  rm_ = r0_;
  r0_ = rp_;
89
  atoms_.sync_positions(r0_);
Francois Gygi committed
90 91
  atoms_.set_positions(r0_);
}