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 : template <class T>
35 66 : struct AdjacencyMatrixData {
36 : T matrixdata;
37 : bool usepbc{true};
38 : bool components{false};
39 : std::size_t nlists{0};
40 : unsigned natoms_per_list{0};
41 : std::vector<std::size_t> nlist_v;
42 : std::size_t *nlist{nullptr};
43 : unsigned natoms_per_three_list{0};
44 : std::vector<std::size_t> nlist_three_v;
45 : std::size_t* nlist_three{nullptr};
46 : void update() {
47 11507 : nlist=nlist_v.data();
48 11507 : nlist_three=nlist_three_v.data();
49 : }
50 : #ifdef __PLUMED_USE_OPENACC
51 : void toACCDevice() const {
52 : #pragma acc enter data copyin(this[0:1],usepbc,components,nlists, \
53 : natoms_per_list,nlist[0:nlist_v.size()], \
54 : natoms_per_three_list, \
55 : nlist_three[0:nlist_three_v.size()])
56 : matrixdata.toACCDevice();
57 : }
58 : void removeFromACCDevice() const {
59 : matrixdata.removeFromACCDevice();
60 : #pragma acc exit data delete(nlist_three[0:nlist_three_v.size()], \
61 : natoms_per_three_list, \
62 : nlist[0:nlist_v.size()], natoms_per_list, \
63 : nlists, components, usepbc, this[0:1])
64 : }
65 : #endif //__PLUMED_USE_OPENACC
66 : };
67 :
68 : struct AdjacencyMatrixInput {
69 : bool noderiv{false};
70 : const Pbc* pbc;
71 : Vector pos;
72 : std::size_t natoms{0};
73 : VectorView extra_positions;
74 : };
75 :
76 : struct MatrixOutput {
77 : View<double,1> val;
78 : View<double> deriv;
79 :
80 : ///doing t= Tensor(v1,v2); deriv[x:x+9]=t with no extra memory allocation
81 : template <typename Iterable1, typename Iterable2>
82 10169006 : void assignOuterProduct(const std::size_t startingIndex,
83 : const Iterable1& v1,
84 : const Iterable2& v2 ) {
85 10169006 : deriv[startingIndex + 0] = v1[0] * v2[0];
86 10169006 : deriv[startingIndex + 1] = v1[0] * v2[1];
87 10169006 : deriv[startingIndex + 2] = v1[0] * v2[2];
88 10169006 : deriv[startingIndex + 3] = v1[1] * v2[0];
89 10169006 : deriv[startingIndex + 4] = v1[1] * v2[1];
90 10169006 : deriv[startingIndex + 5] = v1[1] * v2[2];
91 10169006 : deriv[startingIndex + 6] = v1[2] * v2[0];
92 10169006 : deriv[startingIndex + 7] = v1[2] * v2[1];
93 10169006 : deriv[startingIndex + 8] = v1[2] * v2[2];
94 10169006 : }
95 : };
96 :
97 : template <class T>
98 : class AdjacencyMatrixBase : public ActionWithMatrix {
99 : public:
100 : using input_type = AdjacencyMatrixData<T>;
101 : using PTM = ParallelTaskManager<AdjacencyMatrixBase<T>>;
102 : private:
103 : PTM taskmanager;
104 : bool nopbc, read_one_group;
105 : LinkCells linkcells, threecells;
106 : std::vector<unsigned> ablocks, threeblocks;
107 : double nl_cut, nl_cut2;
108 : unsigned nl_stride;
109 : void setupThirdAtomBlock( const std::vector<AtomNumber>& tc,
110 : std::vector<AtomNumber>& t );
111 : public:
112 : static constexpr size_t virialSize = 9;
113 : static void registerKeywords( Keywords& keys );
114 : explicit AdjacencyMatrixBase(const ActionOptions&);
115 : unsigned getNumberOfDerivatives() override ;
116 : void calculate() override ;
117 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
118 : void getInputData( std::vector<double>& inputdata ) const override;
119 10 : std::string writeInGraph() const override {
120 10 : if( getName()=="CONTACT_MATRIX_PROPER" ) {
121 10 : return "CONTACT_MATRIX";
122 : }
123 0 : return getName();
124 : }
125 : void setLinkCellCutoff( const bool& symmetric, const double& lcut, double tcut=-1.0 );
126 : static void performTask( std::size_t task_index,
127 : const AdjacencyMatrixData<T>& actiondata,
128 : ParallelActionsInput& input,
129 : ParallelActionsOutput output );
130 : static int getNumberOfValuesPerTask( std::size_t task_index,
131 : const AdjacencyMatrixData<T>& actiondata );
132 : static void getForceIndices( std::size_t task_index,
133 : std::size_t colno,
134 : std::size_t ntotal_force,
135 : const AdjacencyMatrixData<T>& actiondata,
136 : const ParallelActionsInput& input,
137 : ForceIndexHolder force_indices );
138 : void getMatrixColumnTitles( std::vector<std::string>& argnames ) const override ;
139 : };
140 :
141 : template <class T>
142 905 : void AdjacencyMatrixBase<T>::registerKeywords( Keywords& keys ) {
143 905 : ActionWithMatrix::registerKeywords( keys );
144 1810 : keys.addInputKeyword("optional","MASK","vector","a vector that is used to used to determine which rows of the adjancency matrix to compute");
145 905 : keys.add("atoms","GROUP","the atoms for which you would like to calculate the adjacency matrix");
146 905 : 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");
147 905 : 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");
148 905 : keys.addDeprecatedKeyword("ATOMS","GROUP");
149 905 : keys.reserve("atoms","GROUPC","a group of atoms that must be summed over when calculating each element of the adjacency matrix");
150 905 : keys.addFlag("COMPONENTS",false,"also calculate the components of the vector connecting the atoms in the contact matrix");
151 905 : keys.addFlag("NOPBC",false,"don't use pbc");
152 905 : 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");
153 905 : keys.add("compulsory","NL_STRIDE","1","The frequency with which we are updating the atoms in the neighbor list");
154 905 : T::registerKeywords( keys );
155 905 : PTM::registerKeywords( keys );
156 1810 : keys.addOutputComponent("w","COMPONENTS","matrix","a matrix containing the weights for the bonds between each pair of atoms");
157 1810 : keys.addOutputComponent("x","COMPONENTS","matrix","the projection of the bond on the x axis");
158 1810 : keys.addOutputComponent("y","COMPONENTS","matrix","the projection of the bond on the y axis");
159 1810 : keys.addOutputComponent("z","COMPONENTS","matrix","the projection of the bond on the z axis");
160 1810 : keys.setValueDescription("matrix","a matrix containing the weights for the bonds between each pair of atoms");
161 905 : keys.addDOI("10.1021/acs.jctc.6b01073");
162 905 : }
163 :
164 : template <class T>
165 276 : AdjacencyMatrixBase<T>::AdjacencyMatrixBase(const ActionOptions& ao):
166 : Action(ao),
167 : ActionWithMatrix(ao),
168 276 : taskmanager(this),
169 276 : read_one_group(false),
170 276 : linkcells(comm),
171 552 : threecells(comm) {
172 276 : std::vector<std::size_t> shape(2);
173 : std::vector<AtomNumber> t;
174 552 : parseAtomList("GROUP", t );
175 276 : if( getName()!="HBOND_MATRIX" ) {
176 274 : if( t.size()==0 ) {
177 148 : parseAtomList("ATOMS", t);
178 74 : if( t.size()>0 ) {
179 0 : warning("using depracated syntax for contact matrix. You are strongly recommended to use GROUP instead of ATOMS");
180 : }
181 : }
182 2 : } else if( t.size()>0 ) {
183 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");
184 : }
185 :
186 276 : if( t.size()==0 ) {
187 : std::vector<AtomNumber> ta;
188 75 : if( getName()=="HBOND_MATRIX") {
189 2 : parseAtomList("DONORS",ta);
190 : } else {
191 148 : parseAtomList("GROUPA",ta);
192 : }
193 : std::vector<AtomNumber> tb;
194 75 : if( getName()=="HBOND_MATRIX") {
195 2 : parseAtomList("ACCEPTORS",tb);
196 : } else {
197 148 : parseAtomList("GROUPB",tb);
198 : }
199 75 : if( ta.size()==0 || tb.size()==0 ) {
200 0 : error("no atoms have been specified in input");
201 : }
202 :
203 : // Create list of tasks
204 75 : log.printf(" atoms in GROUPA ");
205 2215 : for(unsigned i=0; i<ta.size(); ++i) {
206 2140 : log.printf("%d ", ta[i].serial());
207 2140 : t.push_back(ta[i]);
208 : }
209 75 : log.printf("\n");
210 75 : log.printf(" atoms in GROUPB ");
211 75 : ablocks.resize( tb.size() );
212 : unsigned n=0;
213 10335 : for(unsigned i=0; i<tb.size(); ++i) {
214 10260 : log.printf("%d ", tb[i].serial());
215 10260 : ablocks[i]=ta.size()+n;
216 10260 : t.push_back(tb[i]);
217 10260 : n++;
218 : }
219 75 : log.printf("\n");
220 75 : shape[0]=ta.size();
221 75 : shape[1]=tb.size();
222 : } else {
223 : // Create list of tasks
224 201 : log.printf(" atoms in GROUP ");
225 201 : ablocks.resize( t.size() );
226 143642 : for(unsigned i=0; i<t.size(); ++i) {
227 143441 : log.printf("%d ", t[i].serial());
228 143441 : ablocks[i]=i;
229 : }
230 201 : log.printf("\n");
231 201 : shape[0]=shape[1]=t.size();
232 201 : read_one_group=true;
233 : }
234 276 : if( keywords.exists("GROUPC") ) {
235 : std::vector<AtomNumber> tc;
236 12 : parseAtomList("GROUPC",tc);
237 6 : if( tc.size()==0 ) {
238 0 : error("no atoms in GROUPC specified");
239 : }
240 6 : log.printf(" atoms in GROUPC ");
241 6 : setupThirdAtomBlock( tc, t );
242 270 : } else if( keywords.exists("BRIDGING_ATOMS") ) {
243 : std::vector<AtomNumber> tc;
244 16 : parseAtomList("BRIDGING_ATOMS",tc);
245 8 : if( tc.size()==0 ) {
246 0 : error("no BRIDGING_ATOMS specified");
247 : }
248 8 : log.printf(" bridging atoms are ");
249 8 : setupThirdAtomBlock( tc, t );
250 262 : } else if( keywords.exists("HYDROGENS") ) {
251 : std::vector<AtomNumber> tc;
252 4 : parseAtomList("HYDROGENS",tc);
253 2 : if( tc.size()==0 ) {
254 0 : error("no HYDROGEN atoms specified");
255 : }
256 2 : log.printf(" hydrogen atoms are ");
257 2 : setupThirdAtomBlock( tc, t );
258 260 : } else if( keywords.exists("BACKGROUND_ATOMS") ) {
259 : std::vector<AtomNumber> tc;
260 16 : parseAtomList("BACKGROUND_ATOMS",tc);
261 8 : if( tc.size()==0 ) {
262 0 : error("no ATOMS atoms specified");
263 : }
264 8 : log.printf(" atoms for background density are ");
265 8 : setupThirdAtomBlock( tc, t );
266 : }
267 : // Request the atoms from the ActionAtomistic
268 276 : requestAtoms( t, false );
269 : bool components;
270 276 : parseFlag("COMPONENTS",components);
271 276 : parseFlag("NOPBC",nopbc);
272 276 : if( !components ) {
273 192 : addValue( shape );
274 192 : setNotPeriodic();
275 : } else {
276 84 : addComponent( "w", shape );
277 168 : componentIsNotPeriodic("w");
278 : }
279 276 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
280 : // Stuff for neighbor list
281 276 : parse("NL_CUTOFF",nl_cut);
282 276 : nl_cut2=nl_cut*nl_cut;
283 276 : parse("NL_STRIDE",nl_stride);
284 276 : if( nl_cut==0 && nl_stride>1 ) {
285 0 : error("NL_CUTOFF must be set if NL_STRIDE is set greater than 1");
286 : }
287 276 : if( nl_cut>0 ) {
288 0 : log.printf(" using neighbor list with cutoff %f. List is updated every %u steps.\n",nl_cut,nl_stride);
289 : }
290 :
291 276 : if( components ) {
292 84 : addComponent( "x", shape );
293 84 : componentIsNotPeriodic("x");
294 84 : addComponent( "y", shape );
295 84 : componentIsNotPeriodic("y");
296 84 : addComponent( "z", shape );
297 168 : componentIsNotPeriodic("z");
298 : }
299 552 : log<<" Bibliography "<<plumed.cite("Tribello, Giberti, Sosso, Salvalaglio and Parrinello, J. Chem. Theory Comput. 13, 1317 (2017)")<<"\n";
300 243 : AdjacencyMatrixData<T> matdata;
301 276 : matdata.usepbc = !nopbc;
302 276 : matdata.components = components;
303 276 : matdata.nlists = getPntrToComponent(0)->getShape()[0];
304 276 : matdata.matrixdata.parseInput( this );
305 276 : taskmanager.setActionInput( matdata );
306 519 : }
307 :
308 : template <class T>
309 21787 : unsigned AdjacencyMatrixBase<T>::getNumberOfDerivatives() {
310 21787 : return 3*getNumberOfAtoms() + 9;
311 : }
312 :
313 : template <class T>
314 24 : void AdjacencyMatrixBase<T>::setupThirdAtomBlock( const std::vector<AtomNumber>& tc,
315 : std::vector<AtomNumber>& t ) {
316 24 : threeblocks.resize( tc.size() );
317 24 : unsigned base=t.size();
318 53168 : for(unsigned i=0; i<tc.size(); ++i) {
319 53144 : log.printf("%d ", tc[i].serial());
320 53144 : t.push_back(tc[i]);
321 53144 : threeblocks[i]=base+i;
322 : }
323 24 : log.printf("\n");
324 24 : }
325 :
326 : template <class T>
327 276 : void AdjacencyMatrixBase<T>::setLinkCellCutoff( const bool& symmetric,
328 : const double& lcut,
329 : double tcut ) {
330 276 : if( read_one_group && symmetric ) {
331 196 : getPntrToComponent(0)->setSymmetric( true );
332 : }
333 276 : if( nl_cut>0 && lcut>nl_cut ) {
334 0 : error("D_MAX for switching functions should be shorter than neighbor list cutoff");
335 : }
336 :
337 276 : if( tcut<0 ) {
338 268 : tcut=lcut;
339 : }
340 276 : if( nl_cut>0 ) {
341 0 : linkcells.setCutoff( nl_cut );
342 : } else {
343 276 : linkcells.setCutoff( lcut );
344 : }
345 276 : if( linkcells.getCutoff()<std::numeric_limits<double>::max() ) {
346 135 : log.printf(" set link cell cutoff to %f \n", linkcells.getCutoff() );
347 : }
348 276 : threecells.setCutoff( tcut );
349 276 : }
350 :
351 : template <class T>
352 11507 : void AdjacencyMatrixBase<T>::getInputData( std::vector<double>& inputdata ) const {
353 11507 : if( inputdata.size()!=3*getNumberOfAtoms() ) {
354 258 : inputdata.resize( 3*getNumberOfAtoms() );
355 : }
356 :
357 : std::size_t k=0;
358 1359712 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
359 1348205 : Vector mypos( ActionAtomistic::getPosition(i) );
360 1348205 : inputdata[k] = mypos[0];
361 1348205 : k++;
362 1348205 : inputdata[k] = mypos[1];
363 1348205 : k++;
364 1348205 : inputdata[k] = mypos[2];
365 1348205 : k++;
366 : }
367 11507 : }
368 :
369 : template <class T>
370 0 : void AdjacencyMatrixBase<T>::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
371 : std::string num;
372 0 : for(unsigned i=0; i<getConstPntrToComponent(0)->getShape()[1]; ++i) {
373 0 : Tools::convert( i+1, num );
374 0 : argnames.push_back( num );
375 : }
376 0 : }
377 :
378 : template <class T>
379 11507 : void AdjacencyMatrixBase<T>::calculate() {
380 11507 : Value* myval = getPntrToComponent(0);
381 : // Retrieve the task list
382 11507 : std::vector<unsigned> & pTaskList( getListOfActiveTasks(this) );
383 : // Get the number of tasks we have to deal with
384 : unsigned ntasks=myval->getShape()[0];
385 11507 : if( nl_stride==1 ) {
386 11507 : ntasks=pTaskList.size();
387 : } else {
388 0 : error("neighbour list non updates are not actually implemented or tested");
389 : }
390 : unsigned fbsize=0;
391 11507 : unsigned lstart = getConstPntrToComponent(0)->getShape()[0];
392 : // This is triggered if you have GROUPA/GROUPB
393 : // in that case the second index for the matrix is recovered from the neighbour list
394 : // by subtracting the number of atoms in GROUPA as we do here.
395 11507 : if( ablocks[0]>=lstart ) {
396 : fbsize = lstart;
397 : }
398 :
399 : // Get the atoms
400 11507 : std::vector<Vector> ltmp_pos2( ntasks );
401 748752 : for(unsigned i=0; i<ntasks; ++i) {
402 737245 : ltmp_pos2[i] = ActionAtomistic::getPosition(pTaskList[i]);
403 : }
404 11507 : auto & matdata = taskmanager.getActionInput();
405 : // Now update the neighbor lists
406 11507 : if( getStep()%nl_stride==0 ) {
407 : // Build the link cells
408 11507 : std::vector<Vector> ltmp_pos( ablocks.size() );
409 1290139 : for(unsigned i=0; i<ablocks.size(); ++i) {
410 1278632 : ltmp_pos[i]=ActionAtomistic::getPosition( ablocks[i] );
411 : }
412 11507 : linkcells.createNeighborList( getConstPntrToComponent(0)->getShape()[0],
413 : make_const_view(ltmp_pos2),
414 : make_const_view(pTaskList),
415 : make_const_view(pTaskList),
416 : make_const_view(ltmp_pos),
417 : make_const_view(ablocks),
418 : getPbc(),
419 11507 : matdata.natoms_per_list,
420 11507 : matdata.nlist_v );
421 11507 : if( threeblocks.size()>0 ) {
422 89 : std::vector<Vector> ltmp_pos3( threeblocks.size() );
423 66184 : for(unsigned i=0; i<threeblocks.size(); ++i) {
424 66095 : ltmp_pos3[i]=ActionAtomistic::getPosition( threeblocks[i] );
425 : }
426 89 : threecells.createNeighborList( getConstPntrToComponent(0)->getShape()[0],
427 : make_const_view(ltmp_pos2),
428 : make_const_view(pTaskList),
429 : make_const_view(pTaskList),
430 : make_const_view(ltmp_pos3),
431 : make_const_view(threeblocks),
432 : getPbc(),
433 89 : matdata.natoms_per_three_list,
434 89 : matdata.nlist_three_v);
435 : }
436 : matdata.update();
437 : } else {
438 0 : error("neighbour list non updates are not actually implemented or tested");
439 : }
440 : // And finally work out maximum number of columns to use
441 11507 : unsigned maxcol = matdata.nlist[0];
442 1254216 : for(unsigned i=1; i<getConstPntrToComponent(0)->getShape()[0]; ++i) {
443 1242709 : if( matdata.nlist[i]>maxcol ) {
444 1374 : maxcol = matdata.nlist[i];
445 : }
446 : }
447 : // The first element returned by the neighbour list is the central atom. The
448 : // neighbour list thus always has one more atom in it than there are atoms in the
449 : // neighbourhood of the central atom
450 11507 : maxcol = maxcol-1;
451 :
452 : // Reshape the matrix store if the number of columns has changed
453 11507 : if( maxcol!=myval->getNumberOfColumns() ) {
454 1124 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
455 760 : getPntrToComponent(i)->reshapeMatrixStore( maxcol );
456 : }
457 : }
458 : // Clear the bookeeping array
459 1265723 : for(unsigned i=0; i<lstart; ++i) {
460 1254216 : myval->setMatrixBookeepingElement( i*(1+maxcol), 0 );
461 : }
462 : // Transfer neighbor list data to the bookeeping arrays
463 748752 : for(unsigned i=0; i<ntasks; ++i) {
464 737245 : unsigned bstart = pTaskList[i]*(1+maxcol);
465 737245 : unsigned rstart = lstart + pTaskList[i]*(1+matdata.natoms_per_list);
466 737245 : myval->setMatrixBookeepingElement( bstart, matdata.nlist[pTaskList[i]] - 1 );
467 47672136 : for(unsigned j=1; j<matdata.nlist[pTaskList[i]]; ++j) {
468 46934891 : myval->setMatrixBookeepingElement( bstart+j, matdata.nlist[rstart+j] - fbsize );
469 : }
470 : }
471 12266 : for(unsigned i=1; i<getNumberOfComponents(); ++i) {
472 759 : getPntrToComponent(i)->copyBookeepingArrayFromArgument( getPntrToComponent(0) );
473 : }
474 : // We need to setup the task manager here because at this point we know the size of the sparse matrices
475 11507 : unsigned nder = 6 + 3*matdata.natoms_per_three_list + virialSize;
476 : //there was an if on `matdata.components` but both branches were calling the next line
477 11507 : taskmanager.setupParallelTaskManager( nder, getNumberOfDerivatives() );
478 : // Create the workspace that we use in performTask
479 11507 : taskmanager.setWorkspaceSize(3*(maxcol + 2 + matdata.natoms_per_three_list) );
480 : // And run all the tasks
481 11507 : taskmanager.runAllTasks();
482 11507 : }
483 :
484 : template <class T>
485 1044042 : void AdjacencyMatrixBase<T>::performTask( std::size_t task_index,
486 : const AdjacencyMatrixData<T>& actiondata,
487 : ParallelActionsInput& input,
488 : ParallelActionsOutput output ) {
489 1044042 : const unsigned nneigh=actiondata.nlist[task_index];
490 1044042 : const unsigned n3neigh=(actiondata.natoms_per_three_list>0) ?
491 1148 : actiondata.nlist_three[task_index]:0;
492 1044042 : VectorView atoms(output.buffer.data(), nneigh + n3neigh );
493 :
494 1044042 : const unsigned fstart = actiondata.nlists + task_index*(1+actiondata.natoms_per_list);
495 1044042 : const View<double,3> pos0( input.inputdata + 3 * task_index);
496 55631425 : for(unsigned i=0; i<nneigh; ++i) {
497 54587383 : atoms[i][0] = input.inputdata[3*actiondata.nlist[fstart+i]+0] - pos0[0];
498 54587383 : atoms[i][1] = input.inputdata[3*actiondata.nlist[fstart+i]+1] - pos0[1];
499 54587383 : atoms[i][2] = input.inputdata[3*actiondata.nlist[fstart+i]+2] - pos0[2];
500 : }
501 :
502 1044042 : if( actiondata.natoms_per_three_list>0 ) {
503 : // Retrieve the set of third atoms
504 1148 : unsigned fstart3 = actiondata.nlists
505 1148 : + task_index*(1+actiondata.natoms_per_three_list);
506 4182123 : for(unsigned i=1; i<n3neigh; ++i) {
507 4180975 : atoms[nneigh+i-1][0] =
508 4180975 : input.inputdata[3*actiondata.nlist_three[fstart3+i]+0] - pos0[0];
509 4180975 : atoms[nneigh+i-1][1] =
510 4180975 : input.inputdata[3*actiondata.nlist_three[fstart3+i]+1] - pos0[1];
511 4180975 : atoms[nneigh+i-1][2] =
512 4180975 : input.inputdata[3*actiondata.nlist_three[fstart3+i]+2] - pos0[2];
513 : }
514 : }
515 :
516 :
517 : // Apply periodic boundary conditions to all the atoms
518 1044042 : if( actiondata.usepbc ) {
519 1044042 : input.pbc->apply( atoms, atoms.size() );
520 : }
521 1044042 : AdjacencyMatrixInput adjinp {input.noderiv,
522 1044042 : input.pbc,
523 : Vector{0.0,0.0,0.0},
524 : 0,
525 1044042 : VectorView{ atoms.data() + 3*nneigh,n3neigh}};
526 :
527 1044042 : if( n3neigh>1 ) {
528 1148 : adjinp.natoms = n3neigh-1;
529 : }
530 :
531 : // And calculate this row of the matrices
532 1044042 : std::size_t nderiv(6 + 3*adjinp.natoms + virialSize);
533 1044042 : const unsigned ncomponents =(actiondata.components) ? 4 : 1;
534 :
535 : // Must clear the derivatives here as otherwise sparsity pattern
536 : // of derivatives does not match sparsity pattern for forces
537 1044042 : if( !input.noderiv ) {
538 502293030 : for(unsigned i=0; i<output.derivatives.size(); ++i) {
539 501952314 : output.derivatives[i] = 0.0;
540 : }
541 : }
542 :
543 54587383 : for(unsigned i=1; i<nneigh; ++i ) {
544 53543341 : adjinp.pos = Vector(atoms[i][0],atoms[i][1],atoms[i][2]);
545 53543341 : const std::size_t valpos = (i-1)*ncomponents;
546 53543341 : MatrixOutput adjout{View<double,1>{&output.values[ valpos]},
547 53543341 : View{output.derivatives.data() + valpos*nderiv,nderiv}};
548 53543341 : T::calculateWeight( actiondata.matrixdata, adjinp, adjout );
549 53543341 : if( !actiondata.components ) {
550 48208922 : continue ;
551 : }
552 25445714 : output.values[valpos+1] = atoms[i][0];
553 25445714 : output.values[valpos+2] = atoms[i][1];
554 25445714 : output.values[valpos+3] = atoms[i][2];
555 25445714 : if( input.noderiv ) {
556 20111295 : continue ;
557 : }
558 : //sugar for not having to repeat [valpos*nderiv+something]
559 21337676 : for(int ii=1; ii<4; ++ii) {
560 16003257 : auto derivs = output.derivatives.subview_n<5>(valpos*nderiv+ii*nderiv);
561 16003257 : derivs[0] = -1.0;
562 16003257 : derivs[1] = 1.0;
563 16003257 : derivs[2] = -atoms[i][0];
564 16003257 : derivs[3] = -atoms[i][1];
565 16003257 : derivs[4] = -atoms[i][2];
566 : }
567 :
568 : }
569 1044042 : }
570 :
571 : template <class T>
572 10861 : void AdjacencyMatrixBase<T>::applyNonZeroRankForces( std::vector<double>& outforces ) {
573 10861 : taskmanager.applyForces( outforces );
574 10861 : }
575 :
576 : template <class T>
577 340716 : int AdjacencyMatrixBase<T>::getNumberOfValuesPerTask( std::size_t task_index,
578 : const AdjacencyMatrixData<T>& actiondata ) {
579 340716 : return actiondata.nlist[task_index] - 1;
580 : }
581 :
582 : template <class T>
583 12704465 : void AdjacencyMatrixBase<T>::getForceIndices( const std::size_t task_index,
584 : const std::size_t colno,
585 : const std::size_t ntotal_force,
586 : const AdjacencyMatrixData<T>& actiondata,
587 : const ParallelActionsInput& /*input*/,
588 : ForceIndexHolder force_indices ) {
589 :
590 12704465 : const unsigned fstart = actiondata.nlists
591 12704465 : + task_index*(1 + actiondata.natoms_per_list);
592 :
593 12704465 : const auto three_task_index= 3 * task_index;
594 12704465 : force_indices.indices[0][0] = three_task_index;
595 12704465 : force_indices.indices[0][1] = three_task_index + 1;
596 12704465 : force_indices.indices[0][2] = three_task_index + 2;
597 :
598 12704465 : const unsigned myatom = 3 * actiondata.nlist[fstart + 1 + colno];
599 12704465 : force_indices.indices[0][3] = myatom;
600 12704465 : force_indices.indices[0][4] = myatom + 1;
601 12704465 : force_indices.indices[0][5] = myatom + 2;
602 12704465 : if( actiondata.components ) {
603 5334419 : force_indices.indices[1][0] = three_task_index;
604 5334419 : force_indices.indices[1][1] = myatom;
605 5334419 : force_indices.indices[2][0] = three_task_index + 1;
606 5334419 : force_indices.indices[2][1] = myatom + 1;
607 5334419 : force_indices.indices[3][0] = three_task_index + 2;
608 5334419 : force_indices.indices[3][1] = myatom + 2;
609 : }
610 :
611 12704465 : force_indices.threadsafe_derivatives_end[0] = 0;
612 : unsigned n = 6;
613 12704465 : if( actiondata.natoms_per_three_list>0 ) {
614 7754 : const unsigned n3neigh = actiondata.nlist_three[task_index];
615 7754 : const unsigned fstart3 = actiondata.nlists
616 7754 : + task_index*(1 + actiondata.natoms_per_three_list);
617 905962 : for(unsigned j=1; j<n3neigh; ++j) {
618 898208 : unsigned my3atom = 3 * actiondata.nlist_three[fstart3 + j];
619 898208 : force_indices.indices[0][n ] = my3atom;
620 898208 : force_indices.indices[0][n+1] = my3atom + 1;
621 898208 : force_indices.indices[0][n+2] = my3atom + 2;
622 898208 : n += 3;
623 : }
624 : }
625 12704465 : const unsigned virstart = ntotal_force - 9;
626 12704465 : force_indices.indices[0][n ] = virstart + 0;
627 12704465 : force_indices.indices[0][n+1] = virstart + 1;
628 12704465 : force_indices.indices[0][n+2] = virstart + 2;
629 12704465 : force_indices.indices[0][n+3] = virstart + 3;
630 12704465 : force_indices.indices[0][n+4] = virstart + 4;
631 12704465 : force_indices.indices[0][n+5] = virstart + 5;
632 12704465 : force_indices.indices[0][n+6] = virstart + 6;
633 12704465 : force_indices.indices[0][n+7] = virstart + 7;
634 12704465 : force_indices.indices[0][n+8] = virstart + 8;
635 :
636 12704465 : force_indices.tot_indices[0] = n + 9;
637 12704465 : if( actiondata.components ) {
638 5334419 : force_indices.threadsafe_derivatives_end[1] = 0;
639 5334419 : force_indices.threadsafe_derivatives_end[2] = 0;
640 5334419 : force_indices.threadsafe_derivatives_end[3] = 0;
641 5334419 : force_indices.indices[1][2] = virstart;
642 5334419 : force_indices.indices[1][3] = virstart + 1;
643 5334419 : force_indices.indices[1][4] = virstart + 2;
644 5334419 : force_indices.indices[2][2] = virstart + 3;
645 5334419 : force_indices.indices[2][3] = virstart + 4;
646 5334419 : force_indices.indices[2][4] = virstart + 5;
647 5334419 : force_indices.indices[3][2] = virstart + 6;
648 5334419 : force_indices.indices[3][3] = virstart + 7;
649 5334419 : force_indices.indices[3][4] = virstart + 8;
650 5334419 : force_indices.tot_indices[1] = 5;
651 5334419 : force_indices.tot_indices[2] = 5;
652 5334419 : force_indices.tot_indices[3] = 5;
653 : }
654 12704465 : }
655 :
656 : } // namespace adjmat
657 : } // namespace PLMD
658 :
659 : #endif //__PLUMED_adjmat_AdjacencyMatrixBase_h
|
