////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or .
//
////////////////////////////////////////////////////////////////////////////////
//
// CGIonicStepper.C
//
////////////////////////////////////////////////////////////////////////////////
#include "CGIonicStepper.h"
#include "CGOptimizer.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
CGIonicStepper::CGIonicStepper(Sample& s) : IonicStepper(s),
cgopt_(CGOptimizer(3*natoms_))
{
cgopt_.set_alpha_start(1.0);
cgopt_.set_alpha_max(50.0);
cgopt_.set_beta_max(10.0);
#ifdef DEBUG
if ( s.ctxt_.onpe0() )
cgopt_.set_debug_print();
#endif
}
////////////////////////////////////////////////////////////////////////////////
void CGIonicStepper::compute_r(double e0, const vector >& f0)
{
// CG algorithm
valarray x(3*natoms_),xp(3*natoms_),g(3*natoms_);
vector > gvec;
gvec.resize(r0_.size());
for ( int is = 0, i = 0; is < r0_.size(); is++ )
{
gvec[is].resize(r0_[is].size());
for ( int j = 0; j < r0_[is].size(); j++ )
{
x[i] = r0_[is][j];
gvec[is][j] = -f0[is][j];
i++;
}
}
// 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];
cgopt_.compute_xp(x,e0,g,xp);
// check largest displacement
// max_disp: largest acceptable displacement
const double max_disp = 0.2;
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 )
{
if ( s_.ctxt_.onpe0() )
cout << " CGIonicStepper: displacement exceeds limit, rescaling" << endl;
// rescale displacement and reset the CG optimizer
double fac = max_disp/largest_disp;
xp = x + fac * (xp - x);
cgopt_.set_alpha_start(fac*cgopt_.alpha_start());
cgopt_.reset();
}
if ( s_.ctxt_.onpe0() )
{
cout << " CGIonicStepper: alpha = " << cgopt_.alpha() << endl;
}
for ( int is = 0, i = 0; is < r0_.size(); is++ )
for ( int j = 0; j < r0_[is].size(); j++ )
rp_[is][j] = xp[i++];
constraints_.enforce_r(r0_,rp_);
rm_ = r0_;
r0_ = rp_;
atoms_.set_positions(r0_);
atoms_.reset_velocities();
}