CGIonicStepper.C 2.83 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_alpha_start(1.0);
Francois Gygi committed
28
  cgopt_.set_alpha_max(5.0);
29
  cgopt_.set_beta_max(10.0);
30
#ifdef DEBUG
Francois Gygi committed
31 32
  if ( s.ctxt_.onpe0() )
    cgopt_.set_debug_print();
33
#endif
34 35 36 37 38
}
////////////////////////////////////////////////////////////////////////////////
void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
{
  // CG algorithm
39

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

  // 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
61

62
  cgopt_.compute_xp(x,e0,g,xp);
63

64 65 66 67 68 69 70
  // 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
71
  {
72
    if ( s_.ctxt_.onpe0() )
73 74
      cout << "  CGIonicStepper: displacement exceeds limit, rescaling" << endl;
    // rescale displacement and reset the CG optimizer
Francois Gygi committed
75 76
    double fac = max_disp/largest_disp;
    xp = x + fac * (xp - x);
77

Francois Gygi committed
78 79
    // reduce alpha starting value in CG optmizer
    cgopt_.set_alpha_start(fac*cgopt_.alpha_start());
80
    cgopt_.reset();
Francois Gygi committed
81
  }
82 83

  if ( s_.ctxt_.onpe0() )
84 85 86
  {
    cout << "  CGIonicStepper: alpha = " << cgopt_.alpha() << endl;
  }
87

88 89 90
  for ( int is = 0, i = 0; is < r0_.size(); is++ )
    for ( int j = 0; j < r0_[is].size(); j++ )
      rp_[is][j] = xp[i++];
91

Francois Gygi committed
92 93 94 95
  constraints_.enforce_r(r0_,rp_);
  rm_ = r0_;
  r0_ = rp_;
  atoms_.set_positions(r0_);
96
  atoms_.reset_velocities();
Francois Gygi committed
97
}