1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2016-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_gridtools_GridSearch_h
23 #define __PLUMED_gridtools_GridSearch_h
24
25 #include "tools/MinimiseBase.h"
26 #include "GridVessel.h"
27 #include <iostream>
28 #include <memory>
29
30 namespace PLMD {
31 namespace gridtools {
32
33 template <class FCLASS>
34 class GridSearch {
35 private:
36 /// This is the pointer to the member function in the energy
37 /// calculating class that calculates the energy
38 typedef double(FCLASS::*engf_pointer)( const std::vector<double>& p, std::vector<double>& der );
39 FCLASS* myclass_func;
40 std::unique_ptr<GridVessel> mygrid;
41 std::unique_ptr<GridVessel> myfgrid;
42 public:
GridSearch(const std::vector<double> & mmin,const std::vector<double> & mmax,const std::vector<unsigned> & ng,const std::vector<unsigned> & nfg,FCLASS * funcc)43 GridSearch( const std::vector<double>& mmin, const std::vector<double>& mmax, const std::vector<unsigned>& ng, const std::vector<unsigned>& nfg, FCLASS* funcc ) :
44 myclass_func( funcc )
45 {
46 // Create the grid objects
47 std::string nstr, vstring="COMPONENTS=func COORDINATES=x1";
48 for(unsigned i=1; i<mmin.size(); ++i) { Tools::convert(i+1,nstr); vstring += ",x" + nstr; }
49 vstring += " PBC=F"; for(unsigned i=1; i<mmin.size(); ++i) vstring += ",F";
50 vesselbase::VesselOptions da("mygrid","",-1,vstring,NULL);
51 Keywords keys; gridtools::GridVessel::registerKeywords( keys );
52 vesselbase::VesselOptions dar( da, keys );
53 mygrid.reset( new GridVessel(dar) );
54 if( nfg[0]>0 ) myfgrid.reset( new GridVessel(dar) );
55
56 // Now setup the min and max values for the grid
57 std::vector<std::string> gmin( nfg.size() ), gmax( nfg.size() ); std::vector<double> dummy_spacing;
58 for(unsigned i=0; i<nfg.size(); ++i) { Tools::convert(mmin[i],gmin[i]); Tools::convert(mmax[i],gmax[i]); }
59 mygrid->setBounds( gmin, gmax, ng, dummy_spacing ); mygrid->resize();
60 if( myfgrid ) myfgrid->setBounds( gmin, gmax, nfg, dummy_spacing );
61 }
62 bool minimise( std::vector<double>& p, engf_pointer myfunc );
63 };
64
65 template <class FCLASS>
minimise(std::vector<double> & p,engf_pointer myfunc)66 bool GridSearch<FCLASS>::minimise( std::vector<double>& p, engf_pointer myfunc ) {
67 std::vector<double> der( p.size() ); std::vector<double> coords( p.size() );
68 double initial_eng = (myclass_func->*myfunc)( p, der );
69 mygrid->getGridPointCoordinates( 0, coords );
70 double emin=(myclass_func->*myfunc)( coords, der );
71 mygrid->setValueAndDerivatives( 0, 0, emin, der ); unsigned pmin=0;
72 for(unsigned i=1; i<mygrid->getNumberOfPoints(); ++i) {
73 mygrid->getGridPointCoordinates( i, coords );
74 double eng = (myclass_func->*myfunc)( coords, der );
75 mygrid->setValueAndDerivatives( i, 0, eng, der );
76 if( eng<emin ) { emin=eng; pmin=i; }
77 }
78 // This prevents division by zero
79 mygrid->setNorm( 1.0 );
80
81 if( myfgrid ) {
82 myfgrid->getGridPointCoordinates( 0, coords ); pmin=0;
83 double emin=mygrid->getValueAndDerivatives( coords, 0, der );
84 for(unsigned i=1; i<myfgrid->getNumberOfPoints(); ++i) {
85 myfgrid->getGridPointCoordinates( i, coords );
86 double eng = mygrid->getValueAndDerivatives( coords, 0, der );
87 if( eng<emin ) { emin=eng; pmin=i; }
88 }
89 myfgrid->getGridPointCoordinates( pmin, coords );
90 double checkEng = (myclass_func->*myfunc)( coords, der );
91 if( checkEng<initial_eng ) {
92 myfgrid->getGridPointCoordinates( pmin, p );
93 return true;
94 } else {
95 return false;
96 }
97 }
98
99 if( emin<initial_eng ) {
100 mygrid->getGridPointCoordinates( pmin, p );
101 return true;
102 } else {
103 return false;
104 }
105 }
106
107 }
108 }
109 #endif
110
111