Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "MultiColvarBase.h"
23 : #include "ActionVolume.h"
24 : #include "MultiColvarFilter.h"
25 : #include "vesselbase/Vessel.h"
26 : #include "vesselbase/BridgeVessel.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/ActionSet.h"
29 : #include "tools/Pbc.h"
30 : #include "AtomValuePack.h"
31 : #include <vector>
32 : #include <string>
33 : #include <limits>
34 :
35 : namespace PLMD {
36 : namespace multicolvar {
37 :
38 651 : void MultiColvarBase::registerKeywords( Keywords& keys ) {
39 651 : Action::registerKeywords( keys );
40 651 : ActionWithValue::registerKeywords( keys );
41 651 : ActionAtomistic::registerKeywords( keys );
42 1302 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
43 651 : ActionWithVessel::registerKeywords( keys );
44 1302 : keys.add("hidden","NL_STRIDE","the frequency with which the neighbor list should be updated. Between neighbour list update steps all quantities "
45 : "that contributed less than TOL at the previous neighbor list update step are ignored.");
46 651 : keys.setComponentsIntroduction("When the label of this action is used as the input for a second you are not referring to a scalar quantity as you are in "
47 : "regular collective variables. The label is used to reference the full set of quantities calculated by "
48 : "the action. This is usual when using \\ref multicolvarfunction. Generally when doing this the previously calculated "
49 : "multicolvar will be referenced using the DATA keyword rather than ARG.\n\n"
50 : "This Action can be used to calculate the following scalar quantities directly. These quantities are calculated by "
51 : "employing the keywords listed below. "
52 : "These quantities can then be referenced elsewhere in the input file by using this Action's label "
53 : "followed by a dot and the name of the quantity. Some of them can be calculated multiple times "
54 : "with different parameters. In this case the quantities calculated can be referenced elsewhere in the "
55 : "input by using the name of the quantity followed by a numerical identifier "
56 : "e.g. <em>label</em>.lessthan-1, <em>label</em>.lessthan-2 etc. When doing this and, for clarity we have "
57 : "made it so that the user can set a particular label for each of the components. As such by using the LABEL keyword in the description of the keyword "
58 : "input you can customize the component name");
59 1302 : keys.reserve("atoms-3","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
60 : "one coordination number for each of the atoms specified. Each of these coordination numbers specifies how many of the "
61 : "other specified atoms are within a certain cutoff of the central atom. You can specify the atoms here as another multicolvar "
62 : "action or using a MultiColvarFilter or ActionVolume action. When you do so the quantity is calculated for those atoms specified "
63 : "in the previous multicolvar. This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
64 : "coordination number more than four for example");
65 1302 : keys.reserve("atoms-4","SPECIESA","this keyword is used for colvars such as the coordination number. In that context it species that plumed should calculate "
66 : "one coordination number for each of the atoms specified in SPECIESA. Each of these coordination numbers specifies how many "
67 : "of the atoms specifies using SPECIESB is within the specified cutoff. As with the species keyword the input can also be specified "
68 : "using the label of another multicolvar");
69 1302 : keys.reserve("atoms-4","SPECIESB","this keyword is used for colvars such as the coordination number. It must appear with SPECIESA. For a full explanation see "
70 : "the documentation for that keyword");
71 1302 : keys.add("hidden","ALL_INPUT_SAME_TYPE","remove this keyword to remove certain checks in the input on the sanity of your input file. See code for details");
72 651 : }
73 :
74 351 : MultiColvarBase::MultiColvarBase(const ActionOptions&ao):
75 : Action(ao),
76 : ActionAtomistic(ao),
77 : ActionWithValue(ao),
78 : ActionWithVessel(ao),
79 351 : usepbc(false),
80 351 : allthirdblockintasks(false),
81 351 : uselinkforthree(false),
82 351 : linkcells(comm),
83 351 : threecells(comm),
84 351 : setup_completed(false),
85 351 : atomsWereRetrieved(false),
86 351 : matsums(false),
87 351 : usespecies(false),
88 702 : nblock(0) {
89 702 : if( keywords.exists("NOPBC") ) {
90 351 : bool nopbc=!usepbc;
91 351 : parseFlag("NOPBC",nopbc);
92 351 : usepbc=!nopbc;
93 : }
94 702 : if( keywords.exists("SPECIESA") ) {
95 101 : matsums=usespecies=true;
96 : }
97 351 : }
98 :
99 48 : void MultiColvarBase::readAtomsLikeKeyword( const std::string & key, const int& natoms, std::vector<AtomNumber>& all_atoms ) {
100 48 : plumed_assert( !usespecies );
101 48 : if( all_atoms.size()>0 ) {
102 0 : return;
103 : }
104 :
105 : std::vector<AtomNumber> t;
106 48 : for(int i=1;; ++i ) {
107 970 : parseAtomList(key, i, t );
108 970 : if( t.empty() ) {
109 : break;
110 : }
111 :
112 922 : log.printf(" Colvar %d is calculated from atoms : ", i);
113 3450 : for(unsigned j=0; j<t.size(); ++j) {
114 2528 : log.printf("%d ",t[j].serial() );
115 : }
116 922 : log.printf("\n");
117 :
118 922 : if( i==1 && natoms<0 ) {
119 6 : ablocks.resize(t.size());
120 916 : } else if( i==1 ) {
121 42 : ablocks.resize(natoms);
122 : }
123 922 : if( t.size()!=ablocks.size() ) {
124 : std::string ss;
125 0 : Tools::convert(i,ss);
126 0 : error(key + ss + " keyword has the wrong number of atoms");
127 : }
128 3450 : for(unsigned j=0; j<ablocks.size(); ++j) {
129 2528 : ablocks[j].push_back( ablocks.size()*(i-1)+j );
130 2528 : all_atoms.push_back( t[j] );
131 2528 : atom_lab.push_back( std::pair<unsigned,unsigned>( 0, ablocks.size()*(i-1)+j ) );
132 : }
133 922 : t.resize(0);
134 922 : }
135 48 : if( all_atoms.size()>0 ) {
136 48 : nblock=0;
137 970 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
138 922 : addTaskToList( i );
139 : }
140 : }
141 : }
142 :
143 471 : bool MultiColvarBase::parseMultiColvarAtomList(const std::string& key, const int& num, std::vector<AtomNumber>& t) {
144 : std::vector<std::string> mlabs;
145 471 : if( num<0 ) {
146 471 : parseVector(key,mlabs);
147 : } else {
148 0 : parseNumberedVector(key,num,mlabs);
149 : }
150 :
151 471 : if( mlabs.size()==0 ) {
152 : return false;
153 : }
154 :
155 : std::string mname;
156 323 : unsigned found_mcolv=mlabs.size();
157 450 : for(unsigned i=0; i<mlabs.size(); ++i) {
158 328 : MultiColvarBase* mycolv = plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlabs[i]);
159 328 : if(!mycolv) {
160 : found_mcolv=i;
161 201 : break;
162 : }
163 : // Check all base multicolvars are of same type
164 127 : if( i==0 ) {
165 122 : mname = mycolv->getName();
166 244 : if( keywords.exists("ALL_INPUT_SAME_TYPE") && mycolv->isPeriodic() ) {
167 0 : error("multicolvar functions don't work with this multicolvar");
168 : }
169 : } else {
170 10 : if( keywords.exists("ALL_INPUT_SAME_TYPE") && mname!=mycolv->getName() ) {
171 0 : error("All input multicolvars must be of same type");
172 : }
173 : }
174 : // And track which variable stores each colvar
175 4458438 : for(unsigned j=0; j<mycolv->getFullNumberOfTasks(); ++j) {
176 4458311 : atom_lab.push_back( std::pair<unsigned,unsigned>( mybasemulticolvars.size()+1, j ) );
177 : }
178 : // And store the multicolvar base
179 127 : mybasemulticolvars.push_back( mycolv );
180 : // And create a basedata stash
181 127 : mybasedata.push_back( mybasemulticolvars[mybasemulticolvars.size()-1]->buildDataStashes( this ) );
182 : // Check if weight has derivatives
183 127 : if( mybasemulticolvars[ mybasemulticolvars.size()-1 ]->weightHasDerivatives ) {
184 26 : weightHasDerivatives=true;
185 : }
186 127 : plumed_assert( mybasemulticolvars.size()==mybasedata.size() );
187 : }
188 : // Have we conventional atoms to read in
189 323 : if( found_mcolv==0 ) {
190 : std::vector<AtomNumber> tt;
191 201 : ActionAtomistic::interpretAtomList( mlabs, tt );
192 89846 : for(unsigned i=0; i<tt.size(); ++i) {
193 89645 : atom_lab.push_back( std::pair<unsigned,unsigned>( 0, t.size() + i ) );
194 : }
195 402 : log.printf(" keyword %s takes atoms : ", key.c_str() );
196 89846 : for(unsigned i=0; i<tt.size(); ++i) {
197 89645 : t.push_back( tt[i] );
198 89645 : log.printf("%d ",tt[i].serial() );
199 : }
200 201 : log.printf("\n");
201 122 : } else if( found_mcolv==mlabs.size() ) {
202 122 : if( checkNumericalDerivatives() ) {
203 0 : error("cannot use NUMERICAL_DERIVATIVES keyword with dynamic groups as input");
204 : }
205 122 : log.printf(" keyword %s takes dynamic groups of atoms constructed from multicolvars labelled : ", key.c_str() );
206 249 : for(unsigned i=0; i<mlabs.size(); ++i) {
207 127 : log.printf("%s ",mlabs[i].c_str() );
208 : }
209 122 : log.printf("\n");
210 0 : } else if( found_mcolv<mlabs.size() ) {
211 0 : error("cannot mix multicolvar input and atom input in one line");
212 : }
213 : return true;
214 471 : }
215 :
216 76 : void MultiColvarBase::readTwoGroups( const std::string& key0, const std::string& key1, const std::string& key2, std::vector<AtomNumber>& all_atoms ) {
217 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) );
218 76 : ablocks.resize( 2 );
219 :
220 76 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
221 26 : nblock=atom_lab.size();
222 78 : for(unsigned i=0; i<2; ++i) {
223 52 : ablocks[i].resize(nblock);
224 : }
225 26 : resizeBookeepingArray( nblock, nblock );
226 5486 : for(unsigned i=0; i<nblock; ++i) {
227 5460 : ablocks[0][i]=ablocks[1][i]=i;
228 : }
229 5460 : for(unsigned i=1; i<nblock; ++i) {
230 4216329 : for(unsigned j=0; j<i; ++j) {
231 4210895 : bookeeping(i,j).first=getFullNumberOfTasks();
232 4210895 : addTaskToList( i*nblock + j );
233 4210895 : bookeeping(i,j).second=getFullNumberOfTasks();
234 : }
235 : }
236 : } else {
237 50 : parseMultiColvarAtomList(key1,-1,all_atoms);
238 50 : ablocks[0].resize( atom_lab.size() );
239 67 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
240 17 : ablocks[0][i]=i;
241 : }
242 50 : parseMultiColvarAtomList(key2,-1,all_atoms);
243 50 : ablocks[1].resize( atom_lab.size() - ablocks[0].size() );
244 1119 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
245 1069 : ablocks[1][i]=ablocks[0].size() + i;
246 : }
247 :
248 50 : if( ablocks[0].size()>ablocks[1].size() ) {
249 0 : nblock = ablocks[0].size();
250 : } else {
251 50 : nblock=ablocks[1].size();
252 : }
253 :
254 50 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
255 67 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
256 1126 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
257 1109 : bookeeping(i,j).first=getFullNumberOfTasks();
258 1109 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 ) {
259 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
260 0 : atom_lab[ablocks[0][i]].second!=atom_lab[ablocks[1][j]].second ) {
261 0 : addTaskToList( i*nblock + j );
262 : }
263 1109 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] ) {
264 1109 : addTaskToList( i*nblock + j );
265 : }
266 1109 : bookeeping(i,j).second=getFullNumberOfTasks();
267 : }
268 : }
269 : }
270 76 : }
271 :
272 12 : void MultiColvarBase::readGroupKeywords( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& key3,
273 : const bool& no_third_dim_accum, const bool& symmetric, std::vector<AtomNumber>& all_atoms ) {
274 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) );
275 12 : ablocks.resize( 3 );
276 :
277 12 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
278 10 : if( no_third_dim_accum ) {
279 10 : nblock=atom_lab.size();
280 10 : ablocks[0].resize(nblock);
281 10 : ablocks[1].resize( nblock );
282 1481 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
283 1471 : ablocks[0][i]=ablocks[1][i]=i;
284 : }
285 10 : resizeBookeepingArray( nblock, nblock );
286 10 : if( symmetric ) {
287 : // This ensures that later parts of the code don't switch off allthirdblockintasks
288 1343 : for(unsigned i=0; i<nblock; ++i) {
289 1337 : bookeeping(i,i).first=0;
290 1337 : bookeeping(i,i).second=1;
291 : }
292 1337 : for(unsigned i=1; i<nblock; ++i) {
293 222447 : for(unsigned j=0; j<i; ++j) {
294 221116 : bookeeping(j,i).first=bookeeping(i,j).first=getFullNumberOfTasks();
295 221116 : addTaskToList( i*nblock + j );
296 221116 : bookeeping(j,i).second=bookeeping(i,j).second=getFullNumberOfTasks();
297 : }
298 : }
299 : } else {
300 138 : for(unsigned i=0; i<nblock; ++i) {
301 8344 : for(unsigned j=0; j<nblock; ++j) {
302 8210 : if( i==j ) {
303 134 : continue ;
304 : }
305 8076 : bookeeping(i,j).first=getFullNumberOfTasks();
306 8076 : addTaskToList( i*nblock + j );
307 8076 : bookeeping(i,j).second=getFullNumberOfTasks();
308 : }
309 : }
310 : }
311 10 : if( !parseMultiColvarAtomList(key3,-1,all_atoms) ) {
312 0 : error("missing required keyword " + key3 + " in input");
313 : }
314 10 : ablocks[2].resize( atom_lab.size() - ablocks[0].size() );
315 49845 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
316 49835 : ablocks[2][i]=ablocks[0].size() + i;
317 : }
318 : } else {
319 0 : nblock=atom_lab.size();
320 0 : for(unsigned i=0; i<3; ++i) {
321 0 : ablocks[i].resize(nblock);
322 : }
323 0 : resizeBookeepingArray( nblock, nblock );
324 0 : for(unsigned i=0; i<nblock; ++i) {
325 0 : ablocks[0][i]=i;
326 0 : ablocks[1][i]=i;
327 0 : ablocks[2][i]=i;
328 : }
329 0 : if( symmetric ) {
330 0 : for(unsigned i=2; i<nblock; ++i) {
331 0 : for(unsigned j=1; j<i; ++j) {
332 0 : bookeeping(i,j).first=getFullNumberOfTasks();
333 0 : for(unsigned k=0; k<j; ++k) {
334 0 : addTaskToList( i*nblock*nblock + j*nblock + k );
335 : }
336 0 : bookeeping(i,j).second=getFullNumberOfTasks();
337 : }
338 : }
339 : } else {
340 0 : for(unsigned i=0; i<nblock; ++i) {
341 0 : for(unsigned j=0; j<nblock; ++j) {
342 0 : if( i==j ) {
343 0 : continue;
344 : }
345 0 : bookeeping(i,j).first=getFullNumberOfTasks();
346 0 : for(unsigned k=0; k<nblock; ++k) {
347 0 : if( i!=k && j!=k ) {
348 0 : addTaskToList( i*nblock*nblock + j*nblock + k );
349 : }
350 : }
351 0 : bookeeping(i,j).first=getFullNumberOfTasks();
352 : }
353 : }
354 : }
355 : }
356 : } else {
357 2 : readThreeGroups( key1, key2, key3, true, no_third_dim_accum, all_atoms );
358 : }
359 :
360 12 : }
361 :
362 10 : void MultiColvarBase::readThreeGroups( const std::string& key1, const std::string& key2, const std::string& key3,
363 : const bool& allow2, const bool& no_third_dim_accum, std::vector<AtomNumber>& all_atoms ) {
364 : plumed_dbg_assert( keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) );
365 10 : ablocks.resize( 3 );
366 :
367 10 : bool readkey1=parseMultiColvarAtomList(key1,-1,all_atoms);
368 10 : ablocks[0].resize( atom_lab.size() );
369 31 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
370 21 : ablocks[0][i]=i;
371 : }
372 10 : bool readkey2=parseMultiColvarAtomList(key2,-1,all_atoms);
373 10 : if( !readkey1 && !readkey2 ) {
374 : return ;
375 : }
376 10 : ablocks[1].resize( atom_lab.size() - ablocks[0].size() );
377 225 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
378 215 : ablocks[1][i]=ablocks[0].size() + i;
379 : }
380 :
381 10 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
382 10 : bool readkey3=parseMultiColvarAtomList(key3,-1,all_atoms);
383 10 : if( !readkey3 && !allow2 ) {
384 0 : error("missing atom specification " + key3);
385 10 : } else if( !readkey3 ) {
386 2 : if( ablocks[1].size()>ablocks[0].size() ) {
387 2 : nblock=ablocks[1].size();
388 : } else {
389 0 : nblock=ablocks[0].size();
390 : }
391 :
392 2 : ablocks[2].resize( ablocks[1].size() );
393 200 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
394 198 : ablocks[2][i]=ablocks[0].size() + i;
395 : }
396 4 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
397 198 : for(unsigned j=1; j<ablocks[1].size(); ++j) {
398 196 : bookeeping(i,j).first=getFullNumberOfTasks();
399 9898 : for(unsigned k=0; k<j; ++k) {
400 9702 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
401 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
402 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
403 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
404 0 : atom_lab[ablocks[0][i]].second!=atom_lab[ablocks[1][j]].second && atom_lab[ablocks[0][i]].second!=atom_lab[ablocks[2][k]].second &&
405 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) {
406 0 : addTaskToList( nblock*nblock*i + nblock*j + k );
407 : }
408 9702 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
409 9702 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
410 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) {
411 9702 : addTaskToList( nblock*nblock*i + nblock*j + k );
412 : }
413 : }
414 196 : bookeeping(i,j).second=getFullNumberOfTasks();
415 : }
416 : }
417 : } else {
418 8 : ablocks[2].resize( atom_lab.size() - ablocks[1].size() - ablocks[0].size() );
419 3182 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
420 3174 : ablocks[2][i] = ablocks[0].size() + ablocks[1].size() + i;
421 : }
422 :
423 8 : if( ablocks[1].size()>ablocks[0].size() ) {
424 0 : nblock=ablocks[1].size();
425 : } else {
426 8 : nblock=ablocks[0].size();
427 : }
428 8 : if( ablocks[2].size()>nblock ) {
429 5 : nblock=ablocks[2].size();
430 : }
431 :
432 : unsigned kcount;
433 8 : if( no_third_dim_accum ) {
434 : kcount=1;
435 : } else {
436 0 : kcount=ablocks[2].size();
437 : }
438 :
439 27 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
440 128 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
441 109 : bookeeping(i,j).first=getFullNumberOfTasks();
442 218 : for(unsigned k=0; k<kcount; ++k) {
443 109 : if( no_third_dim_accum ) {
444 109 : addTaskToList( nblock*i + j );
445 0 : } else if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
446 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
447 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
448 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
449 0 : atom_lab[ablocks[0][i]].second!=atom_lab[ablocks[1][j]].second && atom_lab[ablocks[0][i]].second!=atom_lab[ablocks[2][k]].second &&
450 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) {
451 0 : addTaskToList( nblock*nblock*i + nblock*j + k );
452 : }
453 0 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
454 0 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
455 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) {
456 0 : addTaskToList( nblock*nblock*i + nblock*j + k );
457 : }
458 : }
459 109 : bookeeping(i,j).second=getFullNumberOfTasks();
460 : }
461 : }
462 : }
463 : }
464 :
465 4 : void MultiColvarBase::buildSets() {
466 : std::vector<AtomNumber> fake_atoms;
467 8 : if( !parseMultiColvarAtomList("DATA",-1,fake_atoms) ) {
468 0 : error("missing DATA keyword");
469 : }
470 4 : if( fake_atoms.size()>0 ) {
471 0 : error("no atoms should appear in the specification for this object. Input should be other multicolvars");
472 : }
473 :
474 4 : nblock = mybasemulticolvars[0]->getFullNumberOfTasks();
475 11 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
476 7 : if( mybasemulticolvars[i]->getFullNumberOfTasks()!=nblock ) {
477 0 : error("mismatch between numbers of tasks in various base multicolvars");
478 : }
479 : }
480 4 : ablocks.resize( mybasemulticolvars.size() );
481 4 : usespecies=false;
482 11 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
483 7 : ablocks[i].resize( nblock );
484 27 : for(unsigned j=0; j<nblock; ++j) {
485 20 : ablocks[i][j]=i*nblock+j;
486 : }
487 : }
488 15 : for(unsigned i=0; i<nblock; ++i) {
489 11 : if( mybasemulticolvars.size()<4 ) {
490 11 : unsigned cvcode=0, tmpc=1;
491 31 : for(unsigned j=0; j<ablocks.size(); ++j) {
492 20 : cvcode += i*tmpc;
493 20 : tmpc *= nblock;
494 : }
495 11 : addTaskToList( cvcode );
496 : } else {
497 0 : addTaskToList( i );
498 : }
499 : }
500 4 : mybasedata[0]->resizeTemporyMultiValues( mybasemulticolvars.size() );
501 4 : setupMultiColvarBase( fake_atoms );
502 4 : }
503 :
504 12512702 : void MultiColvarBase::addTaskToList( const unsigned& taskCode ) {
505 12512702 : plumed_assert( getNumberOfVessels()==0 );
506 12512702 : ActionWithVessel::addTaskToList( taskCode );
507 12512702 : }
508 :
509 96 : void MultiColvarBase::resizeBookeepingArray( const unsigned& num1, const unsigned& num2 ) {
510 96 : bookeeping.resize( num1, num2 );
511 7065 : for(unsigned i=0; i<num1; ++i) {
512 8887414 : for(unsigned j=0; j<num2; ++j) {
513 8880445 : bookeeping(i,j).first=0;
514 8880445 : bookeeping(i,j).second=0;
515 : }
516 : }
517 96 : }
518 :
519 309 : void MultiColvarBase::setupMultiColvarBase( const std::vector<AtomNumber>& atoms ) {
520 309 : if( !matsums && atom_lab.size()==0 ) {
521 0 : error("No atoms have been read in");
522 : }
523 : std::vector<AtomNumber> all_atoms;
524 : // Setup decoder array
525 309 : if( !usespecies && nblock>0 ) {
526 :
527 63 : ncentral=ablocks.size();
528 63 : use_for_central_atom.resize( ablocks.size(), true );
529 63 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
530 63 : if( ablocks.size()==3 ) {
531 20 : allthirdblockintasks=uselinkforthree=true;
532 1512 : for(unsigned i=0; i<bookeeping.nrows(); ++i) {
533 453578 : for(unsigned j=0; j<bookeeping.ncols(); ++j) {
534 452086 : unsigned ntper = bookeeping(i,j).second - bookeeping(i,j).first;
535 452086 : if( i==j && ntper==0 ) {
536 136 : continue;
537 451950 : } else if( ntper == 1 && allthirdblockintasks ) {
538 451756 : allthirdblockintasks=true;
539 194 : } else if( ntper != ablocks[2].size() ) {
540 194 : allthirdblockintasks=uselinkforthree=false;
541 : } else {
542 0 : allthirdblockintasks=false;
543 : }
544 : }
545 : }
546 : }
547 :
548 63 : if( allthirdblockintasks ) {
549 18 : decoder.resize(2);
550 18 : plumed_assert( ablocks.size()==3 );
551 : // Check if number of atoms is too large
552 18 : if( std::pow( double(nblock), 2.0 )>std::numeric_limits<unsigned>::max() ) {
553 0 : error("number of atoms in groups is too big for PLUMED to handle");
554 : }
555 : } else {
556 45 : decoder.resize( ablocks.size() );
557 : // Check if number of atoms is too large
558 45 : if( std::pow( double(nblock), double(ablocks.size()) )>std::numeric_limits<unsigned>::max() ) {
559 0 : error("number of atoms in groups is too big for PLUMED to handle");
560 : }
561 : }
562 : unsigned code=1;
563 190 : for(unsigned i=0; i<decoder.size(); ++i) {
564 127 : decoder[decoder.size()-1-i]=code;
565 127 : code *= nblock;
566 : }
567 246 : } else if( !usespecies ) {
568 87 : ncentral=ablocks.size();
569 87 : use_for_central_atom.resize( ablocks.size(), true );
570 87 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
571 318 : } else if( keywords.exists("SPECIESA") ) {
572 101 : plumed_assert( atom_lab.size()==0 && all_atoms.size()==0 );
573 101 : ablocks.resize( 1 );
574 101 : bool readspecies=parseMultiColvarAtomList("SPECIES", -1, all_atoms);
575 101 : if( readspecies ) {
576 81 : ablocks[0].resize( atom_lab.size() );
577 41493 : for(unsigned i=0; i<atom_lab.size(); ++i) {
578 41412 : addTaskToList(i);
579 41412 : ablocks[0][i]=i;
580 : }
581 : } else {
582 40 : if( !parseMultiColvarAtomList("SPECIESA", -1, all_atoms) ) {
583 0 : error("missing SPECIES/SPECIESA keyword");
584 : }
585 20 : unsigned nat1=atom_lab.size();
586 40 : if( !parseMultiColvarAtomList("SPECIESB", -1, all_atoms) ) {
587 0 : error("missing SPECIESB keyword");
588 : }
589 20 : unsigned nat2=atom_lab.size() - nat1;
590 :
591 759 : for(unsigned i=0; i<nat1; ++i) {
592 739 : addTaskToList(i);
593 : }
594 20 : ablocks[0].resize( nat2 );
595 3726 : for(unsigned i=0; i<nat2; ++i) {
596 : bool found=false;
597 : unsigned inum=0;
598 288729 : for(unsigned j=0; j<nat1; ++j) {
599 285201 : if( atom_lab[nat1+i].first>0 && atom_lab[j].first>0 ) {
600 252720 : if( atom_lab[nat1+i].first==atom_lab[j].first &&
601 0 : mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
602 252720 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) {
603 : found=true;
604 : inum=j;
605 : break;
606 : }
607 32481 : } else if( atom_lab[nat1+i].first>0 ) {
608 0 : if( mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
609 0 : all_atoms[atom_lab[j].second] ) {
610 : found=true;
611 : inum=nat1 + i;
612 : break;
613 : }
614 32481 : } else if( atom_lab[j].first>0 ) {
615 0 : if( all_atoms[atom_lab[nat1+i].second]==
616 0 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) {
617 : found=true;
618 : inum=nat1+i;
619 : break;
620 : }
621 32481 : } else if( all_atoms[atom_lab[nat1+i].second]==all_atoms[atom_lab[j].second] ) {
622 : found=true;
623 : inum=j;
624 : break;
625 : }
626 : }
627 : // This prevents mistakes being made in colvar setup
628 : if( found ) {
629 178 : ablocks[0][i]=inum;
630 : } else {
631 3528 : ablocks[0][i]=nat1 + i;
632 : }
633 : }
634 : }
635 : }
636 309 : if( mybasemulticolvars.size()>0 ) {
637 246 : for(unsigned i=0; i<mybasedata.size(); ++i) {
638 127 : mybasedata[i]->resizeTemporyMultiValues(2);
639 127 : mybasemulticolvars[i]->my_tmp_capacks.resize(2);
640 : }
641 : }
642 :
643 : // Copy lists of atoms involved from base multicolvars
644 : std::vector<AtomNumber> tmp_atoms;
645 436 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
646 127 : tmp_atoms=mybasemulticolvars[i]->getAbsoluteIndexes();
647 416316 : for(unsigned j=0; j<tmp_atoms.size(); ++j) {
648 416189 : all_atoms.push_back( tmp_atoms[j] );
649 : }
650 : }
651 : // Copy atom lists from input
652 59576 : for(unsigned i=0; i<atoms.size(); ++i) {
653 59267 : all_atoms.push_back( atoms[i] );
654 : }
655 :
656 : // Now make sure we get all the atom positions
657 309 : ActionAtomistic::requestAtoms( all_atoms );
658 : // And setup dependencies
659 436 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
660 127 : addDependency( mybasemulticolvars[i] );
661 : }
662 :
663 : // Setup underlying ActionWithVessel
664 309 : readVesselKeywords();
665 309 : }
666 :
667 19 : void MultiColvarBase::setAtomsForCentralAtom( const std::vector<bool>& catom_ind ) {
668 : unsigned nat=0;
669 19 : plumed_assert( catom_ind.size()==ablocks.size() );
670 89 : for(unsigned i=0; i<catom_ind.size(); ++i) {
671 70 : use_for_central_atom[i]=catom_ind[i];
672 70 : if( use_for_central_atom[i] ) {
673 28 : nat++;
674 : }
675 : }
676 : plumed_dbg_assert( nat>0 );
677 19 : ncentral=nat;
678 19 : numberForCentralAtom = 1.0 / static_cast<double>( nat );
679 19 : }
680 :
681 429 : void MultiColvarBase::turnOnDerivatives() {
682 429 : ActionWithValue::turnOnDerivatives();
683 429 : needsDerivatives();
684 429 : forcesToApply.resize( getNumberOfDerivatives() );
685 429 : }
686 :
687 192 : void MultiColvarBase::setLinkCellCutoff( const double& lcut, double tcut ) {
688 192 : plumed_assert( usespecies || ablocks.size()<4 );
689 192 : if( tcut<0 ) {
690 186 : tcut=lcut;
691 : }
692 :
693 192 : if( !linkcells.enabled() ) {
694 192 : linkcells.setCutoff( lcut );
695 192 : threecells.setCutoff( tcut );
696 : } else {
697 0 : if( lcut>linkcells.getCutoff() ) {
698 0 : linkcells.setCutoff( lcut );
699 : }
700 0 : if( tcut>threecells.getCutoff() ) {
701 0 : threecells.setCutoff( tcut );
702 : }
703 : }
704 192 : }
705 :
706 8 : double MultiColvarBase::getLinkCellCutoff() const {
707 8 : return linkcells.getCutoff();
708 : }
709 :
710 1939 : void MultiColvarBase::setupLinkCells() {
711 1939 : if( (!usespecies && nblock==0) || !linkcells.enabled() ) {
712 : return ;
713 : }
714 : // Retrieve any atoms that haven't already been retrieved
715 1319 : for(std::vector<MultiColvarBase*>::iterator p=mybasemulticolvars.begin(); p!=mybasemulticolvars.end(); ++p) {
716 203 : (*p)->retrieveAtoms();
717 : }
718 1116 : retrieveAtoms();
719 :
720 : unsigned iblock;
721 1116 : if( usespecies ) {
722 : iblock=0;
723 217 : } else if( ablocks.size()<4 ) {
724 : iblock=1;
725 : } else {
726 0 : plumed_error();
727 : }
728 :
729 : // Count number of currently active atoms
730 1116 : nactive_atoms=0;
731 387954 : for(unsigned i=0; i<ablocks[iblock].size(); ++i) {
732 386838 : if( isCurrentlyActive( ablocks[iblock][i] ) ) {
733 381822 : nactive_atoms++;
734 : }
735 : }
736 :
737 1116 : if( nactive_atoms>0 ) {
738 1116 : std::vector<Vector> ltmp_pos( nactive_atoms );
739 1116 : std::vector<unsigned> ltmp_ind( nactive_atoms );
740 :
741 1116 : nactive_atoms=0;
742 1116 : if( usespecies ) {
743 366639 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
744 365740 : if( !isCurrentlyActive( ablocks[0][i] ) ) {
745 96 : continue;
746 : }
747 365644 : ltmp_ind[nactive_atoms]=ablocks[0][i];
748 365644 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ltmp_ind[nactive_atoms] );
749 365644 : nactive_atoms++;
750 : }
751 : } else {
752 21315 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
753 21098 : if( !isCurrentlyActive( ablocks[1][i] ) ) {
754 4920 : continue;
755 : }
756 16178 : ltmp_ind[nactive_atoms]=i;
757 16178 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ablocks[1][i] );
758 16178 : nactive_atoms++;
759 : }
760 : }
761 :
762 : // Build the lists for the link cells
763 1116 : linkcells.buildCellLists( ltmp_pos, ltmp_ind, getPbc() );
764 : }
765 : }
766 :
767 4628 : void MultiColvarBase::setupNonUseSpeciesLinkCells( const unsigned& my_always_active ) {
768 4628 : plumed_assert( !usespecies );
769 4628 : if( nblock==0 || !linkcells.enabled() ) {
770 4418 : return ;
771 : }
772 210 : deactivateAllTasks();
773 : std::vector<unsigned> requiredlinkcells;
774 :
775 210 : if( !uselinkforthree && nactive_atoms>0 ) {
776 : // Get some parallel info
777 156 : unsigned stride=comm.Get_size();
778 156 : unsigned rank=comm.Get_rank();
779 156 : if( serialCalculation() ) {
780 : stride=1;
781 : rank=0;
782 : }
783 :
784 : // Ensure we only do tasks where atoms are in appropriate link cells
785 156 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
786 5757 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
787 5601 : if( !isCurrentlyActive( ablocks[0][i] ) ) {
788 2979 : continue;
789 : }
790 2622 : unsigned natomsper=1;
791 2622 : linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
792 2622 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, natomsper, linked_atoms );
793 205705 : for(unsigned j=0; j<natomsper; ++j) {
794 353118 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) {
795 150035 : taskFlags[k]=1;
796 : }
797 : }
798 : }
799 210 : } else if( nactive_atoms>0 ) {
800 : // Get some parallel info
801 54 : unsigned stride=comm.Get_size();
802 54 : unsigned rank=comm.Get_rank();
803 54 : if( serialCalculation() ) {
804 : stride=1;
805 : rank=0;
806 : }
807 :
808 : unsigned nactive_three=0;
809 66599 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
810 66545 : if( isCurrentlyActive( ablocks[2][i] ) ) {
811 66545 : nactive_three++;
812 : }
813 : }
814 :
815 54 : std::vector<Vector> lttmp_pos( nactive_three );
816 54 : std::vector<unsigned> lttmp_ind( nactive_three );
817 :
818 : nactive_three=0;
819 54 : if( allthirdblockintasks ) {
820 66599 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
821 66545 : if( !isCurrentlyActive( ablocks[2][i] ) ) {
822 0 : continue;
823 : }
824 66545 : lttmp_ind[nactive_three]=ablocks[2][i];
825 66545 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
826 66545 : nactive_three++;
827 : }
828 : } else {
829 0 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
830 0 : if( !isCurrentlyActive( ablocks[2][i] ) ) {
831 0 : continue;
832 : }
833 0 : lttmp_ind[nactive_three]=i;
834 0 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
835 0 : nactive_three++;
836 : }
837 : }
838 : // Build the list of the link cells
839 54 : threecells.buildCellLists( lttmp_pos, lttmp_ind, getPbc() );
840 :
841 : // Ensure we only do tasks where atoms are in appropriate link cells
842 54 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
843 54 : std::vector<unsigned> tlinked_atoms( 1+ablocks[2].size() );
844 633 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
845 579 : if( !isCurrentlyActive( ablocks[0][i] ) ) {
846 0 : continue;
847 : }
848 579 : unsigned natomsper=1;
849 579 : linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
850 579 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, natomsper, linked_atoms );
851 579 : if( allthirdblockintasks ) {
852 71832 : for(unsigned j=0; j<natomsper; ++j) {
853 142372 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) {
854 71119 : taskFlags[k]=1;
855 : }
856 : }
857 : } else {
858 0 : unsigned ntatomsper=1;
859 0 : tlinked_atoms[0]=lttmp_ind[0];
860 0 : threecells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, ntatomsper, tlinked_atoms );
861 0 : for(unsigned j=0; j<natomsper; ++j) {
862 0 : for(unsigned k=0; k<ntatomsper; ++k) {
863 0 : taskFlags[bookeeping(i,linked_atoms[j]).first+tlinked_atoms[k]]=1;
864 : }
865 : }
866 : }
867 : }
868 : }
869 210 : if( !serialCalculation() ) {
870 210 : comm.Sum( taskFlags );
871 : }
872 210 : lockContributors();
873 : }
874 :
875 245988 : void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const {
876 : plumed_dbg_assert( !usespecies && nblock>0 );
877 245988 : if( atoms.size()!=decoder.size() ) {
878 24491 : atoms.resize( decoder.size() );
879 : }
880 :
881 245988 : unsigned scode = taskCode;
882 786464 : for(unsigned i=0; i<decoder.size(); ++i) {
883 540476 : unsigned ind=( scode / decoder[i] );
884 540476 : atoms[i] = ablocks[i][ind];
885 540476 : scode -= ind*decoder[i];
886 : }
887 245988 : }
888 :
889 452889 : bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ) const {
890 452889 : if( isDensity() ) {
891 : myatoms.setNumberOfAtoms( 1 );
892 19070 : myatoms.setAtom( 0, taskCode );
893 19070 : return true;
894 433819 : } else if( usespecies ) {
895 182233 : std::vector<unsigned> task_atoms(1);
896 182233 : task_atoms[0]=taskCode;
897 182233 : unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getPositionOfAtomForLinkCells( taskCode ), linkcells );
898 182233 : return natomsper>1;
899 251586 : } else if( matsums ) {
900 : myatoms.setNumberOfAtoms( getNumberOfAtoms() );
901 52798 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
902 52357 : myatoms.setAtom( i, i );
903 : }
904 251145 : } else if( allthirdblockintasks ) {
905 : plumed_dbg_assert( ablocks.size()==3 );
906 39816 : std::vector<unsigned> atoms(2);
907 39816 : decodeIndexToAtoms( taskCode, atoms );
908 39816 : myatoms.setupAtomsFromLinkCells( atoms, getPositionOfAtomForLinkCells( atoms[0] ), threecells );
909 211329 : } else if( nblock>0 ) {
910 179551 : std::vector<unsigned> atoms( ablocks.size() );
911 179551 : decodeIndexToAtoms( taskCode, atoms );
912 179551 : myatoms.setNumberOfAtoms( ablocks.size() );
913 587153 : for(unsigned i=0; i<ablocks.size(); ++i) {
914 407602 : myatoms.setAtom( i, atoms[i] );
915 : }
916 : } else {
917 31778 : myatoms.setNumberOfAtoms( ablocks.size() );
918 124498 : for(unsigned i=0; i<ablocks.size(); ++i) {
919 92720 : myatoms.setAtom( i, ablocks[i][taskCode] );
920 : }
921 : }
922 : return true;
923 : }
924 :
925 9074 : void MultiColvarBase::setupActiveTaskSet( std::vector<unsigned>& active_tasks, const std::string& input_label ) {
926 9074 : if( !setup_completed ) {
927 : bool justVolumes=false;
928 1931 : if( usespecies ) {
929 : justVolumes=true;
930 1438 : for(unsigned i=0; i<getNumberOfVessels(); ++i) {
931 1318 : vesselbase::StoreDataVessel* mys=dynamic_cast<vesselbase::StoreDataVessel*>( getPntrToVessel(i) );
932 1318 : if( mys ) {
933 508 : continue;
934 : }
935 810 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(i) );
936 810 : if( !myb ) {
937 : justVolumes=false;
938 : break;
939 : }
940 38 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
941 38 : if( !myv ) {
942 : justVolumes=false;
943 : break;
944 : }
945 : }
946 : }
947 1931 : deactivateAllTasks();
948 1931 : if( justVolumes && mydata ) {
949 88 : if( mydata->getNumberOfDataUsers()==0 ) {
950 : justVolumes=false;
951 : }
952 :
953 137 : for(unsigned i=0; i<mydata->getNumberOfDataUsers(); ++i) {
954 49 : MultiColvarBase* myu=dynamic_cast<MultiColvarBase*>( mydata->getDataUser(i) );
955 49 : if( myu ) {
956 49 : myu->setupActiveTaskSet( taskFlags, getLabel() );
957 : } else {
958 0 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
959 0 : taskFlags[i]=1;
960 : }
961 : }
962 : }
963 : }
964 1931 : if( justVolumes ) {
965 160 : for(unsigned j=0; j<getNumberOfVessels(); ++j) {
966 80 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(j) );
967 80 : if( !myb ) {
968 48 : continue ;
969 : }
970 32 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
971 32 : if( !myv ) {
972 0 : continue ;
973 : }
974 32 : myv->retrieveAtoms();
975 32 : myv->setupRegions();
976 :
977 148161 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
978 148129 : if( myv->inVolumeOfInterest(i) ) {
979 3438 : taskFlags[i]=1;
980 : }
981 : }
982 : }
983 : } else {
984 4886702 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
985 4884851 : taskFlags[i]=1;
986 : }
987 : }
988 :
989 : // Now activate all this class
990 1931 : lockContributors();
991 : // Setup the link cells
992 1931 : setupLinkCells();
993 : // Ensures that setup is not performed multiple times during one cycle
994 1931 : setup_completed=true;
995 : }
996 :
997 : // And activate the tasks in input action
998 9074 : if( getLabel()!=input_label ) {
999 : int input_code=-1;
1000 61 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
1001 61 : if( mybasemulticolvars[i]->getLabel()==input_label ) {
1002 49 : input_code=i+1;
1003 49 : break;
1004 : }
1005 : }
1006 :
1007 49 : MultiValue my_tvals( getNumberOfQuantities(), getNumberOfDerivatives() );
1008 49 : AtomValuePack mytmp_atoms( my_tvals, this );
1009 148234 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
1010 148185 : if( !taskIsCurrentlyActive(i) ) {
1011 144656 : continue;
1012 : }
1013 3529 : setupCurrentAtomList( getTaskCode(i), mytmp_atoms );
1014 254796 : for(unsigned j=0; j<mytmp_atoms.getNumberOfAtoms(); ++j) {
1015 : unsigned itask=mytmp_atoms.getIndex(j);
1016 251267 : if( atom_lab[itask].first==input_code ) {
1017 250567 : active_tasks[ atom_lab[itask].second ]=1;
1018 : }
1019 : }
1020 : }
1021 49 : }
1022 9074 : }
1023 :
1024 467 : bool MultiColvarBase::filtersUsedAsInput() {
1025 : bool inputAreFilters=false;
1026 728 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
1027 261 : MultiColvarFilter* myfilt=dynamic_cast<MultiColvarFilter*>( mybasemulticolvars[i] );
1028 261 : if( myfilt || mybasemulticolvars[i]->filtersUsedAsInput() ) {
1029 : inputAreFilters=true;
1030 : }
1031 : }
1032 467 : return inputAreFilters;
1033 : }
1034 :
1035 9025 : void MultiColvarBase::calculate() {
1036 : // Recursive function that sets up tasks
1037 9025 : setupActiveTaskSet( taskFlags, getLabel() );
1038 :
1039 : // Check for filters and rerun setup of link cells if there are any
1040 9025 : if( mybasemulticolvars.size()>0 && filtersUsedAsInput() ) {
1041 8 : setupLinkCells();
1042 : }
1043 :
1044 : // Setup the link cells if we are not using species
1045 9025 : if( !usespecies && ablocks.size()>1 ) {
1046 : // This loop finds the first active atom, which is always checked because
1047 : // of a peculiarity in linkcells
1048 4628 : unsigned first_active=std::numeric_limits<unsigned>::max();
1049 4762 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
1050 4762 : if( !isCurrentlyActive( ablocks[1][i] ) ) {
1051 : continue;
1052 : } else {
1053 4628 : first_active=i;
1054 4628 : break;
1055 : }
1056 : }
1057 4628 : setupNonUseSpeciesLinkCells( first_active );
1058 : }
1059 : // And run all tasks
1060 9025 : runAllTasks();
1061 9025 : }
1062 :
1063 213 : void MultiColvarBase::calculateNumericalDerivatives( ActionWithValue* a ) {
1064 213 : if( mybasemulticolvars.size()>0 ) {
1065 0 : plumed_merror("cannot calculate numerical derivatives for this quantity");
1066 : }
1067 213 : calculateAtomicNumericalDerivatives( this, 0 );
1068 213 : }
1069 :
1070 2960 : void MultiColvarBase::prepare() {
1071 2960 : setup_completed=false;
1072 2960 : atomsWereRetrieved=false;
1073 2960 : }
1074 :
1075 6550 : void MultiColvarBase::retrieveAtoms() {
1076 6550 : if( !atomsWereRetrieved ) {
1077 2960 : ActionAtomistic::retrieveAtoms();
1078 2960 : atomsWereRetrieved=true;
1079 : }
1080 6550 : }
1081 :
1082 38245 : void MultiColvarBase::mergeInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
1083 : const unsigned& jatom, const std::vector<double>& der,
1084 : MultiValue& myder, AtomValuePack& myatoms ) const {
1085 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
1086 : plumed_dbg_assert( ival<myatoms.getUnderlyingMultiValue().getNumberOfValues() );
1087 : plumed_dbg_assert( start<myder.getNumberOfValues() && end<=myder.getNumberOfValues() );
1088 : plumed_dbg_assert( der.size()==myder.getNumberOfValues() && jatom<myatoms.getNumberOfAtoms() );
1089 : // Convert input atom to local index
1090 : unsigned katom = myatoms.getIndex( jatom );
1091 : plumed_dbg_assert( katom<atom_lab.size() );
1092 : plumed_dbg_assert( atom_lab[katom].first>0 );
1093 : // Find base colvar
1094 38245 : unsigned mmc=atom_lab[katom].first - 1;
1095 : plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
1096 : // Get start of indices for this atom
1097 : unsigned basen=0;
1098 38855 : for(unsigned i=0; i<mmc; ++i) {
1099 610 : basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
1100 : }
1101 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
1102 38245 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
1103 4805851 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
1104 : unsigned jder=myder.getActiveIndex(j);
1105 4767606 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
1106 4434561 : unsigned kder=basen+jder;
1107 110197260 : for(unsigned icomp=start; icomp<end; ++icomp) {
1108 105762699 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
1109 : }
1110 : } else {
1111 333045 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
1112 7551360 : for(unsigned icomp=start; icomp<end; ++icomp) {
1113 7218315 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
1114 : }
1115 : }
1116 : }
1117 38245 : }
1118 :
1119 106 : void MultiColvarBase::splitInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
1120 : const unsigned& jatom, const std::vector<double>& der,
1121 : MultiValue& myder, AtomValuePack& myatoms ) const {
1122 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
1123 : plumed_dbg_assert( ival<myder.getNumberOfValues() );
1124 : plumed_dbg_assert( start<myvals.getNumberOfValues() && end<=myvals.getNumberOfValues() );
1125 : plumed_dbg_assert( der.size()==myatoms.getUnderlyingMultiValue().getNumberOfValues() && jatom<myatoms.getNumberOfAtoms() );
1126 : // Convert input atom to local index
1127 : unsigned katom = myatoms.getIndex( jatom );
1128 : plumed_dbg_assert( katom<atom_lab.size() );
1129 : plumed_dbg_assert( atom_lab[katom].first>0 );
1130 : // Find base colvar
1131 106 : unsigned mmc=atom_lab[katom].first - 1;
1132 : plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
1133 : // Get start of indices for this atom
1134 : unsigned basen=0;
1135 154 : for(unsigned i=0; i<mmc; ++i) {
1136 48 : basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
1137 : }
1138 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
1139 106 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
1140 19534 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
1141 : unsigned jder=myder.getActiveIndex(j);
1142 19428 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
1143 18474 : unsigned kder=basen+jder;
1144 18474 : plumed_assert( kder<myvals.getNumberOfDerivatives() );
1145 39942 : for(unsigned icomp=start; icomp<end; ++icomp) {
1146 21468 : myvals.addDerivative( icomp, kder, der[icomp]*myder.getDerivative( ival, jder ) );
1147 : }
1148 : } else {
1149 954 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
1150 3942 : for(unsigned icomp=start; icomp<end; ++icomp) {
1151 2988 : myvals.addDerivative( icomp, kder, der[icomp]*myder.getDerivative( ival, jder ) );
1152 : }
1153 : }
1154 : }
1155 106 : }
1156 :
1157 980606 : void MultiColvarBase::addComDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
1158 : plumed_dbg_assert( ival<static_cast<int>(myatoms.getUnderlyingMultiValue().getNumberOfValues()) && iatom<myatoms.getNumberOfAtoms() );
1159 : // Convert input atom to local index
1160 : unsigned katom = myatoms.getIndex( iatom );
1161 : plumed_dbg_assert( atom_lab[katom].first>0 );
1162 : // Find base colvar
1163 980606 : unsigned mmc = atom_lab[katom].first - 1;
1164 : plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
1165 980606 : if( usespecies && iatom==0 ) {
1166 453923 : myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[0] );
1167 453923 : return;
1168 : }
1169 :
1170 : // Get start of indices for this atom
1171 526683 : unsigned basen=0;
1172 721370 : for(unsigned i=0; i<mmc; ++i) {
1173 194687 : basen+=(mybasemulticolvars[i]->getNumberOfDerivatives() - 9) / 3;
1174 : }
1175 526683 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[katom].second, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
1176 526683 : myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
1177 : }
1178 :
1179 1124044 : void MultiColvarBase::getInputData( const unsigned& ind, const bool& normed,
1180 : const multicolvar::AtomValuePack& myatoms,
1181 : std::vector<double>& orient ) const {
1182 : // Converint input atom to local index
1183 : unsigned katom = myatoms.getIndex(ind);
1184 : plumed_dbg_assert( atom_lab[katom].first>0 );
1185 : // Find base colvar
1186 1124044 : unsigned mmc = atom_lab[katom].first - 1;
1187 : plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
1188 : // Check if orient is the correct size
1189 1124044 : if( orient.size()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ) {
1190 0 : orient.resize( mybasemulticolvars[mmc]->getNumberOfQuantities() );
1191 : }
1192 : // Retrieve the value
1193 1124044 : mybasedata[mmc]->retrieveValueWithIndex( atom_lab[katom].second, normed, orient );
1194 1124044 : }
1195 :
1196 54449 : MultiValue& MultiColvarBase::getInputDerivatives( const unsigned& iatom, const bool& normed, const multicolvar::AtomValuePack& myatoms ) const {
1197 : // Converint input atom to local index
1198 : unsigned katom = myatoms.getIndex(iatom);
1199 : plumed_dbg_assert( atom_lab[katom].first>0 );
1200 : // Find base colvar
1201 54449 : unsigned mmc = atom_lab[katom].first - 1;
1202 : plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
1203 54449 : if( usespecies && !normed && iatom==0 ) {
1204 1798 : return mybasedata[mmc]->getTemporyMultiValue(0);
1205 : }
1206 :
1207 52651 : unsigned oval=0;
1208 52651 : if( iatom>0 ) {
1209 36741 : oval=1;
1210 : }
1211 52651 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(oval);
1212 105263 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
1213 52612 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
1214 39 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
1215 : }
1216 52651 : mybasedata[mmc]->retrieveDerivatives( atom_lab[katom].second, normed, myder );
1217 : return myder;
1218 : }
1219 :
1220 65138769 : void MultiColvarBase::accumulateSymmetryFunction( const int& ival, const unsigned& iatom, const double& val, const Vector& der, const Tensor& vir, multicolvar::AtomValuePack& myatoms ) const {
1221 : plumed_dbg_assert( usespecies );
1222 : unsigned katom=myatoms.getIndex(0), jatom=myatoms.getIndex(iatom);
1223 : double weight0=1.0;
1224 65138769 : if( atom_lab[katom].first>0 ) {
1225 794 : weight0=mybasedata[atom_lab[katom].first-1]->retrieveWeightWithIndex( atom_lab[katom].second );
1226 : }
1227 : double weighti=1.0;
1228 65138769 : if( atom_lab[jatom].first>0 ) {
1229 794 : weighti=mybasedata[atom_lab[jatom].first-1]->retrieveWeightWithIndex( atom_lab[jatom].second );
1230 : }
1231 : // Accumulate the value
1232 65138769 : if( ival<0 ) {
1233 2766628 : myatoms.getUnderlyingMultiValue().addTemporyValue( weight0*weighti*val );
1234 : } else {
1235 62372141 : myatoms.addValue( ival, weight0*weighti*val );
1236 : }
1237 :
1238 : // Return if we don't need derivatives
1239 65138769 : if( doNotCalculateDerivatives() ) {
1240 : return ;
1241 : }
1242 : // And virial
1243 58765246 : if( ival<0 ) {
1244 2287875 : myatoms.addTemporyBoxDerivatives( weight0*weighti*vir );
1245 : } else {
1246 56477371 : myatoms.addBoxDerivatives( ival, weight0*weighti*vir );
1247 : }
1248 :
1249 : // Add derivatives of central atom
1250 58765246 : if( atom_lab[katom].first>0 ) {
1251 794 : addComDerivatives( ival, 0, -weight0*weighti*der, myatoms );
1252 794 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
1253 794 : tmpder[0]=weighti*val;
1254 794 : mergeInputDerivatives( ival, 0, 1, 0, tmpder, getInputDerivatives(0, false, myatoms), myatoms );
1255 : } else {
1256 58764452 : if( ival<0 ) {
1257 2287875 : myatoms.addTemporyAtomsDerivatives( 0, -der );
1258 : } else {
1259 56476577 : myatoms.addAtomsDerivatives( ival, 0, -der );
1260 : }
1261 : }
1262 : // Add derivatives of atom in coordination sphere
1263 58765246 : if( atom_lab[jatom].first>0 ) {
1264 794 : addComDerivatives( ival, iatom, weight0*weighti*der, myatoms );
1265 794 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
1266 794 : tmpder[0]=weight0*val;
1267 794 : mergeInputDerivatives( ival, 0, 1, iatom, tmpder, getInputDerivatives(iatom, false, myatoms), myatoms );
1268 : } else {
1269 58764452 : if( ival<0 ) {
1270 2287875 : myatoms.addTemporyAtomsDerivatives( iatom, der );
1271 : } else {
1272 56476577 : myatoms.addAtomsDerivatives( ival, iatom, der );
1273 : }
1274 : }
1275 : }
1276 :
1277 1503315 : void MultiColvarBase::addAtomDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
1278 1503315 : if( doNotCalculateDerivatives() ) {
1279 : return ;
1280 : }
1281 : unsigned jatom=myatoms.getIndex(iatom);
1282 :
1283 1255739 : if( atom_lab[jatom].first>0 ) {
1284 979018 : addComDerivatives( ival, iatom, der, myatoms );
1285 : } else {
1286 276721 : if( ival<0 ) {
1287 0 : myatoms.addTemporyAtomsDerivatives( iatom, der );
1288 : } else {
1289 276721 : myatoms.addAtomsDerivatives( ival, iatom, der );
1290 : }
1291 : }
1292 : }
1293 :
1294 244246 : double MultiColvarBase::calculateWeight( const unsigned& current, const double& weight, AtomValuePack& myvals ) const {
1295 244246 : return 1.0;
1296 : }
1297 :
1298 449360 : void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
1299 449360 : AtomValuePack myatoms( myvals, this );
1300 : // Retrieve the atom list
1301 449360 : if( !setupCurrentAtomList( current, myatoms ) ) {
1302 : return;
1303 : }
1304 : // Get weight due to dynamic groups
1305 449240 : double weight = 1.0;
1306 449240 : if( !matsums ) {
1307 386210280 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
1308 385940036 : if( atom_lab[myatoms.getIndex(i)].first==0 ) {
1309 385696019 : continue;
1310 : }
1311 : // Only need to do first two atoms for things like TopologyMatrix, HbondMatrix, Bridge and so on
1312 244017 : if( allthirdblockintasks && i>1 ) {
1313 : break;
1314 : }
1315 244017 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
1316 244017 : weight *= mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
1317 : }
1318 178996 : } else if( usespecies ) {
1319 178555 : if( atom_lab[myatoms.getIndex(0)].first>0 ) {
1320 9073 : if( mybasedata[atom_lab[myatoms.getIndex(0)].first-1]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(0)].second )<epsilon ) {
1321 78 : weight=0.;
1322 : }
1323 : }
1324 : }
1325 : // Do a quick check on the size of this contribution
1326 449240 : double multweight = calculateWeight( current, weight, myatoms );
1327 449240 : if( weight*multweight<getTolerance() ) {
1328 121536 : updateActiveAtoms( myatoms );
1329 : return;
1330 : }
1331 : myatoms.setValue( 0, weight*multweight );
1332 : // Deal with derivatives of weights due to dynamic groups
1333 327704 : if( !matsums && !doNotCalculateDerivatives() && mybasemulticolvars.size()>0 ) {
1334 : MultiValue& outder=myatoms.getUnderlyingMultiValue();
1335 39416 : MultiValue myder(0,0);
1336 115123 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
1337 : // Neglect any atoms without differentiable weights
1338 75707 : if( atom_lab[myatoms.getIndex(i)].first==0 ) {
1339 0 : continue;
1340 : }
1341 :
1342 : // Retrieve derivatives
1343 75707 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
1344 75707 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() || myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
1345 39416 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
1346 : }
1347 75707 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(i)].second, false, myder );
1348 :
1349 : // Retrieve the prefactor (product of all other weights)
1350 75707 : double prefactor = multweight*weight / mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
1351 : // And accumulate the derivatives
1352 2927141 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
1353 2851434 : unsigned jder=myder.getActiveIndex(j);
1354 2851434 : outder.addDerivative( 0, jder, prefactor*myder.getDerivative(0,jder) );
1355 : }
1356 75707 : myder.clearAll();
1357 : }
1358 39416 : }
1359 : // Retrieve derivative stuff for central atom
1360 327704 : if( !doNotCalculateDerivatives() ) {
1361 208301 : if( usespecies && mybasemulticolvars.size()>0 && atom_lab[myatoms.getIndex(0)].first>0 ) {
1362 1610 : unsigned mmc = atom_lab[0].first - 1;
1363 1610 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(0);
1364 3208 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
1365 1598 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
1366 12 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
1367 : }
1368 1610 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(0)].second, false, myder );
1369 1610 : unsigned basen=0;
1370 1610 : for(unsigned i=0; i<mmc; ++i) {
1371 0 : basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
1372 : }
1373 1610 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[myatoms.getIndex(0)].second, mybasemulticolvars[mmc]->my_tmp_capacks[0] );
1374 : }
1375 : }
1376 : // Compute everything
1377 327704 : double vv=compute( task_index, myatoms );
1378 327704 : updateActiveAtoms( myatoms );
1379 : myatoms.setValue( 1, vv );
1380 327704 : return;
1381 : }
1382 :
1383 662394 : void MultiColvarBase::updateActiveAtoms( AtomValuePack& myatoms ) const {
1384 662394 : if( mybasemulticolvars.size()==0 ) {
1385 520523 : myatoms.updateUsingIndices();
1386 : } else {
1387 141871 : myatoms.updateDynamicList();
1388 : }
1389 662394 : }
1390 :
1391 2558087 : Vector MultiColvarBase::getCentralAtomPos( const unsigned& taskIndex ) {
1392 2558087 : unsigned curr=getTaskCode( taskIndex );
1393 :
1394 2558087 : if( usespecies || isDensity() ) {
1395 1264098 : return getPositionOfAtomForLinkCells(curr);
1396 1293989 : } else if( nblock>0 ) {
1397 : // double factor=1.0/static_cast<double>( ablocks.size() );
1398 2130 : Vector mypos;
1399 2130 : mypos.zero();
1400 2130 : std::vector<unsigned> atoms( ablocks.size() );
1401 2130 : decodeIndexToAtoms( curr, atoms );
1402 6390 : for(unsigned i=0; i<ablocks.size(); ++i) {
1403 4260 : if( use_for_central_atom[i] ) {
1404 4260 : mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(atoms[i]);
1405 : }
1406 : }
1407 2130 : return mypos;
1408 : } else {
1409 1291859 : Vector mypos;
1410 1291859 : mypos.zero();
1411 5138646 : for(unsigned i=0; i<ablocks.size(); ++i) {
1412 3846787 : if( use_for_central_atom[i] ) {
1413 1323451 : mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(ablocks[i][curr]);
1414 : }
1415 : }
1416 1291859 : return mypos;
1417 : }
1418 : }
1419 :
1420 642832 : void MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsigned& taskIndex, CatomPack& mypack ) {
1421 642832 : unsigned curr=getTaskCode( taskIndex );
1422 :
1423 642832 : if(usespecies) {
1424 502418 : if( mypack.getNumberOfAtomsWithDerivatives()!=1 ) {
1425 25571 : mypack.resize(1);
1426 : }
1427 502418 : mypack.setIndex( 0, basn + curr );
1428 502418 : mypack.setDerivative( 0, Tensor::identity() );
1429 140414 : } else if( nblock>0 ) {
1430 0 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) {
1431 0 : mypack.resize(ncentral);
1432 : }
1433 0 : unsigned k=0;
1434 0 : std::vector<unsigned> atoms( ablocks.size() );
1435 0 : decodeIndexToAtoms( curr, atoms );
1436 0 : for(unsigned i=0; i<ablocks.size(); ++i) {
1437 0 : if( use_for_central_atom[i] ) {
1438 0 : mypack.setIndex( k, basn + atoms[i] );
1439 0 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
1440 0 : k++;
1441 : }
1442 : }
1443 : } else {
1444 140414 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) {
1445 89006 : mypack.resize(ncentral);
1446 : }
1447 140414 : unsigned k=0;
1448 316748 : for(unsigned i=0; i<ablocks.size(); ++i) {
1449 176334 : if( use_for_central_atom[i] ) {
1450 173694 : mypack.setIndex( k, basn + ablocks[i][curr] );
1451 173694 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
1452 173694 : k++;
1453 : }
1454 : }
1455 : }
1456 642832 : }
1457 :
1458 121150993 : Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 ) const {
1459 121150993 : if(usepbc) {
1460 121150993 : return pbcDistance( vec1, vec2 );
1461 : } else {
1462 0 : return delta( vec1, vec2 );
1463 : }
1464 : }
1465 :
1466 222049 : void MultiColvarBase::applyPbc(std::vector<Vector>& dlist, unsigned int max_index) const {
1467 222049 : if (usepbc) {
1468 222049 : pbcApply(dlist, max_index);
1469 : }
1470 222049 : }
1471 :
1472 1995 : void MultiColvarBase::apply() {
1473 1995 : if( getForcesFromVessels( forcesToApply ) ) {
1474 86 : setForcesOnAtoms( forcesToApply );
1475 : }
1476 1995 : }
1477 :
1478 : }
1479 : }
|
