Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 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 : using namespace std;
36 :
37 : namespace PLMD {
38 : namespace multicolvar {
39 :
40 399 : void MultiColvarBase::registerKeywords( Keywords& keys ) {
41 399 : Action::registerKeywords( keys );
42 399 : ActionWithValue::registerKeywords( keys );
43 399 : ActionAtomistic::registerKeywords( keys );
44 1197 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
45 399 : ActionWithVessel::registerKeywords( keys );
46 1596 : keys.add("hidden","NL_STRIDE","the frequency with which the neighbor list should be updated. Between neighbour list update steps all quantities "
47 : "that contributed less than TOL at the previous neighbor list update step are ignored.");
48 798 : 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 "
49 : "regular collective variables. The label is used to reference the full set of quantities calculated by "
50 : "the action. This is usual when using \\ref multicolvarfunction. Generally when doing this the previously calculated "
51 : "multicolvar will be referenced using the DATA keyword rather than ARG.\n\n"
52 : "This Action can be used to calculate the following scalar quantities directly. These quantities are calculated by "
53 : "employing the keywords listed below. "
54 : "These quantities can then be referenced elsewhere in the input file by using this Action's label "
55 : "followed by a dot and the name of the quantity. Some amongst them can be calculated multiple times "
56 : "with different parameters. In this case the quantities calculated can be referenced elsewhere in the "
57 : "input by using the name of the quantity followed by a numerical identifier "
58 : "e.g. <em>label</em>.lessthan-1, <em>label</em>.lessthan-2 etc. When doing this and, for clarity we have "
59 : "made the label of the components customizable. As such by using the LABEL keyword in the description of the keyword "
60 : "input you can customize the component name");
61 1596 : keys.reserve("atoms-3","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
62 : "one coordination number for each of the atoms specified. Each of these coordination numbers specifies how many of the "
63 : "other specified atoms are within a certain cutoff of the central atom. You can specify the atoms here as another multicolvar "
64 : "action or using a MultiColvarFilter or ActionVolume action. When you do so the quantity is calculated for those atoms specified "
65 : "in the previous multicolvar. This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
66 : "coordination number more than four for example");
67 1596 : 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 "
68 : "one coordination number for each of the atoms specified in SPECIESA. Each of these cooordination numbers specifies how many "
69 : "of the atoms specifies using SPECIESB is within the specified cutoff. As with the species keyword the input can also be specified "
70 : "using the label of another multicolvar");
71 1596 : 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 "
72 : "the documentation for that keyword");
73 399 : }
74 :
75 328 : MultiColvarBase::MultiColvarBase(const ActionOptions&ao):
76 : Action(ao),
77 : ActionAtomistic(ao),
78 : ActionWithValue(ao),
79 : ActionWithVessel(ao),
80 : usepbc(false),
81 : allthirdblockintasks(false),
82 : uselinkforthree(false),
83 : linkcells(comm),
84 : threecells(comm),
85 : setup_completed(false),
86 : atomsWereRetrieved(false),
87 : matsums(false),
88 : usespecies(false),
89 984 : nblock(0)
90 : {
91 656 : if( keywords.exists("NOPBC") ) {
92 656 : bool nopbc=!usepbc; parseFlag("NOPBC",nopbc);
93 328 : usepbc=!nopbc;
94 : }
95 656 : if( keywords.exists("SPECIESA") ) { matsums=usespecies=true; }
96 328 : }
97 :
98 46 : void MultiColvarBase::readAtomsLikeKeyword( const std::string & key, const int& natoms, std::vector<AtomNumber>& all_atoms ) {
99 46 : plumed_assert( !usespecies );
100 46 : if( all_atoms.size()>0 ) return;
101 :
102 : std::vector<AtomNumber> t;
103 916 : for(int i=1;; ++i ) {
104 962 : parseAtomList(key, i, t );
105 962 : if( t.empty() ) break;
106 :
107 916 : log.printf(" Colvar %d is calculated from atoms : ", i);
108 9344 : for(unsigned j=0; j<t.size(); ++j) log.printf("%d ",t[j].serial() );
109 916 : log.printf("\n");
110 :
111 922 : if( i==1 && natoms<0 ) { ablocks.resize(t.size()); }
112 910 : else if( i==1 ) ablocks.resize(natoms);
113 916 : if( t.size()!=ablocks.size() ) {
114 0 : std::string ss; Tools::convert(i,ss);
115 0 : error(key + ss + " keyword has the wrong number of atoms");
116 : }
117 9344 : for(unsigned j=0; j<ablocks.size(); ++j) {
118 7512 : ablocks[j].push_back( ablocks.size()*(i-1)+j ); all_atoms.push_back( t[j] );
119 7512 : atom_lab.push_back( std::pair<unsigned,unsigned>( 0, ablocks.size()*(i-1)+j ) );
120 : }
121 916 : t.resize(0);
122 916 : }
123 46 : if( all_atoms.size()>0 ) {
124 46 : nblock=0;
125 1008 : for(unsigned i=0; i<ablocks[0].size(); ++i) addTaskToList( i );
126 : }
127 : }
128 :
129 440 : bool MultiColvarBase::parseMultiColvarAtomList(const std::string& key, const int& num, std::vector<AtomNumber>& t) {
130 440 : std::vector<std::string> mlabs;
131 440 : if( num<0 ) parseVector(key,mlabs);
132 0 : else parseNumberedVector(key,num,mlabs);
133 :
134 440 : if( mlabs.size()==0 ) return false;
135 :
136 296 : std::string mname; unsigned found_mcolv=mlabs.size();
137 958 : for(unsigned i=0; i<mlabs.size(); ++i) {
138 600 : MultiColvarBase* mycolv = plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlabs[i]);
139 300 : if(!mycolv) { found_mcolv=i; break; }
140 : // Check all base multicolvars are of same type
141 122 : if( i==0 ) {
142 118 : mname = mycolv->getName();
143 118 : if( mycolv->isPeriodic() ) error("multicolvar functions don't work with this multicolvar");
144 : } else {
145 8 : if( mname!=mycolv->getName() ) error("All input multicolvars must be of same type");
146 : }
147 : // And track which variable stores each colvar
148 17801200 : for(unsigned j=0; j<mycolv->getFullNumberOfTasks(); ++j) atom_lab.push_back( std::pair<unsigned,unsigned>( mybasemulticolvars.size()+1, j ) );
149 : // And store the multicolvar base
150 122 : mybasemulticolvars.push_back( mycolv );
151 : // And create a basedata stash
152 488 : mybasedata.push_back( mybasemulticolvars[mybasemulticolvars.size()-1]->buildDataStashes( this ) );
153 : // Check if weight has derivatives
154 244 : if( mybasemulticolvars[ mybasemulticolvars.size()-1 ]->weightHasDerivatives ) weightHasDerivatives=true;
155 122 : plumed_assert( mybasemulticolvars.size()==mybasedata.size() );
156 : }
157 : // Have we conventional atoms to read in
158 296 : if( found_mcolv==0 ) {
159 178 : std::vector<AtomNumber> tt; ActionAtomistic::interpretAtomList( mlabs, tt );
160 344772 : for(unsigned i=0; i<tt.size(); ++i) { atom_lab.push_back( std::pair<unsigned,unsigned>( 0, t.size() + i ) ); }
161 356 : log.printf(" keyword %s takes atoms : ", key.c_str() );
162 344772 : for(unsigned i=0; i<tt.size(); ++i) { t.push_back( tt[i] ); log.printf("%d ",tt[i].serial() ); }
163 178 : log.printf("\n");
164 118 : } else if( found_mcolv==mlabs.size() ) {
165 118 : if( checkNumericalDerivatives() ) error("cannot use NUMERICAL_DERIVATIVES keyword with dynamic groups as input");
166 236 : log.printf(" keyword %s takes dynamic groups of atoms constructed from multicolvars labelled : ", key.c_str() );
167 602 : for(unsigned i=0; i<mlabs.size(); ++i) log.printf("%s ",mlabs[i].c_str() );
168 118 : log.printf("\n");
169 0 : } else if( found_mcolv<mlabs.size() ) {
170 0 : error("cannot mix multicolvar input and atom input in one line");
171 : }
172 : return true;
173 : }
174 :
175 76 : void MultiColvarBase::readTwoGroups( const std::string& key0, const std::string& key1, const std::string& key2, std::vector<AtomNumber>& all_atoms ) {
176 76 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) ); ablocks.resize( 2 );
177 :
178 76 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
179 78 : nblock=atom_lab.size(); for(unsigned i=0; i<2; ++i) ablocks[i].resize(nblock);
180 26 : resizeBookeepingArray( nblock, nblock );
181 16406 : for(unsigned i=0; i<nblock; ++i) ablocks[0][i]=ablocks[1][i]=i;
182 10894 : for(unsigned i=1; i<nblock; ++i) {
183 8427224 : for(unsigned j=0; j<i; ++j) {
184 4210895 : bookeeping(i,j).first=getFullNumberOfTasks();
185 4210895 : addTaskToList( i*nblock + j );
186 4210895 : bookeeping(i,j).second=getFullNumberOfTasks();
187 : }
188 : }
189 : } else {
190 50 : parseMultiColvarAtomList(key1,-1,all_atoms);
191 134 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=i;
192 50 : parseMultiColvarAtomList(key2,-1,all_atoms);
193 3307 : ablocks[1].resize( atom_lab.size() - ablocks[0].size() ); for(unsigned i=0; i<ablocks[1].size(); ++i) ablocks[1][i]=ablocks[0].size() + i;
194 :
195 50 : if( ablocks[0].size()>ablocks[1].size() ) nblock = ablocks[0].size();
196 50 : else nblock=ablocks[1].size();
197 :
198 100 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
199 151 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
200 3361 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
201 1109 : bookeeping(i,j).first=getFullNumberOfTasks();
202 2218 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 ) {
203 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
204 0 : atom_lab[ablocks[0][i]].second!=atom_lab[ablocks[1][j]].second ) addTaskToList( i*nblock + j );
205 5545 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] ) addTaskToList( i*nblock + j );
206 1109 : bookeeping(i,j).second=getFullNumberOfTasks();
207 : }
208 : }
209 : }
210 76 : }
211 :
212 10 : void MultiColvarBase::readGroupKeywords( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& key3,
213 : const bool& no_third_dim_accum, const bool& symmetric, std::vector<AtomNumber>& all_atoms ) {
214 10 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) ); ablocks.resize( 3 );
215 :
216 10 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
217 8 : if( no_third_dim_accum ) {
218 24 : nblock=atom_lab.size(); ablocks[0].resize(nblock); ablocks[1].resize( nblock );
219 5632 : for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=ablocks[1][i]=i;
220 8 : resizeBookeepingArray( nblock, nblock );
221 8 : if( symmetric ) {
222 : // This ensures that later parts of the code don't switch off allthirdblockintasks
223 4017 : for(unsigned i=0; i<nblock; ++i) { bookeeping(i,i).first=0; bookeeping(i,i).second=1; }
224 2668 : for(unsigned i=1; i<nblock; ++i) {
225 443563 : for(unsigned j=0; j<i; ++j) {
226 221116 : bookeeping(j,i).first=bookeeping(i,j).first=getFullNumberOfTasks();
227 221116 : addTaskToList( i*nblock + j );
228 442232 : bookeeping(j,i).second=bookeeping(i,j).second=getFullNumberOfTasks();
229 : }
230 : }
231 : } else {
232 136 : for(unsigned i=0; i<nblock; ++i) {
233 8277 : for(unsigned j=0; j<nblock; ++j) {
234 4105 : if( i==j ) continue ;
235 4038 : bookeeping(i,j).first=getFullNumberOfTasks();
236 4038 : addTaskToList( i*nblock + j );
237 4038 : bookeeping(i,j).second=getFullNumberOfTasks();
238 : }
239 : }
240 : }
241 8 : if( !parseMultiColvarAtomList(key3,-1,all_atoms) ) error("missing required keyword " + key3 + " in input");
242 8 : ablocks[2].resize( atom_lab.size() - ablocks[0].size() );
243 198820 : for(unsigned i=0; i<ablocks[2].size(); ++i) ablocks[2][i]=ablocks[0].size() + i;
244 : } else {
245 0 : nblock=atom_lab.size(); for(unsigned i=0; i<3; ++i) ablocks[i].resize(nblock);
246 0 : resizeBookeepingArray( nblock, nblock );
247 0 : for(unsigned i=0; i<nblock; ++i) { ablocks[0][i]=i; ablocks[1][i]=i; ablocks[2][i]=i; }
248 0 : if( symmetric ) {
249 0 : for(unsigned i=2; i<nblock; ++i) {
250 0 : for(unsigned j=1; j<i; ++j) {
251 0 : bookeeping(i,j).first=getFullNumberOfTasks();
252 0 : for(unsigned k=0; k<j; ++k) addTaskToList( i*nblock*nblock + j*nblock + k );
253 0 : bookeeping(i,j).second=getFullNumberOfTasks();
254 : }
255 : }
256 : } else {
257 0 : for(unsigned i=0; i<nblock; ++i) {
258 0 : for(unsigned j=0; j<nblock; ++j) {
259 0 : if( i==j ) continue;
260 0 : bookeeping(i,j).first=getFullNumberOfTasks();
261 0 : for(unsigned k=0; k<nblock; ++k) {
262 0 : if( i!=k && j!=k ) addTaskToList( i*nblock*nblock + j*nblock + k );
263 : }
264 0 : bookeeping(i,j).first=getFullNumberOfTasks();
265 : }
266 : }
267 : }
268 : }
269 : } else {
270 2 : readThreeGroups( key1, key2, key3, true, no_third_dim_accum, all_atoms );
271 : }
272 :
273 10 : }
274 :
275 11 : void MultiColvarBase::readThreeGroups( const std::string& key1, const std::string& key2, const std::string& key3,
276 : const bool& allow2, const bool& no_third_dim_accum, std::vector<AtomNumber>& all_atoms ) {
277 11 : plumed_dbg_assert( keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) ); ablocks.resize( 3 );
278 :
279 11 : bool readkey1=parseMultiColvarAtomList(key1,-1,all_atoms);
280 84 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=i;
281 11 : bool readkey2=parseMultiColvarAtomList(key2,-1,all_atoms);
282 11 : if( !readkey1 && !readkey2 ) return ;
283 697 : ablocks[1].resize( atom_lab.size() - ablocks[0].size() ); for(unsigned i=0; i<ablocks[1].size(); ++i) ablocks[1][i]=ablocks[0].size() + i;
284 :
285 22 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
286 11 : bool readkey3=parseMultiColvarAtomList(key3,-1,all_atoms);
287 11 : if( !readkey3 && !allow2 ) {
288 0 : error("missing atom specification " + key3);
289 11 : } else if( !readkey3 ) {
290 2 : if( ablocks[1].size()>ablocks[0].size() ) nblock=ablocks[1].size();
291 0 : else nblock=ablocks[0].size();
292 :
293 2 : ablocks[2].resize( ablocks[1].size() );
294 796 : for(unsigned i=0; i<ablocks[1].size(); ++i) ablocks[2][i]=ablocks[0].size() + i;
295 10 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
296 592 : for(unsigned j=1; j<ablocks[1].size(); ++j) {
297 196 : bookeeping(i,j).first=getFullNumberOfTasks();
298 19600 : for(unsigned k=0; k<j; ++k) {
299 19404 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
300 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
301 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
302 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
303 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 &&
304 0 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) addTaskToList( nblock*nblock*i + nblock*j + k );
305 38808 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
306 48510 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
307 9702 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) addTaskToList( nblock*nblock*i + nblock*j + k );
308 : }
309 196 : bookeeping(i,j).second=getFullNumberOfTasks();
310 : }
311 : }
312 : } else {
313 18 : ablocks[2].resize( atom_lab.size() - ablocks[1].size() - ablocks[0].size() );
314 15938 : for(unsigned i=0; i<ablocks[2].size(); ++i) ablocks[2][i] = ablocks[0].size() + ablocks[1].size() + i;
315 :
316 9 : if( ablocks[1].size()>ablocks[0].size() ) nblock=ablocks[1].size();
317 9 : else nblock=ablocks[0].size();
318 9 : if( ablocks[2].size()>nblock ) nblock=ablocks[2].size();
319 :
320 9 : unsigned kcount; if( no_third_dim_accum ) kcount=1; else kcount=ablocks[2].size();
321 :
322 105 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
323 685 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
324 209 : bookeeping(i,j).first=getFullNumberOfTasks();
325 627 : for(unsigned k=0; k<kcount; ++k) {
326 209 : if( no_third_dim_accum ) addTaskToList( nblock*i + j );
327 0 : else if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
328 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
329 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
330 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
331 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 &&
332 0 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) addTaskToList( nblock*nblock*i + nblock*j + k );
333 0 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
334 0 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
335 0 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) addTaskToList( nblock*nblock*i + nblock*j + k );
336 : }
337 209 : bookeeping(i,j).second=getFullNumberOfTasks();
338 : }
339 : }
340 : }
341 : }
342 :
343 2 : void MultiColvarBase::buildSets() {
344 : std::vector<AtomNumber> fake_atoms;
345 4 : if( !parseMultiColvarAtomList("DATA",-1,fake_atoms) ) error("missing DATA keyword");
346 2 : if( fake_atoms.size()>0 ) error("no atoms should appear in the specification for this object. Input should be other multicolvars");
347 :
348 4 : nblock = mybasemulticolvars[0]->getFullNumberOfTasks();
349 16 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
350 8 : if( mybasemulticolvars[i]->getFullNumberOfTasks()!=nblock ) {
351 0 : error("mismatch between numbers of tasks in various base multicolvars");
352 : }
353 : }
354 2 : ablocks.resize( mybasemulticolvars.size() ); usespecies=false;
355 16 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
356 4 : ablocks[i].resize( nblock );
357 28 : for(unsigned j=0; j<nblock; ++j) ablocks[i][j]=i*nblock+j;
358 : }
359 8 : for(unsigned i=0; i<nblock; ++i) {
360 6 : if( mybasemulticolvars.size()<4 ) {
361 6 : unsigned cvcode=0, tmpc=1;
362 36 : for(unsigned j=0; j<ablocks.size(); ++j) { cvcode += i*tmpc; tmpc *= nblock; }
363 6 : addTaskToList( cvcode );
364 : } else {
365 0 : addTaskToList( i );
366 : }
367 : }
368 4 : mybasedata[0]->resizeTemporyMultiValues( mybasemulticolvars.size() ); setupMultiColvarBase( fake_atoms );
369 2 : }
370 :
371 12506830 : void MultiColvarBase::addTaskToList( const unsigned& taskCode ) {
372 12506830 : plumed_assert( getNumberOfVessels()==0 );
373 12506830 : ActionWithVessel::addTaskToList( taskCode );
374 12506830 : }
375 :
376 95 : void MultiColvarBase::resizeBookeepingArray( const unsigned& num1, const unsigned& num2 ) {
377 95 : bookeeping.resize( num1, num2 );
378 13919 : for(unsigned i=0; i<num1; ++i) {
379 26636232 : for(unsigned j=0; j<num2; ++j) { bookeeping(i,j).first=0; bookeeping(i,j).second=0; }
380 : }
381 95 : }
382 :
383 284 : void MultiColvarBase::setupMultiColvarBase( const std::vector<AtomNumber>& atoms ) {
384 402 : if( !matsums && atom_lab.size()==0 ) error("No atoms have been read in");
385 : std::vector<AtomNumber> all_atoms;
386 : // Setup decoder array
387 284 : if( !usespecies && nblock>0 ) {
388 :
389 60 : ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
390 60 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
391 60 : if( ablocks.size()==3 ) {
392 19 : allthirdblockintasks=uselinkforthree=true;
393 2889 : for(unsigned i=0; i<bookeeping.nrows(); ++i) {
394 897597 : for(unsigned j=0; j<bookeeping.ncols(); ++j) {
395 448081 : unsigned ntper = bookeeping(i,j).second - bookeeping(i,j).first;
396 448081 : if( i==j && ntper==0 ) {
397 69 : continue;
398 448012 : } else if( ntper == 1 && allthirdblockintasks ) {
399 447818 : allthirdblockintasks=true;
400 388 : } else if( ntper != ablocks[2].size() ) {
401 194 : allthirdblockintasks=uselinkforthree=false;
402 : } else {
403 0 : allthirdblockintasks=false;
404 : }
405 : }
406 : }
407 : }
408 :
409 60 : if( allthirdblockintasks ) {
410 34 : decoder.resize(2); plumed_assert( ablocks.size()==3 );
411 : // Check if number of atoms is too large
412 17 : if( pow( double(nblock), 2.0 )>std::numeric_limits<unsigned>::max() ) error("number of atoms in groups is too big for PLUMED to handle");
413 : } else {
414 43 : decoder.resize( ablocks.size() );
415 : // Check if number of atoms is too large
416 43 : if( pow( double(nblock), double(ablocks.size()) )>std::numeric_limits<unsigned>::max() ) error("number of atoms in groups is too big for PLUMED to handle");
417 : }
418 486 : unsigned code=1; for(unsigned i=0; i<decoder.size(); ++i) { decoder[decoder.size()-1-i]=code; code *= nblock; }
419 224 : } else if( !usespecies ) {
420 83 : ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
421 83 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
422 282 : } else if( keywords.exists("SPECIESA") ) {
423 174 : plumed_assert( atom_lab.size()==0 && all_atoms.size()==0 );
424 261 : ablocks.resize( 1 ); bool readspecies=parseMultiColvarAtomList("SPECIES", -1, all_atoms);
425 87 : if( readspecies ) {
426 79880 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<atom_lab.size(); ++i) { addTaskToList(i); ablocks[0][i]=i; }
427 : } else {
428 32 : if( !parseMultiColvarAtomList("SPECIESA", -1, all_atoms) ) error("missing SPECIES/SPECIESA keyword");
429 16 : unsigned nat1=atom_lab.size();
430 32 : if( !parseMultiColvarAtomList("SPECIESB", -1, all_atoms) ) error("missing SPECIESB keyword");
431 16 : unsigned nat2=atom_lab.size() - nat1;
432 :
433 16 : for(unsigned i=0; i<nat1; ++i) addTaskToList(i);
434 16 : ablocks[0].resize( nat2 );
435 4412 : for(unsigned i=0; i<nat2; ++i) {
436 : bool found=false; unsigned inum;
437 521066 : for(unsigned j=0; j<nat1; ++j) {
438 771748 : if( atom_lab[nat1+i].first>0 && atom_lab[j].first>0 ) {
439 252720 : if( atom_lab[nat1+i].first==atom_lab[j].first &&
440 252720 : mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
441 252720 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) { found=true; inum=j; break; }
442 6794 : } else if( atom_lab[nat1+i].first>0 ) {
443 0 : if( mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
444 0 : all_atoms[atom_lab[j].second] ) { found=true; inum=nat1 + i; break; }
445 13588 : } else if( atom_lab[j].first>0 ) {
446 0 : if( all_atoms[atom_lab[nat1+i].second]==
447 0 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) { found=true; inum=nat1+i; break; }
448 20382 : } else if( all_atoms[atom_lab[nat1+i].second]==all_atoms[atom_lab[j].second] ) { found=true; inum=j; break; }
449 : }
450 : // This prevents mistakes being made in colvar setup
451 160 : if( found ) { ablocks[0][i]=inum; }
452 4236 : else { ablocks[0][i]=nat1 + i; }
453 : }
454 : }
455 : }
456 284 : if( mybasemulticolvars.size()>0 ) {
457 596 : for(unsigned i=0; i<mybasedata.size(); ++i) {
458 244 : mybasedata[i]->resizeTemporyMultiValues(2); mybasemulticolvars[i]->my_tmp_capacks.resize(2);
459 : }
460 : }
461 :
462 : // Copy lists of atoms involved from base multicolvars
463 : std::vector<AtomNumber> tmp_atoms;
464 934 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
465 122 : tmp_atoms=mybasemulticolvars[i]->getAbsoluteIndexes();
466 1243981 : for(unsigned j=0; j<tmp_atoms.size(); ++j) all_atoms.push_back( tmp_atoms[j] );
467 : }
468 : // Copy atom lists from input
469 284 : if( !usespecies ) {
470 176296 : for(unsigned i=0; i<atoms.size(); ++i) all_atoms.push_back( atoms[i] );
471 : }
472 :
473 : // Now make sure we get all the atom positions
474 284 : ActionAtomistic::requestAtoms( all_atoms );
475 : // And setup dependencies
476 934 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) addDependency( mybasemulticolvars[i] );
477 :
478 : // Setup underlying ActionWithVessel
479 284 : readVesselKeywords();
480 284 : }
481 :
482 18 : void MultiColvarBase::setAtomsForCentralAtom( const std::vector<bool>& catom_ind ) {
483 18 : unsigned nat=0; plumed_assert( catom_ind.size()==ablocks.size() );
484 231 : for(unsigned i=0; i<catom_ind.size(); ++i) {
485 : use_for_central_atom[i]=catom_ind[i];
486 65 : if( use_for_central_atom[i] ) nat++;
487 : }
488 18 : plumed_dbg_assert( nat>0 ); ncentral=nat;
489 18 : numberForCentralAtom = 1.0 / static_cast<double>( nat );
490 18 : }
491 :
492 399 : void MultiColvarBase::turnOnDerivatives() {
493 399 : ActionWithValue::turnOnDerivatives();
494 399 : needsDerivatives();
495 399 : forcesToApply.resize( getNumberOfDerivatives() );
496 399 : }
497 :
498 175 : void MultiColvarBase::setLinkCellCutoff( const double& lcut, double tcut ) {
499 263 : plumed_assert( usespecies || ablocks.size()<4 );
500 344 : if( tcut<0 ) tcut=lcut;
501 175 : linkcells.setCutoff( lcut );
502 175 : threecells.setCutoff( tcut );
503 175 : }
504 :
505 8 : double MultiColvarBase::getLinkCellCutoff() const {
506 8 : return linkcells.getCutoff();
507 : }
508 :
509 1688 : void MultiColvarBase::setupLinkCells() {
510 2599 : if( (!usespecies && nblock==0) || !linkcells.enabled() ) return ;
511 : // Retrieve any atoms that haven't already been retrieved
512 1093 : for(std::vector<MultiColvarBase*>::iterator p=mybasemulticolvars.begin(); p!=mybasemulticolvars.end(); ++p) {
513 202 : (*p)->retrieveAtoms();
514 : }
515 891 : retrieveAtoms();
516 :
517 : unsigned iblock;
518 891 : if( usespecies ) {
519 : iblock=0;
520 220 : } else if( ablocks.size()<4 ) {
521 : iblock=1;
522 : } else {
523 0 : plumed_error();
524 : }
525 :
526 : // Count number of currently active atoms
527 891 : nactive_atoms=0;
528 1134057 : for(unsigned i=0; i<ablocks[iblock].size(); ++i) {
529 754850 : if( isCurrentlyActive( ablocks[iblock][i] ) ) nactive_atoms++;
530 : }
531 :
532 891 : if( nactive_atoms>0 ) {
533 891 : std::vector<Vector> ltmp_pos( nactive_atoms );
534 891 : std::vector<unsigned> ltmp_ind( nactive_atoms );
535 :
536 891 : nactive_atoms=0;
537 891 : if( usespecies ) {
538 1070374 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
539 712688 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
540 712496 : ltmp_ind[nactive_atoms]=ablocks[0][i];
541 712496 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ltmp_ind[nactive_atoms] );
542 356248 : nactive_atoms++;
543 : }
544 : } else {
545 63683 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
546 42162 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
547 32322 : ltmp_ind[nactive_atoms]=i;
548 48483 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ablocks[1][i] );
549 16161 : nactive_atoms++;
550 : }
551 : }
552 :
553 : // Build the lists for the link cells
554 1782 : linkcells.buildCellLists( ltmp_pos, ltmp_ind, getPbc() );
555 : }
556 : }
557 :
558 5108 : void MultiColvarBase::setupNonUseSpeciesLinkCells( const unsigned& my_always_active ) {
559 5108 : plumed_assert( !usespecies );
560 6201 : if( nblock==0 || !linkcells.enabled() ) return ;
561 708 : deactivateAllTasks();
562 : std::vector<unsigned> requiredlinkcells;
563 :
564 708 : if( !uselinkforthree && nactive_atoms>0 ) {
565 : // Get some parallel info
566 156 : unsigned stride=comm.Get_size();
567 156 : unsigned rank=comm.Get_rank();
568 156 : if( serialCalculation() ) { stride=1; rank=0; }
569 :
570 : // Ensure we only do tasks where atoms are in appropriate link cells
571 156 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
572 17115 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
573 11202 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
574 5244 : unsigned natomsper=1; linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
575 5244 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, natomsper, linked_atoms );
576 408788 : for(unsigned j=0; j<natomsper; ++j) {
577 909319 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) taskFlags[k]=1;
578 : }
579 156 : }
580 552 : } else if( nactive_atoms>0 ) {
581 : // Get some parallel info
582 552 : unsigned stride=comm.Get_size();
583 552 : unsigned rank=comm.Get_rank();
584 552 : if( serialCalculation() ) { stride=1; rank=0; }
585 :
586 : unsigned nactive_three=0;
587 215337 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
588 142822 : if( isCurrentlyActive( ablocks[2][i] ) ) nactive_three++;
589 : }
590 :
591 552 : std::vector<Vector> lttmp_pos( nactive_three );
592 552 : std::vector<unsigned> lttmp_ind( nactive_three );
593 :
594 : nactive_three=0;
595 552 : if( allthirdblockintasks ) {
596 215337 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
597 142822 : if( !isCurrentlyActive( ablocks[2][i] ) ) continue;
598 142822 : lttmp_ind[nactive_three]=ablocks[2][i];
599 142822 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
600 71411 : nactive_three++;
601 : }
602 : } else {
603 0 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
604 0 : if( !isCurrentlyActive( ablocks[2][i] ) ) continue;
605 0 : lttmp_ind[nactive_three]=i;
606 0 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
607 0 : nactive_three++;
608 : }
609 : }
610 : // Build the list of the link cells
611 1104 : threecells.buildCellLists( lttmp_pos, lttmp_ind, getPbc() );
612 :
613 : // Ensure we only do tasks where atoms are in appropriate link cells
614 552 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
615 552 : std::vector<unsigned> tlinked_atoms( 1+ablocks[2].size() );
616 17640 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
617 11024 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
618 11024 : unsigned natomsper=1; linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
619 11024 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, natomsper, linked_atoms );
620 5512 : if( allthirdblockintasks ) {
621 239808 : for(unsigned j=0; j<natomsper; ++j) {
622 585606 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) taskFlags[k]=1;
623 : }
624 : } else {
625 0 : unsigned ntatomsper=1; tlinked_atoms[0]=lttmp_ind[0];
626 0 : threecells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, ntatomsper, tlinked_atoms );
627 0 : for(unsigned j=0; j<natomsper; ++j) {
628 0 : for(unsigned k=0; k<ntatomsper; ++k) taskFlags[bookeeping(i,linked_atoms[j]).first+tlinked_atoms[k]]=1;
629 : }
630 : }
631 : }
632 : }
633 708 : if( !serialCalculation() ) comm.Sum( taskFlags );
634 708 : lockContributors();
635 : }
636 :
637 291922 : void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const {
638 : plumed_dbg_assert( !usespecies && nblock>0 );
639 291922 : if( atoms.size()!=decoder.size() ) atoms.resize( decoder.size() );
640 :
641 291922 : unsigned scode = taskCode;
642 2480906 : for(unsigned i=0; i<decoder.size(); ++i) {
643 632354 : unsigned ind=( scode / decoder[i] );
644 1264708 : atoms[i] = ablocks[i][ind];
645 632354 : scode -= ind*decoder[i];
646 : }
647 291922 : }
648 :
649 465455 : bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ) const {
650 465455 : if( isDensity() ) {
651 19065 : myatoms.setNumberOfAtoms( 1 ); myatoms.setAtom( 0, taskCode ); return true;
652 446390 : } else if( usespecies ) {
653 298068 : std::vector<unsigned> task_atoms(1); task_atoms[0]=taskCode;
654 149034 : unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getPositionOfAtomForLinkCells( taskCode ), linkcells );
655 149034 : return natomsper>1;
656 297356 : } else if( matsums ) {
657 : myatoms.setNumberOfAtoms( getNumberOfAtoms() );
658 28407 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) myatoms.setAtom( i, i );
659 297043 : } else if( allthirdblockintasks ) {
660 85778 : plumed_dbg_assert( ablocks.size()==3 ); std::vector<unsigned> atoms(2); decodeIndexToAtoms( taskCode, atoms );
661 171556 : myatoms.setupAtomsFromLinkCells( atoms, getPositionOfAtomForLinkCells( atoms[0] ), threecells );
662 211265 : } else if( nblock>0 ) {
663 179523 : std::vector<unsigned> atoms( ablocks.size() );
664 359046 : decodeIndexToAtoms( taskCode, atoms ); myatoms.setNumberOfAtoms( ablocks.size() );
665 1174158 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, atoms[i] );
666 : } else {
667 31742 : myatoms.setNumberOfAtoms( ablocks.size() );
668 248636 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, ablocks[i][taskCode] );
669 : }
670 : return true;
671 : }
672 :
673 9021 : void MultiColvarBase::setupActiveTaskSet( std::vector<unsigned>& active_tasks, const std::string& input_label ) {
674 9021 : if( !setup_completed ) {
675 : bool justVolumes=false;
676 1680 : if( usespecies ) {
677 : justVolumes=true;
678 1328 : for(unsigned i=0; i<getNumberOfVessels(); ++i) {
679 670 : vesselbase::StoreDataVessel* mys=dynamic_cast<vesselbase::StoreDataVessel*>( getPntrToVessel(i) );
680 670 : if( mys ) continue;
681 373 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(i) );
682 373 : if( !myb ) { justVolumes=false; break; }
683 38 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
684 38 : if( !myv ) { justVolumes=false; break; }
685 : }
686 : }
687 1680 : deactivateAllTasks();
688 1680 : if( justVolumes && mydata ) {
689 297 : if( mydata->getNumberOfDataUsers()==0 ) justVolumes=false;
690 :
691 741 : for(unsigned i=0; i<mydata->getNumberOfDataUsers(); ++i) {
692 49 : MultiColvarBase* myu=dynamic_cast<MultiColvarBase*>( mydata->getDataUser(i) );
693 49 : if( myu ) {
694 98 : myu->setupActiveTaskSet( taskFlags, getLabel() );
695 : } else {
696 0 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
697 : }
698 : }
699 : }
700 1680 : if( justVolumes ) {
701 240 : for(unsigned j=0; j<getNumberOfVessels(); ++j) {
702 80 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(j) );
703 80 : if( !myb ) continue ;
704 32 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
705 32 : if( !myv ) continue ;
706 32 : myv->retrieveAtoms(); myv->setupRegions();
707 :
708 296322 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
709 151567 : if( myv->inVolumeOfInterest(i) ) taskFlags[i]=1;
710 : }
711 : }
712 : } else {
713 9757270 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
714 : }
715 :
716 : // Now activate all this class
717 1680 : lockContributors();
718 : // Setup the link cells
719 1680 : setupLinkCells();
720 : // Ensures that setup is not performed multiple times during one cycle
721 1680 : setup_completed=true;
722 : }
723 :
724 : // And activate the tasks in input action
725 18042 : if( getLabel()!=input_label ) {
726 : int input_code=-1;
727 134 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
728 122 : if( mybasemulticolvars[i]->getLabel()==input_label ) { input_code=i+1; break; }
729 : }
730 :
731 98 : MultiValue my_tvals( getNumberOfQuantities(), getNumberOfDerivatives() );
732 49 : AtomValuePack mytmp_atoms( my_tvals, this );
733 296468 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
734 148185 : if( !taskIsCurrentlyActive(i) ) continue;
735 3529 : setupCurrentAtomList( getTaskCode(i), mytmp_atoms );
736 506063 : for(unsigned j=0; j<mytmp_atoms.getNumberOfAtoms(); ++j) {
737 : unsigned itask=mytmp_atoms.getIndex(j);
738 753101 : if( atom_lab[itask].first==input_code ) active_tasks[ atom_lab[itask].second ]=1;
739 : }
740 : }
741 : }
742 9021 : }
743 :
744 433 : bool MultiColvarBase::filtersUsedAsInput() {
745 : bool inputAreFilters=false;
746 1589 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
747 241 : MultiColvarFilter* myfilt=dynamic_cast<MultiColvarFilter*>( mybasemulticolvars[i] );
748 241 : if( myfilt || mybasemulticolvars[i]->filtersUsedAsInput() ) inputAreFilters=true;
749 : }
750 433 : return inputAreFilters;
751 : }
752 :
753 8972 : void MultiColvarBase::calculate() {
754 : // Recursive function that sets up tasks
755 17944 : setupActiveTaskSet( taskFlags, getLabel() );
756 :
757 : // Check for filters and rerun setup of link cells if there are any
758 8972 : if( mybasemulticolvars.size()>0 && filtersUsedAsInput() ) setupLinkCells();
759 :
760 : // Setup the link cells if we are not using species
761 16089 : if( !usespecies && ablocks.size()>1 ) {
762 : // This loop finds the first active atom, which is always checked because
763 : // of a peculiarity in linkcells
764 5108 : unsigned first_active=std::numeric_limits<unsigned>::max();
765 10618 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
766 10484 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
767 : else {
768 5108 : first_active=i; break;
769 : }
770 : }
771 5108 : setupNonUseSpeciesLinkCells( first_active );
772 : }
773 : // And run all tasks
774 8972 : runAllTasks();
775 8972 : }
776 :
777 217 : void MultiColvarBase::calculateNumericalDerivatives( ActionWithValue* a ) {
778 217 : if( mybasemulticolvars.size()>0 ) plumed_merror("cannot calculate numerical derivatives for this quantity");
779 217 : calculateAtomicNumericalDerivatives( this, 0 );
780 217 : }
781 :
782 2709 : void MultiColvarBase::prepare() {
783 2709 : setup_completed=false; atomsWereRetrieved=false;
784 2709 : }
785 :
786 6073 : void MultiColvarBase::retrieveAtoms() {
787 6073 : if( !atomsWereRetrieved ) { ActionAtomistic::retrieveAtoms(); atomsWereRetrieved=true; }
788 6073 : }
789 :
790 38205 : void MultiColvarBase::mergeInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
791 : const unsigned& jatom, const std::vector<double>& der,
792 : MultiValue& myder, AtomValuePack& myatoms ) const {
793 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
794 : plumed_dbg_assert( ival<myatoms.getUnderlyingMultiValue().getNumberOfValues() );
795 : plumed_dbg_assert( start<myder.getNumberOfValues() && end<=myder.getNumberOfValues() );
796 : plumed_dbg_assert( der.size()==myder.getNumberOfValues() && jatom<myatoms.getNumberOfAtoms() );
797 : // Convert input atom to local index
798 : unsigned katom = myatoms.getIndex( jatom ); plumed_dbg_assert( katom<atom_lab.size() ); plumed_dbg_assert( atom_lab[katom].first>0 );
799 : // Find base colvar
800 76410 : unsigned mmc=atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
801 : // Get start of indices for this atom
802 39425 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
803 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
804 :
805 38205 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
806 9573417 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
807 : unsigned jder=myder.getActiveIndex(j);
808 9535212 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
809 4434561 : unsigned kder=basen+jder;
810 110197260 : for(unsigned icomp=start; icomp<end; ++icomp) {
811 317288097 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
812 : }
813 : } else {
814 333045 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
815 7551360 : for(unsigned icomp=start; icomp<end; ++icomp) {
816 14436630 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
817 : }
818 : }
819 : }
820 38205 : }
821 :
822 60 : void MultiColvarBase::splitInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
823 : const unsigned& jatom, const std::vector<double>& der,
824 : MultiValue& myder, AtomValuePack& myatoms ) const {
825 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
826 : plumed_dbg_assert( ival<myder.getNumberOfValues() );
827 : plumed_dbg_assert( start<myvals.getNumberOfValues() && end<=myvals.getNumberOfValues() );
828 : plumed_dbg_assert( der.size()==myatoms.getUnderlyingMultiValue().getNumberOfValues() && jatom<myatoms.getNumberOfAtoms() );
829 : // Convert input atom to local index
830 : unsigned katom = myatoms.getIndex( jatom ); plumed_dbg_assert( katom<atom_lab.size() ); plumed_dbg_assert( atom_lab[katom].first>0 );
831 : // Find base colvar
832 120 : unsigned mmc=atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
833 : // Get start of indices for this atom
834 120 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
835 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
836 60 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
837 36420 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
838 : unsigned jder=myder.getActiveIndex(j);
839 36360 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
840 35280 : unsigned kder=basen+jder; plumed_assert( kder<myvals.getNumberOfDerivatives() );
841 35280 : for(unsigned icomp=start; icomp<end; ++icomp) {
842 52920 : myvals.addDerivative( icomp, kder, der[icomp]*myder.getDerivative( ival, jder ) );
843 : }
844 : } else {
845 540 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
846 1080 : for(unsigned icomp=start; icomp<end; ++icomp) {
847 1620 : myvals.addDerivative( icomp, kder, der[icomp]*myder.getDerivative( ival, jder ) );
848 : }
849 : }
850 : }
851 60 : }
852 :
853 905722 : void MultiColvarBase::addComDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
854 : plumed_dbg_assert( ival<static_cast<int>(myatoms.getUnderlyingMultiValue().getNumberOfValues()) && iatom<myatoms.getNumberOfAtoms() );
855 : // Convert input atom to local index
856 : unsigned katom = myatoms.getIndex( iatom ); plumed_dbg_assert( atom_lab[katom].first>0 );
857 : // Find base colvar
858 1811444 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
859 1738684 : if( usespecies && iatom==0 ) { myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[0] ); return; }
860 :
861 : // Get start of indices for this atom
862 646506 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=(mybasemulticolvars[i]->getNumberOfDerivatives() - 9) / 3;
863 1467723 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[katom].second, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
864 978482 : myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
865 : }
866 :
867 1122612 : void MultiColvarBase::getInputData( const unsigned& ind, const bool& normed,
868 : const multicolvar::AtomValuePack& myatoms,
869 : std::vector<double>& orient ) const {
870 : // Converint input atom to local index
871 : unsigned katom = myatoms.getIndex(ind); plumed_dbg_assert( atom_lab[katom].first>0 );
872 : // Find base colvar
873 2245224 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
874 : // Check if orient is the correct size
875 2245224 : if( orient.size()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ) orient.resize( mybasemulticolvars[mmc]->getNumberOfQuantities() );
876 : // Retrieve the value
877 2245224 : mybasedata[mmc]->retrieveValueWithIndex( atom_lab[katom].second, normed, orient );
878 1122612 : }
879 :
880 52977 : MultiValue& MultiColvarBase::getInputDerivatives( const unsigned& iatom, const bool& normed, const multicolvar::AtomValuePack& myatoms ) const {
881 : // Converint input atom to local index
882 : unsigned katom = myatoms.getIndex(iatom); plumed_dbg_assert( atom_lab[katom].first>0 );
883 : // Find base colvar
884 105954 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
885 54755 : if( usespecies && !normed && iatom==0 ) return mybasedata[mmc]->getTemporyMultiValue(0);
886 :
887 51199 : unsigned oval=0; if( iatom>0 ) oval=1;
888 102398 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(oval);
889 102362 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
890 51163 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
891 72 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
892 : }
893 102398 : mybasedata[mmc]->retrieveDerivatives( atom_lab[katom].second, normed, myder );
894 : return myder;
895 : }
896 :
897 60207197 : void MultiColvarBase::accumulateSymmetryFunction( const int& ival, const unsigned& iatom, const double& val, const Vector& der, const Tensor& vir, multicolvar::AtomValuePack& myatoms ) const {
898 : plumed_dbg_assert( usespecies ); unsigned katom=myatoms.getIndex(0), jatom=myatoms.getIndex(iatom);
899 120415168 : double weight0=1.0; if( atom_lab[katom].first>0 ) weight0=mybasedata[atom_lab[katom].first-1]->retrieveWeightWithIndex( atom_lab[katom].second );
900 120415168 : double weighti=1.0; if( atom_lab[jatom].first>0 ) weighti=mybasedata[atom_lab[jatom].first-1]->retrieveWeightWithIndex( atom_lab[jatom].second );
901 : // Accumulate the value
902 62943725 : if( ival<0 ) myatoms.getUnderlyingMultiValue().addTemporyValue( weight0*weighti*val );
903 57470669 : else myatoms.addValue( ival, weight0*weighti*val );
904 :
905 : // Return if we don't need derivatives
906 60207197 : if( doNotCalculateDerivatives() ) return ;
907 : // And virial
908 54044650 : if( ival<0 ) myatoms.addTemporyBoxDerivatives( weight0*weighti*vir );
909 51786875 : else myatoms.addBoxDerivatives( ival, weight0*weighti*vir );
910 :
911 : // Add derivatives of central atom
912 54044650 : if( atom_lab[katom].first>0 ) {
913 774 : addComDerivatives( ival, 0, -weight0*weighti*der, myatoms );
914 2322 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
915 1548 : tmpder[0]=weighti*val; mergeInputDerivatives( ival, 0, 1, 0, tmpder, getInputDerivatives(0, false, myatoms), myatoms );
916 : } else {
917 54043876 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( 0, -der );
918 51786101 : else myatoms.addAtomsDerivatives( ival, 0, -der );
919 : }
920 : // Add derivatives of atom in coordination sphere
921 54044650 : if( atom_lab[jatom].first>0 ) {
922 774 : addComDerivatives( ival, iatom, weight0*weighti*der, myatoms );
923 2322 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
924 1548 : tmpder[0]=weight0*val; mergeInputDerivatives( ival, 0, 1, iatom, tmpder, getInputDerivatives(iatom, false, myatoms), myatoms );
925 : } else {
926 54043876 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
927 51786101 : else myatoms.addAtomsDerivatives( ival, iatom, der );
928 : }
929 : }
930 :
931 2927751 : void MultiColvarBase::addAtomDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
932 2927751 : if( doNotCalculateDerivatives() ) return ;
933 : unsigned jatom=myatoms.getIndex(iatom);
934 :
935 5360350 : if( atom_lab[jatom].first>0 ) {
936 904174 : addComDerivatives( ival, iatom, der, myatoms );
937 : } else {
938 1776001 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
939 1776001 : else myatoms.addAtomsDerivatives( ival, iatom, der );
940 : }
941 : }
942 :
943 256812 : double MultiColvarBase::calculateWeight( const unsigned& current, const double& weight, AtomValuePack& myvals ) const {
944 256812 : return 1.0;
945 : }
946 :
947 461926 : void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
948 461926 : AtomValuePack myatoms( myvals, this );
949 : // Retrieve the atom list
950 461926 : if( !setupCurrentAtomList( current, myatoms ) ) return;
951 : // Get weight due to dynamic groups
952 461806 : double weight = 1.0;
953 461806 : if( !matsums ) {
954 772330581 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
955 772014578 : if( atom_lab[myatoms.getIndex(i)].first==0 ) continue;
956 : // Only need to do first two atoms for thigns like TopologyMatrix, HbondMatrix, Bridge and so on
957 243971 : if( allthirdblockintasks && i>1 ) break;
958 243971 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
959 487942 : weight *= mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
960 : }
961 145803 : } else if( usespecies ) {
962 290980 : if( atom_lab[myatoms.getIndex(0)].first>0 ) {
963 18136 : if( mybasedata[atom_lab[myatoms.getIndex(0)].first-1]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(0)].second )<epsilon ) weight=0.;
964 : }
965 : }
966 : // Do a quick check on the size of this contribution
967 461806 : double multweight = calculateWeight( current, weight, myatoms );
968 923612 : if( weight*multweight<getTolerance() ) {
969 121536 : updateActiveAtoms( myatoms );
970 : return;
971 : }
972 : myatoms.setValue( 0, weight*multweight );
973 : // Deal with derivatives of weights due to dynamic groups
974 480460 : if( !matsums && !doNotCalculateDerivatives() && mybasemulticolvars.size()>0 ) {
975 78776 : MultiValue& outder=myatoms.getUnderlyingMultiValue(); MultiValue myder(0,0);
976 190710 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
977 : // Neglect any atoms without differentiable weights
978 151322 : if( atom_lab[myatoms.getIndex(i)].first==0 ) continue;
979 :
980 : // Retrieve derivatives
981 75661 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
982 187595 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() || myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
983 78776 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
984 : }
985 226983 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(i)].second, false, myder );
986 :
987 : // Retrieve the prefactor (product of all other weights)
988 302644 : double prefactor = multweight*weight / mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
989 : // And accumulate the derivatives
990 8626219 : for(unsigned j=0; j<myder.getNumberActive(); ++j) { unsigned jder=myder.getActiveIndex(j); outder.addDerivative( 0, jder, prefactor*myder.getDerivative(0,jder) ); }
991 75661 : myder.clearAll();
992 : }
993 : }
994 : // Retrieve derivative stuff for central atom
995 340270 : if( !doNotCalculateDerivatives() ) {
996 318525 : if( usespecies && mybasemulticolvars.size()>0 && atom_lab[myatoms.getIndex(0)].first>0 ) {
997 1605 : unsigned mmc = atom_lab[0].first - 1;
998 3210 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(0);
999 3198 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
1000 1593 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
1001 24 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
1002 : }
1003 4815 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(0)].second, false, myder );
1004 1605 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
1005 4815 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[myatoms.getIndex(0)].second, mybasemulticolvars[mmc]->my_tmp_capacks[0] );
1006 : }
1007 : }
1008 : // Compute everything
1009 340270 : double vv=compute( task_index, myatoms ); updateActiveAtoms( myatoms );
1010 : myatoms.setValue( 1, vv );
1011 340270 : return;
1012 : }
1013 :
1014 672158 : void MultiColvarBase::updateActiveAtoms( AtomValuePack& myatoms ) const {
1015 672158 : if( mybasemulticolvars.size()==0 ) myatoms.updateUsingIndices();
1016 141710 : else myatoms.updateDynamicList();
1017 672158 : }
1018 :
1019 2753562 : Vector MultiColvarBase::getCentralAtomPos( const unsigned& taskIndex ) {
1020 2753562 : unsigned curr=getTaskCode( taskIndex );
1021 :
1022 2753562 : if( usespecies || isDensity() ) {
1023 1459773 : return getPositionOfAtomForLinkCells(curr);
1024 1293789 : } else if( nblock>0 ) {
1025 : // double factor=1.0/static_cast<double>( ablocks.size() );
1026 2130 : Vector mypos; mypos.zero();
1027 2130 : std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
1028 17040 : for(unsigned i=0; i<ablocks.size(); ++i) {
1029 8520 : if( use_for_central_atom[i] ) mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(atoms[i]);
1030 : }
1031 2130 : return mypos;
1032 : } else {
1033 1291659 : Vector mypos; mypos.zero();
1034 14122431 : for(unsigned i=0; i<ablocks.size(); ++i) {
1035 5169550 : if( use_for_central_atom[i] ) mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(ablocks[i][curr]);
1036 : }
1037 1291659 : return mypos;
1038 : }
1039 : }
1040 :
1041 605385 : void MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsigned& taskIndex, CatomPack& mypack ) {
1042 605385 : unsigned curr=getTaskCode( taskIndex );
1043 :
1044 605385 : if(usespecies) {
1045 464996 : if( mypack.getNumberOfAtomsWithDerivatives()!=1 ) mypack.resize(1);
1046 464996 : mypack.setIndex( 0, basn + curr );
1047 929992 : mypack.setDerivative( 0, Tensor::identity() );
1048 140389 : } else if( nblock>0 ) {
1049 0 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral);
1050 : unsigned k=0;
1051 0 : std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
1052 0 : for(unsigned i=0; i<ablocks.size(); ++i) {
1053 0 : if( use_for_central_atom[i] ) {
1054 0 : mypack.setIndex( k, basn + atoms[i] );
1055 0 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
1056 0 : k++;
1057 : }
1058 : }
1059 : } else {
1060 140389 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral);
1061 : unsigned k=0;
1062 809705 : for(unsigned i=0; i<ablocks.size(); ++i) {
1063 176309 : if( use_for_central_atom[i] ) {
1064 347338 : mypack.setIndex( k, basn + ablocks[i][curr] );
1065 347338 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
1066 173669 : k++;
1067 : }
1068 : }
1069 : }
1070 605385 : }
1071 :
1072 120056683 : Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 ) const {
1073 120056683 : if(usepbc) { return pbcDistance( vec1, vec2 ); }
1074 0 : else { return delta( vec1, vec2 ); }
1075 : }
1076 :
1077 234812 : void MultiColvarBase::applyPbc(std::vector<Vector>& dlist, unsigned int max_index) const {
1078 234812 : if (usepbc) pbcApply(dlist, max_index);
1079 234812 : }
1080 :
1081 1744 : void MultiColvarBase::apply() {
1082 1744 : if( getForcesFromVessels( forcesToApply ) ) setForcesOnAtoms( forcesToApply );
1083 1744 : }
1084 :
1085 : }
1086 4839 : }
|
