Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 : #ifndef __PLUMED_adjmat_AdjacencyMatrixBase_h
23 : #define __PLUMED_adjmat_AdjacencyMatrixBase_h
24 :
25 : #include <vector>
26 : #include "core/ParallelTaskManager.h"
27 : #include "core/ActionWithMatrix.h"
28 : #include "core/PlumedMain.h"
29 : #include "tools/LinkCells.h"
30 :
31 : namespace PLMD {
32 : namespace adjmat {
33 :
34 : struct ContactMatrix;
35 :
36 : template <class T>
37 92 : struct AdjacencyMatrixData {
38 : T matrixdata;
39 : bool usepbc{true};
40 : bool components{false};
41 : bool no_thread_gather{false};
42 : bool gatherForceOnColumns{false};
43 : unsigned fbsize{0};
44 : RequiredMatrixElements outmat;
45 : std::size_t nlists{0};
46 : unsigned natoms_per_three_list{0};
47 : std::vector<std::size_t> nlist_three_v;
48 : std::size_t* nlist_three{nullptr};
49 : void update() {
50 11828 : nlist_three=nlist_three_v.data();
51 : }
52 : #ifdef __PLUMED_HAS_OPENACC
53 : void toACCDevice() const {
54 : #pragma acc enter data copyin(this[0:1],usepbc,components,no_thread_gather,gatherForceOnColumns,fbsize, \
55 : nlists,natoms_per_three_list, \
56 : nlist_three[0:nlist_three_v.size()])
57 : matrixdata.toACCDevice();
58 : outmat.toACCDevice();
59 : }
60 : void removeFromACCDevice() const {
61 : matrixdata.removeFromACCDevice();
62 : outmat.removeFromACCDevice();
63 : #pragma acc exit data delete(nlist_three[0:nlist_three_v.size()], \
64 : natoms_per_three_list,nlists, \
65 : fbsize,gatherForceOnColumns,no_thread_gather,components,usepbc,this[0:1])
66 : }
67 : #endif //__PLUMED_HAS_OPENACC
68 : };
69 :
70 : struct AdjacencyMatrixInput {
71 : bool noderiv{false};
72 : const Pbc* pbc;
73 : Vector pos;
74 : std::size_t natoms{0};
75 : VectorView extra_positions;
76 : };
77 :
78 : struct MatrixOutput {
79 : View<double,1> val;
80 : View<double> deriv;
81 :
82 : ///doing t= Tensor(v1,v2); deriv[x:x+9]=t with no extra memory allocation
83 : template <typename Iterable1, typename Iterable2>
84 11490856 : void assignOuterProduct(const std::size_t startingIndex,
85 : const Iterable1& v1,
86 : const Iterable2& v2 ) {
87 11490856 : deriv[startingIndex + 0] = v1[0] * v2[0];
88 11490856 : deriv[startingIndex + 1] = v1[0] * v2[1];
89 11490856 : deriv[startingIndex + 2] = v1[0] * v2[2];
90 11490856 : deriv[startingIndex + 3] = v1[1] * v2[0];
91 11490856 : deriv[startingIndex + 4] = v1[1] * v2[1];
92 11490856 : deriv[startingIndex + 5] = v1[1] * v2[2];
93 11490856 : deriv[startingIndex + 6] = v1[2] * v2[0];
94 11490856 : deriv[startingIndex + 7] = v1[2] * v2[1];
95 11490856 : deriv[startingIndex + 8] = v1[2] * v2[2];
96 11490856 : }
97 : };
98 :
99 : template <class CV, typename myPTM=defaultPTM>
100 : class AdjacencyMatrixBase : public ActionWithMatrix {
101 : public:
102 : using input_type = AdjacencyMatrixData<CV>;
103 : using PTM = typename myPTM::template PTM<AdjacencyMatrixBase<CV, myPTM>>;
104 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
105 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
106 : private:
107 : PTM taskmanager;
108 : bool nopbc, read_one_group;
109 : unsigned natoms_per_list;
110 : std::vector<std::size_t> nlist_v;
111 : LinkCells linkcells, threecells;
112 : /// True if the matrix elements are exactly zero beyond the link cell cutoff: the
113 : /// neighbour lists then hold only the atoms within the cutoff
114 : bool exact_cutoff{false};
115 : std::vector<unsigned> ablocks, threeblocks;
116 : double nl_cut, nl_cut2;
117 : unsigned nl_stride;
118 : void setupThirdAtomBlock( const std::vector<AtomNumber>& tc,
119 : std::vector<AtomNumber>& t );
120 : public:
121 : static constexpr size_t virialSize = 9;
122 : static void registerKeywords( Keywords& keys );
123 : explicit AdjacencyMatrixBase(const ActionOptions&);
124 : unsigned getNumberOfDerivatives() override ;
125 : void calculate() override ;
126 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
127 : void getAtomicInputData( std::vector<double>& inputdata ) const override;
128 : void getAtomicInputData( std::vector<float>& inputdata ) const override;
129 10 : std::string writeInGraph() const override {
130 : if constexpr ( std::is_same_v<ContactMatrix,CV> ) {
131 : //TODO:this will change to std::is_same_v<ContactMatrix<precision>,CV> when adding the mixed precision to the contact matrix
132 10 : return "CONTACT_MATRIX";
133 : }
134 0 : return getName();
135 : }
136 : /// Set the cutoff of the link cells. With exact=true the matrix elements are
137 : /// known to be zero when the atoms are further apart than lcut, so the rows of
138 : /// the sparse matrix hold only the atoms within the cutoff instead of all the
139 : /// atoms in the neighbouring link cells (many more at short cutoffs)
140 : void setLinkCellCutoff( const bool& symmetric, const double& lcut, double tcut=-1.0, const bool exact=false );
141 : static void performTask( std::size_t task_index,
142 : const AdjacencyMatrixData<CV>& actiondata,
143 : ParallelActionsInput& input,
144 : ParallelActionsOutput output );
145 : static int getNumberOfValuesPerTask( std::size_t task_index,
146 : const AdjacencyMatrixData<CV>& actiondata );
147 : static void getForceIndices( std::size_t task_index,
148 : std::size_t colno,
149 : std::size_t ntotal_force,
150 : const AdjacencyMatrixData<CV>& actiondata,
151 : const ParallelActionsInput& input,
152 : ForceIndexHolder force_indices );
153 : void getMatrixColumnTitles( std::vector<std::string>& argnames ) const override ;
154 : };
155 :
156 : template <class CV, typename myPTM>
157 1130 : void AdjacencyMatrixBase<CV, myPTM>::registerKeywords( Keywords& keys ) {
158 1130 : ActionWithMatrix::registerKeywords( keys );
159 2260 : keys.addInputKeyword("optional","MASK","vector","a vector that is used to used to determine which rows of the adjancency matrix to compute");
160 1130 : keys.add("atoms","GROUP","the atoms for which you would like to calculate the adjacency matrix");
161 1130 : keys.add("atoms","GROUPA","when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPB");
162 1130 : keys.add("atoms","GROUPB","when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPA");
163 1130 : keys.addDeprecatedKeyword("ATOMS","GROUP");
164 1130 : keys.reserve("atoms","GROUPC","a group of atoms that must be summed over when calculating each element of the adjacency matrix");
165 1130 : keys.addFlag("COMPONENTS",false,"also calculate the components of the vector connecting the atoms in the contact matrix");
166 1130 : keys.addFlag("NOPBC",false,"don't use pbc");
167 1130 : keys.add("compulsory","NL_CUTOFF","0.0","The cutoff for the neighbor list. A value of 0 means we are not using a neighbor list");
168 1130 : keys.add("compulsory","NL_STRIDE","1","The frequency with which we are updating the atoms in the neighbor list");
169 1130 : CV::registerKeywords( keys );
170 1130 : PTM::registerKeywords( keys );
171 2260 : keys.addOutputComponent("w","COMPONENTS","matrix","a matrix containing the weights for the bonds between each pair of atoms");
172 2260 : keys.addOutputComponent("x","COMPONENTS","matrix","the projection of the bond on the x axis");
173 2260 : keys.addOutputComponent("y","COMPONENTS","matrix","the projection of the bond on the y axis");
174 2260 : keys.addOutputComponent("z","COMPONENTS","matrix","the projection of the bond on the z axis");
175 2260 : keys.setValueDescription("matrix","a matrix containing the weights for the bonds between each pair of atoms");
176 1130 : keys.addDOI("10.1021/acs.jctc.6b01073");
177 1130 : }
178 :
179 : template <class CV, typename myPTM>
180 340 : AdjacencyMatrixBase<CV, myPTM>::AdjacencyMatrixBase(const ActionOptions& ao):
181 : Action(ao),
182 : ActionWithMatrix(ao),
183 340 : taskmanager(this),
184 340 : read_one_group(false),
185 340 : natoms_per_list(0),
186 340 : linkcells(comm),
187 680 : threecells(comm) {
188 340 : std::vector<std::size_t> shape(2);
189 340 : std::vector<AtomNumber> t;
190 680 : parseAtomList("GROUP", t );
191 340 : if( getName()!="HBOND_MATRIX" ) {
192 338 : if( t.size()==0 ) {
193 202 : parseAtomList("ATOMS", t);
194 101 : if( t.size()>0 ) {
195 0 : warning("using deprecated syntax for contact matrix. You are strongly recommended to use GROUP instead of ATOMS");
196 : }
197 : }
198 2 : } else if( t.size()>0 ) {
199 2 : warning("GROUP keyword has been deprecated for HBOND_MATRIX as it may lead users to wrongly assume that the matrices calculated by this action are symmetric. We strongly recommend using DONORS/ACCEPTORS instead");
200 : }
201 :
202 340 : if( t.size()==0 ) {
203 102 : std::vector<AtomNumber> ta;
204 102 : if( getName()=="HBOND_MATRIX") {
205 2 : parseAtomList("DONORS",ta);
206 : } else {
207 202 : parseAtomList("GROUPA",ta);
208 : }
209 102 : std::vector<AtomNumber> tb;
210 102 : if( getName()=="HBOND_MATRIX") {
211 2 : parseAtomList("ACCEPTORS",tb);
212 : } else {
213 202 : parseAtomList("GROUPB",tb);
214 : }
215 102 : if( ta.size()==0 || tb.size()==0 ) {
216 0 : error("no atoms have been specified in input");
217 : }
218 :
219 : // Create list of tasks
220 102 : log.printf(" atoms in GROUPA ");
221 5906 : for(unsigned i=0; i<ta.size(); ++i) {
222 5804 : log.printf("%d ", ta[i].serial());
223 5804 : t.push_back(ta[i]);
224 : }
225 102 : log.printf("\n");
226 102 : log.printf(" atoms in GROUPB ");
227 102 : ablocks.resize( tb.size() );
228 : unsigned n=0;
229 15014 : for(unsigned i=0; i<tb.size(); ++i) {
230 14912 : log.printf("%d ", tb[i].serial());
231 14912 : ablocks[i]=ta.size()+n;
232 14912 : t.push_back(tb[i]);
233 14912 : n++;
234 : }
235 102 : log.printf("\n");
236 102 : shape[0]=ta.size();
237 102 : shape[1]=tb.size();
238 102 : } else {
239 : // Create list of tasks
240 238 : log.printf(" atoms in GROUP ");
241 238 : ablocks.resize( t.size() );
242 145777 : for(unsigned i=0; i<t.size(); ++i) {
243 145539 : log.printf("%d ", t[i].serial());
244 145539 : ablocks[i]=i;
245 : }
246 238 : log.printf("\n");
247 238 : shape[0]=shape[1]=t.size();
248 238 : read_one_group=true;
249 : }
250 340 : if( keywords.exists("GROUPC") ) {
251 6 : std::vector<AtomNumber> tc;
252 12 : parseAtomList("GROUPC",tc);
253 6 : if( tc.size()==0 ) {
254 0 : error("no atoms in GROUPC specified");
255 : }
256 6 : log.printf(" atoms in GROUPC ");
257 6 : setupThirdAtomBlock( tc, t );
258 340 : } else if( keywords.exists("BRIDGING_ATOMS") ) {
259 0 : std::vector<AtomNumber> tc;
260 0 : parseAtomList("BRIDGING_ATOMS",tc);
261 0 : if( tc.size()==0 ) {
262 0 : error("no BRIDGING_ATOMS specified");
263 : }
264 0 : log.printf(" bridging atoms are ");
265 0 : setupThirdAtomBlock( tc, t );
266 334 : } else if( keywords.exists("HYDROGENS") ) {
267 2 : std::vector<AtomNumber> tc;
268 4 : parseAtomList("HYDROGENS",tc);
269 2 : if( tc.size()==0 ) {
270 0 : error("no HYDROGEN atoms specified");
271 : }
272 2 : log.printf(" hydrogen atoms are ");
273 2 : setupThirdAtomBlock( tc, t );
274 334 : } else if( keywords.exists("BACKGROUND_ATOMS") ) {
275 8 : std::vector<AtomNumber> tc;
276 16 : parseAtomList("BACKGROUND_ATOMS",tc);
277 8 : if( tc.size()==0 ) {
278 0 : error("no ATOMS atoms specified");
279 : }
280 8 : log.printf(" atoms for background density are ");
281 8 : setupThirdAtomBlock( tc, t );
282 8 : }
283 : // Request the atoms from the ActionAtomistic
284 340 : requestAtoms( t, false );
285 : bool components;
286 340 : parseFlag("COMPONENTS",components);
287 340 : parseFlag("NOPBC",nopbc);
288 340 : if( !components ) {
289 250 : addValue( shape );
290 250 : setNotPeriodic();
291 : } else {
292 90 : addComponent( "w", shape );
293 180 : componentIsNotPeriodic("w");
294 : }
295 340 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
296 : // Stuff for neighbor list
297 340 : parse("NL_CUTOFF",nl_cut);
298 340 : nl_cut2=nl_cut*nl_cut;
299 340 : parse("NL_STRIDE",nl_stride);
300 : // calculate() rebuilds the neighbor list at every step and cannot reuse it:
301 : // reject NL_STRIDE>1 here, so that the input fails at setup and not at the first step
302 340 : if( nl_stride!=1 ) {
303 0 : error("NL_STRIDE different from 1 is not implemented: the neighbor list is rebuilt at every step with link cells (cutoff from D_MAX, or NL_CUTOFF if set). Remove NL_STRIDE from the input");
304 : }
305 340 : if( nl_cut>0 ) {
306 0 : log.printf(" using neighbor list with cutoff %f. List is updated every %u steps.\n",nl_cut,nl_stride);
307 : }
308 :
309 340 : if( components ) {
310 90 : addComponent( "x", shape );
311 90 : componentIsNotPeriodic("x");
312 90 : addComponent( "y", shape );
313 90 : componentIsNotPeriodic("y");
314 90 : addComponent( "z", shape );
315 180 : componentIsNotPeriodic("z");
316 : }
317 386 : log<<" Bibliography "<<plumed.cite("Tribello, Giberti, Sosso, Salvalaglio and Parrinello, J. Chem. Theory Comput. 13, 1317 (2017)")<<"\n";
318 294 : AdjacencyMatrixData<CV> matdata;
319 340 : matdata.usepbc = !nopbc;
320 340 : matdata.components = components;
321 340 : matdata.nlists = getPntrToComponent(0)->getShape()[0];
322 340 : matdata.fbsize = 0;
323 340 : if( ablocks[0]>=getConstPntrToComponent(0)->getShape()[0] ) {
324 102 : matdata.fbsize = getConstPntrToComponent(0)->getShape()[0];
325 : }
326 340 : matdata.matrixdata.parseInput( this );
327 340 : matdata.gatherForceOnColumns = false;
328 340 : matdata.no_thread_gather = no_thread_gather;
329 340 : if( no_thread_gather && threeblocks.size()>0 ) {
330 0 : error("AVOID_THREAD_GATHER flag cannot be used with this action");
331 : }
332 340 : taskmanager.setActionInput( matdata );
333 340 : }
334 :
335 : template <class CV, typename myPTM>
336 22247 : unsigned AdjacencyMatrixBase<CV, myPTM>::getNumberOfDerivatives() {
337 22247 : return 3*getNumberOfAtoms() + 9;
338 : }
339 :
340 : template <class CV, typename myPTM>
341 16 : void AdjacencyMatrixBase<CV, myPTM>::setupThirdAtomBlock( const std::vector<AtomNumber>& tc,
342 : std::vector<AtomNumber>& t ) {
343 16 : threeblocks.resize( tc.size() );
344 16 : unsigned base=t.size();
345 49986 : for(unsigned i=0; i<tc.size(); ++i) {
346 49970 : log.printf("%d ", tc[i].serial());
347 49970 : t.push_back(tc[i]);
348 49970 : threeblocks[i]=base+i;
349 : }
350 16 : log.printf("\n");
351 16 : }
352 :
353 : template <class CV, typename myPTM>
354 340 : void AdjacencyMatrixBase<CV, myPTM>::setLinkCellCutoff( const bool& symmetric,
355 : const double& lcut,
356 : double tcut,
357 : const bool exact ) {
358 340 : if( read_one_group && symmetric ) {
359 233 : getPntrToComponent(0)->setSymmetric( true );
360 : }
361 340 : exact_cutoff = exact;
362 340 : if( nl_cut>0 && lcut>nl_cut ) {
363 0 : error("D_MAX for switching functions should be shorter than neighbor list cutoff");
364 : }
365 :
366 340 : if( tcut<0 ) {
367 332 : tcut=lcut;
368 : }
369 340 : if( nl_cut>0 ) {
370 0 : linkcells.setCutoff( nl_cut );
371 : } else {
372 340 : linkcells.setCutoff( lcut );
373 : }
374 340 : if( linkcells.getCutoff()<std::numeric_limits<double>::max() ) {
375 164 : log.printf(" set link cell cutoff to %f \n", linkcells.getCutoff() );
376 164 : if( exact_cutoff ) {
377 124 : log.printf(" only atoms within the cutoff are stored in the rows of the matrix \n");
378 : }
379 : }
380 340 : threecells.setCutoff( tcut );
381 340 : }
382 :
383 : template <class CV, typename myPTM>
384 11828 : void AdjacencyMatrixBase<CV, myPTM>::getAtomicInputData( std::vector<double>& inputdata ) const {
385 11828 : if( inputdata.size()!=3*getNumberOfAtoms() ) {
386 315 : inputdata.resize( 3*getNumberOfAtoms() );
387 : }
388 :
389 : std::size_t k=0;
390 1391660 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
391 1379832 : Vector mypos( ActionAtomistic::getPosition(i) );
392 1379832 : inputdata[k] = mypos[0];
393 1379832 : k++;
394 1379832 : inputdata[k] = mypos[1];
395 1379832 : k++;
396 1379832 : inputdata[k] = mypos[2];
397 1379832 : k++;
398 : }
399 11828 : }
400 :
401 : template <class CV, typename myPTM>
402 0 : void AdjacencyMatrixBase<CV, myPTM>::getAtomicInputData( std::vector<float>& inputdata ) const {
403 0 : if( inputdata.size()!=3*getNumberOfAtoms() ) {
404 0 : inputdata.resize( 3*getNumberOfAtoms() );
405 : }
406 :
407 : std::size_t k=0;
408 0 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
409 0 : Vector mypos( ActionAtomistic::getPosition(i) );
410 0 : inputdata[k] = mypos[0];
411 0 : k++;
412 0 : inputdata[k] = mypos[1];
413 0 : k++;
414 0 : inputdata[k] = mypos[2];
415 0 : k++;
416 : }
417 0 : }
418 :
419 : template <class CV, typename myPTM>
420 0 : void AdjacencyMatrixBase<CV, myPTM>::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
421 : std::string num;
422 0 : for(unsigned i=0; i<getConstPntrToComponent(0)->getShape()[1]; ++i) {
423 0 : Tools::convert( i+1, num );
424 0 : argnames.push_back( num );
425 : }
426 0 : }
427 :
428 : template <class CV, typename myPTM>
429 11828 : void AdjacencyMatrixBase<CV, myPTM>::calculate() {
430 11828 : Value* myval = getPntrToComponent(0);
431 : // Retrieve the task list
432 11828 : std::vector<unsigned> & pTaskList( getListOfActiveTasks() );
433 : // Get the number of tasks we have to deal with
434 : unsigned ntasks=myval->getShape()[0];
435 11828 : if( nl_stride==1 ) {
436 11828 : ntasks=pTaskList.size();
437 : } else {
438 0 : error("neighbour list non updates are not actually implemented or tested");
439 : }
440 :
441 : // Get the atoms
442 11828 : std::vector<Vector> ltmp_pos2( ntasks );
443 776352 : for(unsigned i=0; i<ntasks; ++i) {
444 764524 : ltmp_pos2[i] = ActionAtomistic::getPosition(pTaskList[i]);
445 : }
446 11828 : auto & matdata = taskmanager.getActionInput();
447 : // Now update the neighbor lists
448 11828 : if( getStep()%nl_stride==0 ) {
449 : // Build the link cells
450 11828 : std::vector<Vector> ltmp_pos( ablocks.size() );
451 1320172 : for(unsigned i=0; i<ablocks.size(); ++i) {
452 1308344 : ltmp_pos[i]=ActionAtomistic::getPosition( ablocks[i] );
453 : }
454 11828 : linkcells.createNeighborList( getConstPntrToComponent(0)->getShape()[0],
455 : make_const_view(ltmp_pos2),
456 : make_const_view(pTaskList),
457 : make_const_view(pTaskList),
458 : make_const_view(ltmp_pos),
459 : make_const_view(ablocks),
460 : getPbc(),
461 11828 : natoms_per_list,
462 11828 : nlist_v,
463 11828 : exact_cutoff );
464 11828 : if( threeblocks.size()>0 ) {
465 53 : std::vector<Vector> ltmp_pos3( threeblocks.size() );
466 50282 : for(unsigned i=0; i<threeblocks.size(); ++i) {
467 50229 : ltmp_pos3[i]=ActionAtomistic::getPosition( threeblocks[i] );
468 : }
469 53 : threecells.createNeighborList( getConstPntrToComponent(0)->getShape()[0],
470 : make_const_view(ltmp_pos2),
471 : make_const_view(pTaskList),
472 : make_const_view(pTaskList),
473 : make_const_view(ltmp_pos3),
474 : make_const_view(threeblocks),
475 : getPbc(),
476 53 : matdata.natoms_per_three_list,
477 53 : matdata.nlist_three_v);
478 : }
479 : matdata.update();
480 : } else {
481 0 : error("neighbour list non updates are not actually implemented or tested");
482 : }
483 : // And finally work out maximum number of columns to use
484 11828 : unsigned maxcol = nlist_v[0];
485 1281495 : for(unsigned i=1; i<getConstPntrToComponent(0)->getShape()[0]; ++i) {
486 1269667 : if( nlist_v[i]>maxcol ) {
487 1241 : maxcol = nlist_v[i];
488 : }
489 : }
490 : // The first element returned by the neighbour list is the central atom. The
491 : // neighbour list thus always has one more atom in it than there are atoms in the
492 : // neighbourhood of the central atom
493 11828 : maxcol = maxcol-1;
494 : // A matrix with no elements at all (no atom within the cutoff of any other)
495 : // still needs one column of storage: zero columns is the flag of the
496 : // lower-triangular storage in Value and leaves the data array empty
497 11828 : if( maxcol<1 ) {
498 1 : maxcol = 1;
499 : }
500 :
501 : // Reshape the matrix store if the number of columns has changed
502 11828 : if( maxcol!=myval->getNumberOfColumns() ) {
503 1346 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
504 889 : getPntrToComponent(i)->reshapeMatrixStore( maxcol );
505 : }
506 : }
507 : // Transfer neighbor list data to the bookeeping arrays
508 11828 : unsigned lstart = getConstPntrToComponent(0)->getShape()[0];
509 11828 : if( ntasks<lstart ) {
510 519941 : for(unsigned i=0; i<lstart; ++i) {
511 519814 : myval->setMatrixBookeepingElement( i*(1+maxcol), 0 );
512 : }
513 : }
514 776352 : for(unsigned i=0; i<ntasks; ++i) {
515 764524 : unsigned bstart = pTaskList[i]*(1+maxcol);
516 764524 : unsigned rstart = lstart + pTaskList[i]*(1+natoms_per_list);
517 764524 : myval->setMatrixBookeepingElement( bstart, nlist_v[pTaskList[i]] - 1 );
518 20722551 : for(unsigned j=1; j<nlist_v[pTaskList[i]]; ++j) {
519 19958027 : myval->setMatrixBookeepingElement( bstart+j, nlist_v[rstart+j] - matdata.fbsize );
520 : }
521 : }
522 11828 : copyMatrixBookeepingFromFirstComponent( taskmanager.getActionInput().outmat );
523 12620 : for(unsigned i=1; i<getNumberOfComponents(); ++i) {
524 792 : getPntrToComponent(i)->copyBookeepingArrayFromArgument( getPntrToComponent(0) );
525 : }
526 : // We need to setup the task manager here because at this point we know the size of the sparse matrices
527 11828 : if( no_thread_gather ) {
528 25 : unsigned maxatoms = maxcol;
529 25 : if( maxcolsize>maxatoms ) {
530 : maxatoms = maxcolsize;
531 : }
532 25 : taskmanager.setupParallelTaskManager( 3 + 3*matdata.natoms_per_three_list + virialSize, 3*threeblocks.size()+virialSize );
533 25 : taskmanager.setWorkspaceSize(3*(maxatoms + 1 + matdata.natoms_per_three_list) + 3*(6+matdata.natoms_per_three_list)+virialSize );
534 : } else {
535 11803 : taskmanager.setupParallelTaskManager( 6 + 3*matdata.natoms_per_three_list + virialSize, getNumberOfDerivatives() );
536 11803 : taskmanager.setWorkspaceSize(3*(maxcol + 1 + matdata.natoms_per_three_list) );
537 : }
538 : // And run all the tasks
539 11828 : taskmanager.runAllTasks();
540 11828 : }
541 :
542 : template <class CV, typename myPTM>
543 1074691 : void AdjacencyMatrixBase<CV, myPTM>::performTask( std::size_t task_index,
544 : const AdjacencyMatrixData<CV>& actiondata,
545 : ParallelActionsInput& input,
546 : ParallelActionsOutput output ) {
547 1074691 : const unsigned fstart=task_index*(1+actiondata.outmat.ncols);
548 1074691 : const unsigned nneigh=actiondata.outmat[fstart];
549 1074691 : const unsigned n3neigh=(actiondata.natoms_per_three_list>0) ?
550 1011 : actiondata.nlist_three[task_index]:0;
551 1074691 : VectorView atoms(output.buffer.data(), nneigh + n3neigh );
552 :
553 1074691 : if( actiondata.gatherForceOnColumns ) {
554 479 : const View<double,3> pos0( input.inputdata + 3 * (task_index+actiondata.fbsize) );
555 5659 : for(unsigned i=0; i<nneigh; ++i) {
556 5180 : atoms[i][0] = pos0[0] - input.inputdata[3*actiondata.outmat[fstart+1+i]+0];
557 5180 : atoms[i][1] = pos0[1] - input.inputdata[3*actiondata.outmat[fstart+1+i]+1];
558 5180 : atoms[i][2] = pos0[2] - input.inputdata[3*actiondata.outmat[fstart+1+i]+2];
559 : }
560 : } else {
561 1074212 : const View<double,3> pos0( input.inputdata + 3 * task_index);
562 22733349 : for(unsigned i=0; i<nneigh; ++i) {
563 21659137 : atoms[i][0] = input.inputdata[3*(actiondata.outmat[fstart+1+i]+actiondata.fbsize)+0] - pos0[0];
564 21659137 : atoms[i][1] = input.inputdata[3*(actiondata.outmat[fstart+1+i]+actiondata.fbsize)+1] - pos0[1];
565 21659137 : atoms[i][2] = input.inputdata[3*(actiondata.outmat[fstart+1+i]+actiondata.fbsize)+2] - pos0[2];
566 : }
567 :
568 1074212 : if( actiondata.natoms_per_three_list>0 ) {
569 : // Retrieve the set of third atoms
570 1011 : unsigned fstart3 = actiondata.nlists
571 1011 : + task_index*(1+actiondata.natoms_per_three_list);
572 4165164 : for(unsigned i=1; i<n3neigh; ++i) {
573 4164153 : atoms[nneigh+i-1][0] =
574 4164153 : input.inputdata[3*actiondata.nlist_three[fstart3+i]+0] - pos0[0];
575 4164153 : atoms[nneigh+i-1][1] =
576 4164153 : input.inputdata[3*actiondata.nlist_three[fstart3+i]+1] - pos0[1];
577 4164153 : atoms[nneigh+i-1][2] =
578 4164153 : input.inputdata[3*actiondata.nlist_three[fstart3+i]+2] - pos0[2];
579 : }
580 : }
581 : }
582 :
583 :
584 : // Apply periodic boundary conditions to all the atoms
585 1074691 : if( actiondata.usepbc ) {
586 1074691 : input.pbc->apply( atoms, atoms.size() );
587 : }
588 1074691 : AdjacencyMatrixInput adjinp {input.noderiv,
589 1074691 : input.pbc,
590 : Vector{0.0,0.0,0.0},
591 : 0,
592 1074691 : VectorView{ atoms.data() + 3*nneigh,n3neigh}};
593 :
594 1074691 : if( n3neigh>1 ) {
595 1011 : adjinp.natoms = n3neigh-1;
596 : }
597 :
598 : // And calculate this row of the matrices
599 1074691 : std::size_t nderiv(6 + 3*adjinp.natoms + virialSize);
600 1074691 : if( actiondata.no_thread_gather ) {
601 887 : nderiv = 3 + + 3*adjinp.natoms + virialSize;
602 : }
603 1074691 : const unsigned ncomponents =(actiondata.components) ? 4 : 1;
604 :
605 : // Must clear the derivatives here as otherwise sparsity pattern
606 : // of derivatives does not match sparsity pattern for forces
607 1074691 : if( !input.noderiv ) {
608 209733748 : for(unsigned i=0; i<output.derivatives.size(); ++i) {
609 209389662 : output.derivatives[i] = 0.0;
610 : }
611 : }
612 :
613 1074691 : if( actiondata.no_thread_gather ) {
614 16427 : for(unsigned i=0; i<nneigh; ++i ) {
615 15540 : adjinp.pos = Vector(atoms[i][0],atoms[i][1],atoms[i][2]);
616 15540 : const std::size_t valpos = static_cast<std::size_t>(i) * ncomponents;
617 15540 : MatrixOutput adjout{View<double,1>{&output.values[valpos]},
618 15540 : View{output.buffer.data() + 3*(nneigh+n3neigh),nderiv+3}};
619 15540 : CV::calculateWeight( actiondata.matrixdata, adjinp, adjout );
620 15540 : if( actiondata.components ) {
621 12096 : output.values[valpos+1] = atoms[i][0];
622 12096 : output.values[valpos+2] = atoms[i][1];
623 12096 : output.values[valpos+3] = atoms[i][2];
624 : }
625 15540 : if( input.noderiv ) {
626 7476 : continue;
627 : }
628 10360 : if( fabs(adjout.val[0])>epsilon ) {
629 10318 : if( actiondata.gatherForceOnColumns ) {
630 5159 : auto derivs = output.derivatives.template subview_n<3>(valpos*nderiv);
631 5159 : derivs[0] = adjout.deriv[3];
632 5159 : derivs[1] = adjout.deriv[4];
633 5159 : derivs[2] = adjout.deriv[5];
634 : } else {
635 5159 : auto derivs = output.derivatives.template subview_n<3>(valpos*nderiv);
636 5159 : derivs[0] = adjout.deriv[0];
637 5159 : derivs[1] = adjout.deriv[1];
638 5159 : derivs[2] = adjout.deriv[2];
639 51590 : for(unsigned j=6; j<nderiv+3; ++j) {
640 46431 : output.derivatives[valpos*nderiv + j - 3] = adjout.deriv[j];
641 : }
642 : }
643 : }
644 10360 : if( !actiondata.components ) {
645 2296 : continue;
646 : }
647 8064 : if( actiondata.gatherForceOnColumns ) {
648 16128 : for(int ii=1; ii<4; ++ii) {
649 12096 : auto derivs = output.derivatives.template subview_n<1>(valpos*nderiv+ii*nderiv);
650 12096 : derivs[0] = 1.0;
651 : }
652 : } else {
653 16128 : for(int ii=1; ii<4; ++ii) {
654 12096 : auto derivs = output.derivatives.template subview_n<4>(valpos*nderiv+ii*nderiv);
655 12096 : derivs[0] = -1.0;
656 12096 : derivs[1] = -atoms[i][0];
657 12096 : derivs[2] = -atoms[i][1];
658 12096 : derivs[3] = -atoms[i][2];
659 : }
660 : }
661 : }
662 : } else {
663 22722581 : for(unsigned i=0; i<nneigh; ++i ) {
664 21648777 : adjinp.pos = Vector(atoms[i][0],atoms[i][1],atoms[i][2]);
665 21648777 : const std::size_t valpos = static_cast<std::size_t>(i) * ncomponents;
666 21648777 : MatrixOutput adjout{View<double,1>{&output.values[valpos]},
667 21648777 : View{output.derivatives.data() + valpos*nderiv,nderiv}};
668 21648777 : CV::calculateWeight( actiondata.matrixdata, adjinp, adjout );
669 21648777 : if( !actiondata.components ) {
670 19738267 : continue ;
671 : }
672 6023474 : output.values[valpos+1] = atoms[i][0];
673 6023474 : output.values[valpos+2] = atoms[i][1];
674 6023474 : output.values[valpos+3] = atoms[i][2];
675 6023474 : if( input.noderiv ) {
676 4112964 : continue ;
677 : }
678 : //sugar for not having to repeat [valpos*nderiv+something]
679 7642040 : for(int ii=1; ii<4; ++ii) {
680 5731530 : auto derivs = output.derivatives.template subview_n<5>(valpos*nderiv+ii*nderiv);
681 5731530 : derivs[0] = -1.0;
682 5731530 : derivs[1] = 1.0;
683 5731530 : derivs[2] = -atoms[i][0];
684 5731530 : derivs[3] = -atoms[i][1];
685 5731530 : derivs[4] = -atoms[i][2];
686 : }
687 : }
688 : }
689 1074691 : }
690 :
691 : template <class CV, typename myPTM>
692 11085 : void AdjacencyMatrixBase<CV, myPTM>::applyNonZeroRankForces( std::vector<double>& outforces ) {
693 11085 : taskmanager.applyForces( outforces );
694 11085 : if( no_thread_gather ) {
695 25 : getColumnBookeepingArrays( taskmanager.getActionInput().outmat );
696 25 : taskmanager.getActionInput().gatherForceOnColumns = gatherForceOnColumns = true;
697 25 : taskmanager.setNForceScalars( getNumberOfComponents()*maxcolsize );
698 25 : taskmanager.applyForces( outforces, false );
699 25 : taskmanager.setNForceScalars( getNumberOfComponents()*getPntrToComponent(0)->getNumberOfColumns() );
700 25 : taskmanager.getActionInput().gatherForceOnColumns = gatherForceOnColumns = false;
701 : }
702 11085 : }
703 :
704 : template <class CV, typename myPTM>
705 344086 : int AdjacencyMatrixBase<CV, myPTM>::getNumberOfValuesPerTask( std::size_t task_index,
706 : const AdjacencyMatrixData<CV>& actiondata ) {
707 344086 : return actiondata.outmat[task_index*(1+actiondata.outmat.ncols)];
708 : }
709 :
710 : template <class CV, typename myPTM>
711 5524898 : void AdjacencyMatrixBase<CV, myPTM>::getForceIndices( const std::size_t task_index,
712 : const std::size_t colno,
713 : const std::size_t ntotal_force,
714 : const AdjacencyMatrixData<CV>& actiondata,
715 : const ParallelActionsInput& /*input*/,
716 : ForceIndexHolder force_indices ) {
717 :
718 : // Map a (task, column) pair onto the contiguous force-index ranges used by
719 : // the parallel force accumulation code.
720 : //
721 : // The mapping depends on how partial thread results are gathered:
722 : // 1) standard gather (!no_thread_gather): use the 3-list decomposition;
723 : // 2) no-thread-gather + gather-on-columns: use per-column contiguous blocks;
724 : // 3) no-thread-gather + default layout: use the regular component-major layout.
725 : //
726 : // `ntotal_force` is the global scalar-force extent for the current layout and
727 : // is used to offset into the correct component/column segment.
728 5524898 : if( !actiondata.no_thread_gather ) {
729 : // Standard threaded gather path: each logical task is expanded into three
730 : // sub-tasks (three-list bookkeeping), hence the 3*task_index conversion.
731 5514580 : const auto three_task_index= 3 * task_index;
732 5514580 : force_indices.indices[0][0] = three_task_index;
733 5514580 : force_indices.indices[0][1] = three_task_index + 1;
734 5514580 : force_indices.indices[0][2] = three_task_index + 2;
735 :
736 5514580 : const unsigned myatom = 3 * (actiondata.outmat[ task_index*(1+actiondata.outmat.ncols) + 1 + colno] + actiondata.fbsize);
737 5514580 : force_indices.indices[0][3] = myatom;
738 5514580 : force_indices.indices[0][4] = myatom + 1;
739 5514580 : force_indices.indices[0][5] = myatom + 2;
740 5514580 : if( actiondata.components ) {
741 1910510 : force_indices.indices[1][0] = three_task_index;
742 1910510 : force_indices.indices[1][1] = myatom;
743 1910510 : force_indices.indices[2][0] = three_task_index + 1;
744 1910510 : force_indices.indices[2][1] = myatom + 1;
745 1910510 : force_indices.indices[3][0] = three_task_index + 2;
746 1910510 : force_indices.indices[3][1] = myatom + 2;
747 : }
748 :
749 5514580 : force_indices.threadsafe_derivatives_end[0] = 0;
750 : unsigned n = 6;
751 5514580 : if( actiondata.natoms_per_three_list>0 ) {
752 4065 : const unsigned n3neigh = actiondata.nlist_three[task_index];
753 4065 : const unsigned fstart3 = actiondata.nlists
754 4065 : + task_index*(1 + actiondata.natoms_per_three_list);
755 500948 : for(unsigned j=1; j<n3neigh; ++j) {
756 496883 : unsigned my3atom = 3 * actiondata.nlist_three[fstart3 + j];
757 496883 : force_indices.indices[0][n ] = my3atom;
758 496883 : force_indices.indices[0][n+1] = my3atom + 1;
759 496883 : force_indices.indices[0][n+2] = my3atom + 2;
760 496883 : n += 3;
761 : }
762 : }
763 5514580 : const unsigned virstart = ntotal_force - 9;
764 5514580 : force_indices.indices[0][n ] = virstart + 0;
765 5514580 : force_indices.indices[0][n+1] = virstart + 1;
766 5514580 : force_indices.indices[0][n+2] = virstart + 2;
767 5514580 : force_indices.indices[0][n+3] = virstart + 3;
768 5514580 : force_indices.indices[0][n+4] = virstart + 4;
769 5514580 : force_indices.indices[0][n+5] = virstart + 5;
770 5514580 : force_indices.indices[0][n+6] = virstart + 6;
771 5514580 : force_indices.indices[0][n+7] = virstart + 7;
772 5514580 : force_indices.indices[0][n+8] = virstart + 8;
773 :
774 5514580 : force_indices.tot_indices[0] = n + 9;
775 5514580 : if( actiondata.components ) {
776 1910510 : force_indices.threadsafe_derivatives_end[1] = 0;
777 1910510 : force_indices.threadsafe_derivatives_end[2] = 0;
778 1910510 : force_indices.threadsafe_derivatives_end[3] = 0;
779 1910510 : force_indices.indices[1][2] = virstart;
780 1910510 : force_indices.indices[1][3] = virstart + 1;
781 1910510 : force_indices.indices[1][4] = virstart + 2;
782 1910510 : force_indices.indices[2][2] = virstart + 3;
783 1910510 : force_indices.indices[2][3] = virstart + 4;
784 1910510 : force_indices.indices[2][4] = virstart + 5;
785 1910510 : force_indices.indices[3][2] = virstart + 6;
786 1910510 : force_indices.indices[3][3] = virstart + 7;
787 1910510 : force_indices.indices[3][4] = virstart + 8;
788 1910510 : force_indices.tot_indices[1] = 5;
789 1910510 : force_indices.tot_indices[2] = 5;
790 1910510 : force_indices.tot_indices[3] = 5;
791 : }
792 10318 : } else if( actiondata.gatherForceOnColumns ) {
793 : // Alternative no-thread-gather layout: forces are grouped by matrix columns.
794 : // This branch computes indices into column-contiguous storage to reduce
795 : // gather/scatter overhead when operating on sparse column subsets.
796 5159 : const auto three_task_index = 3 * (task_index+actiondata.fbsize);
797 5159 : force_indices.indices[0][0] = three_task_index;
798 5159 : force_indices.indices[0][1] = three_task_index + 1;
799 5159 : force_indices.indices[0][2] = three_task_index + 2;
800 5159 : force_indices.threadsafe_derivatives_end[0] = 3;
801 5159 : force_indices.tot_indices[0] = 3;
802 5159 : if( actiondata.components ) {
803 4032 : force_indices.indices[1][0] = three_task_index;
804 4032 : force_indices.indices[2][0] = three_task_index + 1;
805 4032 : force_indices.indices[3][0] = three_task_index + 2;
806 4032 : force_indices.threadsafe_derivatives_end[1] = 1;
807 4032 : force_indices.threadsafe_derivatives_end[2] = 1;
808 4032 : force_indices.threadsafe_derivatives_end[3] = 1;
809 4032 : force_indices.tot_indices[1] = 1;
810 4032 : force_indices.tot_indices[2] = 1;
811 4032 : force_indices.tot_indices[3] = 1;
812 : }
813 : } else {
814 5159 : const auto three_task_index= 3 * task_index;
815 5159 : force_indices.indices[0][0] = three_task_index;
816 5159 : force_indices.indices[0][1] = three_task_index + 1;
817 5159 : force_indices.indices[0][2] = three_task_index + 2;
818 5159 : force_indices.threadsafe_derivatives_end[0] = 3;
819 5159 : const unsigned virstart = ntotal_force - 9;
820 5159 : force_indices.indices[0][3] = virstart + 0;
821 5159 : force_indices.indices[0][4] = virstart + 1;
822 5159 : force_indices.indices[0][5] = virstart + 2;
823 5159 : force_indices.indices[0][6] = virstart + 3;
824 5159 : force_indices.indices[0][7] = virstart + 4;
825 5159 : force_indices.indices[0][8] = virstart + 5;
826 5159 : force_indices.indices[0][9] = virstart + 6;
827 5159 : force_indices.indices[0][10] = virstart + 7;
828 5159 : force_indices.indices[0][11] = virstart + 8;
829 5159 : force_indices.tot_indices[0] = 12;
830 5159 : if( actiondata.components ) {
831 4032 : force_indices.indices[1][0] = three_task_index;
832 4032 : force_indices.indices[2][0] = three_task_index + 1;
833 4032 : force_indices.indices[3][0] = three_task_index + 2;
834 4032 : force_indices.threadsafe_derivatives_end[1] = 1;
835 4032 : force_indices.threadsafe_derivatives_end[2] = 1;
836 4032 : force_indices.threadsafe_derivatives_end[3] = 1;
837 4032 : force_indices.indices[1][1] = virstart;
838 4032 : force_indices.indices[1][2] = virstart + 1;
839 4032 : force_indices.indices[1][3] = virstart + 2;
840 4032 : force_indices.indices[2][1] = virstart + 3;
841 4032 : force_indices.indices[2][2] = virstart + 4;
842 4032 : force_indices.indices[2][3] = virstart + 5;
843 4032 : force_indices.indices[3][1] = virstart + 6;
844 4032 : force_indices.indices[3][2] = virstart + 7;
845 4032 : force_indices.indices[3][3] = virstart + 8;
846 4032 : force_indices.tot_indices[1] = 4;
847 4032 : force_indices.tot_indices[2] = 4;
848 4032 : force_indices.tot_indices[3] = 4;
849 : }
850 : }
851 5524898 : }
852 :
853 : } // namespace adjmat
854 : } // namespace PLMD
855 :
856 : #endif //__PLUMED_adjmat_AdjacencyMatrixBase_h
|
