1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2015-2021 The plumed team
3 (see the PEOPLE file at the root of the distribution for a list of names)
4
5 See http://www.plumed.org for more information.
6
7 This file is part of plumed, version 2.
8
9 plumed is free software: you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13
14 plumed is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 #ifndef __PLUMED_tools_ConjugateGradient_h
23 #define __PLUMED_tools_ConjugateGradient_h
24
25 #include "MinimiseBase.h"
26
27 namespace PLMD {
28
29 template <class FCLASS>
30 class ConjugateGradient : public MinimiseBase<FCLASS> {
31 private:
32 /// This is the pointer to the member function in the energy
33 /// calculating class that calculates the energy
34 typedef double(FCLASS::*engf_pointer)( const std::vector<double>& p, std::vector<double>& der );
35 const unsigned ITMAX;
36 const double EPS;
37 public:
ConjugateGradient(FCLASS * funcc)38 explicit ConjugateGradient( FCLASS* funcc ) : MinimiseBase<FCLASS>(funcc), ITMAX(200), EPS(1E-10) {}
39 void minimise( const double& ftol, std::vector<double>& p, engf_pointer myfunc );
40 };
41
42 template <class FCLASS>
minimise(const double & ftol,std::vector<double> & p,engf_pointer myfunc)43 void ConjugateGradient<FCLASS>::minimise( const double& ftol, std::vector<double>& p, engf_pointer myfunc ) {
44 std::vector<double> xi( p.size() ), g( p.size() ), h( p.size() );
45 double fp = this->calcDerivatives( p, xi, myfunc );
46 for(unsigned j=0; j<p.size(); ++j) { g[j] = -xi[j]; xi[j]=h[j]=g[j]; }
47
48 for(unsigned its=0; its<ITMAX; ++its) {
49 double fret=this->linemin( xi, p, myfunc );
50 // The exit condition
51 if( 2.0*fabs(fret-fp) <= ftol*(fabs(fret)+fabs(fp)+EPS)) { return; }
52 fp = fret; this->calcDerivatives( p, xi, myfunc );
53 double ddg=0., gg=0.;
54 for(unsigned j=0; j<p.size(); ++j) { gg += g[j]*g[j]; ddg += (xi[j]+g[j])*xi[j]; }
55
56 if( gg==0.0 ) return;
57
58 double gam=ddg/gg;
59 for(unsigned j=0; j<p.size(); ++j) { g[j] = -xi[j]; xi[j]=h[j]=g[j]+gam*h[j]; }
60 }
61 plumed_merror("Too many interactions in conjugate gradient");
62 }
63
64 }
65
66 #endif
67