Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2020 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 : #ifdef __PLUMED_HAS_OPENACC
23 : #define __PLUMED_USE_OPENACC 1
24 : #endif //__PLUMED_HAS_OPENACC
25 : #include "core/ActionRegister.h"
26 : #include "tools/Pbc.h"
27 : #include "tools/HistogramBead.h"
28 : #include "ActionVolume.h"
29 : #include "VolumeShortcut.h"
30 :
31 : //+PLUMEDOC VOLUMES AROUND
32 : /*
33 : This quantity can be used to calculate functions of the distribution of collective variables for the atoms that lie in a particular, user-specified part of of the cell.
34 :
35 : This action can be used to calculate whether each of the atoms are within a particular part of the simulation box or not as illustrated by the following example:
36 :
37 : ```plumed
38 : f: FIXEDATOM AT=0,0,0
39 : a: AROUND ...
40 : ATOMS=1-100 ORIGIN=f
41 : SIGMA=0.2 KERNEL=gaussian
42 : XLOWER=-1.0 XUPPER=1.0
43 : YLOWER=-1.0 YUPPER=1.0
44 : ZLOWER=-1.0 ZUPPER=1.0
45 : ...
46 : PRINT ARG=a FILE=colvar
47 : ```
48 :
49 : The 100 elements of the vector `a` that is returned from the AROUND action in the above input are calculated using:
50 :
51 : $$
52 : w(x_i,y_i,z_i) = \int_{xl}^{xu} \int_{yl}^{yu} \int_{zl}^{zu} \textrm{d}x\textrm{d}y\textrm{d}z K\left( \frac{x - x_i}{\sigma} \right)K\left( \frac{y - y_i}{\sigma} \right)K\left( \frac{z - z_i}{\sigma} \right)
53 : $$
54 :
55 : where $K$ is one of the kernel functions described in the documentation for the function [BETWEEN](BETWEEN.md), $\sigma$ is a bandwidth parameter and the limits
56 : for the integrals are the values specified using the keywords XLOWER, XUPPER, YLOWER, YUPPER, YUPPER, ZLOWER and ZUPPER. $x_i$, $y_i$ and $z_i$ are then the components
57 : of the vector that connects the $i$th atom that was specified using the ATOMS keyword to the atom that was specified using the ORIGIN keyword. In other words,
58 : $w(x_i,y_i,z_i)$ is 1 if the atom is within a rectangular box that is centered on the atom that is specified as the origin and zero otherwise.
59 :
60 : If instead of calculating if the atoms are inside this box you want to calculate if they are outside this box you can use the following input:
61 :
62 : ```plumed
63 : f: FIXEDATOM AT=0,0,0
64 : a: AROUND ...
65 : ATOMS=1-100 ORIGIN=f
66 : SIGMA=0.2 KERNEL=gaussian
67 : XLOWER=-1.0 XUPPER=1.0
68 : YLOWER=-1.0 YUPPER=1.0
69 : ZLOWER=-1.0 ZUPPER=1.0
70 : OUTSIDE
71 : ...
72 : PRINT ARG=a FILE=colvar
73 : ```
74 :
75 : The 100 elements of the vector `a` that is returned from the AROUND action in the above input are calculated using:
76 :
77 : $$
78 : v(x_i,y_i,z_i) = 1 - w(x_i,y_i,z_i)
79 : $$
80 :
81 : ## Calculating the number of atoms in a particular part of the box
82 :
83 : Lets suppose that you want to calculate how many atoms are in have an $x$ coordinate that is between -1.0 and 1.0. You can do this using the following PLUMED input:
84 :
85 : ```plumed
86 : f: FIXEDATOM AT=0,0,0
87 : a: AROUND ATOMS=1-100 ORIGIN=f SIGMA=0.2 XLOWER=-1.0 XUPPER=1.0
88 : s: SUM ARG=a PERIODIC=NO
89 : PRINT ARG=s FILE=colvar
90 : ```
91 :
92 : In this example the components of `a` are calculated as:
93 :
94 : $$
95 : w(x_i,y_i,z_i) = \int_{xl}^{xu} \textrm{d}x K\left( \frac{x - x_i}{\sigma} \right)
96 : $$
97 :
98 : as the YLOWER, YUPPER, YUPPER, ZLOWER and ZUPPER flags have not been included. The [SUM](SUM.md) command then adds together all the elements of the vector `a` to calculate the total number of atoms in the region
99 : of the box that is of interest.
100 :
101 : ## Calculating the average value for an order parameter in a particular part of the box
102 :
103 : Suppose that you have calculated a vector of order parameters that can be assigned to a particular point in three dimensional space.
104 : The symmetry functions in the [symfunc](module_symfunc.md) module are examples of order parameters that satisfy this criteria. You can use
105 : the AROUND command to calculate the average value of the symmetry function in a particular part of the box as follows:
106 :
107 : ```plumed
108 : f: FIXEDATOM AT=0,0,0
109 : a: AROUND ...
110 : ATOMS=1-100 ORIGIN=f SIGMA=0.2
111 : XLOWER=-1.0 XUPPER=1.0
112 : YLOWER=-1.0 YUPPER=1.0
113 : ...
114 :
115 : c: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1.0} MASK=a
116 : p: CUSTOM ARG=c,a FUNC=x*y PERIODIC=NO
117 : n: SUM ARG=p PERIODIC=NO
118 : d: SUM ARG=a PERIODIC=NO
119 : av: CUSTOM ARG=n,d FUNC=x/y PERIODIC=NO
120 : PRINT ARG=av FILE=colvar
121 : ```
122 :
123 : The final quantity `av` here is:
124 :
125 : $$
126 : \overline{s}_{\tau} = \frac{ \sum_i c_i w(x_i,y_i,z_i) }{ \sum_i w(x_i,y_i,z_i) }
127 : $$
128 :
129 : where $c_i$ are the coordination numbers and $w_i$ is:
130 :
131 : $$
132 : w(x_i,y_i,z_i) = \int_{xl}^{xu} \int_{yl}^{yu} \textrm{d}x \textrm{d}y K\left( \frac{x - x_i}{\sigma} \right) K\left( \frac{y - y_i}{\sigma} \right)
133 : $$
134 :
135 : ## Old syntax
136 :
137 : In earlier versions of PLUMED the syntax for the calculation in the previous section is as follows:
138 :
139 : ```plumed
140 : c: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1.0}
141 : f: FIXEDATOM AT=0,0,0
142 : a: AROUND ...
143 : DATA=c ORIGIN=f SIGMA=0.2
144 : XLOWER=-1.0 XUPPER=1.0
145 : YLOWER=-1.0 YUPPER=1.0
146 : MEAN
147 : ...
148 : PRINT ARG=a.mean FILE=colvar
149 : ```
150 :
151 : This old syntax still works but __we highly recommend you use the newer syntax__ as it is easlier to understand,
152 : more flexible and calculations with this newer input will run faster. You will notice that AROUND in the input above
153 : is a shortcut that expands to a longer input
154 : that is similar to that given in the previous section. The main difference is that the order of the AROUND
155 : and [COORDINATIONNUMBER](COORDINATIONNUMBER.md) actions is reversed in the new syntax. The reason this reversal is necessary
156 : is that the vector output by AROUND must be passed as as a MASK action to the [COORDINATIONNUMBER](COORDINATIONNUMBER.md)
157 : action in order to optimize code performance. Passing the vector from AROUND as a MASK in coordination number ensures that
158 : we only calculate the coordination numbers for those atomms that are in the region of interest. We thus avoid a lot of computational
159 : expense that would otherwise be associated with calculating coordination numbers for atoms that are not within the region of
160 : interest and would thus make no difference to the final average that we are calculating.
161 :
162 : The old syntax also allowed you to compute the sum of the coordination numbers in the region of interest using an input like the one shown below:
163 :
164 : ```plumed
165 : c: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1.0}
166 : f: FIXEDATOM AT=0,0,0
167 : a: AROUND ...
168 : DATA=c ORIGIN=f SIGMA=0.2
169 : XLOWER=-1.0 XUPPER=1.0
170 : YLOWER=-1.0 YUPPER=1.0
171 : SUM
172 : ...
173 : PRINT ARG=a.sum FILE=colvar
174 : ```
175 :
176 : The final CV that is computed here is:
177 :
178 : $$
179 : \overline{s}_{\tau} = \sum_i c_i w(x_i,y_i,z_i)
180 : $$
181 :
182 : the equivalent input with the new syntax is thus:
183 :
184 : ```plumed
185 : f: FIXEDATOM AT=0,0,0
186 : a: AROUND ...
187 : ATOMS=1-100 ORIGIN=f SIGMA=0.2
188 : XLOWER=-1.0 XUPPER=1.0
189 : YLOWER=-1.0 YUPPER=1.0
190 : ...
191 :
192 : c: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1.0} MASK=a
193 : p: CUSTOM ARG=c,a FUNC=x*y PERIODIC=NO
194 : n: SUM ARG=p PERIODIC=NO
195 : PRINT ARG=n FILE=colvar
196 : ```
197 :
198 : That old syntax also allowed you to accumulate quantities such as:
199 :
200 : $$
201 : \overline{s}_{\tau} = \sum_i f(c_i) w(x_i,y_i,z_i)
202 : $$
203 :
204 : where $f$ could be one of the switching functions discussed in the documentation for [LESS_THAN](LESS_THAN.md),
205 : one of the reverse switching functions discussed in the documentation for [MORE_THAN](MORE_THAN.md) or one of the
206 : two sided switching functions discussed in the documentation for [BETWEEN](BETWEEN.md). An example of an old input
207 : that computes all three of three types of sum is shown below:
208 :
209 : ```plumed
210 : c: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1.0}
211 : f: FIXEDATOM AT=0,0,0
212 : a: AROUND ...
213 : DATA=c ORIGIN=f SIGMA=0.2
214 : XLOWER=-1.0 XUPPER=1.0
215 : YLOWER=-1.0 YUPPER=1.0
216 : LESS_THAN={RATIONAL R_0=3}
217 : MORE_THAN={RATIONAL R_0=6}
218 : BETWEEN={GAUSSIAN LOWER=3 UPPER=6 SMEAR=0.5}
219 : ...
220 : PRINT ARG=a.lessthan,a.between,a.morethan FILE=colvar
221 : ```
222 :
223 : With the new syntax we can achieve the same result using the following input:
224 :
225 : ```plumed
226 : f: FIXEDATOM AT=0,0,0
227 : a: AROUND ...
228 : ATOMS=1-100 ORIGIN=f SIGMA=0.2
229 : XLOWER=-1.0 XUPPER=1.0
230 : YLOWER=-1.0 YUPPER=1.0
231 : ...
232 :
233 : c: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1.0} MASK=a
234 : # This part does the LESS_THAN={RATIONAL R_0=3}
235 : lt: LESS_THAN ARG=c SWITCH={RATIONAL R_0=3}
236 : wlt: CUSTOM ARG=a,lt FUNC=x*y PERIODIC=NO
237 : lessthan: SUM ARG=wlt PERIODIC=NO
238 : # This part does the BETWEEN={GAUSSIAN LOWER=3 UPPER=6 SMEAR=0.5}
239 : bt: BETWEEN ARG=c SWITCH={GAUSSIAN LOWER=3 UPPER=6 SMEAR=0.5}
240 : wbt: CUSTOM ARG=a,bt FUNC=x*y PERIODIC=NO
241 : between: SUM ARG=wbt PERIODIC=NO
242 : # This part does the MORE_THAN={RATIONAL R_0=6}
243 : mt: MORE_THAN ARG=c SWITCH={RATIONAL R_0=6}
244 : wmt: CUSTOM ARG=a,mt FUNC=x*y PERIODIC=NO
245 : morethan: SUM ARG=wmt PERIODIC=NO
246 : PRINT ARG=lessthan,between,morethan FILE=colvar
247 : ```
248 :
249 : */
250 : //+ENDPLUMEDOC
251 :
252 : namespace PLMD {
253 : namespace volumes {
254 :
255 98 : class VolumeAround {
256 : public:
257 : bool dox{true}, doy{true}, doz{true};
258 : double sigma;
259 : double xlow{0.0}, xhigh{0.0};
260 : double ylow{0.0}, yhigh{0.0};
261 : double zlow{0.0}, zhigh{0.0};
262 : HistogramBead::KernelType kerneltype;
263 : static void registerKeywords( Keywords& keys );
264 : void parseInput( ActionVolume<VolumeAround>* action );
265 : void setupRegions( ActionVolume<VolumeAround>* action,
266 : const Pbc& pbc,
267 : const std::vector<Vector>& positions ) {}
268 : static void parseAtoms( ActionVolume<VolumeAround>* action,
269 : std::vector<AtomNumber>& atom );
270 : static void calculateNumberInside( const VolumeInput& input,
271 : const VolumeAround& actioninput,
272 : VolumeOutput& output );
273 : #ifdef __PLUMED_USE_OPENACC
274 : void toACCDevice() const {
275 : #pragma acc enter data copyin(this[0:1],dox,doy,doz,sigma,\
276 : xlow,xhigh,ylow,yhigh,zlow,zhigh,kerneltype)
277 : }
278 : void removeFromACCDevice() const {
279 : #pragma acc exit data delete(kerneltype,zhigh,zlow,yhigh,ylow,xhigh,xlow,\
280 : sigma,doz,doy,dox,this[0:1])
281 : }
282 : #endif //__PLUMED_USE_OPENACC
283 : };
284 :
285 : typedef ActionVolume<VolumeAround> Vola;
286 : PLUMED_REGISTER_ACTION(Vola,"AROUND_CALC")
287 : char glob_around[] = "AROUND";
288 : typedef VolumeShortcut<glob_around> VolumeAroundShortcut;
289 : PLUMED_REGISTER_ACTION(VolumeAroundShortcut,"AROUND")
290 :
291 165 : void VolumeAround::registerKeywords( Keywords& keys ) {
292 330 : keys.setDisplayName("AROUND");
293 165 : keys.add("atoms","ORIGIN","the atom whose vicinity we are interested in examining");
294 165 : keys.addDeprecatedKeyword("ATOM","ORIGIN");
295 165 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
296 165 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
297 165 : keys.add("compulsory","XLOWER","0.0","the lower boundary in x relative to the x coordinate of the atom (0 indicates use full extent of box).");
298 165 : keys.add("compulsory","XUPPER","0.0","the upper boundary in x relative to the x coordinate of the atom (0 indicates use full extent of box).");
299 165 : keys.add("compulsory","YLOWER","0.0","the lower boundary in y relative to the y coordinate of the atom (0 indicates use full extent of box).");
300 165 : keys.add("compulsory","YUPPER","0.0","the upper boundary in y relative to the y coordinate of the atom (0 indicates use full extent of box).");
301 165 : keys.add("compulsory","ZLOWER","0.0","the lower boundary in z relative to the z coordinate of the atom (0 indicates use full extent of box).");
302 165 : keys.add("compulsory","ZUPPER","0.0","the upper boundary in z relative to the z coordinate of the atom (0 indicates use full extent of box).");
303 165 : }
304 :
305 49 : void VolumeAround::parseAtoms( ActionVolume<VolumeAround>* action, std::vector<AtomNumber>& atom ) {
306 98 : action->parseAtomList("ORIGIN",atom);
307 49 : if( atom.size()==0 ) {
308 0 : action->parseAtomList("ATOM",atom);
309 : }
310 49 : if( atom.size()!=1 ) {
311 0 : action->error("should only be one atom specified");
312 : }
313 49 : action->log.printf(" boundaries for region are calculated based on positions of atom : %d\n",atom[0].serial() );
314 49 : }
315 :
316 49 : void VolumeAround::parseInput( ActionVolume<VolumeAround>* action ) {
317 98 : action->parse("SIGMA",sigma);
318 : std::string mykerneltype;
319 49 : action->parse("KERNEL",mykerneltype);
320 49 : kerneltype=HistogramBead::getKernelType(mykerneltype);
321 49 : dox=true;
322 49 : action->parse("XLOWER",xlow);
323 49 : action->parse("XUPPER",xhigh);
324 49 : doy=true;
325 49 : action->parse("YLOWER",ylow);
326 49 : action->parse("YUPPER",yhigh);
327 49 : doz=true;
328 49 : action->parse("ZLOWER",zlow);
329 49 : action->parse("ZUPPER",zhigh);
330 49 : if( xlow==0.0 && xhigh==0.0 ) {
331 8 : dox=false;
332 : }
333 49 : if( ylow==0.0 && yhigh==0.0 ) {
334 16 : doy=false;
335 : }
336 49 : if( zlow==0.0 && zhigh==0.0 ) {
337 16 : doz=false;
338 : }
339 49 : if( !dox && !doy && !doz ) {
340 0 : action->error("no subregion defined use XLOWER, XUPPER, YLOWER, YUPPER, ZLOWER, ZUPPER");
341 : }
342 49 : action->log.printf(" boundaries for region (region of interest about atom) : x %f %f, y %f %f, z %f %f \n",xlow,xhigh,ylow,yhigh,zlow,zhigh);
343 49 : }
344 :
345 140716 : void VolumeAround::calculateNumberInside( const VolumeInput& input,
346 : const VolumeAround& actioninput,
347 : VolumeOutput& output ) {
348 : // Setup the histogram bead
349 140716 : HistogramBead bead(actioninput.kerneltype, actioninput.xlow, actioninput.xhigh, actioninput.sigma );
350 :
351 : // Calculate position of atom wrt to origin
352 140716 : Vector fpos=input.pbc.distance( Vector(input.refpos[0][0],input.refpos[0][1],input.refpos[0][2]),
353 140716 : Vector(input.cpos[0],input.cpos[1],input.cpos[2]) );
354 : double xcontr=1.0;
355 140716 : double xder=0.0;
356 140716 : if( actioninput.dox ) {
357 : //bead parameters set in the constructor
358 96716 : xcontr=bead.calculate( fpos[0], xder );
359 : }
360 : double ycontr=1.0;
361 140716 : double yder=0.0;
362 140716 : if( actioninput.doy ) {
363 92288 : bead.set( actioninput.ylow, actioninput.yhigh, actioninput.sigma );
364 92288 : ycontr=bead.calculate( fpos[1], yder );
365 : }
366 : double zcontr=1.0;
367 140716 : double zder=0.0;
368 140716 : if( actioninput.doz ) {
369 92288 : bead.set( actioninput.zlow, actioninput.zhigh, actioninput.sigma );
370 92288 : zcontr=bead.calculate( fpos[2], zder );
371 : }
372 :
373 140716 : output.derivatives[0]=xder*ycontr*zcontr;
374 140716 : output.derivatives[1]=xcontr*yder*zcontr;
375 140716 : output.derivatives[2]=xcontr*ycontr*zder;
376 : // Add derivatives wrt to position of origin atom
377 140716 : output.refders[0][0] = -output.derivatives[0];
378 140716 : output.refders[0][1] = -output.derivatives[1];
379 140716 : output.refders[0][2] = -output.derivatives[2];
380 : // Add virial contribution
381 140716 : output.virial.set( 0, -Tensor(fpos,
382 140716 : Vector(output.derivatives[0], output.derivatives[1], output.derivatives[2])) );
383 140716 : output.values[0] = xcontr*ycontr*zcontr;
384 140716 : }
385 :
386 : }
387 : }
|
