Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 : #include "ActionWithGrid.h"
23 : #include "core/ActionRegister.h"
24 : #include "lepton/Lepton.h"
25 : #include "tools/IFile.h"
26 :
27 : //+PLUMEDOC GRIDCALC REFERENCE_GRID
28 : /*
29 : Setup a constant grid by either reading values from a file or definining a function in input
30 :
31 : \par Examples
32 :
33 : */
34 : //+ENDPLUMEDOC
35 :
36 : namespace PLMD {
37 : namespace gridtools {
38 :
39 : static std::map<std::string, double> leptonConstants= {
40 : {"e", std::exp(1.0)},
41 : {"log2e", 1.0/std::log(2.0)},
42 : {"log10e", 1.0/std::log(10.0)},
43 : {"ln2", std::log(2.0)},
44 : {"ln10", std::log(10.0)},
45 : {"pi", pi},
46 : {"pi_2", pi*0.5},
47 : {"pi_4", pi*0.25},
48 : // {"1_pi", 1.0/pi},
49 : // {"2_pi", 2.0/pi},
50 : // {"2_sqrtpi", 2.0/std::sqrt(pi)},
51 : {"sqrt2", std::sqrt(2.0)},
52 : {"sqrt1_2", std::sqrt(0.5)}
53 : };
54 :
55 : class ReadGridInSetup : public ActionWithGrid {
56 : private:
57 : GridCoordinatesObject gridobject;
58 : std::vector<std::string> dernames;
59 : void createGridAndValue( const std::string& gtype, const std::vector<bool>& ipbc, const unsigned& nfermi,
60 : const std::vector<std::string>& gmin, const std::vector<std::string>& gmax,
61 : const std::vector<unsigned>& gbin );
62 : public:
63 : static void registerKeywords( Keywords& keys );
64 : explicit ReadGridInSetup(const ActionOptions&ao);
65 : unsigned getNumberOfDerivatives() override ;
66 0 : void setupOnFirstStep( const bool incalc ) override { plumed_merror("should not be in ReadGridInSetup setupOnFirstStep" ); }
67 0 : void performTask( const unsigned& current, MultiValue& myvals ) const override { plumed_merror("should not be in readGridInSetup performTask"); }
68 : std::vector<std::string> getGridCoordinateNames() const override ;
69 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
70 : };
71 :
72 : PLUMED_REGISTER_ACTION(ReadGridInSetup,"REFERENCE_GRID")
73 :
74 21 : void ReadGridInSetup::registerKeywords( Keywords& keys ) {
75 21 : ActionWithGrid::registerKeywords(keys); keys.remove("SERIAL");
76 42 : keys.add("optional","FUNC","the function to compute on the grid");
77 42 : keys.add("compulsory","GRID_MIN","auto","the lower bounds for the grid");
78 42 : keys.add("compulsory","GRID_MAX","auto","the upper bounds for the grid");
79 42 : keys.add("compulsory","PERIODIC","are the grid directions periodic");
80 42 : keys.add("optional","GRID_BIN","the number of bins for the grid");
81 42 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
82 42 : keys.add("optional","VAR","the names to give each of the grid directions in the function. If you have up to three grid coordinates in your function you can use x, y and z to refer to them. Otherwise you must use this flag to give your variables names.");
83 42 : keys.add("compulsory","FILE","the name of the file that contains the reference data");
84 42 : keys.add("compulsory","VALUE","the name of the value that should be read from the grid");
85 21 : keys.setValueDescription("the constant function on the grid that was specified in input");
86 21 : }
87 :
88 11 : ReadGridInSetup::ReadGridInSetup(const ActionOptions&ao):
89 : Action(ao),
90 11 : ActionWithGrid(ao)
91 : {
92 22 : std::string func; parse("FUNC",func);
93 11 : if( func.length()>0 ) {
94 : // Read in stuff for grid
95 10 : std::vector<std::string> gmin; parseVector("GRID_MIN",gmin);
96 10 : std::vector<std::string> gmax(gmin.size()); parseVector("GRID_MAX",gmax);
97 10 : std::vector<unsigned> gbin(gmin.size()); parseVector("GRID_BIN",gbin);
98 10 : std::vector<std::string> pbc(gmin.size()); parseVector("PERIODIC",pbc);
99 5 : std::vector<bool> ipbc( pbc.size() );
100 10 : for(unsigned i=0; i<ipbc.size(); ++i) {
101 5 : if( pbc[i]=="YES" ) ipbc[i]=true;
102 5 : else if( pbc[i]=="NO" ) ipbc[i]=false;
103 0 : else error( pbc[i] + " is not a valid instruction to the PERIODIC keyword");
104 : }
105 :
106 : // Read in the variables
107 10 : parseVector("VAR",dernames);
108 5 : if(dernames.size()==0) {
109 3 : dernames.resize(gmin.size());
110 3 : if(gmin.size()>3)
111 0 : error("Using more than 3 arguments you should explicitly write their names with VAR");
112 3 : if(dernames.size()>0) dernames[0]="x";
113 3 : if(dernames.size()>1) dernames[1]="y";
114 3 : if(dernames.size()>2) dernames[2]="z";
115 : }
116 5 : if(dernames.size()!=gmin.size()) error("Size of VAR array should be the same as number of grid dimensions");
117 :
118 : // Create the grid and the value of the grid
119 5 : createGridAndValue( "flat", ipbc, 0, gmin, gmax, gbin );
120 :
121 : // Read in stuff for function
122 5 : log.printf(" evaluating function : %s\n",func.c_str());
123 5 : log.printf(" with variables :");
124 10 : for(unsigned i=0; i<dernames.size(); i++) log.printf(" %s",dernames[i].c_str());
125 5 : log.printf("\n");
126 5 : log.printf(" on %d", gbin[0]);
127 5 : for(unsigned i=1; i<gbin.size(); ++i) log.printf(" by %d \n", gbin[i]);
128 5 : log.printf(" grid of points between (%s", gmin[0].c_str() );
129 5 : for(unsigned i=1; i<gmin.size(); ++i) log.printf(", %s", gmin[i].c_str() );
130 5 : log.printf(") and (%s", gmax[0].c_str() );
131 5 : for(unsigned i=1; i<gmax.size(); ++i) log.printf(", %s", gmax[i].c_str() );
132 5 : log.printf(")\n");
133 :
134 5 : lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(leptonConstants);
135 5 : log<<" function as parsed by lepton: "<<pe<<"\n";
136 5 : lepton::CompiledExpression expression=pe.createCompiledExpression();
137 10 : for(auto &p: expression.getVariables()) {
138 5 : if(std::find(dernames.begin(),dernames.end(),p)==dernames.end()) {
139 0 : error("variable " + p + " is not defined");
140 : }
141 : }
142 5 : log<<" derivatives as computed by lepton:\n";
143 5 : std::vector<lepton::CompiledExpression> expression_deriv( dernames.size() );
144 10 : for(unsigned i=0; i<dernames.size(); i++) {
145 10 : lepton::ParsedExpression pe=lepton::Parser::parse(func).differentiate(dernames[i]).optimize(leptonConstants);
146 5 : log<<" "<<pe<<"\n";
147 5 : expression_deriv[i]=pe.createCompiledExpression();
148 : }
149 : // And finally calculate all the grid points
150 5 : std::vector<double> dder( dernames.size() ), xx( dernames.size() ); Value* valout=getPntrToComponent(0);
151 115 : for(unsigned index=0; index<valout->getNumberOfValues(); ++index) {
152 110 : gridobject.getGridPointCoordinates( index, xx );
153 220 : for(unsigned j=0; j<xx.size(); ++j) {
154 : try {
155 110 : expression.getVariableReference(dernames[j])=xx[j];
156 0 : } catch(PLMD::lepton::Exception& exc) {
157 : // this is necessary since in some cases lepton things a variable is not present even though it is present
158 : // e.g. func=0*x
159 0 : }
160 : }
161 110 : valout->set( index, expression.evaluate() );
162 220 : for(unsigned k=0; k<xx.size(); ++k) {
163 220 : for(unsigned j=0; j<xx.size(); ++j) {
164 : try {
165 110 : expression_deriv[k].getVariableReference(dernames[j])=xx[j];
166 0 : } catch(PLMD::lepton::Exception& exc) {
167 : // this is necessary since in some cases lepton things a variable is not present even though it is present
168 : // e.g. func=0*x
169 0 : }
170 : }
171 110 : valout->addGridDerivatives( index, k, expression_deriv[k].evaluate() );
172 : }
173 : }
174 15 : } else {
175 12 : std::string valuestr; parse("VALUE",valuestr);
176 12 : std::string tstyle, filen; parse("FILE",filen);
177 6 : if( filen.length()>0 ) {
178 6 : std::size_t dot=filen.find_first_of(".");
179 12 : if( dot!=std::string::npos ) tstyle=filen.substr(dot+1);
180 6 : if( tstyle!="grid" ) error("can only read in grids using read value in setup");
181 6 : log.printf(" reading function %s on grid from file %s \n", valuestr.c_str(), filen.c_str() );
182 : }
183 6 : IFile ifile; ifile.open(filen);
184 6 : if( !ifile.FieldExist( valuestr ) ) error("could not find grid value in input file");
185 6 : std::vector<std::string> fieldnames; ifile.scanFieldList( fieldnames );
186 :
187 : // Retrieve the names of the variables the grid is computed over
188 : bool flatgrid=false;
189 54 : for(unsigned i=0; i<fieldnames.size(); ++i) {
190 48 : if( fieldnames[i].find("min_")!=std::string::npos ) flatgrid=true;
191 96 : std::size_t dot = fieldnames[i].find_first_of("d" + valuestr + "_" );
192 114 : if( fieldnames[i].find("d" + valuestr + "_")!=std::string::npos ) dernames.push_back( fieldnames[i].substr(dot+2+valuestr.length()) );
193 : }
194 6 : if( flatgrid && dernames.size()==0 ) error("could not find any derivatives for value " + valuestr + " in input file. Header should contain at least columns with a name starting d" + valuestr + "_");
195 : // Now get all the header data for the grid
196 6 : std::vector<std::string> gmin( dernames.size() ), gmax( dernames.size() ); std::string pstring;
197 6 : int gbin1; std::vector<unsigned> gbin( dernames.size() ); std::vector<bool> ipbc( dernames.size() );
198 6 : if( !flatgrid ) {
199 12 : ifile.scanField( "nbins", gbin1); gbin[0]=gbin1;
200 12 : createGridAndValue( "fibonacci", ipbc, gbin[0], gmin, gmax, gbin );
201 : } else {
202 0 : for(unsigned i=0; i<dernames.size(); ++i) {
203 0 : ifile.scanField( "min_" + dernames[i], gmin[i]);
204 0 : ifile.scanField( "max_" + dernames[i], gmax[i]);
205 0 : ifile.scanField( "periodic_" + dernames[i], pstring );
206 0 : ifile.scanField( "nbins_" + dernames[i], gbin1); gbin[i]=gbin1;
207 0 : if( pstring=="true" ) {
208 0 : log.printf(" for periodic coordinate %s minimum is %s maximum is %s and number of bins is %d \n",dernames[i].c_str(),gmin[i].c_str(),gmax[i].c_str(),gbin[i]);
209 : ipbc[i]=true;
210 0 : } else if( pstring=="false" ) {
211 0 : log.printf(" for coordinate %s minimum is %s maximum is %s and number of bins is %d \n",dernames[i].c_str(),gmin[i].c_str(),gmax[i].c_str(),gbin[i]);
212 : ipbc[i]=false;
213 0 : } else error("do not understand periodidicy of " + dernames[i] );
214 :
215 0 : bool hasder=ifile.FieldExist( "d" + valuestr + "_" + dernames[i] );
216 0 : if( !hasder ) plumed_merror("missing derivatives from grid file");
217 : }
218 0 : createGridAndValue( "flat", ipbc, 0, gmin, gmax, gbin );
219 : }
220 : // And finally read all the grid points
221 6 : Value* valout=getPntrToComponent(0);
222 6 : std::vector<double> dder( dernames.size() ), xx( dernames.size() );
223 2106 : for(unsigned i=0; i<valout->getNumberOfValues(); ++i) {
224 2100 : double x, val; ifile.scanField( valuestr, val );
225 8400 : for(unsigned j=0; j<dernames.size(); ++j) {
226 6300 : ifile.scanField(dernames[j],x);
227 6300 : if( !flatgrid ) {
228 12600 : ifile.scanField("nbins", gbin1);
229 : } else {
230 0 : xx[j]=x+gridobject.getGridSpacing()[j]/2.0;
231 0 : ifile.scanField( "min_" + dernames[j], gmin[j]);
232 0 : ifile.scanField( "max_" + dernames[j], gmax[j]);
233 0 : ifile.scanField( "nbins_" + dernames[j], gbin1);
234 0 : ifile.scanField( "periodic_" + dernames[j], pstring );
235 : }
236 : }
237 8400 : for(unsigned j=0; j<dernames.size(); ++j) ifile.scanField( "d" + valuestr + "_" + dernames[j], dder[j] );
238 :
239 2100 : unsigned index=gridobject.getIndex(xx); if( !flatgrid ) index=i;
240 2100 : valout->set( index, val );
241 8400 : for(unsigned j=0; j<dernames.size(); ++j) valout->addGridDerivatives( index, j, dder[j] );
242 2100 : ifile.scanField();
243 : }
244 6 : ifile.close();
245 6 : }
246 11 : }
247 :
248 11 : void ReadGridInSetup::createGridAndValue( const std::string& gtype, const std::vector<bool>& ipbc, const unsigned& nfermi,
249 : const std::vector<std::string>& gmin, const std::vector<std::string>& gmax,
250 : const std::vector<unsigned>& gbin ) {
251 11 : gridobject.setup( gtype, ipbc, nfermi, 0.0 ); std::vector<double> gspacing;
252 11 : if( gtype=="flat" ) {
253 5 : gridobject.setBounds( gmin, gmax, gbin, gspacing );
254 : // Now create the value
255 5 : std::vector<unsigned> shape( gridobject.getNbin(true) );
256 5 : ActionWithValue::addValueWithDerivatives( shape ); setNotPeriodic();
257 : } else {
258 6 : std::vector<unsigned> shape( 3 ); shape[0]=gbin[0]; shape[1]=shape[2]=1;
259 6 : ActionWithValue::addValueWithDerivatives( shape ); setNotPeriodic();
260 : }
261 22 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
262 11 : getPntrToComponent(i)->setConstant();
263 11 : getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
264 : }
265 : // This ensures we set the flag to never active the action. We can say we have atoms here as we don't need them
266 : // to calculate the CV
267 11 : setupConstantValues( true );
268 11 : }
269 :
270 34 : unsigned ReadGridInSetup::getNumberOfDerivatives() {
271 34 : return dernames.size();
272 : }
273 :
274 3 : std::vector<std::string> ReadGridInSetup::getGridCoordinateNames() const {
275 3 : return dernames;
276 : }
277 :
278 12 : const GridCoordinatesObject& ReadGridInSetup::getGridCoordinatesObject() const {
279 12 : return gridobject;
280 : }
281 :
282 : }
283 : }
284 :
|
