1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // SDAIonicStepper.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18
19 #include "MPIdata.h"
20 #include "SDAIonicStepper.h"
21 using namespace std;
22
23 ////////////////////////////////////////////////////////////////////////////////
compute_r(double e0,const vector<vector<double>> & f0)24 void SDAIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
25 {
26 double fp0;
27 bool wolfe1, wolfe2;
28
29 // check if the largest component of the force f0 is larger than max_force
30 const double max_force = 0.1;
31 double largest_force = 0.0;
32 for ( int is = 0; is < r0_.size(); is++ )
33 {
34 for ( int i = 0; i < r0_[is].size(); i++ )
35 {
36 largest_force = max(largest_force,fabs(f0[is][i]));
37 }
38 }
39
40 if ( largest_force > max_force )
41 {
42 if ( MPIdata::onpe0() )
43 cout << " SDAIonicStepper: force exceeds limit, taking SD step " << endl;
44 // take a steepest descent step with limited displacement and exit
45 const double alpha_sd = max_force/largest_force;
46 // SD step
47 for ( int is = 0; is < r0_.size(); is++ )
48 {
49 for ( int i = 0; i < r0_[is].size(); i++ )
50 {
51 rp_[is][i] = r0_[is][i] + alpha_sd * f0[is][i];
52 }
53 }
54 constraints_.enforce_r(r0_,rp_);
55 rm_ = r0_;
56 r0_ = rp_;
57 atoms_.set_positions(r0_);
58 // reset the SDA algorithm
59 first_step_ = true;
60 return;
61 }
62
63 // SDA algorithm
64
65 if ( !first_step_ )
66 {
67 wolfe1 = e0 < ec_ + fpc_ * sigma1_ * alpha_;
68 // fp0 = -proj(f0,pc)
69 fp0 = 0.0;
70 for ( int is = 0; is < r0_.size(); is++ )
71 {
72 for ( int i = 0; i < r0_[is].size(); i++ )
73 {
74 fp0 -= f0[is][i] * pc_[is][i];
75 }
76 }
77 wolfe2 = fabs(fp0) < sigma2_ * fabs(fpc_);
78 if ( MPIdata::onpe0() )
79 {
80 cout << " SDAIonicStepper: fpc = " << fpc_ << endl;
81 cout << " SDAIonicStepper: fp0 = " << fp0 << endl;
82 cout << " SDAIonicStepper: ec = " << ec_ << " e0 = " << e0 << endl;
83 cout << " SDAIonicStepper: ec_ + fpc_ * sigma1_ * alpha_ ="
84 << ec_ + fpc_ * sigma1_ * alpha_ << endl;
85 cout << " SDAIonicStepper: wolfe1/wolfe2 = "
86 << wolfe1 << "/" << wolfe2 << endl;
87 }
88 }
89
90 if ( first_step_ || (wolfe1 && wolfe2) || linmin_.done() )
91 {
92 // set new descent direction
93 // pc = f0
94 fc_ = f0;
95 pc_ = fc_;
96 // fpc = d_e / d_alpha in direction pc
97 fpc_ = 0.0;
98 for ( int is = 0; is < r0_.size(); is++ )
99 {
100 for ( int i = 0; i < r0_[is].size(); i++ )
101 {
102 fpc_ -= fc_[is][i] * pc_[is][i];
103 }
104 }
105 ec_ = e0;
106 rc_ = r0_;
107 fp0 = fpc_;
108 // reset line minimizer
109 linmin_.reset();
110 }
111
112 alpha_ = linmin_.next_alpha(alpha_,e0,fp0);
113
114 if ( MPIdata::onpe0() )
115 cout << " SDAIonicStepper: alpha = " << alpha_ << endl;
116
117 // rp = rc + alpha * pc
118 for ( int is = 0; is < r0_.size(); is++ )
119 {
120 for ( int i = 0; i < r0_[is].size(); i++ )
121 {
122 rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i];
123 }
124 }
125
126 constraints_.enforce_r(r0_,rp_);
127 rm_ = r0_;
128 r0_ = rp_;
129 atoms_.set_positions(r0_);
130
131 first_step_ = false;
132 }
133
134