Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2017 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 "core/ActionWithVector.h"
23 : #include "core/ParallelTaskManager.h"
24 : #include "core/ActionShortcut.h"
25 : #include "core/ActionRegister.h"
26 :
27 : //+PLUMEDOC MCOLVAR NEIGHBORS
28 : /*
29 : Build a matrix with ones in for the N nearest neighbours of an atom
30 :
31 : The following input illustrates how to use this action in tandem with [DISTANCE_MATRIX](DISTANCE_MATRIX.md) to find the six
32 : nearest atoms to each of the first 100 atoms in the input file:
33 :
34 : ```plumed
35 : d1: DISTANCE_MATRIX GROUP=1-100
36 : n: NEIGHBORS ARG=d1 NLOWEST=6
37 : ```
38 :
39 : Alternatively, if you would like to use a [CONTACT_MATRIX](CONTACT_MATRIX.md) to do something similar you would do the following:
40 :
41 : ```plumed
42 : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.5}
43 : n: NEIGHBORS ARG=c1 NHIGHEST=6
44 : ```
45 :
46 : This command is useful for implementing alternatives to the symmatry functions that are defined by the shortcuts
47 : in the module symfunc. For example, suppose that you want to calculate a variant on the [TETRAHEDRAL](TETRAHEDRAL.md) symmetry function.
48 : In this variant on the CV the coordination sphere around each central atom is not defined using a switching function. Instead
49 : this coordination sphere contains only the four nearest atoms. You can implement this CV by using the following input:
50 :
51 : ```plumed
52 : d1: DISTANCE_MATRIX GROUP=1-100 COMPONENTS
53 : n: NEIGHBORS ARG=d1.w NLOWEST=4
54 : f: CUSTOM ARG=n,d1.x,d1.y,d1.z,d1.w VAR=w,x,y,z,r PERIODIC=NO FUNC=w*(((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3)
55 : ones: ONES SIZE=100
56 : ucv: MATRIX_VECTOR_PRODUCT ARG=f,ones
57 : cv: CUSTOM ARG=ucv PERIODIC=NO FUNC=x/4
58 : ```
59 :
60 : You can even use these ideas with the functionality that is in the [volumes module](module_volumes.md) as shown below:
61 :
62 : ```plumed
63 : # The atoms that are of interest
64 : ow: GROUP ATOMS=1-16500
65 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
66 : center: FIXEDATOM AT=2.5,2.5,2.5
67 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
68 : sphere: INSPHERE ATOMS=ow CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
69 : # The distance matrix
70 : dmap: DISTANCE_MATRIX COMPONENTS GROUP=ow LINKCELL_CUTOFF=1.0 MASK=sphere
71 : # Find the four nearest neighbors
72 : acv_neigh: NEIGHBORS ARG=dmap.w NLOWEST=4 MASK=sphere
73 : # Compute a function for the atoms that are in the first coordination sphere
74 : acv_g8: GSYMFUNC_THREEBODY ...
75 : WEIGHT=acv_neigh ARG=dmap.x,dmap.y,dmap.z
76 : FUNCTION1={FUNC=(cos(ajik)+1/3)^2 LABEL=g8}
77 : MASK=sphere
78 : ...
79 : # Now compute the value of the function above for those atoms that are in the
80 : # sphere of interest
81 : acv: CUSTOM ARG=acv_g8.g8,sphere FUNC=y*(1-(3*x/8)) PERIODIC=NO
82 : # And now compute the final average
83 : acv_sum: SUM ARG=acv PERIODIC=NO
84 : acv_norm: SUM ARG=sphere PERIODIC=NO
85 : mean: CUSTOM ARG=acv_sum,acv_norm FUNC=x/y PERIODIC=NO
86 : PRINT ARG=mean FILE=colvar
87 : ```
88 :
89 : */
90 : //+ENDPLUMEDOC
91 :
92 : namespace PLMD {
93 : namespace adjmat {
94 :
95 : class NeighborsShortcut : public ActionShortcut {
96 : public:
97 : static void registerKeywords( Keywords& keys );
98 : explicit NeighborsShortcut(const ActionOptions&);
99 : };
100 :
101 : PLUMED_REGISTER_ACTION(NeighborsShortcut,"NEIGHBORS")
102 :
103 21 : void NeighborsShortcut::registerKeywords( Keywords& keys ) {
104 21 : ActionShortcut::registerKeywords( keys );
105 42 : keys.addInputKeyword("optional","MASK","vector","a vector that is used to used to determine which rows of the neighbors matrix to compute");
106 42 : keys.addInputKeyword("compulsory","ARG","matrix","the label of an adjacency/distance matrix that will be used to find the nearest neighbors");
107 21 : keys.add("compulsory","NLOWEST","0","in each row of the output matrix set the elements that correspond to the n lowest elements in each row of the input matrix equal to one");
108 21 : keys.add("compulsory","NHIGHEST","0","in each row of the output matrix set the elements that correspond to the n highest elements in each row of the input matrix equal to one");
109 42 : keys.setValueDescription("matrix","a matrix in which the ij element is one if the ij-element of the input matrix is one of the NLOWEST/NHIGHEST elements on that row of the input matrix and zero otherwise");
110 21 : keys.addActionNameSuffix("_1LOW");
111 21 : keys.addActionNameSuffix("_NLOW");
112 21 : keys.addActionNameSuffix("_1HIGH");
113 21 : keys.addActionNameSuffix("_NHIGH");
114 21 : }
115 :
116 9 : NeighborsShortcut::NeighborsShortcut(const ActionOptions& ao):
117 : Action(ao),
118 9 : ActionShortcut(ao) {
119 : unsigned nlow;
120 9 : parse("NLOWEST",nlow);
121 : unsigned nhigh;
122 9 : parse("NHIGHEST",nhigh);
123 9 : if( nlow==0 && nhigh==0 ) {
124 0 : error("missing NLOWEST or NHIGHEST keyword one of these two keywords must be set in input");
125 : }
126 9 : if( nlow>0 && nhigh>0 ) {
127 0 : error("should only be one of NLOWEST or NHIGHEST set in input");
128 : }
129 9 : if( nlow==1 ) {
130 6 : log.printf(" output matrix will have non-zero values for elements that correspond to the %d lowest elements in each row of the input matrix\n",nlow);
131 12 : readInputLine( getShortcutLabel() + ": NEIGHBORS_1LOW N=1 " + convertInputLineToString() );
132 3 : } else if( nlow>1 ) {
133 : std::string str_n;
134 3 : Tools::convert( nlow, str_n );
135 3 : log.printf(" output matrix will have non-zero values for elements that correspond to the %d lowest elements in each row of the input matrix\n",nlow);
136 9 : readInputLine( getShortcutLabel() + ": NEIGHBORS_NLOW N=" + str_n + " " + convertInputLineToString() );
137 0 : } else if( nhigh==1 ) {
138 0 : log.printf(" output matrix will have non-zero values for elements that correspond to the %d highest elements in each row of the input matrix\n",nhigh);
139 0 : readInputLine( getShortcutLabel() + ": NEIGHBORS_1HIGH N=1 " + convertInputLineToString() );
140 0 : } else if( nhigh>1 ) {
141 : std::string str_n;
142 0 : Tools::convert( nlow, str_n );
143 0 : log.printf(" output matrix will have non-zero values for elements that correspond to the %d highest elements in each row of the input matrix\n",nhigh);
144 0 : readInputLine( getShortcutLabel() + ": NEIGHBORS_NHIGH N=" + str_n + " " + convertInputLineToString() );
145 : } else {
146 0 : error("do not know what I am supposed to do");
147 : }
148 9 : }
149 :
150 : class NeighborCalcInput {
151 : public:
152 : unsigned number;
153 : unsigned nind;
154 : unsigned ncols;
155 : unsigned nbonds;
156 : View<const double> matrow;
157 : View<const std::size_t> bookrow;
158 3548 : NeighborCalcInput( unsigned task_index,
159 : unsigned n,
160 : const ArgumentBookeepingHolder& arg,
161 3548 : double* d ):
162 3548 : number(n),
163 3548 : nind(0),
164 3548 : ncols(arg.ncols),
165 3548 : nbonds(arg.bookeeping[task_index*(1+ncols)]),
166 3548 : matrow( d + task_index*ncols, nbonds ),
167 3548 : bookrow( arg.bookeeping.data() + task_index*(1+ncols)+1, nbonds ) {
168 431553 : for(unsigned i=0; i<nbonds; ++i) {
169 428005 : if( matrow[i]<epsilon ) {
170 517 : continue ;
171 : }
172 427488 : nind++;
173 : }
174 3548 : if( number>nind ) {
175 0 : plumed_merror("not enough matrix elements were stored");
176 : }
177 3548 : }
178 : void getSortedData( std::vector<std::pair<double,unsigned> >& rows ) const ;
179 : };
180 :
181 2512 : void NeighborCalcInput::getSortedData( std::vector<std::pair<double,unsigned> >& rows ) const {
182 : unsigned n=0;
183 180378 : for(unsigned i=0; i<nbonds; ++i) {
184 177866 : if( matrow[i]<epsilon ) {
185 0 : continue ;
186 : }
187 177866 : rows[n].first=matrow[i];
188 177866 : rows[n].second=bookrow[i];
189 177866 : n++;
190 : }
191 2512 : std::sort( rows.begin(), rows.end() );
192 2512 : }
193 :
194 : class OneLowInput {
195 : public:
196 : unsigned number;
197 : static void calculate( const NeighborCalcInput& input, View<double>& output );
198 : };
199 :
200 1036 : void OneLowInput::calculate( const NeighborCalcInput& input, View<double>& output ) {
201 1036 : unsigned nv = input.bookrow[0];
202 1036 : double min = input.matrow[0];
203 250139 : for(unsigned i=1; i<input.nbonds; ++i) {
204 249103 : if( input.matrow[i]<min ) {
205 : min = input.matrow[i];
206 5047 : nv = input.bookrow[i];
207 : }
208 : }
209 1036 : output[0] = nv;
210 1036 : }
211 :
212 : class NLowInput {
213 : public:
214 : unsigned number;
215 : static void calculate( const NeighborCalcInput& input, View<double>& output );
216 : };
217 :
218 2512 : void NLowInput::calculate( const NeighborCalcInput& input, View<double>& output ) {
219 2512 : std::vector<std::pair<double,unsigned> > rows( input.nind );
220 2512 : input.getSortedData( rows );
221 :
222 16560 : for(unsigned i=0; i<input.number; ++i) {
223 14048 : output[i] = rows[i].second;
224 : }
225 2512 : }
226 :
227 : class OneHighInput {
228 : public:
229 : unsigned number;
230 : static void calculate( const NeighborCalcInput& input, View<double>& output );
231 : };
232 :
233 0 : void OneHighInput::calculate( const NeighborCalcInput& input, View<double>& output ) {
234 0 : unsigned nv = input.bookrow[0];
235 0 : double max = input.matrow[0];
236 0 : for(unsigned i=1; i<input.nbonds; ++i) {
237 0 : if( input.matrow[i]>max ) {
238 : max = input.matrow[i];
239 0 : nv = input.bookrow[i];
240 : }
241 : }
242 0 : output[0] = nv;
243 0 : }
244 :
245 : class NHighInput {
246 : public:
247 : unsigned number;
248 : static void calculate( const NeighborCalcInput& input, View<double>& output );
249 : };
250 :
251 0 : void NHighInput::calculate( const NeighborCalcInput& input, View<double>& output ) {
252 0 : std::vector<std::pair<double,unsigned> > rows( input.nind );
253 0 : input.getSortedData( rows );
254 :
255 0 : for(unsigned i=0; i<input.number; ++i) {
256 0 : output[i] = rows[input.nind-1-i].second;
257 : }
258 0 : }
259 :
260 : template <class T>
261 : class Neighbors : public ActionWithVector {
262 : public:
263 : using input_type = T;
264 : using PTM = ParallelTaskManager<Neighbors<T>>;
265 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
266 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
267 : private:
268 : /// The parallel task manager
269 : PTM taskmanager;
270 : public:
271 : static void registerKeywords( Keywords& keys );
272 : explicit Neighbors(const ActionOptions&);
273 : unsigned getNumberOfDerivatives() override;
274 : void turnOnDerivatives() override ;
275 : void transferStashToValues( const std::vector<unsigned>& partialTaskList, const std::vector<double>& stash ) override ;
276 : void prepare() override ;
277 : void calculate() override ;
278 : static void performTask( std::size_t task_index,
279 : const T& actiondata,
280 : ParallelActionsInput& input,
281 : ParallelActionsOutput& output );
282 : };
283 :
284 : typedef Neighbors<OneLowInput> olinp;
285 : PLUMED_REGISTER_ACTION(olinp,"NEIGHBORS_1LOW")
286 : typedef Neighbors<NLowInput> nlinp;
287 : PLUMED_REGISTER_ACTION(nlinp,"NEIGHBORS_NLOW")
288 : typedef Neighbors<OneHighInput> ohinp;
289 : PLUMED_REGISTER_ACTION(ohinp,"NEIGHBORS_1HIGH")
290 : typedef Neighbors<NHighInput> nhinp;
291 : PLUMED_REGISTER_ACTION(nhinp,"NEIGHBORS_NHIGH")
292 :
293 : template <class T>
294 34 : void Neighbors<T>::registerKeywords( Keywords& keys ) {
295 34 : ActionWithVector::registerKeywords( keys );
296 34 : keys.setDisplayName("NEIGHBORS");
297 68 : keys.addInputKeyword("optional","MASK","vector","a vector that is used to used to determine which rows of the neighbors matrix to compute");
298 68 : keys.addInputKeyword("compulsory","ARG","matrix","the label of an adjacency/distance matrix that will be used to find the nearest neighbors");
299 34 : keys.add("compulsory","N","the number of non-zero elements in each row of the output matrix");
300 68 : keys.setValueDescription("matrix","a matrix in which the ij element is one if the ij-element of the input matrix is one of the NLOWEST/NHIGHEST elements on that row of the input matrix and zero otherwise");
301 34 : PTM::registerKeywords( keys );
302 34 : }
303 :
304 : template <class T>
305 9 : Neighbors<T>::Neighbors(const ActionOptions&ao):
306 : Action(ao),
307 : ActionWithVector(ao),
308 9 : taskmanager(this) {
309 : unsigned nargs = getNumberOfArguments();
310 9 : if( getNumberOfMasks()>0 ) {
311 0 : nargs = nargs - getNumberOfMasks();
312 : }
313 9 : if( nargs!=1 ) {
314 0 : error("found wrong number of arguments in input");
315 : }
316 9 : if( getPntrToArgument(0)->getRank()!=2 ) {
317 0 : error("input argument should be a matrix");
318 : }
319 :
320 : T myinp;
321 18 : parse("N",myinp.number);
322 : // And get the shape
323 9 : std::vector<std::size_t> shape( getPntrToArgument(0)->getShape() );
324 9 : addValue( shape );
325 9 : setNotPeriodic();
326 9 : getPntrToComponent(0)->reshapeMatrixStore( myinp.number );
327 9 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
328 : // Setup the parallel task manager
329 9 : taskmanager.setupParallelTaskManager( 0, 0 );
330 9 : taskmanager.setActionInput( myinp );
331 9 : }
332 :
333 : template <class T>
334 0 : void Neighbors<T>::turnOnDerivatives() {
335 0 : error("If the symmetry function can be calculated from distances only then the derivatives of a function that is computed with NEIGHBORS will be continuous. If you calculate angles between vectors or use the vectors directly then the symmetry function computed using NEIGHBORS is not continuous. Out of an abundance of caution we thus forbid forces for this action. If you are interested on applying forces with this action email gareth.tribello@gmail.com");
336 : }
337 :
338 : template <class T>
339 0 : unsigned Neighbors<T>::getNumberOfDerivatives() {
340 0 : return 0;
341 : }
342 :
343 : template <class T>
344 19 : void Neighbors<T>::prepare() {
345 19 : ActionWithVector::prepare();
346 19 : Value* myval = getPntrToComponent(0);
347 19 : if( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0]
348 19 : && myval->getShape()[1]==getPntrToArgument(0)->getShape()[1] ) {
349 13 : return;
350 : }
351 6 : std::vector<std::size_t> shape( getPntrToArgument(0)->getShape() );
352 6 : myval->setShape(shape);
353 6 : myval->reshapeMatrixStore( taskmanager.getActionInput().number );
354 6 : }
355 :
356 : template <class T>
357 19 : void Neighbors<T>::calculate() {
358 19 : taskmanager.runAllTasks();
359 19 : }
360 :
361 : template <class T>
362 3548 : void Neighbors<T>::performTask( std::size_t task_index,
363 : const T& actiondata,
364 : ParallelActionsInput& input,
365 : ParallelActionsOutput& output ) {
366 3548 : T::calculate( NeighborCalcInput( task_index,
367 3548 : actiondata.number,
368 7096 : ArgumentBookeepingHolder::create( 0, input ),
369 : input.inputdata ),
370 3548 : output.values );
371 3548 : }
372 :
373 : template <class T>
374 19 : void Neighbors<T>::transferStashToValues( const std::vector<unsigned>& partialTaskList, const std::vector<double>& stash ) {
375 19 : Value* myval = getPntrToComponent(0);
376 : // All values are set equal to one
377 15103 : for(unsigned i=0; i<myval->getNumberOfStoredValues(); ++i) {
378 15084 : myval->set( i, 1 );
379 : }
380 :
381 : // And we set the bookeeping data from the data in the stash
382 : unsigned k=0;
383 19 : std::vector<std::size_t> rowindices( taskmanager.getActionInput().number );
384 3567 : for(unsigned i=0; i<myval->getShape()[0]; ++i) {
385 18632 : for(unsigned j=0; j<rowindices.size(); ++j) {
386 15084 : rowindices[j] = static_cast<std::size_t>( stash[k] );
387 15084 : ++k;
388 : }
389 3548 : myval->setRowIndices( i, rowindices );
390 : }
391 19 : }
392 :
393 : }
394 : }
|
