Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2023 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 "KDE.h"
23 : #include "core/ActionShortcut.h"
24 : #include "core/ActionRegister.h"
25 :
26 : //+PLUMEDOC ANALYSIS KDE
27 : /*
28 : Create a histogram from the input scalar/vector/matrix using KDE
29 :
30 : This action can be used to construct instantaneous distributions for quantities by using [kernel density esstimation](https://en.wikipedia.org/wiki/Kernel_density_estimation).
31 : The input arguments must all have the same rank and size but you can use a scalar, vector or matrix in input. The distribution
32 : of this quantity on a grid is then computed using kernel density estimation.
33 :
34 : The following example demonstrates how this action can be used with a scalar as input:
35 :
36 : ```plumed
37 : d1: DISTANCE ATOMS=1,2
38 : kde: KDE ARG=d1 GRID_MIN=0.0 GRID_MAX=1.0 GRID_BIN=100 BANDWIDTH=0.2
39 : DUMPGRID ARG=kde STRIDE=1 FILE=kde.grid
40 : ```
41 :
42 : This input outputs a different file on every time step. These files contain a function stored on a grid. The function output in this case
43 : consists of a single Gaussian with $\sigma=0.2$ that is centered on the instantaneous value of the distance between atoms 1 and 2. Obviously,
44 : you are unlikely to use an input like the one above. The more usual thing to do would be to accumulate the histogram over the course of a
45 : few trajectory frames using the [ACCUMULATE](ACCUMULATE.md) command as has been done in the input below, which estimates a histogram as a function
46 : of two collective variables:
47 :
48 : ```plumed
49 : d1: DISTANCE ATOMS=1,2
50 : d2: DISTANCE ATOMS=1,2
51 : kde: KDE ARG=d1,d2 GRID_MIN=0.0,0.0 GRID_MAX=1.0,1.0 GRID_BIN=100,100 BANDWIDTH=0.2,0.2
52 : histo: ACCUMULATE ARG=kde STRIDE=1
53 : DUMPGRID ARG=histo FILE=histo.grid STRIDE=10000
54 : ```
55 :
56 : Notice, that you can also achieve something similar by using the [HISTOGRAM](HISTOGRAM.md) shortcut.
57 :
58 : ## Controlloing the grid
59 :
60 : If you prefer to specify the grid spacing rather than the number of bins you can do so using the GRID_SPACING keyword as shown below:
61 :
62 : ```plumed
63 : d1: DISTANCE ATOMS=1,2
64 : kde: KDE ARG=d1 GRID_MIN=0.0 GRID_MAX=1.0 GRID_SPACING=0.01 BANDWIDTH=0.2
65 : DUMPGRID ARG=kde STRIDE=1 FILE=kde.grid
66 : ```
67 :
68 : If $x$ is one of the input arguments to the KDE action and $x<g_{min}$ or $x>g_{max}$, where $g_{min}$ and $g_{max}$ are the minimum
69 : and maximum values on the grid for that argument that were specified using GRID_MIN and GRID_MAX, then by PLUMED will crash.
70 :
71 : Notice also that when you use Gaussian kernels to accumulate a denisty as in the input above you need to define a cutoff beyond, which the
72 : Gaussian (which is a function with infinite support) is assumed not to contribute to the accumulated density. When setting this cutoff you
73 : set the value of $x$ in the following expression $\sigma \sqrt{2*x}$, where $\sigma$ is the bandwidth. By default $x$ is set equal to 6.25 but
74 : you can change this value by using the CUTOFF keyword as shown below:
75 :
76 : ```plumed
77 : d1: DISTANCE ATOMS=1,2
78 : kde: KDE ARG=d1 GRID_MIN=0.0 GRID_MAX=1.0 GRID_SPACING=0.01 BANDWIDTH=0.2 CUTOFF=6.25
79 : DUMPGRID ARG=kde STRIDE=1 FILE=kde.grid
80 : ```
81 :
82 : ## Constructing the density
83 :
84 : If you are performing a simulation in the NVT ensemble and wish to look at the density as a function of position in the cell you can use an input like the one shown below:
85 :
86 : ```plumed
87 : a: FIXEDATOM AT=0,0,0
88 : dens: DISTANCES ATOMS=1-100 ORIGIN=a COMPONENTS
89 : kde: KDE ARG=dens.x,dens.y,dens.z GRID_BIN=100,100,100 BANDWIDTH=0.05,0.05,0.05
90 : DUMPGRID ARG=kde STRIDE=1 FILE=density
91 : ```
92 :
93 : Notice that you do not need to specify GRID_MIN and GRID_MAX values with this input. In this case PLUMED gets the extent of the grid from the cell vectors during the first
94 : step of the simulation.
95 :
96 : ## Specifying a non diagonal bandwidth
97 :
98 : If for any reason you want to use a bandwidth that is not diagonal when doing kensity density estimation you can do by using an input similar to the one shown below:
99 :
100 : ```plumed
101 : d1: DISTANCE ATOMS=1,2
102 : d2: DISTANCE ATOMS=1,2
103 : kde: KDE ...
104 : ARG=d1,d2 GRID_MIN=0.0,0.0
105 : GRID_MAX=1.0,1.0 GRID_BIN=100,100
106 : BANDWIDTH=0.2,0.1,0.1,0.2 HEIGHTS=1
107 : ...
108 : histo: ACCUMULATE ARG=kde STRIDE=1
109 : DUMPGRID ARG=histo FILE=histo.grid STRIDE=10000
110 : ```
111 :
112 : As there are two arguments for this KDE action the four numbers passed in the bandwdith parameter are interepretted as a $2\times 2$ matrix.
113 : Notice that you can also pass the information for the bandwidth in from another argument as has been done here:
114 :
115 : ```plumed
116 : m: CONSTANT VALUES=0.2,0.1,0.1,0.2 NROWS=2 NCOLS=2
117 :
118 : d1: DISTANCE ATOMS=1,2
119 : d2: DISTANCE ATOMS=1,2
120 : kde: KDE ...
121 : ARG=d1,d2 GRID_MIN=0.0,0.0
122 : GRID_MAX=1.0,1.0 GRID_BIN=100,100
123 : BANDWIDTH=m HEIGHTS=1
124 : ...
125 : histo: ACCUMULATE ARG=kde STRIDE=1
126 : DUMPGRID ARG=histo FILE=histo.grid STRIDE=10000
127 : ```
128 :
129 : In this case the input is equivalent to the first input above and the bandwidth is a constant. You could, however, also use a non-constant value as input to the BANDWIDTH keyword.
130 :
131 : ## Working with vectors and scalars
132 :
133 : If the input to your KDE action is a set of scalars it appears odd to separate the process of computing the KDE from the process of accumulating the histogram. However, if
134 : you are using vectors as in the example below, this division can be helpful.
135 :
136 : ```plumed
137 : d1: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8 ATOMS5=9,10
138 : kde: KDE ARG=d1 GRID_MIN=0.0 GRID_MAX=1.0 GRID_BIN=100 BANDWIDTH=0.2
139 : ```
140 :
141 : In the papea cited in the bibliography below, the [KL_ENTROPY](KL_ENTROPY.md) between the instantaneous distribution of CVs and a reference distribution was introduced
142 : as a collective variable. As is detailed in the documentation for that action, the ability to calculate the instaneous histogram from an input vector is essential to
143 : reproducing these calculations.
144 :
145 : Notice that you can also use a one or multiple matrices in the input for a KDE object. The example below uses the angles between the z axis and set of bonds aroud two
146 : atoms:
147 :
148 : ```plumed
149 : d1: DISTANCE_MATRIX GROUPA=1,2 GROUPB=3-10 COMPONENTS
150 : phi: CUSTOM ARG=d1.z,d1.w FUNC=acos(x/y) PERIODIC=NO
151 : kde: KDE ARG=phi GRID_MIN=0 GRID_MAX=pi GRID_BIN=200 BANDWIDTH=0.1
152 : ```
153 :
154 : ## Using different weights
155 :
156 : In all the inputs above the kernels that are added to the grid on each step are Gaussians with that are normalised so that their integral over all space is one. If you want your
157 : Gaussians to have a particular height you can use the HEIGHT keyword as illustrated below:
158 :
159 : ```plumed
160 : d1: CONTACT_MATRIX GROUPA=1,2 GROUPB=3-10 SWITCH={RATIONAL R_0=0.1} COMPONENTS
161 : mag: CUSTOM ARG=d1.x,d1.y,d1.z FUNC=x*x+y*y+z*z PERIODIC=NO
162 : phi: CUSTOM ARG=d1.z,mag FUNC=acos(x/sqrt(y)) PERIODIC=NO
163 : kde: KDE ARG=phi GRID_MIN=0 GRID_MAX=pi HEIGHTS=d1.w GRID_BIN=200 BANDWIDTH=0.1
164 : ```
165 :
166 : As indicated above, the HEIGHTS keyword should be passed a Value that has the same rank and size as the arguments that are passed using the ARG keyword. Each of the Gaussian kernels
167 : that are added to the grid in this case have a value equal to the weight at the maximum of the function.
168 :
169 : Notice that you can also use the VOLUMES keyword in a similar way as shown below:
170 :
171 : ```plumed
172 : d1: CONTACT_MATRIX GROUPA=1,2 GROUPB=3-10 SWITCH={RATIONAL R_0=0.1} COMPONENTS
173 : mag: CUSTOM ARG=d1.x,d1.y,d1.z FUNC=x*x+y*y+z*z PERIODIC=NO
174 : phi: CUSTOM ARG=d1.z,mag FUNC=acos(x/sqrt(y)) PERIODIC=NO
175 : kde: KDE ARG=phi GRID_MIN=0 GRID_MAX=pi VOLUMES=d1.w GRID_BIN=200 BANDWIDTH=0.1
176 : ```
177 :
178 : Now, however, the integral of the Gaussians over all space are equal to the elements of d1.w.
179 :
180 : */
181 : //+ENDPLUMEDOC
182 :
183 : //+PLUMEDOC ANALYSIS SPHERICAL_KDE
184 : /*
185 : Create a histogram from the input scalar/vector/matrix using SPHERICAL_KDE
186 :
187 : This action operates similarly to [KDE](KDE.md) but it is designed to be used for investigating [directional statistics]().
188 : It is particularly useful if you are looking at the distribution of bond vectors as illustrated in the input below:
189 :
190 : ```plumed
191 : # Calculate all the bond vectors
192 : d1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.1} COMPONENTS
193 : # Normalise the bond vectors
194 : mag: CUSTOM ARG=d1.x,d1.y,d1.z FUNC=sqrt(x*x+y*y+z*z) PERIODIC=NO
195 : d1x: CUSTOM ARG=d1.x,mag FUNC=x/y PERIODIC=NO
196 : d1y: CUSTOM ARG=d1.y,mag FUNC=x/y PERIODIC=NO
197 : d1z: CUSTOM ARG=d1.z,mag FUNC=x/y PERIODIC=NO
198 : # And construct the KDE
199 : kde: SPHERICAL_KDE ARG=d1x,d1y,d1z HEIGHTS=d1.w CONCENTRATION=100 GRID_BIN=144
200 : ```
201 :
202 : Each bond vector here contributes a [Fisher von-Mises kernel](https://en.wikipedia.org/wiki/Von_Mises–Fisher_distribution) to the spherical grid. This spherical grid is constructed
203 : using a [Fibonnacci sphere algorithm](https://stackoverflow.com/questions/9600801/evenly-distributing-n-points-on-a-sphere) so the number of specified using the GRID_BIN keyword must be a Fibonacci number.
204 :
205 : */
206 : //+ENDPLUMEDOC
207 :
208 : namespace PLMD {
209 : namespace gridtools {
210 :
211 : template <class K, class P>
212 : class KDEGridTools {
213 : public:
214 : double dp2cutoff;
215 : std::vector<double> gspacing;
216 : std::vector<std::size_t> nbin;
217 : std::vector<std::string> gmin, gmax;
218 : static void registerKeywords( Keywords& keys );
219 : static void readHeightKeyword( bool canusevol, std::size_t nargs, const std::vector<std::string>& bw, ActionWithArguments* action );
220 : static void readBandwidth( std::size_t nargs, ActionWithArguments* action, std::vector<std::string>& bw );
221 : static void readBandwidthKeyword( std::size_t nargs, ActionWithArguments* action, std::vector<std::string>& bw, std::vector<Value*>& bwargs );
222 : static void readBandwidthAndHeight( const P& params, ActionWithArguments* action );
223 : static void convertHeightsToVolumes( const std::size_t& nargs, const std::vector<std::string>& bw, const std::string& volstr, ActionWithArguments* action );
224 : static void readGridParameters( KDEGridTools<K,P>& g, ActionWithArguments* action, GridCoordinatesObject& gridobject, std::vector<std::size_t>& shape );
225 : static void setupGridBounds( KDEGridTools<K,P>& g, const Tensor& box, GridCoordinatesObject& gridobject, const std::vector<Value*>& args, Value* myval );
226 : static void getDiscreteSupport( const KDEGridTools<K,P>& g, P& p, const K& kp, std::vector<unsigned>& nneigh, GridCoordinatesObject& gridobject );
227 : static void getNeighbors( const P& p, K& kp, const GridCoordinatesObject& gridobject, const std::vector<unsigned>& nneigh, unsigned& num_neighbors, std::vector<unsigned>& neighbors );
228 : };
229 :
230 : template <class K, class P>
231 311 : void KDEGridTools<K,P>::registerKeywords( Keywords& keys ) {
232 311 : keys.add("optional","BANDWIDTH","the bandwidths for kernel density esimtation");
233 311 : keys.add("optional","VOLUMES","this keyword take the label of an action that calculates a vector of values. The elements of this vector "
234 : "divided by the volume of the Gaussian are used as weights for the Gaussians");
235 311 : keys.add("optional","HEIGHTS","this keyword takes the label of an action that calculates a vector of values. The elements of this vector "
236 : "are used as weights for the Gaussians.");
237 311 : keys.add("compulsory","GRID_MIN","auto","the lower bounds for the grid");
238 311 : keys.add("compulsory","GRID_MAX","auto","the upper bounds for the grid");
239 311 : keys.add("compulsory","CUTOFF","6.25","the cutoff at which to stop evaluating the kernel functions is set equal to sqrt(2*x)*bandwidth in each direction where x is this number");
240 311 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
241 311 : keys.add("optional","GRID_BIN","the number of bins for the grid");
242 311 : }
243 :
244 : template <class K, class P>
245 65 : void KDEGridTools<K,P>::readHeightKeyword( bool canusevol, std::size_t nargs, const std::vector<std::string>& bw, ActionWithArguments* action ) {
246 : std::string weight_str;
247 130 : action->parse("HEIGHTS",weight_str);
248 : std::string str_nvals;
249 65 : if( (action->getPntrToArgument(0))->getRank()==2 ) {
250 : std::string nr, nc;
251 3 : Tools::convert( (action->getPntrToArgument(0))->getShape()[0], nr );
252 3 : Tools::convert( (action->getPntrToArgument(0))->getShape()[1], nc );
253 6 : str_nvals = nr + "," + nc;
254 : } else {
255 62 : Tools::convert( (action->getPntrToArgument(0))->getNumberOfValues(), str_nvals );
256 : }
257 65 : if( weight_str.length()>0 ) {
258 34 : KDEHelper<K,P,KDEGridTools<K,P>>::readKernelParameters( weight_str, action, "_heights", true );
259 48 : } else if( canusevol ) {
260 29 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_volumes: ONES SIZE=" + str_nvals ), false );
261 58 : KDEGridTools<K,P>::convertHeightsToVolumes(nargs,bw,action->getLabel() + "_volumes",action);
262 : } else {
263 19 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: ONES SIZE=" + str_nvals ), false );
264 38 : KDEHelper<K,P,KDEGridTools<K,P>>::addArgument( action->getLabel() + "_heights", action );
265 : }
266 65 : }
267 :
268 : template <>
269 13 : void KDEGridTools<DiagonalKernelParams,DiscreteKernel>::readBandwidthAndHeight( const DiscreteKernel& params, ActionWithArguments* action ) {
270 13 : std::size_t nargs = action->getNumberOfArguments();
271 13 : KDEGridTools<DiagonalKernelParams,DiscreteKernel>::readHeightKeyword( false, nargs, std::vector<std::string>(), action );
272 13 : }
273 :
274 : template<class K, class P>
275 66 : void KDEGridTools<K, P>::readBandwidthKeyword( std::size_t nargs, ActionWithArguments* action, std::vector<std::string>& bw, std::vector<Value*>& bwargs ) {
276 66 : action->parseVector("BANDWIDTH",bw);
277 66 : if( nargs>1 && bw.size()==1 ) {
278 0 : ActionWithArguments::interpretArgumentList( bw, action->plumed.getActionSet(), action, bwargs );
279 0 : if( bwargs.size()!=1 ) {
280 0 : action->error("invalid bandwidth found");
281 : }
282 : // Create a vector of ones with the right size
283 : std::string nvals;
284 0 : if( (action->getPntrToArgument(0))->getRank()==2 ) {
285 : std::string nr, nc;
286 0 : Tools::convert( (action->getPntrToArgument(0))->getShape()[0], nr );
287 0 : Tools::convert( (action->getPntrToArgument(0))->getShape()[1], nc );
288 0 : nvals = nr + "," + nc;
289 : } else {
290 0 : Tools::convert( (action->getPntrToArgument(0))->getNumberOfValues(), nvals );
291 : }
292 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_bwones: ONES SIZE=" + nvals ), false );
293 : }
294 66 : }
295 :
296 : template <class K, class P>
297 66 : void KDEGridTools<K,P>::readBandwidth( std::size_t nargs, ActionWithArguments* action, std::vector<std::string>& bw ) {
298 66 : plumed_assert( typeid(K)==typeid(DiagonalKernelParams) );
299 : std::vector<Value*> bwargs;
300 66 : readBandwidthKeyword( nargs, action, bw, bwargs );
301 66 : if( nargs>1 && bw.size()==1 ) {
302 0 : if( bwargs[0]->getRank()!=1 || bwargs[0]->getNumberOfValues()!=nargs ) {
303 0 : action->error("invalid input for bandwidth parameter");
304 : }
305 : std::string str_i;
306 0 : bw.resize( nargs );
307 0 : for(unsigned i=0; i<nargs; ++i) {
308 0 : Tools::convert( i+1, str_i );
309 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_scalar_bw" + str_i + ": SELECT_COMPONENTS ARG=" + bwargs[0]->getName() + " COMPONENTS=" + str_i ), false );
310 0 : bw[i] = action->getLabel() + "_bw" + str_i;
311 0 : action->plumed.readInputWords( Tools::getWords( bw[i] + ": CUSTOM ARG=" + action->getLabel() + "_bwones," + action->getLabel() + "_scalar_bw" + str_i + " FUNC=x*y PERIODIC=NO"), false );
312 : }
313 66 : } else if( bw.size()==nargs ) {
314 : double bwval;
315 : std::string str_i;
316 164 : for(unsigned i=0; i<nargs; ++i) {
317 98 : Tools::convert( i+1, str_i );
318 98 : if( Tools::convertNoexcept( bw[i], bwval ) && fabs(bwval)<epsilon ) {
319 4 : KDEHelper<K,P,KDEGridTools<K,P>>::readKernelParameters( bw[i], action, "_bwz" + str_i, true );
320 : } else {
321 192 : KDEHelper<K,P,KDEGridTools<K,P>>::readKernelParameters( bw[i], action, "_bw" + str_i, true );
322 : }
323 : }
324 : } else {
325 0 : action->error("wrong number of arguments specified in input to bandwidth parameter");
326 : }
327 66 : }
328 :
329 : template <>
330 18 : void KDEGridTools<DiagonalKernelParams,HistogramBeadKernel>::readBandwidthAndHeight( const HistogramBeadKernel& params, ActionWithArguments* action ) {
331 : std::vector<std::string> bw;
332 18 : std::size_t nargs = action->getNumberOfArguments();
333 18 : readBandwidth( nargs, action, bw );
334 18 : KDEGridTools<DiagonalKernelParams,HistogramBeadKernel>::readHeightKeyword( false, nargs, bw, action );
335 18 : }
336 :
337 : template <>
338 48 : void KDEGridTools<DiagonalKernelParams,RegularKernel<DiagonalKernelParams>>::readBandwidthAndHeight( const RegularKernel<DiagonalKernelParams>& params, ActionWithArguments* action ) {
339 : std::vector<Value*> bwargs;
340 : std::vector<std::string> bw;
341 48 : std::size_t nargs = action->getNumberOfArguments();
342 48 : readBandwidth( nargs, action, bw );
343 : std::string volstr;
344 96 : action->parse("VOLUMES",volstr);
345 48 : if( volstr.length()>0 ) {
346 14 : if( !params.canusevol ) {
347 0 : action->error("cannot use normalized kernels with selected kernel type");
348 : }
349 : // Check if we are using Gaussian kernels
350 14 : KDEHelper<DiagonalKernelParams,RegularKernel<DiagonalKernelParams>,KDEGridTools<DiagonalKernelParams,RegularKernel<DiagonalKernelParams>>>::readKernelParameters( volstr, action, "_volumes", false );
351 14 : convertHeightsToVolumes(nargs, bw, volstr, action);
352 : } else {
353 34 : KDEGridTools<DiagonalKernelParams,RegularKernel<DiagonalKernelParams>>::readHeightKeyword( params.canusevol, nargs, bw, action );
354 : }
355 96 : }
356 :
357 : template <>
358 0 : void KDEGridTools<NonDiagonalKernelParams,RegularKernel<NonDiagonalKernelParams>>::readBandwidthAndHeight( const RegularKernel<NonDiagonalKernelParams>& params, ActionWithArguments* action ) {
359 : std::vector<Value*> bwargs;
360 : std::vector<std::string> bw;
361 0 : std::size_t nargs = action->getNumberOfArguments();
362 0 : readBandwidthKeyword( nargs, action, bw, bwargs );
363 0 : if( nargs>1 && bw.size()==1 ) {
364 0 : if( bwargs[0]->getRank()!=2 || bwargs[0]->getShape()[0]!=nargs || bwargs[0]->getShape()[1]!=nargs ) {
365 0 : action->error("invalid input for bandwidth parameter");
366 : }
367 : std::string str_i, str_j;
368 0 : bw.resize( nargs*nargs );
369 0 : for(unsigned i=0; i<nargs; ++i) {
370 0 : Tools::convert( i+1, str_i );
371 0 : for(unsigned j=0; j<nargs; ++j) {
372 0 : Tools::convert( j+1, str_j );
373 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_scalar_bw" + str_i + "_" + str_j + ": SELECT_COMPONENTS ARG=" + bwargs[0]->getName() + " COMPONENTS=" + str_i + "." + str_j ), false );
374 0 : bw[i*nargs+j] = action->getLabel() + "_bw" + str_i + "_" + str_j;
375 0 : action->plumed.readInputWords( Tools::getWords( bw[i*nargs+j] + ": CUSTOM ARG=" + action->getLabel() + "_bwones," + action->getLabel() + "_scalar_bw" + str_i + "_" + str_j + " FUNC=x*y PERIODIC=NO"), false );
376 : }
377 : }
378 0 : } else if( bw.size()==nargs*nargs ) {
379 : std::string str_i, str_j;
380 0 : for(unsigned i=0; i<nargs; ++i) {
381 0 : Tools::convert( i+1, str_i );
382 0 : for(unsigned j=0; j<nargs; ++j) {
383 0 : Tools::convert( j+1, str_j );
384 0 : KDEHelper<NonDiagonalKernelParams,RegularKernel<NonDiagonalKernelParams>,KDEGridTools<NonDiagonalKernelParams,RegularKernel<NonDiagonalKernelParams>>>::readKernelParameters( bw[i*nargs+j], action, "_bw" + str_i + "_" + str_j, true );
385 : }
386 : }
387 : } else {
388 0 : action->error("wrong number of arguments specified in input to bandwidth parameter");
389 : }
390 : std::string volstr;
391 0 : action->parse("VOLUMES",volstr);
392 0 : if( volstr.length()>0 ) {
393 0 : if( !params.canusevol ) {
394 0 : action->error("cannot use normalized kernels with selected kernel type");
395 : }
396 : // Check if we are using Gaussian kernels
397 0 : KDEHelper<NonDiagonalKernelParams,RegularKernel<NonDiagonalKernelParams>,KDEGridTools<DiagonalKernelParams,RegularKernel<DiagonalKernelParams>>>::readKernelParameters( volstr, action, "_volumes", false );
398 0 : convertHeightsToVolumes(nargs, bw, volstr, action);
399 : } else {
400 0 : KDEGridTools<NonDiagonalKernelParams,RegularKernel<NonDiagonalKernelParams>>::readHeightKeyword( params.canusevol, nargs, bw, action );
401 : }
402 0 : }
403 :
404 : template <class K, class P>
405 43 : void KDEGridTools<K, P>::convertHeightsToVolumes( const std::size_t& nargs, const std::vector<std::string>& bw, const std::string& volstr, ActionWithArguments* action ) {
406 43 : if( bw.size()==nargs ) {
407 43 : unsigned nonzeroargs=0;
408 110 : for(unsigned i=0; i<nargs; ++i) {
409 67 : if( bw[i].find("_bwz")==std::string::npos ) {
410 65 : nonzeroargs++;
411 : }
412 : }
413 : std::string str_i, nargs_str;
414 43 : Tools::convert( nonzeroargs, nargs_str );
415 86 : std::string varstr = "VAR=h", funcstr = "FUNC=h/(sqrt((2*pi)^" + nargs_str + ")", argstr = "ARG=" + volstr;
416 110 : for(unsigned i=0; i<nargs; ++i) {
417 67 : if( bw[i].find("_bwz")!=std::string::npos ) {
418 2 : continue;
419 : }
420 65 : Tools::convert( i+1, str_i );
421 65 : varstr += ",b" + str_i;
422 130 : funcstr += "*b" + str_i;
423 130 : argstr += "," + bw[i];
424 : }
425 : funcstr += ")";
426 43 : if( (action->getPntrToArgument(0))->getNumberOfValues()==1 ) {
427 27 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: CUSTOM PERIODIC=NO " + argstr + " " + varstr + " " + funcstr), false );
428 : } else {
429 16 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: CUSTOM PERIODIC=NO " + argstr + " " + varstr + " " + funcstr + " MASK=" + volstr), false );
430 : }
431 86 : KDEHelper<K,P,KDEGridTools<K,P>>::addArgument( action->getLabel() + "_heights", action );
432 0 : } else if( bw.size()==nargs*nargs ) {
433 0 : action->error("have not implemented normalization parameter for non-diagonal kernels");
434 : }
435 :
436 43 : }
437 :
438 : template <class K, class P>
439 79 : void KDEGridTools<K,P>::readGridParameters( KDEGridTools<K,P>& g, ActionWithArguments* action, GridCoordinatesObject& gridobject, std::vector<std::size_t>& shape ) {
440 79 : g.gmin.resize( shape.size() );
441 79 : g.gmax.resize( shape.size() );
442 79 : action->parseVector("GRID_MIN",g.gmin);
443 79 : action->parseVector("GRID_MAX",g.gmax);
444 191 : for(unsigned i=0; i<g.gmin.size(); ++i) {
445 112 : if( g.gmin[i]=="auto" ) {
446 52 : action->log.printf(" for %dth coordinate min and max are set automatically \n", (i+1) );
447 52 : if( g.gmax[i]!="auto" ) {
448 0 : action->error("if gmin is set automatically gmax must also be set automatically");
449 : }
450 : Value* myarg = action->getPntrToArgument(i);
451 52 : if( myarg->isPeriodic() ) {
452 0 : if( g.gmin[i]=="auto" ) {
453 0 : myarg->getDomain( g.gmin[i], g.gmax[i] );
454 : } else {
455 : std::string str_min, str_max;
456 0 : myarg->getDomain( str_min, str_max );
457 0 : if( str_min!=g.gmin[i] || str_max!=g.gmax[i] ) {
458 0 : action->error("all periodic arguments should have the same domain");
459 : }
460 : }
461 52 : } else if( myarg->getName().find(".")!=std::string::npos ) {
462 52 : std::size_t dot = myarg->getName().find_first_of(".");
463 52 : std::string name = myarg->getName().substr(dot+1);
464 90 : if( name!="x" && name!="y" && name!="z" ) {
465 0 : action->error("cannot set GRID_MIN and GRID_MAX automatically if input argument is not component of distance");
466 : }
467 : } else {
468 0 : action->error("cannot set GRID_MIN and GRID_MAX automatically if input argument is not component of distance");
469 : }
470 : } else {
471 60 : action->log.printf(" for %dth coordinate min is set to %s and max is set to %s \n", (i+1), g.gmin[i].c_str(), g.gmax[i].c_str() );
472 : }
473 : }
474 :
475 79 : action->parseVector("GRID_BIN",g.nbin);
476 79 : action->parseVector("GRID_SPACING",g.gspacing);
477 158 : action->parse("CUTOFF",g.dp2cutoff);
478 :
479 79 : if( g.nbin.size()!=shape.size() && g.gspacing.size()!=shape.size() ) {
480 0 : action->error("GRID_BIN or GRID_SPACING must be set");
481 : }
482 : // Create a value
483 79 : std::vector<bool> ipbc( shape.size() );
484 191 : for(unsigned i=0; i<shape.size(); ++i) {
485 210 : if( (action->getPntrToArgument( i ))->isPeriodic() || g.gmin[i]=="auto" ) {
486 : ipbc[i]=true;
487 66 : if( g.nbin.size()==shape.size() ) {
488 66 : shape[i] = g.nbin[i];
489 : }
490 : } else {
491 : ipbc[i]=false;
492 46 : if( g.nbin.size()==shape.size() ) {
493 46 : shape[i] = g.nbin[i]+1;
494 : }
495 : }
496 : }
497 158 : gridobject.setup( "flat", ipbc, 0, 0.0 );
498 79 : }
499 :
500 : template <class K, class P>
501 76 : void KDEGridTools<K,P>::setupGridBounds( KDEGridTools<K,P>& g, const Tensor& box, GridCoordinatesObject& gridobject, const std::vector<Value*>& args, Value* myval ) {
502 181 : for(unsigned i=0; i<gridobject.getDimension(); ++i) {
503 105 : if( g.gmin[i]=="auto" ) {
504 : double lcoord, ucoord;
505 46 : std::size_t dot = args[i]->getName().find_first_of(".");
506 46 : std::string name = args[i]->getName().substr(dot+1);
507 46 : if( name=="x" ) {
508 24 : lcoord=-0.5*box(0,0);
509 24 : ucoord=0.5*box(0,0);
510 22 : } else if( name=="y" ) {
511 12 : lcoord=-0.5*box(1,1);
512 12 : ucoord=0.5*box(1,1);
513 10 : } else if( name=="z" ) {
514 10 : lcoord=-0.5*box(2,2);
515 10 : ucoord=0.5*box(2,2);
516 : } else {
517 0 : plumed_error();
518 : }
519 : // And convert to strings for bin and bmax
520 46 : Tools::convert( lcoord, g.gmin[i] );
521 46 : Tools::convert( ucoord, g.gmax[i] );
522 : }
523 : }
524 : // And setup the grid object
525 76 : gridobject.setBounds( g.gmin, g.gmax, g.nbin, g.gspacing );
526 76 : myval->setShape( gridobject.getNbin(true) );
527 76 : }
528 :
529 : template <class K, class P>
530 64 : void KDEGridTools<K, P>::getDiscreteSupport( const KDEGridTools<K,P>& g, P& p, const K& kp, std::vector<unsigned>& nneigh, GridCoordinatesObject& gridobject ) {
531 64 : std::size_t ng = gridobject.getDimension();
532 64 : plumed_assert( nneigh.size()==ng );
533 64 : std::vector<double> support( ng );
534 64 : P::getSupport( p, kp, g.dp2cutoff, support );
535 156 : for(unsigned i=0; i<ng; ++i) {
536 92 : nneigh[i] = static_cast<unsigned>( ceil( support[i]/gridobject.getGridSpacing()[i] ));
537 : }
538 64 : }
539 :
540 : template <>
541 12 : void KDEGridTools<DiagonalKernelParams,DiscreteKernel>::getDiscreteSupport( const KDEGridTools<DiagonalKernelParams,DiscreteKernel>& g, DiscreteKernel& p, const DiagonalKernelParams& kp, std::vector<unsigned>& nneigh, GridCoordinatesObject& gridobject ) {
542 12 : return;
543 : }
544 :
545 : template <class K, class P>
546 191342 : void KDEGridTools<K,P>::getNeighbors( const P& p, K& kp, const GridCoordinatesObject& gridobject, const std::vector<unsigned>& nneigh, unsigned& num_neighbors, std::vector<unsigned>& neighbors ) {
547 191342 : gridobject.getNeighbors( kp.at, nneigh, num_neighbors, neighbors );
548 191342 : }
549 :
550 : template <>
551 32290 : void KDEGridTools<DiagonalKernelParams,DiscreteKernel>::getNeighbors( const DiscreteKernel& p, DiagonalKernelParams& kp, const GridCoordinatesObject& gridobject, const std::vector<unsigned>& nneigh, unsigned& num_neighbors, std::vector<unsigned>& neighbors ) {
552 32290 : num_neighbors=1;
553 32290 : neighbors.resize(1);
554 65364 : for(unsigned i=0; i<kp.at.size(); ++i) {
555 33074 : kp.at[i] += 0.5*gridobject.getGridSpacing()[i];
556 : }
557 32290 : neighbors[0]=gridobject.getIndex( kp.at );
558 32290 : }
559 :
560 : class SphericalKDEGridTools {
561 : public:
562 : std::size_t nbins;
563 : static void registerKeywords( Keywords& keys );
564 : static void readBandwidthAndHeight( const UniversalVonMisses& params, ActionWithArguments* action );
565 : static void readGridParameters( SphericalKDEGridTools& g, ActionWithArguments* action, GridCoordinatesObject& gridobject, std::vector<std::size_t>& shape );
566 : static void setupGridBounds( SphericalKDEGridTools& g, const Tensor& box, GridCoordinatesObject& gridobject, const std::vector<Value*>& args, Value* myval ) {}
567 : static void getDiscreteSupport( const SphericalKDEGridTools& g, const UniversalVonMisses& p, const VonMissesKernelParams& kp, std::vector<unsigned>& nneigh, GridCoordinatesObject& gridobject );
568 : static void getNeighbors( const UniversalVonMisses& p, const VonMissesKernelParams& kp, const GridCoordinatesObject& gridobject, const std::vector<unsigned>& nneigh, unsigned& num_neighbors, std::vector<unsigned>& neighbors );
569 : };
570 :
571 12 : void SphericalKDEGridTools::registerKeywords( Keywords& keys ) {
572 12 : keys.add("compulsory","CONCENTRATION","the concentration parameter for the Von Mises-Fisher distributions");
573 12 : keys.add("compulsory","HEIGHTS","1.0","this keyword takes the label of an action that calculates a vector of values. The elements of this vector "
574 : "are used as weights for the Gaussians.");
575 12 : keys.add("compulsory","GRID_BIN","the number of points on the fibonacci sphere at which the density should be evaluated");
576 12 : }
577 :
578 10 : void SphericalKDEGridTools::readBandwidthAndHeight( const UniversalVonMisses& params, ActionWithArguments* action ) {
579 : // Read in the concentration parameters
580 : std::string von_misses_concentration;
581 10 : action->parse("CONCENTRATION",von_misses_concentration);
582 10 : KDEHelper<VonMissesKernelParams,UniversalVonMisses,SphericalKDEGridTools>::readKernelParameters( von_misses_concentration, action, "_vmconcentration", true );
583 10 : action->log.printf(" getting concentration parameters from %s \n", von_misses_concentration.c_str() );
584 : // Read in the heights
585 : std::string weight_str;
586 10 : action->parse("HEIGHTS",weight_str);
587 20 : KDEHelper<VonMissesKernelParams,UniversalVonMisses,SphericalKDEGridTools>::readKernelParameters( weight_str, action, "_volumes", false );
588 10 : if( (action->getPntrToArgument(0))->getNumberOfValues()==1 ) {
589 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: CUSTOM PERIODIC=NO ARG=" + weight_str + "," + action->getLabel() + "_vmconcentration FUNC=x*y/(4*pi*sinh(y))" ), false );
590 : } else {
591 10 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: CUSTOM PERIODIC=NO ARG=" + weight_str + "," + action->getLabel() + "_vmconcentration FUNC=x*y/(4*pi*sinh(y)) MASK=" + weight_str ), false );
592 : }
593 10 : KDEHelper<VonMissesKernelParams,UniversalVonMisses,SphericalKDEGridTools>::addArgument( action->getLabel() + "_heights", action );
594 10 : action->log.printf(" getting heights from %s \n", weight_str.c_str() );
595 10 : }
596 :
597 10 : void SphericalKDEGridTools::readGridParameters( SphericalKDEGridTools& g, ActionWithArguments* action, GridCoordinatesObject& gridobject, std::vector<std::size_t>& shape ) {
598 10 : if( shape.size()!=3 ) {
599 0 : action->error("should have three coordinates in input to this action");
600 : }
601 10 : action->parse("GRID_BIN",g.nbins);
602 10 : action->log.printf(" setting number of bins to %zu \n", g.nbins );
603 10 : std::vector<bool> ipbc( 3, false );
604 10 : gridobject.setup( "fibonacci", ipbc, g.nbins, 0 );
605 10 : shape[0]=g.nbins;
606 10 : shape[1]=shape[2]=1;
607 10 : }
608 :
609 10 : void SphericalKDEGridTools::getDiscreteSupport( const SphericalKDEGridTools& g, const UniversalVonMisses& p, const VonMissesKernelParams& kp, std::vector<unsigned>& nneigh, GridCoordinatesObject& gridobject ) {
610 10 : plumed_assert( nneigh.size()==gridobject.getDimension() );
611 10 : std::vector<bool> ipbc( 3, false );
612 10 : double fib_cutoff = std::log( epsilon / (kp.concentration/(4*pi*sinh(kp.concentration))) ) / kp.concentration;
613 20 : gridobject.setup( "fibonacci", ipbc, gridobject.getNumberOfPoints(), fib_cutoff );
614 10 : }
615 :
616 205624 : void SphericalKDEGridTools::getNeighbors( const UniversalVonMisses& p, const VonMissesKernelParams& kp, const GridCoordinatesObject& gridobject, const std::vector<unsigned>& nneigh, unsigned& num_neighbors, std::vector<unsigned>& neighbors ) {
617 205624 : gridobject.getNeighbors( kp.at, nneigh, num_neighbors, neighbors );
618 205624 : }
619 :
620 : typedef KDE<DiagonalKernelParams,DiscreteKernel,KDEGridTools<DiagonalKernelParams,DiscreteKernel>> discretekde;
621 : PLUMED_REGISTER_ACTION(discretekde,"KDE_DISCRETE")
622 : typedef KDE<DiagonalKernelParams,HistogramBeadKernel,KDEGridTools<DiagonalKernelParams,HistogramBeadKernel>> beadkde;
623 : PLUMED_REGISTER_ACTION(beadkde,"KDE_BEADS")
624 : typedef KDE<DiagonalKernelParams,RegularKernel<DiagonalKernelParams>,KDEGridTools<DiagonalKernelParams,RegularKernel<DiagonalKernelParams>>> flatkde;
625 : PLUMED_REGISTER_ACTION(flatkde,"KDE_KERNELS")
626 : typedef KDE<VonMissesKernelParams,UniversalVonMisses,SphericalKDEGridTools> sphericalkde;
627 : PLUMED_REGISTER_ACTION(sphericalkde,"SPHERICAL_KDE")
628 : typedef KDE<NonDiagonalKernelParams,RegularKernel<NonDiagonalKernelParams>,KDEGridTools<NonDiagonalKernelParams,RegularKernel<NonDiagonalKernelParams>>> flatfkde;
629 : PLUMED_REGISTER_ACTION(flatfkde,"KDE_FULLCOVAR")
630 :
631 :
632 : class KDEShortcut : public ActionShortcut {
633 : public:
634 : static void registerKeywords(Keywords& keys);
635 : explicit KDEShortcut(const ActionOptions&);
636 : };
637 :
638 : PLUMED_REGISTER_ACTION(KDEShortcut,"KDE")
639 :
640 142 : void KDEShortcut::registerKeywords(Keywords& keys) {
641 142 : KDE<DiagonalKernelParams,RegularKernel<DiagonalKernelParams>,KDEGridTools<DiagonalKernelParams,RegularKernel<DiagonalKernelParams>>>::registerKeywords( keys );
642 142 : keys.addActionNameSuffix("_DISCRETE");
643 142 : keys.addActionNameSuffix("_KERNELS");
644 142 : keys.addActionNameSuffix("_BEADS");
645 142 : keys.addActionNameSuffix("_FULLCOVAR");
646 142 : }
647 :
648 79 : KDEShortcut::KDEShortcut(const ActionOptions&ao):
649 : Action(ao),
650 79 : ActionShortcut(ao) {
651 : std::string kerneltype;
652 158 : parse("KERNEL",kerneltype);
653 79 : if( kerneltype=="DISCRETE" ) {
654 26 : readInputLine( getShortcutLabel() + ": KDE_DISCRETE " + convertInputLineToString() );
655 : return;
656 : }
657 : std::vector<std::string> args;
658 132 : parseVector("ARG", args );
659 : std::vector<Value*> argvals;
660 66 : ActionWithArguments::interpretArgumentList( args, plumed.getActionSet(), this, argvals );
661 66 : std::string argstr = " ARG=" + argvals[0]->getName();
662 98 : for(unsigned i=1; i<argvals.size(); ++i) {
663 64 : argstr += "," + argvals[i]->getName();
664 : }
665 : std::vector<std::string> bw;
666 132 : parseVector("BANDWIDTH",bw);
667 66 : std::string bwstr = " BANDWIDTH=" + bw[0];
668 98 : for(unsigned i=1; i<bw.size(); ++i) {
669 64 : bwstr += "," + bw[i];
670 : }
671 66 : if( bw.size() == 1 && argvals.size()>1 ) {
672 : std::vector<Value*> bwargs;
673 0 : ActionWithArguments::interpretArgumentList( bw, plumed.getActionSet(), this, bwargs );
674 0 : if( bwargs.size()!=1 ) {
675 0 : error("invalid input for bandwidth parameter");
676 0 : } else if( bwargs[0]->getRank()<=1 ) {
677 0 : if( kerneltype.find("bin")==std::string::npos ) {
678 0 : readInputLine( getShortcutLabel() + ": KDE_KERNELS " + argstr + " " + bwstr + " KERNEL=" + kerneltype + " " + convertInputLineToString() );
679 : } else {
680 0 : std::size_t dd = kerneltype.find("-bin");
681 0 : readInputLine( getShortcutLabel() + ": KDE_BEADS " + argstr + " " + bwstr + " KERNEL=" + kerneltype.substr(0,dd) + " " + convertInputLineToString() );
682 : }
683 0 : } else if( bwargs[0]->getRank()==2 ) {
684 0 : readInputLine( getShortcutLabel() + ": KDE_FULLCOVAR" + argstr + " " + bwstr + " KERNEL=" + kerneltype + " " + convertInputLineToString() );
685 : } else {
686 0 : error("found strange rank for bandwidth parameter");
687 : }
688 66 : } else if( bw.size()==argvals.size() ) {
689 66 : if( kerneltype.find("bin")==std::string::npos ) {
690 96 : readInputLine( getShortcutLabel() + ": KDE_KERNELS " + argstr + " " + bwstr + " KERNEL=" + kerneltype + " " + convertInputLineToString() );
691 : } else {
692 18 : std::size_t dd = kerneltype.find("-bin");
693 36 : readInputLine( getShortcutLabel() + ": KDE_BEADS " + argstr + " " + bwstr + " KERNEL=" + kerneltype.substr(0,dd) + " " + convertInputLineToString() );
694 : }
695 0 : } else if( bw.size()==argvals.size()*argvals.size() ) {
696 0 : readInputLine( getShortcutLabel() + ": KDE_FULLCOVAR" + argstr + " " + bwstr + " KERNEL=" + kerneltype + " " + convertInputLineToString() );
697 : } else {
698 0 : error("invalid input for bandwidth");
699 : }
700 132 : }
701 :
702 : }
703 : }
|
