Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2024 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 :
23 : #include "CLTool.h"
24 : #include "core/CLToolRegister.h"
25 : #include "tools/Tools.h"
26 : #include "tools/SwitchingFunction.h"
27 : #include <string>
28 : #include <iostream>
29 : #include <iomanip>
30 : #include <sstream>
31 : #include <cmath>
32 : #include <cstdio>
33 : #include <limits>
34 :
35 : //+PLUMEDOC TOOLS plotswitch
36 : /*
37 : plotswitch is a tool that takes a the input of a switching function and tabulates the output on the terminal
38 :
39 : The tabulated data is compatible with gnuplot and numpy.loadtxt so you can use this tool to plot any functions that you plan to use with commands like
40 : [LESS_THAN](LESS_THAN.md) or [MORE_THAN](MORE_THAN.md).
41 :
42 : Without options plotswitch will tabulate 50 points between 0 and `R_0`, and then continue in tabulating points with the same step until the value of `D_MAX` is reached. If
43 : `D_MAX is unset then a value of `2R_0` is used in place of `D_MAX`.
44 :
45 : Without options plotswitch will tabulate data calling calculateSqr, since should be the most used option within the various colvars
46 :
47 : Note that if `R_0` happens to be between "from" and "to" the number of steps may not be exacly the number requested as this command forces the value of the function at r0 to be calculated.
48 :
49 : The various `--rational**` options use the special set option for the rational, like in COORDINATION.
50 :
51 : ## Examples
52 :
53 : If this command is called without options like this:
54 :
55 : ```plumed
56 : plumed plotswitch > plot.dat
57 : ```
58 :
59 : The function calculated is:
60 :
61 : $$
62 : f(x) = \frac{1 - x^6}{1 - x^{12}}
63 : $$
64 :
65 : If you want to use a function that is differnt to this one you use the `--switch` keyword as shown below:
66 :
67 : ```plumed
68 : plumed plotswitch --switch="RATIONAL NN=5 MM=9 R_0=1.3" --from=1.29999 --to=1.30001 --steps=100 > plot.dat
69 : ```
70 :
71 : The `--switch` keyword here takes the input for switching functions that is discussed in the documentation for
72 : [LESS_THAN](LESS_THAN.md). Notice also that if you use this command with an older plumed version you will see a discontinuity in dfunc at around 1.3
73 : (if you use gnuplot with "p 'plot.dat' u 1:3 w l t 'dfunc', 'plot.dat' u 1:2 w l axis x1y2 t 'res'")
74 :
75 : The following example shows another way of generating the `plot.dat` that is output by the command above:
76 :
77 : ```plumed
78 : plumed plotswitch --rationalR_0=1.3 --rationalNN=5 --rationalMM=9 --rationalD_0=0 --from=1.29999 --to=1.30001 --steps=100 > plot.dat
79 : ```
80 :
81 : As with [LESS_THAN](LESS_THAN.md), there is a alternative to `--switch` that can be used for sepcifying the parameters of RATIONAL switching function.
82 :
83 : */
84 : //+ENDPLUMEDOC
85 :
86 : namespace PLMD {
87 : namespace cltools {
88 : class SwitchingPlotter : public CLTool {
89 : public:
90 : explicit SwitchingPlotter(const CLToolOptions&);
91 : static void registerKeywords( Keywords& );
92 :
93 : int main( FILE*, FILE*, Communicator& ) override;
94 :
95 : };
96 21976 : PLUMED_REGISTER_CLTOOL(SwitchingPlotter,"plotswitch")
97 :
98 7321 : void SwitchingPlotter::registerKeywords( Keywords& keys ) {
99 7321 : CLTool::registerKeywords( keys );
100 7321 : keys.add("compulsory","--switch",
101 : "RATIONAL NN=6 R_0=1.0","the input to pass to the switching function,"
102 : " please remeber the quotes");
103 7321 : keys.add("compulsory","--steps","50",
104 : "the number of steps between 0 and R_O, or in the specified interval");
105 7321 : keys.add("compulsory","--from","-1",
106 : "the start of the interval, if negative will be set to 0");
107 7321 : keys.add("compulsory","--to","-1",
108 : "the end of the interval, will be D_MAX or 2*R_0 if D_MAX is not set");
109 7321 : keys.add("compulsory","--plotprecision","8",
110 : "the precision to use for the tabulated results");
111 7321 : keys.add("compulsory","--rationalR_0","-1",
112 : "The r_0 parameter of the switching function, this will activate the "
113 : "--rational options, note that this will ignore the --switch option silently");
114 7321 : keys.add("compulsory","--rationalNN","6",
115 : "The n parameter of the switching function");
116 7321 : keys.add("compulsory","--rationalMM","0",
117 : "The m parameter of the switching function; 0 implies 2*NN");
118 7321 : keys.add("compulsory","--rationalD_0","0.0",
119 : "The d_0 parameter of the switching function");
120 7321 : keys.addFlag("--nosquare",false,"use calculate instead of calculateSqr");
121 7321 : keys.add("compulsory","--centerrange","-1",
122 : "centers the visualization in R_0 in a range given epsilons times r_0"
123 : ", note that specifying this will overide all the other range options");
124 7321 : }
125 :
126 13 : SwitchingPlotter::SwitchingPlotter(const CLToolOptions& co ):
127 13 : CLTool(co) {
128 13 : inputdata=inputType::commandline;
129 13 : }
130 6 : int SwitchingPlotter::main( FILE*, FILE* out, Communicator& ) {
131 : //collecting options:
132 : std::string swInput;
133 6 : parse("--switch",swInput);
134 : bool dontOptimize;
135 6 : parseFlag("--nosquare",dontOptimize);
136 : int Nsteps;
137 6 : parse("--steps",Nsteps);
138 : double lowerLimit;
139 6 : parse("--from",lowerLimit);
140 : double upperLimit;
141 6 : parse("--to",upperLimit);
142 : unsigned plotPrecision;
143 6 : parse("--plotprecision",plotPrecision);
144 : int rationalNN;
145 6 : parse("--rationalNN",rationalNN);
146 : int rationalMM;
147 6 : parse("--rationalMM",rationalMM);
148 : double rationalD_0;
149 6 : parse("--rationalD_0",rationalD_0);
150 : double rationalR_0;
151 6 : parse("--rationalR_0",rationalR_0);
152 : //this works only because we use lepton to parse the numbers
153 : double centerrange;
154 8 : parse("--centerrange",centerrange);
155 : //setting up the switching function
156 6 : PLMD::SwitchingFunction switchingFunction;
157 6 : if (rationalR_0>0) {
158 1 : switchingFunction.set(rationalNN,rationalMM,rationalR_0,rationalD_0);
159 : } else {
160 : std::string errors;
161 5 : switchingFunction.set(swInput,errors);
162 5 : if( errors.length()!=0 ) {
163 0 : error("problem reading SWITCH keyword : " + errors );
164 : }
165 : }
166 :
167 : //setting up the limits:
168 6 : const double r0 = switchingFunction.get_r0();
169 6 : const double dmax = switchingFunction.get_dmax();
170 :
171 6 : if (lowerLimit <0) {
172 4 : lowerLimit=0.0;
173 : }
174 6 : if (upperLimit < 0) {
175 4 : upperLimit = dmax;
176 4 : if (! (upperLimit < std::numeric_limits<double>::max())) {
177 3 : upperLimit = 2*r0;
178 : }
179 : }
180 6 : if(centerrange>0) {
181 1 : upperLimit=(1.0+centerrange*PLMD::epsilon)*r0;
182 1 : lowerLimit=(1.0-centerrange*PLMD::epsilon)*r0;
183 : }
184 6 : if (Nsteps<=0) {
185 1 : error("--steps should be a positive integer");
186 : }
187 5 : if (!(upperLimit>lowerLimit)) {
188 1 : error("the end of the interval should be larger than its start");
189 : }
190 12 : const double step = [=]() {
191 4 : if(r0 > lowerLimit && r0< upperLimit) {
192 : //this will make the step pass trough r0:
193 : //r0 is the point number nr0 (rounded fraction of Nsteps)
194 4 : long nr0=std::lround((r0-lowerLimit)/(upperLimit-lowerLimit)*Nsteps);
195 4 : if(nr0<1) {
196 : nr0=1;
197 : }
198 4 : return (r0-lowerLimit)/nr0;
199 : }
200 0 : return (upperLimit-lowerLimit)/double(Nsteps);
201 : }
202 4 : ();
203 4 : if (!(step>0.0)) {
204 0 : error("I calculated a non positive step");
205 : }
206 :
207 : //finally doing the job
208 : //descriptions starts with the values of "r_0"
209 4 : std::ostringstream ostr;
210 12 : ostr <<"#r val dfunc ( r_0="<<switchingFunction.description()<<")\n";
211 : // the points are computed as lowerLimit+i*step rather than by accumulating the step:
212 : // this avoids accumulating rounding errors, and avoids an infinite loop
213 : // if the step is smaller than the spacing of the floating point numbers around x
214 : // (e.g. with a small --centerrange)
215 32 : for(long i=0;; i++) {
216 36 : const double x=lowerLimit+i*step;
217 36 : if(!(x<upperLimit)) {
218 : break;
219 : }
220 32 : double dfunc=0.0;
221 : double res;
222 32 : if(dontOptimize) {
223 8 : res=switchingFunction.calculate(x,dfunc);
224 : } else {
225 24 : res=switchingFunction.calculateSqr(x*x,dfunc);
226 : }
227 32 : ostr << std::setprecision(plotPrecision) << x << "\t"
228 32 : << std::setprecision(plotPrecision) << res << "\t"
229 32 : << std::setprecision(plotPrecision) << dfunc << '\n';
230 32 : }
231 4 : std::fprintf(out,"%s",ostr.str().c_str());
232 4 : return 0;
233 6 : }
234 :
235 : } //namespace cltools
236 : } // namespace PLMD
|
