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 414 : void MultiColvarBase::registerKeywords( Keywords& keys ) {
41 414 : Action::registerKeywords( keys );
42 414 : ActionWithValue::registerKeywords( keys );
43 414 : ActionAtomistic::registerKeywords( keys );
44 1242 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
45 414 : ActionWithVessel::registerKeywords( keys );
46 : keys.add("hidden","NL_STRIDE","the frequency with which the neighbor list should be updated. Between neighbour list update steps all quantities "
47 1656 : "that contributed less than TOL at the previous neighbor list update step are ignored.");
48 : 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 of 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 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 "
60 828 : "input you can customize the component name");
61 : 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 1656 : "coordination number more than four for example");
67 : 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 coordination 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 1656 : "using the label of another multicolvar");
71 : 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 1656 : "the documentation for that keyword");
73 1656 : 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");
74 414 : }
75 :
76 340 : MultiColvarBase::MultiColvarBase(const ActionOptions&ao):
77 : Action(ao),
78 : ActionAtomistic(ao),
79 : ActionWithValue(ao),
80 : ActionWithVessel(ao),
81 : usepbc(false),
82 : allthirdblockintasks(false),
83 : uselinkforthree(false),
84 : linkcells(comm),
85 : threecells(comm),
86 : setup_completed(false),
87 : atomsWereRetrieved(false),
88 : matsums(false),
89 : usespecies(false),
90 1360 : nblock(0)
91 : {
92 680 : if( keywords.exists("NOPBC") ) {
93 680 : bool nopbc=!usepbc; parseFlag("NOPBC",nopbc);
94 340 : usepbc=!nopbc;
95 : }
96 680 : if( keywords.exists("SPECIESA") ) { matsums=usespecies=true; }
97 340 : }
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 ) return;
102 :
103 : std::vector<AtomNumber> t;
104 922 : for(int i=1;; ++i ) {
105 970 : parseAtomList(key, i, t );
106 970 : if( t.empty() ) break;
107 :
108 922 : log.printf(" Colvar %d is calculated from atoms : ", i);
109 5978 : for(unsigned j=0; j<t.size(); ++j) log.printf("%d ",t[j].serial() );
110 922 : log.printf("\n");
111 :
112 928 : if( i==1 && natoms<0 ) { ablocks.resize(t.size()); }
113 916 : else if( i==1 ) ablocks.resize(natoms);
114 922 : if( t.size()!=ablocks.size() ) {
115 0 : std::string ss; Tools::convert(i,ss);
116 0 : error(key + ss + " keyword has the wrong number of atoms");
117 : }
118 5978 : for(unsigned j=0; j<ablocks.size(); ++j) {
119 7584 : ablocks[j].push_back( ablocks.size()*(i-1)+j ); all_atoms.push_back( t[j] );
120 7584 : atom_lab.push_back( std::pair<unsigned,unsigned>( 0, ablocks.size()*(i-1)+j ) );
121 : }
122 922 : t.resize(0);
123 922 : }
124 48 : if( all_atoms.size()>0 ) {
125 48 : nblock=0;
126 1018 : for(unsigned i=0; i<ablocks[0].size(); ++i) addTaskToList( i );
127 : }
128 : }
129 :
130 456 : bool MultiColvarBase::parseMultiColvarAtomList(const std::string& key, const int& num, std::vector<AtomNumber>& t) {
131 : std::vector<std::string> mlabs;
132 456 : if( num<0 ) parseVector(key,mlabs);
133 0 : else parseNumberedVector(key,num,mlabs);
134 :
135 456 : if( mlabs.size()==0 ) return false;
136 :
137 311 : std::string mname; unsigned found_mcolv=mlabs.size();
138 876 : for(unsigned i=0; i<mlabs.size(); ++i) {
139 632 : MultiColvarBase* mycolv = plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlabs[i]);
140 316 : if(!mycolv) { found_mcolv=i; break; }
141 : // Check all base multicolvars are of same type
142 127 : if( i==0 ) {
143 122 : mname = mycolv->getName();
144 244 : if( keywords.exists("ALL_INPUT_SAME_TYPE") && mycolv->isPeriodic() ) error("multicolvar functions don't work with this multicolvar");
145 : } else {
146 14 : if( keywords.exists("ALL_INPUT_SAME_TYPE") && mname!=mycolv->getName() ) error("All input multicolvars must be of same type");
147 : }
148 : // And track which variable stores each colvar
149 17833371 : for(unsigned j=0; j<mycolv->getFullNumberOfTasks(); ++j) atom_lab.push_back( std::pair<unsigned,unsigned>( mybasemulticolvars.size()+1, j ) );
150 : // And store the multicolvar base
151 127 : mybasemulticolvars.push_back( mycolv );
152 : // And create a basedata stash
153 508 : mybasedata.push_back( mybasemulticolvars[mybasemulticolvars.size()-1]->buildDataStashes( this ) );
154 : // Check if weight has derivatives
155 254 : if( mybasemulticolvars[ mybasemulticolvars.size()-1 ]->weightHasDerivatives ) weightHasDerivatives=true;
156 127 : plumed_assert( mybasemulticolvars.size()==mybasedata.size() );
157 : }
158 : // Have we conventional atoms to read in
159 311 : if( found_mcolv==0 ) {
160 189 : std::vector<AtomNumber> tt; ActionAtomistic::interpretAtomList( mlabs, tt );
161 351049 : for(unsigned i=0; i<tt.size(); ++i) { atom_lab.push_back( std::pair<unsigned,unsigned>( 0, t.size() + i ) ); }
162 189 : log.printf(" keyword %s takes atoms : ", key.c_str() );
163 263334 : for(unsigned i=0; i<tt.size(); ++i) { t.push_back( tt[i] ); log.printf("%d ",tt[i].serial() ); }
164 189 : log.printf("\n");
165 122 : } else if( found_mcolv==mlabs.size() ) {
166 122 : if( checkNumericalDerivatives() ) error("cannot use NUMERICAL_DERIVATIVES keyword with dynamic groups as input");
167 122 : log.printf(" keyword %s takes dynamic groups of atoms constructed from multicolvars labelled : ", key.c_str() );
168 376 : for(unsigned i=0; i<mlabs.size(); ++i) log.printf("%s ",mlabs[i].c_str() );
169 122 : log.printf("\n");
170 0 : } else if( found_mcolv<mlabs.size() ) {
171 0 : error("cannot mix multicolvar input and atom input in one line");
172 : }
173 456 : return true;
174 : }
175 :
176 76 : void MultiColvarBase::readTwoGroups( const std::string& key0, const std::string& key1, const std::string& key2, std::vector<AtomNumber>& all_atoms ) {
177 76 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) ); ablocks.resize( 2 );
178 :
179 76 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
180 78 : nblock=atom_lab.size(); for(unsigned i=0; i<2; ++i) ablocks[i].resize(nblock);
181 26 : resizeBookeepingArray( nblock, nblock );
182 5486 : for(unsigned i=0; i<nblock; ++i) ablocks[0][i]=ablocks[1][i]=i;
183 5434 : for(unsigned i=1; i<nblock; ++i) {
184 4210895 : for(unsigned j=0; j<i; ++j) {
185 4210895 : bookeeping(i,j).first=getFullNumberOfTasks();
186 4210895 : addTaskToList( i*nblock + j );
187 4210895 : bookeeping(i,j).second=getFullNumberOfTasks();
188 : }
189 : }
190 : } else {
191 50 : parseMultiColvarAtomList(key1,-1,all_atoms);
192 134 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=i;
193 50 : parseMultiColvarAtomList(key2,-1,all_atoms);
194 2238 : 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;
195 :
196 50 : if( ablocks[0].size()>ablocks[1].size() ) nblock = ablocks[0].size();
197 50 : else nblock=ablocks[1].size();
198 :
199 100 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
200 134 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
201 2235 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
202 1109 : bookeeping(i,j).first=getFullNumberOfTasks();
203 2218 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 ) {
204 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
205 0 : atom_lab[ablocks[0][i]].second!=atom_lab[ablocks[1][j]].second ) addTaskToList( i*nblock + j );
206 4436 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] ) addTaskToList( i*nblock + j );
207 1109 : bookeeping(i,j).second=getFullNumberOfTasks();
208 : }
209 : }
210 : }
211 76 : }
212 :
213 12 : void MultiColvarBase::readGroupKeywords( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& key3,
214 : const bool& no_third_dim_accum, const bool& symmetric, std::vector<AtomNumber>& all_atoms ) {
215 12 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) ); ablocks.resize( 3 );
216 :
217 12 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
218 10 : if( no_third_dim_accum ) {
219 30 : nblock=atom_lab.size(); ablocks[0].resize(nblock); ablocks[1].resize( nblock );
220 2962 : for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=ablocks[1][i]=i;
221 10 : resizeBookeepingArray( nblock, nblock );
222 10 : if( symmetric ) {
223 : // This ensures that later parts of the code don't switch off allthirdblockintasks
224 2674 : for(unsigned i=0; i<nblock; ++i) { bookeeping(i,i).first=0; bookeeping(i,i).second=1; }
225 1331 : for(unsigned i=1; i<nblock; ++i) {
226 221116 : for(unsigned j=0; j<i; ++j) {
227 221116 : bookeeping(j,i).first=bookeeping(i,j).first=getFullNumberOfTasks();
228 221116 : addTaskToList( i*nblock + j );
229 221116 : bookeeping(j,i).second=bookeeping(i,j).second=getFullNumberOfTasks();
230 : }
231 : }
232 : } else {
233 134 : for(unsigned i=0; i<nblock; ++i) {
234 8210 : for(unsigned j=0; j<nblock; ++j) {
235 8210 : if( i==j ) continue ;
236 8076 : bookeeping(i,j).first=getFullNumberOfTasks();
237 8076 : addTaskToList( i*nblock + j );
238 8076 : bookeeping(i,j).second=getFullNumberOfTasks();
239 : }
240 : }
241 : }
242 10 : if( !parseMultiColvarAtomList(key3,-1,all_atoms) ) error("missing required keyword " + key3 + " in input");
243 10 : ablocks[2].resize( atom_lab.size() - ablocks[0].size() );
244 99690 : for(unsigned i=0; i<ablocks[2].size(); ++i) ablocks[2][i]=ablocks[0].size() + i;
245 : } else {
246 0 : nblock=atom_lab.size(); for(unsigned i=0; i<3; ++i) ablocks[i].resize(nblock);
247 0 : resizeBookeepingArray( nblock, nblock );
248 0 : for(unsigned i=0; i<nblock; ++i) { ablocks[0][i]=i; ablocks[1][i]=i; ablocks[2][i]=i; }
249 0 : if( symmetric ) {
250 0 : for(unsigned i=2; i<nblock; ++i) {
251 0 : for(unsigned j=1; j<i; ++j) {
252 0 : bookeeping(i,j).first=getFullNumberOfTasks();
253 0 : for(unsigned k=0; k<j; ++k) addTaskToList( i*nblock*nblock + j*nblock + k );
254 0 : bookeeping(i,j).second=getFullNumberOfTasks();
255 : }
256 : }
257 : } else {
258 0 : for(unsigned i=0; i<nblock; ++i) {
259 0 : for(unsigned j=0; j<nblock; ++j) {
260 0 : if( i==j ) continue;
261 0 : bookeeping(i,j).first=getFullNumberOfTasks();
262 0 : for(unsigned k=0; k<nblock; ++k) {
263 0 : if( i!=k && j!=k ) addTaskToList( i*nblock*nblock + j*nblock + k );
264 : }
265 0 : bookeeping(i,j).first=getFullNumberOfTasks();
266 : }
267 : }
268 : }
269 : }
270 : } else {
271 2 : readThreeGroups( key1, key2, key3, true, no_third_dim_accum, all_atoms );
272 : }
273 :
274 12 : }
275 :
276 11 : void MultiColvarBase::readThreeGroups( const std::string& key1, const std::string& key2, const std::string& key3,
277 : const bool& allow2, const bool& no_third_dim_accum, std::vector<AtomNumber>& all_atoms ) {
278 11 : plumed_dbg_assert( keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) ); ablocks.resize( 3 );
279 :
280 11 : bool readkey1=parseMultiColvarAtomList(key1,-1,all_atoms);
281 84 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=i;
282 11 : bool readkey2=parseMultiColvarAtomList(key2,-1,all_atoms);
283 22 : if( !readkey1 && !readkey2 ) return ;
284 472 : 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;
285 :
286 22 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
287 11 : bool readkey3=parseMultiColvarAtomList(key3,-1,all_atoms);
288 11 : if( !readkey3 && !allow2 ) {
289 0 : error("missing atom specification " + key3);
290 11 : } else if( !readkey3 ) {
291 2 : if( ablocks[1].size()>ablocks[0].size() ) nblock=ablocks[1].size();
292 0 : else nblock=ablocks[0].size();
293 :
294 2 : ablocks[2].resize( ablocks[1].size() );
295 400 : for(unsigned i=0; i<ablocks[1].size(); ++i) ablocks[2][i]=ablocks[0].size() + i;
296 6 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
297 394 : for(unsigned j=1; j<ablocks[1].size(); ++j) {
298 196 : bookeeping(i,j).first=getFullNumberOfTasks();
299 9898 : for(unsigned k=0; k<j; ++k) {
300 19404 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
301 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
302 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
303 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
304 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 &&
305 0 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) addTaskToList( nblock*nblock*i + nblock*j + k );
306 48510 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
307 38808 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
308 9702 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) addTaskToList( nblock*nblock*i + nblock*j + k );
309 : }
310 196 : bookeeping(i,j).second=getFullNumberOfTasks();
311 : }
312 : }
313 : } else {
314 18 : ablocks[2].resize( atom_lab.size() - ablocks[1].size() - ablocks[0].size() );
315 9570 : for(unsigned i=0; i<ablocks[2].size(); ++i) ablocks[2][i] = ablocks[0].size() + ablocks[1].size() + i;
316 :
317 9 : if( ablocks[1].size()>ablocks[0].size() ) nblock=ablocks[1].size();
318 9 : else nblock=ablocks[0].size();
319 9 : if( ablocks[2].size()>nblock ) nblock=ablocks[2].size();
320 :
321 9 : unsigned kcount; if( no_third_dim_accum ) kcount=1; else kcount=ablocks[2].size();
322 :
323 76 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
324 447 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
325 209 : bookeeping(i,j).first=getFullNumberOfTasks();
326 418 : for(unsigned k=0; k<kcount; ++k) {
327 209 : if( no_third_dim_accum ) addTaskToList( nblock*i + j );
328 0 : else if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
329 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
330 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
331 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
332 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 &&
333 0 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) addTaskToList( nblock*nblock*i + nblock*j + k );
334 0 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
335 0 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
336 0 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) addTaskToList( nblock*nblock*i + nblock*j + k );
337 : }
338 209 : bookeeping(i,j).second=getFullNumberOfTasks();
339 : }
340 : }
341 : }
342 : }
343 :
344 4 : void MultiColvarBase::buildSets() {
345 : std::vector<AtomNumber> fake_atoms;
346 8 : if( !parseMultiColvarAtomList("DATA",-1,fake_atoms) ) error("missing DATA keyword");
347 4 : if( fake_atoms.size()>0 ) error("no atoms should appear in the specification for this object. Input should be other multicolvars");
348 :
349 8 : nblock = mybasemulticolvars[0]->getFullNumberOfTasks();
350 22 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
351 14 : if( mybasemulticolvars[i]->getFullNumberOfTasks()!=nblock ) {
352 0 : error("mismatch between numbers of tasks in various base multicolvars");
353 : }
354 : }
355 4 : ablocks.resize( mybasemulticolvars.size() ); usespecies=false;
356 22 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
357 14 : ablocks[i].resize( nblock );
358 40 : for(unsigned j=0; j<nblock; ++j) ablocks[i][j]=i*nblock+j;
359 : }
360 15 : for(unsigned i=0; i<nblock; ++i) {
361 11 : if( mybasemulticolvars.size()<4 ) {
362 11 : unsigned cvcode=0, tmpc=1;
363 42 : for(unsigned j=0; j<ablocks.size(); ++j) { cvcode += i*tmpc; tmpc *= nblock; }
364 11 : addTaskToList( cvcode );
365 : } else {
366 0 : addTaskToList( i );
367 : }
368 : }
369 8 : mybasedata[0]->resizeTemporyMultiValues( mybasemulticolvars.size() ); setupMultiColvarBase( fake_atoms );
370 4 : }
371 :
372 12511150 : void MultiColvarBase::addTaskToList( const unsigned& taskCode ) {
373 12511150 : plumed_assert( getNumberOfVessels()==0 );
374 12511150 : ActionWithVessel::addTaskToList( taskCode );
375 12511150 : }
376 :
377 97 : void MultiColvarBase::resizeBookeepingArray( const unsigned& num1, const unsigned& num2 ) {
378 97 : bookeeping.resize( num1, num2 );
379 7076 : for(unsigned i=0; i<num1; ++i) {
380 17761090 : for(unsigned j=0; j<num2; ++j) { bookeeping(i,j).first=0; bookeeping(i,j).second=0; }
381 : }
382 97 : }
383 :
384 298 : void MultiColvarBase::setupMultiColvarBase( const std::vector<AtomNumber>& atoms ) {
385 426 : if( !matsums && atom_lab.size()==0 ) error("No atoms have been read in");
386 : std::vector<AtomNumber> all_atoms;
387 : // Setup decoder array
388 298 : if( !usespecies && nblock>0 ) {
389 :
390 64 : ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
391 64 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
392 64 : if( ablocks.size()==3 ) {
393 21 : allthirdblockintasks=uselinkforthree=true;
394 3046 : for(unsigned i=0; i<bookeeping.nrows(); ++i) {
395 905874 : for(unsigned j=0; j<bookeeping.ncols(); ++j) {
396 452186 : unsigned ntper = bookeeping(i,j).second - bookeeping(i,j).first;
397 452186 : if( i==j && ntper==0 ) {
398 : continue;
399 452050 : } else if( ntper == 1 && allthirdblockintasks ) {
400 451856 : allthirdblockintasks=true;
401 388 : } else if( ntper != ablocks[2].size() ) {
402 194 : allthirdblockintasks=uselinkforthree=false;
403 : } else {
404 0 : allthirdblockintasks=false;
405 : }
406 : }
407 : }
408 : }
409 :
410 64 : if( allthirdblockintasks ) {
411 38 : decoder.resize(2); plumed_assert( ablocks.size()==3 );
412 : // Check if number of atoms is too large
413 19 : if( pow( double(nblock), 2.0 )>std::numeric_limits<unsigned>::max() ) error("number of atoms in groups is too big for PLUMED to handle");
414 : } else {
415 45 : decoder.resize( ablocks.size() );
416 : // Check if number of atoms is too large
417 45 : 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");
418 : }
419 451 : unsigned code=1; for(unsigned i=0; i<decoder.size(); ++i) { decoder[decoder.size()-1-i]=code; code *= nblock; }
420 234 : } else if( !usespecies ) {
421 87 : ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
422 87 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
423 294 : } else if( keywords.exists("SPECIESA") ) {
424 178 : plumed_assert( atom_lab.size()==0 && all_atoms.size()==0 );
425 178 : ablocks.resize( 1 ); bool readspecies=parseMultiColvarAtomList("SPECIES", -1, all_atoms);
426 89 : if( readspecies ) {
427 79896 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<atom_lab.size(); ++i) { addTaskToList(i); ablocks[0][i]=i; }
428 : } else {
429 34 : if( !parseMultiColvarAtomList("SPECIESA", -1, all_atoms) ) error("missing SPECIES/SPECIESA keyword");
430 17 : unsigned nat1=atom_lab.size();
431 34 : if( !parseMultiColvarAtomList("SPECIESB", -1, all_atoms) ) error("missing SPECIESB keyword");
432 17 : unsigned nat2=atom_lab.size() - nat1;
433 :
434 17 : for(unsigned i=0; i<nat1; ++i) addTaskToList(i);
435 34 : ablocks[0].resize( nat2 );
436 3398 : for(unsigned i=0; i<nat2; ++i) {
437 : bool found=false; unsigned inum;
438 261831 : for(unsigned j=0; j<nat1; ++j) {
439 776546 : if( atom_lab[nat1+i].first>0 && atom_lab[j].first>0 ) {
440 1010880 : if( atom_lab[nat1+i].first==atom_lab[j].first &&
441 252720 : mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
442 0 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) { found=true; inum=j; break; }
443 9193 : } else if( atom_lab[nat1+i].first>0 ) {
444 0 : if( mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
445 0 : all_atoms[atom_lab[j].second] ) { found=true; inum=nat1 + i; break; }
446 18386 : } else if( atom_lab[j].first>0 ) {
447 0 : if( all_atoms[atom_lab[nat1+i].second]==
448 0 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) { found=true; inum=nat1+i; break; }
449 27579 : } else if( all_atoms[atom_lab[nat1+i].second]==all_atoms[atom_lab[j].second] ) { found=true; inum=j; break; }
450 : }
451 : // This prevents mistakes being made in colvar setup
452 3480 : if( found ) { ablocks[0][i]=inum; }
453 6632 : else { ablocks[0][i]=nat1 + i; }
454 : }
455 : }
456 : }
457 298 : if( mybasemulticolvars.size()>0 ) {
458 373 : for(unsigned i=0; i<mybasedata.size(); ++i) {
459 381 : mybasedata[i]->resizeTemporyMultiValues(2); mybasemulticolvars[i]->my_tmp_capacks.resize(2);
460 : }
461 : }
462 :
463 : // Copy lists of atoms involved from base multicolvars
464 : std::vector<AtomNumber> tmp_atoms;
465 850 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
466 127 : tmp_atoms=mybasemulticolvars[i]->getAbsoluteIndexes();
467 832505 : for(unsigned j=0; j<tmp_atoms.size(); ++j) all_atoms.push_back( tmp_atoms[j] );
468 : }
469 : // Copy atom lists from input
470 118892 : for(unsigned i=0; i<atoms.size(); ++i) all_atoms.push_back( atoms[i] );
471 :
472 : // Now make sure we get all the atom positions
473 298 : ActionAtomistic::requestAtoms( all_atoms );
474 : // And setup dependencies
475 552 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) addDependency( mybasemulticolvars[i] );
476 :
477 : // Setup underlying ActionWithVessel
478 298 : readVesselKeywords();
479 298 : }
480 :
481 20 : void MultiColvarBase::setAtomsForCentralAtom( const std::vector<bool>& catom_ind ) {
482 20 : unsigned nat=0; plumed_assert( catom_ind.size()==ablocks.size() );
483 166 : for(unsigned i=0; i<catom_ind.size(); ++i) {
484 : use_for_central_atom[i]=catom_ind[i];
485 73 : if( use_for_central_atom[i] ) nat++;
486 : }
487 20 : plumed_dbg_assert( nat>0 ); ncentral=nat;
488 20 : numberForCentralAtom = 1.0 / static_cast<double>( nat );
489 20 : }
490 :
491 427 : void MultiColvarBase::turnOnDerivatives() {
492 427 : ActionWithValue::turnOnDerivatives();
493 427 : needsDerivatives();
494 427 : forcesToApply.resize( getNumberOfDerivatives() );
495 427 : }
496 :
497 181 : void MultiColvarBase::setLinkCellCutoff( const double& lcut, double tcut ) {
498 271 : plumed_assert( usespecies || ablocks.size()<4 );
499 181 : if( tcut<0 ) tcut=lcut;
500 :
501 181 : if( !linkcells.enabled() ) {
502 181 : linkcells.setCutoff( lcut );
503 181 : threecells.setCutoff( tcut );
504 : } else {
505 0 : if( lcut>linkcells.getCutoff() ) linkcells.setCutoff( lcut );
506 0 : if( tcut>threecells.getCutoff() ) threecells.setCutoff( tcut );
507 : }
508 181 : }
509 :
510 8 : double MultiColvarBase::getLinkCellCutoff() const {
511 8 : return linkcells.getCutoff();
512 : }
513 :
514 1924 : void MultiColvarBase::setupLinkCells() {
515 4980 : if( (!usespecies && nblock==0) || !linkcells.enabled() ) return ;
516 : // Retrieve any atoms that haven't already been retrieved
517 2405 : for(std::vector<MultiColvarBase*>::iterator p=mybasemulticolvars.begin(); p!=mybasemulticolvars.end(); ++p) {
518 203 : (*p)->retrieveAtoms();
519 : }
520 1101 : retrieveAtoms();
521 :
522 : unsigned iblock;
523 1101 : if( usespecies ) {
524 : iblock=0;
525 222 : } else if( ablocks.size()<4 ) {
526 : iblock=1;
527 : } else {
528 0 : plumed_error();
529 : }
530 :
531 : // Count number of currently active atoms
532 1101 : nactive_atoms=0;
533 772130 : for(unsigned i=0; i<ablocks[iblock].size(); ++i) {
534 769928 : if( isCurrentlyActive( ablocks[iblock][i] ) ) nactive_atoms++;
535 : }
536 :
537 1101 : if( nactive_atoms>0 ) {
538 1101 : std::vector<Vector> ltmp_pos( nactive_atoms );
539 1101 : std::vector<unsigned> ltmp_ind( nactive_atoms );
540 :
541 1101 : nactive_atoms=0;
542 1101 : if( usespecies ) {
543 728511 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
544 727632 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
545 727440 : ltmp_ind[nactive_atoms]=ablocks[0][i];
546 1091160 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ltmp_ind[nactive_atoms] );
547 363720 : nactive_atoms++;
548 : }
549 : } else {
550 42518 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
551 42296 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
552 32456 : ltmp_ind[nactive_atoms]=i;
553 48684 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ablocks[1][i] );
554 16228 : nactive_atoms++;
555 : }
556 : }
557 :
558 : // Build the lists for the link cells
559 1101 : linkcells.buildCellLists( ltmp_pos, ltmp_ind, getPbc() );
560 : }
561 : }
562 :
563 5128 : void MultiColvarBase::setupNonUseSpeciesLinkCells( const unsigned& my_always_active ) {
564 5128 : plumed_assert( !usespecies );
565 10647 : if( nblock==0 || !linkcells.enabled() ) return ;
566 710 : deactivateAllTasks();
567 : std::vector<unsigned> requiredlinkcells;
568 :
569 710 : if( !uselinkforthree && nactive_atoms>0 ) {
570 : // Get some parallel info
571 156 : unsigned stride=comm.Get_size();
572 156 : unsigned rank=comm.Get_rank();
573 156 : if( serialCalculation() ) { stride=1; rank=0; }
574 :
575 : // Ensure we only do tasks where atoms are in appropriate link cells
576 156 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
577 11514 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
578 14181 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
579 5244 : unsigned natomsper=1; linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
580 5244 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, natomsper, linked_atoms );
581 205705 : for(unsigned j=0; j<natomsper; ++j) {
582 909319 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) taskFlags[k]=1;
583 : }
584 156 : }
585 554 : } else if( nactive_atoms>0 ) {
586 : // Get some parallel info
587 554 : unsigned stride=comm.Get_size();
588 554 : unsigned rank=comm.Get_rank();
589 554 : if( serialCalculation() ) { stride=1; rank=0; }
590 :
591 : unsigned nactive_three=0;
592 144198 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
593 143090 : if( isCurrentlyActive( ablocks[2][i] ) ) nactive_three++;
594 : }
595 :
596 554 : std::vector<Vector> lttmp_pos( nactive_three );
597 554 : std::vector<unsigned> lttmp_ind( nactive_three );
598 :
599 : nactive_three=0;
600 554 : if( allthirdblockintasks ) {
601 143644 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
602 143090 : if( !isCurrentlyActive( ablocks[2][i] ) ) continue;
603 143090 : lttmp_ind[nactive_three]=ablocks[2][i];
604 143090 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
605 71545 : nactive_three++;
606 : }
607 : } else {
608 0 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
609 0 : if( !isCurrentlyActive( ablocks[2][i] ) ) continue;
610 0 : lttmp_ind[nactive_three]=i;
611 0 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
612 0 : nactive_three++;
613 : }
614 : }
615 : // Build the list of the link cells
616 554 : threecells.buildCellLists( lttmp_pos, lttmp_ind, getPbc() );
617 :
618 : // Ensure we only do tasks where atoms are in appropriate link cells
619 554 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
620 554 : std::vector<unsigned> tlinked_atoms( 1+ablocks[2].size() );
621 12266 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
622 11158 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
623 11158 : unsigned natomsper=1; linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
624 11158 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, natomsper, linked_atoms );
625 5579 : if( allthirdblockintasks ) {
626 121253 : for(unsigned j=0; j<natomsper; ++j) {
627 605997 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) taskFlags[k]=1;
628 : }
629 : } else {
630 0 : unsigned ntatomsper=1; tlinked_atoms[0]=lttmp_ind[0];
631 0 : threecells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, ntatomsper, tlinked_atoms );
632 0 : for(unsigned j=0; j<natomsper; ++j) {
633 0 : for(unsigned k=0; k<ntatomsper; ++k) taskFlags[bookeeping(i,linked_atoms[j]).first+tlinked_atoms[k]]=1;
634 : }
635 : }
636 : }
637 : }
638 710 : if( !serialCalculation() ) comm.Sum( taskFlags );
639 710 : lockContributors();
640 : }
641 :
642 295976 : void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const {
643 : plumed_dbg_assert( !usespecies && nblock>0 );
644 320467 : if( atoms.size()!=decoder.size() ) atoms.resize( decoder.size() );
645 :
646 295965 : unsigned scode = taskCode;
647 1872590 : for(unsigned i=0; i<decoder.size(); ++i) {
648 640330 : unsigned ind=( scode / decoder[i] );
649 1280660 : atoms[i] = ablocks[i][ind];
650 640330 : scode -= ind*decoder[i];
651 : }
652 295965 : }
653 :
654 471197 : bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ) const {
655 471197 : if( isDensity() ) {
656 19047 : myatoms.setNumberOfAtoms( 1 ); myatoms.setAtom( 0, taskCode ); return true;
657 452085 : } else if( usespecies ) {
658 301195 : std::vector<unsigned> task_atoms(1); task_atoms[0]=taskCode;
659 150598 : unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getPositionOfAtomForLinkCells( taskCode ), linkcells );
660 150598 : return natomsper>1;
661 301488 : } else if( matsums ) {
662 : myatoms.setNumberOfAtoms( getNumberOfAtoms() );
663 53239 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) myatoms.setAtom( i, i );
664 301047 : } else if( allthirdblockintasks ) {
665 89798 : plumed_dbg_assert( ablocks.size()==3 ); std::vector<unsigned> atoms(2); decodeIndexToAtoms( taskCode, atoms );
666 179626 : myatoms.setupAtomsFromLinkCells( atoms, getPositionOfAtomForLinkCells( atoms[0] ), threecells );
667 211249 : } else if( nblock>0 ) {
668 179544 : std::vector<unsigned> atoms( ablocks.size() );
669 359068 : decodeIndexToAtoms( taskCode, atoms ); myatoms.setNumberOfAtoms( ablocks.size() );
670 1174092 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, atoms[i] );
671 : } else {
672 31705 : myatoms.setNumberOfAtoms( ablocks.size() );
673 248616 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, ablocks[i][taskCode] );
674 : }
675 : return true;
676 : }
677 :
678 9257 : void MultiColvarBase::setupActiveTaskSet( std::vector<unsigned>& active_tasks, const std::string& input_label ) {
679 9257 : if( !setup_completed ) {
680 : bool justVolumes=false;
681 1916 : if( usespecies ) {
682 : justVolumes=true;
683 1938 : for(unsigned i=0; i<getNumberOfVessels(); ++i) {
684 1298 : vesselbase::StoreDataVessel* mys=dynamic_cast<vesselbase::StoreDataVessel*>( getPntrToVessel(i) );
685 1298 : if( mys ) continue;
686 800 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(i) );
687 800 : if( !myb ) { justVolumes=false; break; }
688 38 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
689 38 : if( !myv ) { justVolumes=false; break; }
690 : }
691 : }
692 1916 : deactivateAllTasks();
693 1916 : if( justVolumes && mydata ) {
694 78 : if( mydata->getNumberOfDataUsers()==0 ) justVolumes=false;
695 :
696 254 : for(unsigned i=0; i<mydata->getNumberOfDataUsers(); ++i) {
697 49 : MultiColvarBase* myu=dynamic_cast<MultiColvarBase*>( mydata->getDataUser(i) );
698 49 : if( myu ) {
699 49 : myu->setupActiveTaskSet( taskFlags, getLabel() );
700 : } else {
701 0 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
702 : }
703 : }
704 : }
705 1916 : if( justVolumes ) {
706 240 : for(unsigned j=0; j<getNumberOfVessels(); ++j) {
707 80 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(j) );
708 80 : if( !myb ) continue ;
709 32 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
710 32 : if( !myv ) continue ;
711 32 : myv->retrieveAtoms(); myv->setupRegions();
712 :
713 296322 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
714 151567 : if( myv->inVolumeOfInterest(i) ) taskFlags[i]=1;
715 : }
716 : }
717 : } else {
718 14652693 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
719 : }
720 :
721 : // Now activate all this class
722 1916 : lockContributors();
723 : // Setup the link cells
724 1916 : setupLinkCells();
725 : // Ensures that setup is not performed multiple times during one cycle
726 1916 : setup_completed=true;
727 : }
728 :
729 : // And activate the tasks in input action
730 18514 : if( getLabel()!=input_label ) {
731 : int input_code=-1;
732 73 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
733 61 : if( mybasemulticolvars[i]->getLabel()==input_label ) { input_code=i+1; break; }
734 : }
735 :
736 49 : MultiValue my_tvals( getNumberOfQuantities(), getNumberOfDerivatives() );
737 49 : AtomValuePack mytmp_atoms( my_tvals, this );
738 296468 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
739 148185 : if( !taskIsCurrentlyActive(i) ) continue;
740 3529 : setupCurrentAtomList( getTaskCode(i), mytmp_atoms );
741 509592 : for(unsigned j=0; j<mytmp_atoms.getNumberOfAtoms(); ++j) {
742 : unsigned itask=mytmp_atoms.getIndex(j);
743 753101 : if( atom_lab[itask].first==input_code ) active_tasks[ atom_lab[itask].second ]=1;
744 : }
745 49 : }
746 : }
747 9257 : }
748 :
749 467 : bool MultiColvarBase::filtersUsedAsInput() {
750 : bool inputAreFilters=false;
751 1456 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
752 261 : MultiColvarFilter* myfilt=dynamic_cast<MultiColvarFilter*>( mybasemulticolvars[i] );
753 261 : if( myfilt || mybasemulticolvars[i]->filtersUsedAsInput() ) inputAreFilters=true;
754 : }
755 467 : return inputAreFilters;
756 : }
757 :
758 9208 : void MultiColvarBase::calculate() {
759 : // Recursive function that sets up tasks
760 9208 : setupActiveTaskSet( taskFlags, getLabel() );
761 :
762 : // Check for filters and rerun setup of link cells if there are any
763 9208 : if( mybasemulticolvars.size()>0 && filtersUsedAsInput() ) setupLinkCells();
764 :
765 : // Setup the link cells if we are not using species
766 16353 : if( !usespecies && ablocks.size()>1 ) {
767 : // This loop finds the first active atom, which is always checked because
768 : // of a peculiarity in linkcells
769 5128 : unsigned first_active=std::numeric_limits<unsigned>::max();
770 10524 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
771 5262 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
772 : else {
773 5128 : first_active=i; break;
774 : }
775 : }
776 5128 : setupNonUseSpeciesLinkCells( first_active );
777 : }
778 : // And run all tasks
779 9208 : runAllTasks();
780 9208 : }
781 :
782 217 : void MultiColvarBase::calculateNumericalDerivatives( ActionWithValue* a ) {
783 217 : if( mybasemulticolvars.size()>0 ) plumed_merror("cannot calculate numerical derivatives for this quantity");
784 217 : calculateAtomicNumericalDerivatives( this, 0 );
785 217 : }
786 :
787 2945 : void MultiColvarBase::prepare() {
788 2945 : setup_completed=false; atomsWereRetrieved=false;
789 2945 : }
790 :
791 6520 : void MultiColvarBase::retrieveAtoms() {
792 6520 : if( !atomsWereRetrieved ) { ActionAtomistic::retrieveAtoms(); atomsWereRetrieved=true; }
793 6520 : }
794 :
795 38245 : void MultiColvarBase::mergeInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
796 : const unsigned& jatom, const std::vector<double>& der,
797 : MultiValue& myder, AtomValuePack& myatoms ) const {
798 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
799 : plumed_dbg_assert( ival<myatoms.getUnderlyingMultiValue().getNumberOfValues() );
800 : plumed_dbg_assert( start<myder.getNumberOfValues() && end<=myder.getNumberOfValues() );
801 : plumed_dbg_assert( der.size()==myder.getNumberOfValues() && jatom<myatoms.getNumberOfAtoms() );
802 : // Convert input atom to local index
803 : unsigned katom = myatoms.getIndex( jatom ); plumed_dbg_assert( katom<atom_lab.size() ); plumed_dbg_assert( atom_lab[katom].first>0 );
804 : // Find base colvar
805 76490 : unsigned mmc=atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
806 : // Get start of indices for this atom
807 38855 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
808 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
809 38245 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
810 9611702 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
811 : unsigned jder=myder.getActiveIndex(j);
812 9535212 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
813 4434561 : unsigned kder=basen+jder;
814 110197260 : for(unsigned icomp=start; icomp<end; ++icomp) {
815 317288097 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
816 : }
817 : } else {
818 333045 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
819 7551360 : for(unsigned icomp=start; icomp<end; ++icomp) {
820 21654945 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
821 : }
822 : }
823 : }
824 38245 : }
825 :
826 106 : void MultiColvarBase::splitInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
827 : const unsigned& jatom, const std::vector<double>& der,
828 : MultiValue& myder, AtomValuePack& myatoms ) const {
829 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
830 : plumed_dbg_assert( ival<myder.getNumberOfValues() );
831 : plumed_dbg_assert( start<myvals.getNumberOfValues() && end<=myvals.getNumberOfValues() );
832 : plumed_dbg_assert( der.size()==myatoms.getUnderlyingMultiValue().getNumberOfValues() && jatom<myatoms.getNumberOfAtoms() );
833 : // Convert input atom to local index
834 : unsigned katom = myatoms.getIndex( jatom ); plumed_dbg_assert( katom<atom_lab.size() ); plumed_dbg_assert( atom_lab[katom].first>0 );
835 : // Find base colvar
836 212 : unsigned mmc=atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
837 : // Get start of indices for this atom
838 154 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
839 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
840 106 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
841 39068 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
842 : unsigned jder=myder.getActiveIndex(j);
843 38856 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
844 36948 : unsigned kder=basen+jder; plumed_assert( kder<myvals.getNumberOfDerivatives() );
845 39942 : for(unsigned icomp=start; icomp<end; ++icomp) {
846 64404 : myvals.addDerivative( icomp, kder, der[icomp]*myder.getDerivative( ival, jder ) );
847 : }
848 : } else {
849 954 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
850 3942 : for(unsigned icomp=start; icomp<end; ++icomp) {
851 8964 : myvals.addDerivative( icomp, kder, der[icomp]*myder.getDerivative( ival, jder ) );
852 : }
853 : }
854 : }
855 106 : }
856 :
857 905762 : void MultiColvarBase::addComDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
858 : plumed_dbg_assert( ival<static_cast<int>(myatoms.getUnderlyingMultiValue().getNumberOfValues()) && iatom<myatoms.getNumberOfAtoms() );
859 : // Convert input atom to local index
860 : unsigned katom = myatoms.getIndex( iatom ); plumed_dbg_assert( atom_lab[katom].first>0 );
861 : // Find base colvar
862 1811524 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
863 2644526 : if( usespecies && iatom==0 ) { myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[0] ); return; }
864 :
865 : // Get start of indices for this atom
866 646526 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=(mybasemulticolvars[i]->getNumberOfDerivatives() - 9) / 3;
867 1957044 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[katom].second, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
868 978522 : myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
869 : }
870 :
871 1122658 : void MultiColvarBase::getInputData( const unsigned& ind, const bool& normed,
872 : const multicolvar::AtomValuePack& myatoms,
873 : std::vector<double>& orient ) const {
874 : // Converint input atom to local index
875 : unsigned katom = myatoms.getIndex(ind); plumed_dbg_assert( atom_lab[katom].first>0 );
876 : // Find base colvar
877 2245316 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
878 : // Check if orient is the correct size
879 2245316 : if( orient.size()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ) orient.resize( mybasemulticolvars[mmc]->getNumberOfQuantities() );
880 : // Retrieve the value
881 2245316 : mybasedata[mmc]->retrieveValueWithIndex( atom_lab[katom].second, normed, orient );
882 1122658 : }
883 :
884 53063 : MultiValue& MultiColvarBase::getInputDerivatives( const unsigned& iatom, const bool& normed, const multicolvar::AtomValuePack& myatoms ) const {
885 : // Converint input atom to local index
886 : unsigned katom = myatoms.getIndex(iatom); plumed_dbg_assert( atom_lab[katom].first>0 );
887 : // Find base colvar
888 106126 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
889 54861 : if( usespecies && !normed && iatom==0 ) return mybasedata[mmc]->getTemporyMultiValue(0);
890 :
891 51265 : unsigned oval=0; if( iatom>0 ) oval=1;
892 102530 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(oval);
893 102491 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
894 51226 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
895 78 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
896 : }
897 102530 : mybasedata[mmc]->retrieveDerivatives( atom_lab[katom].second, normed, myder );
898 : return myder;
899 : }
900 :
901 60104530 : void MultiColvarBase::accumulateSymmetryFunction( const int& ival, const unsigned& iatom, const double& val, const Vector& der, const Tensor& vir, multicolvar::AtomValuePack& myatoms ) const {
902 : plumed_dbg_assert( usespecies ); unsigned katom=myatoms.getIndex(0), jatom=myatoms.getIndex(iatom);
903 120209854 : double weight0=1.0; if( atom_lab[katom].first>0 ) weight0=mybasedata[atom_lab[katom].first-1]->retrieveWeightWithIndex( atom_lab[katom].second );
904 120338542 : double weighti=1.0; if( atom_lab[jatom].first>0 ) weighti=mybasedata[atom_lab[jatom].first-1]->retrieveWeightWithIndex( atom_lab[jatom].second );
905 : // Accumulate the value
906 60168874 : if( ival<0 ) myatoms.getUnderlyingMultiValue().addTemporyValue( weight0*weighti*val );
907 57431983 : else myatoms.addValue( ival, weight0*weighti*val );
908 :
909 : // Return if we don't need derivatives
910 120372271 : if( doNotCalculateDerivatives() ) return ;
911 : // And virial
912 54045327 : if( ival<0 ) myatoms.addTemporyBoxDerivatives( weight0*weighti*vir );
913 51787180 : else myatoms.addBoxDerivatives( ival, weight0*weighti*vir );
914 :
915 : // Add derivatives of central atom
916 54045021 : if( atom_lab[katom].first>0 ) {
917 794 : addComDerivatives( ival, 0, -weight0*weighti*der, myatoms );
918 2382 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
919 794 : tmpder[0]=weighti*val; mergeInputDerivatives( ival, 0, 1, 0, tmpder, getInputDerivatives(0, false, myatoms), myatoms );
920 : } else {
921 54044227 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( 0, -der );
922 51786083 : else myatoms.addAtomsDerivatives( ival, 0, -der );
923 : }
924 : // Add derivatives of atom in coordination sphere
925 54045016 : if( atom_lab[jatom].first>0 ) {
926 794 : addComDerivatives( ival, iatom, weight0*weighti*der, myatoms );
927 2382 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
928 794 : tmpder[0]=weight0*val; mergeInputDerivatives( ival, 0, 1, iatom, tmpder, getInputDerivatives(iatom, false, myatoms), myatoms );
929 : } else {
930 54044222 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
931 51786078 : else myatoms.addAtomsDerivatives( ival, iatom, der );
932 : }
933 : }
934 :
935 2925334 : void MultiColvarBase::addAtomDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
936 5852968 : if( doNotCalculateDerivatives() ) return ;
937 : unsigned jatom=myatoms.getIndex(iatom);
938 :
939 5360230 : if( atom_lab[jatom].first>0 ) {
940 904174 : addComDerivatives( ival, iatom, der, myatoms );
941 : } else {
942 1775941 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
943 1775941 : else myatoms.addAtomsDerivatives( ival, iatom, der );
944 : }
945 : }
946 :
947 262596 : double MultiColvarBase::calculateWeight( const unsigned& current, const double& weight, AtomValuePack& myvals ) const {
948 262596 : return 1.0;
949 : }
950 :
951 467683 : void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
952 467683 : AtomValuePack myatoms( myvals, this );
953 : // Retrieve the atom list
954 467693 : if( !setupCurrentAtomList( current, myatoms ) ) return;
955 : // Get weight due to dynamic groups
956 467579 : double weight = 1.0;
957 467579 : if( !matsums ) {
958 772154973 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
959 772161296 : if( atom_lab[myatoms.getIndex(i)].first==0 ) continue;
960 : // Only need to do first two atoms for thigns like TopologyMatrix, HbondMatrix, Bridge and so on
961 407286 : if( allthirdblockintasks && i>1 ) break;
962 407286 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
963 814572 : weight *= mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
964 : }
965 147364 : } else if( usespecies ) {
966 293840 : if( atom_lab[myatoms.getIndex(0)].first>0 ) {
967 18146 : if( mybasedata[atom_lab[myatoms.getIndex(0)].first-1]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(0)].second )<epsilon ) weight=0.;
968 : }
969 : }
970 : // Do a quick check on the size of this contribution
971 304310 : double multweight = calculateWeight( current, weight, myatoms );
972 935146 : if( weight*multweight<getTolerance() ) {
973 121526 : updateActiveAtoms( myatoms );
974 : return;
975 : }
976 : myatoms.setValue( 0, weight*multweight );
977 : // Deal with derivatives of weights due to dynamic groups
978 486313 : if( !matsums && !doNotCalculateDerivatives() && mybasemulticolvars.size()>0 ) {
979 39416 : MultiValue& outder=myatoms.getUnderlyingMultiValue(); MultiValue myder(0,0);
980 230246 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
981 : // Neglect any atoms without differentiable weights
982 151414 : if( atom_lab[myatoms.getIndex(i)].first==0 ) continue;
983 :
984 : // Retrieve derivatives
985 75707 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
986 187705 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() || myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
987 78832 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
988 : }
989 302828 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(i)].second, false, myder );
990 :
991 : // Retrieve the prefactor (product of all other weights)
992 302828 : double prefactor = multweight*weight / mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
993 : // And accumulate the derivatives
994 8705716 : for(unsigned j=0; j<myder.getNumberActive(); ++j) { unsigned jder=myder.getActiveIndex(j); outder.addDerivative( 0, jder, prefactor*myder.getDerivative(0,jder) ); }
995 75707 : myder.clearAll();
996 39416 : }
997 : }
998 : // Retrieve derivative stuff for central atom
999 346056 : if( !doNotCalculateDerivatives() ) {
1000 318676 : if( usespecies && mybasemulticolvars.size()>0 && atom_lab[myatoms.getIndex(0)].first>0 ) {
1001 1610 : unsigned mmc = atom_lab[0].first - 1;
1002 3220 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(0);
1003 3208 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
1004 1598 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
1005 24 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
1006 : }
1007 6440 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(0)].second, false, myder );
1008 1610 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
1009 4830 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[myatoms.getIndex(0)].second, mybasemulticolvars[mmc]->my_tmp_capacks[0] );
1010 : }
1011 : }
1012 : // Compute everything
1013 346041 : double vv=compute( task_index, myatoms ); updateActiveAtoms( myatoms );
1014 : myatoms.setValue( 1, vv );
1015 346043 : return;
1016 : }
1017 :
1018 677866 : void MultiColvarBase::updateActiveAtoms( AtomValuePack& myatoms ) const {
1019 677866 : if( mybasemulticolvars.size()==0 ) myatoms.updateUsingIndices();
1020 141871 : else myatoms.updateDynamicList();
1021 677922 : }
1022 :
1023 2753917 : Vector MultiColvarBase::getCentralAtomPos( const unsigned& taskIndex ) {
1024 2753917 : unsigned curr=getTaskCode( taskIndex );
1025 :
1026 2753917 : if( usespecies || isDensity() ) {
1027 1459934 : return getPositionOfAtomForLinkCells(curr);
1028 1293989 : } else if( nblock>0 ) {
1029 : // double factor=1.0/static_cast<double>( ablocks.size() );
1030 2130 : Vector mypos; mypos.zero();
1031 2130 : std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
1032 10650 : for(unsigned i=0; i<ablocks.size(); ++i) {
1033 8520 : if( use_for_central_atom[i] ) mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(atoms[i]);
1034 : }
1035 2130 : return mypos;
1036 : } else {
1037 1291859 : Vector mypos; mypos.zero();
1038 10277264 : for(unsigned i=0; i<ablocks.size(); ++i) {
1039 5170210 : if( use_for_central_atom[i] ) mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(ablocks[i][curr]);
1040 : }
1041 1291859 : return mypos;
1042 : }
1043 : }
1044 :
1045 605381 : void MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsigned& taskIndex, CatomPack& mypack ) {
1046 605381 : unsigned curr=getTaskCode( taskIndex );
1047 :
1048 605381 : if(usespecies) {
1049 464995 : if( mypack.getNumberOfAtomsWithDerivatives()!=1 ) mypack.resize(1);
1050 464996 : mypack.setIndex( 0, basn + curr );
1051 464996 : mypack.setDerivative( 0, Tensor::identity() );
1052 140386 : } else if( nblock>0 ) {
1053 0 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral);
1054 0 : unsigned k=0;
1055 0 : std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
1056 0 : for(unsigned i=0; i<ablocks.size(); ++i) {
1057 0 : if( use_for_central_atom[i] ) {
1058 0 : mypack.setIndex( k, basn + atoms[i] );
1059 0 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
1060 0 : k++;
1061 : }
1062 : }
1063 : } else {
1064 140386 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral);
1065 140336 : unsigned k=0;
1066 633310 : for(unsigned i=0; i<ablocks.size(); ++i) {
1067 176246 : if( use_for_central_atom[i] ) {
1068 347208 : mypack.setIndex( k, basn + ablocks[i][curr] );
1069 173633 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
1070 173677 : k++;
1071 : }
1072 : }
1073 : }
1074 605404 : }
1075 :
1076 122146561 : Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 ) const {
1077 122146561 : if(usepbc) { return pbcDistance( vec1, vec2 ); }
1078 0 : else { return delta( vec1, vec2 ); }
1079 : }
1080 :
1081 240395 : void MultiColvarBase::applyPbc(std::vector<Vector>& dlist, unsigned int max_index) const {
1082 240395 : if (usepbc) pbcApply(dlist, max_index);
1083 240404 : }
1084 :
1085 1980 : void MultiColvarBase::apply() {
1086 1980 : if( getForcesFromVessels( forcesToApply ) ) setForcesOnAtoms( forcesToApply );
1087 1980 : }
1088 :
1089 : }
1090 5874 : }
|
