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 {
67 0 : plumed_merror("should not be in ReadGridInSetup setupOnFirstStep" );
68 : }
69 0 : void performTask( const unsigned& current, MultiValue& myvals ) const override {
70 0 : plumed_merror("should not be in readGridInSetup performTask");
71 : }
72 : std::vector<std::string> getGridCoordinateNames() const override ;
73 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
74 : };
75 :
76 : PLUMED_REGISTER_ACTION(ReadGridInSetup,"REFERENCE_GRID")
77 :
78 21 : void ReadGridInSetup::registerKeywords( Keywords& keys ) {
79 21 : ActionWithGrid::registerKeywords(keys);
80 42 : keys.remove("SERIAL");
81 21 : keys.add("optional","FUNC","the function to compute on the grid");
82 21 : keys.add("compulsory","GRID_MIN","auto","the lower bounds for the grid");
83 21 : keys.add("compulsory","GRID_MAX","auto","the upper bounds for the grid");
84 21 : keys.add("compulsory","PERIODIC","are the grid directions periodic");
85 21 : keys.add("optional","GRID_BIN","the number of bins for the grid");
86 21 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
87 21 : 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.");
88 21 : keys.add("compulsory","FILE","the name of the file that contains the reference data");
89 21 : keys.add("compulsory","VALUE","the name of the value that should be read from the grid");
90 42 : keys.setValueDescription("grid","the constant function on the grid that was specified in input");
91 21 : }
92 :
93 11 : ReadGridInSetup::ReadGridInSetup(const ActionOptions&ao):
94 : Action(ao),
95 11 : ActionWithGrid(ao) {
96 : std::string func;
97 22 : parse("FUNC",func);
98 11 : if( func.length()>0 ) {
99 : // Read in stuff for grid
100 : std::vector<std::string> gmin;
101 10 : parseVector("GRID_MIN",gmin);
102 5 : std::vector<std::string> gmax(gmin.size());
103 10 : parseVector("GRID_MAX",gmax);
104 5 : std::vector<unsigned> gbin(gmin.size());
105 10 : parseVector("GRID_BIN",gbin);
106 5 : std::vector<std::string> pbc(gmin.size());
107 10 : parseVector("PERIODIC",pbc);
108 5 : std::vector<bool> ipbc( pbc.size() );
109 10 : for(unsigned i=0; i<ipbc.size(); ++i) {
110 5 : if( pbc[i]=="YES" ) {
111 : ipbc[i]=true;
112 5 : } else if( pbc[i]=="NO" ) {
113 : ipbc[i]=false;
114 : } else {
115 0 : error( pbc[i] + " is not a valid instruction to the PERIODIC keyword");
116 : }
117 : }
118 :
119 : // Read in the variables
120 10 : parseVector("VAR",dernames);
121 5 : if(dernames.size()==0) {
122 3 : dernames.resize(gmin.size());
123 3 : if(gmin.size()>3) {
124 0 : error("Using more than 3 arguments you should explicitly write their names with VAR");
125 : }
126 3 : if(dernames.size()>0) {
127 : dernames[0]="x";
128 : }
129 3 : if(dernames.size()>1) {
130 : dernames[1]="y";
131 : }
132 3 : if(dernames.size()>2) {
133 : dernames[2]="z";
134 : }
135 : }
136 5 : if(dernames.size()!=gmin.size()) {
137 0 : error("Size of VAR array should be the same as number of grid dimensions");
138 : }
139 :
140 : // Create the grid and the value of the grid
141 5 : createGridAndValue( "flat", ipbc, 0, gmin, gmax, gbin );
142 :
143 : // Read in stuff for function
144 5 : log.printf(" evaluating function : %s\n",func.c_str());
145 5 : log.printf(" with variables :");
146 10 : for(unsigned i=0; i<dernames.size(); i++) {
147 5 : log.printf(" %s",dernames[i].c_str());
148 : }
149 5 : log.printf("\n");
150 5 : log.printf(" on %d", gbin[0]);
151 5 : for(unsigned i=1; i<gbin.size(); ++i) {
152 0 : log.printf(" by %d \n", gbin[i]);
153 : }
154 5 : log.printf(" grid of points between (%s", gmin[0].c_str() );
155 5 : for(unsigned i=1; i<gmin.size(); ++i) {
156 0 : log.printf(", %s", gmin[i].c_str() );
157 : }
158 5 : log.printf(") and (%s", gmax[0].c_str() );
159 5 : for(unsigned i=1; i<gmax.size(); ++i) {
160 0 : log.printf(", %s", gmax[i].c_str() );
161 : }
162 5 : log.printf(")\n");
163 :
164 5 : lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(leptonConstants);
165 5 : log<<" function as parsed by lepton: "<<pe<<"\n";
166 5 : lepton::CompiledExpression expression=pe.createCompiledExpression();
167 10 : for(auto &p: expression.getVariables()) {
168 5 : if(std::find(dernames.begin(),dernames.end(),p)==dernames.end()) {
169 0 : error("variable " + p + " is not defined");
170 : }
171 : }
172 5 : log<<" derivatives as computed by lepton:\n";
173 5 : std::vector<lepton::CompiledExpression> expression_deriv( dernames.size() );
174 10 : for(unsigned i=0; i<dernames.size(); i++) {
175 10 : lepton::ParsedExpression pe=lepton::Parser::parse(func).differentiate(dernames[i]).optimize(leptonConstants);
176 5 : log<<" "<<pe<<"\n";
177 5 : expression_deriv[i]=pe.createCompiledExpression();
178 : }
179 : // And finally calculate all the grid points
180 5 : std::vector<double> dder( dernames.size() ), xx( dernames.size() );
181 5 : Value* valout=getPntrToComponent(0);
182 115 : for(unsigned index=0; index<valout->getNumberOfValues(); ++index) {
183 110 : gridobject.getGridPointCoordinates( index, xx );
184 220 : for(unsigned j=0; j<xx.size(); ++j) {
185 : try {
186 110 : expression.getVariableReference(dernames[j])=xx[j];
187 0 : } catch(PLMD::lepton::Exception& exc) {
188 : // this is necessary since in some cases lepton things a variable is not present even though it is present
189 : // e.g. func=0*x
190 0 : }
191 : }
192 110 : valout->set( index, expression.evaluate() );
193 220 : for(unsigned k=0; k<xx.size(); ++k) {
194 220 : for(unsigned j=0; j<xx.size(); ++j) {
195 : try {
196 110 : expression_deriv[k].getVariableReference(dernames[j])=xx[j];
197 0 : } catch(PLMD::lepton::Exception& exc) {
198 : // this is necessary since in some cases lepton things a variable is not present even though it is present
199 : // e.g. func=0*x
200 0 : }
201 : }
202 110 : valout->addGridDerivatives( index, k, expression_deriv[k].evaluate() );
203 : }
204 : }
205 15 : } else {
206 : std::string valuestr;
207 12 : parse("VALUE",valuestr);
208 : std::string tstyle, filen;
209 12 : parse("FILE",filen);
210 6 : if( filen.length()>0 ) {
211 6 : std::size_t dot=filen.find_first_of(".");
212 6 : if( dot!=std::string::npos ) {
213 12 : tstyle=filen.substr(dot+1);
214 : }
215 6 : if( tstyle!="grid" ) {
216 0 : error("can only read in grids using read value in setup");
217 : }
218 6 : log.printf(" reading function %s on grid from file %s \n", valuestr.c_str(), filen.c_str() );
219 : }
220 6 : IFile ifile;
221 6 : ifile.open(filen);
222 6 : if( !ifile.FieldExist( valuestr ) ) {
223 0 : error("could not find grid value in input file");
224 : }
225 : std::vector<std::string> fieldnames;
226 6 : ifile.scanFieldList( fieldnames );
227 :
228 : // Retrieve the names of the variables the grid is computed over
229 : bool flatgrid=false;
230 54 : for(unsigned i=0; i<fieldnames.size(); ++i) {
231 48 : if( fieldnames[i].find("min_")!=std::string::npos ) {
232 : flatgrid=true;
233 : }
234 96 : std::size_t dot = fieldnames[i].find_first_of("d" + valuestr + "_" );
235 96 : if( fieldnames[i].find("d" + valuestr + "_")!=std::string::npos ) {
236 36 : dernames.push_back( fieldnames[i].substr(dot+2+valuestr.length()) );
237 : }
238 : }
239 6 : if( flatgrid && dernames.size()==0 ) {
240 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 + "_");
241 : }
242 : // Now get all the header data for the grid
243 6 : std::vector<std::string> gmin( dernames.size() ), gmax( dernames.size() );
244 : std::string pstring;
245 : int gbin1;
246 6 : std::vector<unsigned> gbin( dernames.size() );
247 0 : std::vector<bool> ipbc( dernames.size() );
248 6 : if( !flatgrid ) {
249 6 : ifile.scanField( "nbins", gbin1);
250 6 : gbin[0]=gbin1;
251 12 : createGridAndValue( "fibonacci", ipbc, gbin[0], gmin, gmax, gbin );
252 : } else {
253 0 : for(unsigned i=0; i<dernames.size(); ++i) {
254 0 : ifile.scanField( "min_" + dernames[i], gmin[i]);
255 0 : ifile.scanField( "max_" + dernames[i], gmax[i]);
256 0 : ifile.scanField( "periodic_" + dernames[i], pstring );
257 0 : ifile.scanField( "nbins_" + dernames[i], gbin1);
258 0 : gbin[i]=gbin1;
259 0 : if( pstring=="true" ) {
260 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]);
261 : ipbc[i]=true;
262 0 : } else if( pstring=="false" ) {
263 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]);
264 : ipbc[i]=false;
265 : } else {
266 0 : error("do not understand periodidicy of " + dernames[i] );
267 : }
268 :
269 0 : bool hasder=ifile.FieldExist( "d" + valuestr + "_" + dernames[i] );
270 0 : if( !hasder ) {
271 0 : plumed_merror("missing derivatives from grid file");
272 : }
273 : }
274 0 : createGridAndValue( "flat", ipbc, 0, gmin, gmax, gbin );
275 : }
276 : // And finally read all the grid points
277 6 : Value* valout=getPntrToComponent(0);
278 6 : std::vector<double> dder( dernames.size() ), xx( dernames.size() );
279 2106 : for(unsigned i=0; i<valout->getNumberOfValues(); ++i) {
280 : double x, val;
281 2100 : ifile.scanField( valuestr, val );
282 8400 : for(unsigned j=0; j<dernames.size(); ++j) {
283 6300 : ifile.scanField(dernames[j],x);
284 6300 : if( !flatgrid ) {
285 12600 : ifile.scanField("nbins", gbin1);
286 : } else {
287 0 : xx[j]=x+gridobject.getGridSpacing()[j]/2.0;
288 0 : ifile.scanField( "min_" + dernames[j], gmin[j]);
289 0 : ifile.scanField( "max_" + dernames[j], gmax[j]);
290 0 : ifile.scanField( "nbins_" + dernames[j], gbin1);
291 0 : ifile.scanField( "periodic_" + dernames[j], pstring );
292 : }
293 : }
294 8400 : for(unsigned j=0; j<dernames.size(); ++j) {
295 12600 : ifile.scanField( "d" + valuestr + "_" + dernames[j], dder[j] );
296 : }
297 :
298 2100 : unsigned index=gridobject.getIndex(xx);
299 2100 : if( !flatgrid ) {
300 2100 : index=i;
301 : }
302 2100 : valout->set( index, val );
303 8400 : for(unsigned j=0; j<dernames.size(); ++j) {
304 6300 : valout->addGridDerivatives( index, j, dder[j] );
305 : }
306 2100 : ifile.scanField();
307 : }
308 6 : ifile.close();
309 6 : }
310 11 : }
311 :
312 11 : void ReadGridInSetup::createGridAndValue( const std::string& gtype, const std::vector<bool>& ipbc, const unsigned& nfermi,
313 : const std::vector<std::string>& gmin, const std::vector<std::string>& gmax,
314 : const std::vector<unsigned>& gbin ) {
315 11 : gridobject.setup( gtype, ipbc, nfermi, 0.0 );
316 : std::vector<double> gspacing;
317 11 : if( gtype=="flat" ) {
318 5 : gridobject.setBounds( gmin, gmax, gbin, gspacing );
319 : // Now create the value
320 5 : std::vector<unsigned> shape( gridobject.getNbin(true) );
321 5 : ActionWithValue::addValueWithDerivatives( shape );
322 5 : setNotPeriodic();
323 : } else {
324 6 : std::vector<unsigned> shape( 3 );
325 6 : shape[0]=gbin[0];
326 6 : shape[1]=shape[2]=1;
327 6 : ActionWithValue::addValueWithDerivatives( shape );
328 6 : setNotPeriodic();
329 : }
330 22 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
331 11 : getPntrToComponent(i)->setConstant();
332 11 : getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
333 : }
334 : // This ensures we set the flag to never active the action. We can say we have atoms here as we don't need them
335 : // to calculate the CV
336 11 : setupConstantValues( true );
337 11 : }
338 :
339 34 : unsigned ReadGridInSetup::getNumberOfDerivatives() {
340 34 : return dernames.size();
341 : }
342 :
343 3 : std::vector<std::string> ReadGridInSetup::getGridCoordinateNames() const {
344 3 : return dernames;
345 : }
346 :
347 12 : const GridCoordinatesObject& ReadGridInSetup::getGridCoordinatesObject() const {
348 12 : return gridobject;
349 : }
350 :
351 : }
352 : }
353 :
|
