Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2020 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #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 415 : void MultiColvarBase::registerKeywords( Keywords& keys ) {
41 415 : Action::registerKeywords( keys );
42 415 : ActionWithValue::registerKeywords( keys );
43 415 : ActionAtomistic::registerKeywords( keys );
44 1245 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
45 415 : ActionWithVessel::registerKeywords( keys );
46 1660 : 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 830 : 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 : "input you can customize the component name");
61 1660 : 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 1660 : 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 : "using the label of another multicolvar");
71 1660 : 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 1660 : 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 415 : }
75 :
76 341 : 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 2046 : nblock(0)
91 : {
92 682 : if( keywords.exists("NOPBC") ) {
93 682 : bool nopbc=!usepbc; parseFlag("NOPBC",nopbc);
94 341 : usepbc=!nopbc;
95 : }
96 682 : if( keywords.exists("SPECIESA") ) { matsums=usespecies=true; }
97 341 : }
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 9428 : 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 9428 : 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 459 : bool MultiColvarBase::parseMultiColvarAtomList(const std::string& key, const int& num, std::vector<AtomNumber>& t) {
131 459 : std::vector<std::string> mlabs;
132 459 : if( num<0 ) parseVector(key,mlabs);
133 0 : else parseNumberedVector(key,num,mlabs);
134 :
135 459 : if( mlabs.size()==0 ) return false;
136 :
137 312 : std::string mname; unsigned found_mcolv=mlabs.size();
138 1005 : for(unsigned i=0; i<mlabs.size(); ++i) {
139 634 : MultiColvarBase* mycolv = plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlabs[i]);
140 317 : 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 17833498 : 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 312 : if( found_mcolv==0 ) {
160 190 : std::vector<AtomNumber> tt; ActionAtomistic::interpretAtomList( mlabs, tt );
161 351280 : for(unsigned i=0; i<tt.size(); ++i) { atom_lab.push_back( std::pair<unsigned,unsigned>( 0, t.size() + i ) ); }
162 380 : log.printf(" keyword %s takes atoms : ", key.c_str() );
163 351280 : for(unsigned i=0; i<tt.size(); ++i) { t.push_back( tt[i] ); log.printf("%d ",tt[i].serial() ); }
164 190 : 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 244 : log.printf(" keyword %s takes dynamic groups of atoms constructed from multicolvars labelled : ", key.c_str() );
168 625 : 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 : 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 16406 : for(unsigned i=0; i<nblock; ++i) ablocks[0][i]=ablocks[1][i]=i;
183 10894 : for(unsigned i=1; i<nblock; ++i) {
184 8427224 : 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 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;
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 151 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
201 3361 : 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 5545 : } 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 5904 : 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 4017 : for(unsigned i=0; i<nblock; ++i) { bookeeping(i,i).first=0; bookeeping(i,i).second=1; }
225 2668 : for(unsigned i=1; i<nblock; ++i) {
226 443563 : for(unsigned j=0; j<i; ++j) {
227 442232 : bookeeping(j,i).first=bookeeping(i,j).first=getFullNumberOfTasks();
228 221116 : addTaskToList( i*nblock + j );
229 442232 : bookeeping(j,i).second=bookeeping(i,j).second=getFullNumberOfTasks();
230 : }
231 : }
232 : } else {
233 272 : for(unsigned i=0; i<nblock; ++i) {
234 16554 : 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 199360 : 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 10 : 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 10 : plumed_dbg_assert( keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) ); ablocks.resize( 3 );
279 :
280 10 : bool readkey1=parseMultiColvarAtomList(key1,-1,all_atoms);
281 62 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=i;
282 10 : bool readkey2=parseMultiColvarAtomList(key2,-1,all_atoms);
283 10 : if( !readkey1 && !readkey2 ) return ;
284 665 : 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 20 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
287 10 : bool readkey3=parseMultiColvarAtomList(key3,-1,all_atoms);
288 10 : if( !readkey3 && !allow2 ) {
289 0 : error("missing atom specification " + key3);
290 10 : } 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 796 : for(unsigned i=0; i<ablocks[1].size(); ++i) ablocks[2][i]=ablocks[0].size() + i;
296 10 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
297 592 : for(unsigned j=1; j<ablocks[1].size(); ++j) {
298 196 : bookeeping(i,j).first=getFullNumberOfTasks();
299 19600 : 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 38808 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
307 48510 : 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 16 : ablocks[2].resize( atom_lab.size() - ablocks[1].size() - ablocks[0].size() );
315 15886 : for(unsigned i=0; i<ablocks[2].size(); ++i) ablocks[2][i] = ablocks[0].size() + ablocks[1].size() + i;
316 :
317 8 : if( ablocks[1].size()>ablocks[0].size() ) nblock=ablocks[1].size();
318 8 : else nblock=ablocks[0].size();
319 8 : if( ablocks[2].size()>nblock ) nblock=ablocks[2].size();
320 :
321 8 : unsigned kcount; if( no_third_dim_accum ) kcount=1; else kcount=ablocks[2].size();
322 :
323 73 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
324 365 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
325 109 : bookeeping(i,j).first=getFullNumberOfTasks();
326 327 : for(unsigned k=0; k<kcount; ++k) {
327 109 : 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 109 : 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 29 : 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 29 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
357 7 : ablocks[i].resize( nblock );
358 47 : 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 62 : 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 12511070 : void MultiColvarBase::addTaskToList( const unsigned& taskCode ) {
373 12511070 : plumed_assert( getNumberOfVessels()==0 );
374 12511070 : ActionWithVessel::addTaskToList( taskCode );
375 12511070 : }
376 :
377 96 : void MultiColvarBase::resizeBookeepingArray( const unsigned& num1, const unsigned& num2 ) {
378 96 : bookeeping.resize( num1, num2 );
379 14034 : for(unsigned i=0; i<num1; ++i) {
380 26648304 : for(unsigned j=0; j<num2; ++j) { bookeeping(i,j).first=0; bookeeping(i,j).second=0; }
381 : }
382 96 : }
383 :
384 299 : 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 299 : if( !usespecies && nblock>0 ) {
389 :
390 63 : ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
391 63 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
392 63 : if( ablocks.size()==3 ) {
393 20 : allthirdblockintasks=uselinkforthree=true;
394 3004 : for(unsigned i=0; i<bookeeping.nrows(); ++i) {
395 905664 : for(unsigned j=0; j<bookeeping.ncols(); ++j) {
396 452086 : unsigned ntper = bookeeping(i,j).second - bookeeping(i,j).first;
397 452086 : if( i==j && ntper==0 ) {
398 136 : continue;
399 451950 : } else if( ntper == 1 && allthirdblockintasks ) {
400 451756 : 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 63 : if( allthirdblockintasks ) {
411 36 : decoder.resize(2); plumed_assert( ablocks.size()==3 );
412 : // Check if number of atoms is too large
413 18 : 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 507 : unsigned code=1; for(unsigned i=0; i<decoder.size(); ++i) { decoder[decoder.size()-1-i]=code; code *= nblock; }
420 236 : } 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 298 : } else if( keywords.exists("SPECIESA") ) {
424 182 : plumed_assert( atom_lab.size()==0 && all_atoms.size()==0 );
425 273 : ablocks.resize( 1 ); bool readspecies=parseMultiColvarAtomList("SPECIES", -1, all_atoms);
426 91 : 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 38 : if( !parseMultiColvarAtomList("SPECIESA", -1, all_atoms) ) error("missing SPECIES/SPECIESA keyword");
430 19 : unsigned nat1=atom_lab.size();
431 38 : if( !parseMultiColvarAtomList("SPECIESB", -1, all_atoms) ) error("missing SPECIESB keyword");
432 19 : unsigned nat2=atom_lab.size() - nat1;
433 :
434 19 : for(unsigned i=0; i<nat1; ++i) addTaskToList(i);
435 19 : ablocks[0].resize( nat2 );
436 6855 : for(unsigned i=0; i<nat2; ++i) {
437 : bool found=false; unsigned inum;
438 527480 : for(unsigned j=0; j<nat1; ++j) {
439 776946 : if( atom_lab[nat1+i].first>0 && atom_lab[j].first>0 ) {
440 252720 : 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 252720 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) { found=true; inum=j; break; }
443 9393 : } 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 18786 : } 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 28179 : } 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 164 : if( found ) { ablocks[0][i]=inum; }
453 6672 : else { ablocks[0][i]=nat1 + i; }
454 : }
455 : }
456 : }
457 299 : if( mybasemulticolvars.size()>0 ) {
458 619 : for(unsigned i=0; i<mybasedata.size(); ++i) {
459 254 : 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 979 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
466 127 : tmp_atoms=mybasemulticolvars[i]->getAbsoluteIndexes();
467 1248821 : for(unsigned j=0; j<tmp_atoms.size(); ++j) all_atoms.push_back( tmp_atoms[j] );
468 : }
469 : // Copy atom lists from input
470 178399 : 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 299 : ActionAtomistic::requestAtoms( all_atoms );
474 : // And setup dependencies
475 979 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) addDependency( mybasemulticolvars[i] );
476 :
477 : // Setup underlying ActionWithVessel
478 299 : readVesselKeywords();
479 299 : }
480 :
481 19 : void MultiColvarBase::setAtomsForCentralAtom( const std::vector<bool>& catom_ind ) {
482 19 : unsigned nat=0; plumed_assert( catom_ind.size()==ablocks.size() );
483 248 : for(unsigned i=0; i<catom_ind.size(); ++i) {
484 : use_for_central_atom[i]=catom_ind[i];
485 70 : if( use_for_central_atom[i] ) nat++;
486 : }
487 19 : plumed_dbg_assert( nat>0 ); ncentral=nat;
488 19 : numberForCentralAtom = 1.0 / static_cast<double>( nat );
489 19 : }
490 :
491 426 : void MultiColvarBase::turnOnDerivatives() {
492 426 : ActionWithValue::turnOnDerivatives();
493 426 : needsDerivatives();
494 426 : forcesToApply.resize( getNumberOfDerivatives() );
495 426 : }
496 :
497 182 : void MultiColvarBase::setLinkCellCutoff( const double& lcut, double tcut ) {
498 271 : plumed_assert( usespecies || ablocks.size()<4 );
499 358 : if( tcut<0 ) tcut=lcut;
500 :
501 182 : if( !linkcells.enabled() ) {
502 182 : linkcells.setCutoff( lcut );
503 182 : 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 182 : }
509 :
510 8 : double MultiColvarBase::getLinkCellCutoff() const {
511 8 : return linkcells.getCutoff();
512 : }
513 :
514 1927 : void MultiColvarBase::setupLinkCells() {
515 3063 : if( (!usespecies && nblock==0) || !linkcells.enabled() ) return ;
516 : // Retrieve any atoms that haven't already been retrieved
517 1307 : for(std::vector<MultiColvarBase*>::iterator p=mybasemulticolvars.begin(); p!=mybasemulticolvars.end(); ++p) {
518 202 : (*p)->retrieveAtoms();
519 : }
520 1105 : retrieveAtoms();
521 :
522 : unsigned iblock;
523 1105 : if( usespecies ) {
524 : iblock=0;
525 217 : } else if( ablocks.size()<4 ) {
526 : iblock=1;
527 : } else {
528 0 : plumed_error();
529 : }
530 :
531 : // Count number of currently active atoms
532 1105 : nactive_atoms=0;
533 1157237 : for(unsigned i=0; i<ablocks[iblock].size(); ++i) {
534 770018 : if( isCurrentlyActive( ablocks[iblock][i] ) ) nactive_atoms++;
535 : }
536 :
537 1105 : if( nactive_atoms>0 ) {
538 1105 : std::vector<Vector> ltmp_pos( nactive_atoms );
539 1105 : std::vector<unsigned> ltmp_ind( nactive_atoms );
540 :
541 1105 : nactive_atoms=0;
542 1105 : if( usespecies ) {
543 1093509 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
544 727822 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
545 727630 : ltmp_ind[nactive_atoms]=ablocks[0][i];
546 727630 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ltmp_ind[nactive_atoms] );
547 363815 : nactive_atoms++;
548 : }
549 : } else {
550 63728 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
551 42196 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
552 32356 : ltmp_ind[nactive_atoms]=i;
553 48534 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ablocks[1][i] );
554 16178 : nactive_atoms++;
555 : }
556 : }
557 :
558 : // Build the lists for the link cells
559 2210 : linkcells.buildCellLists( ltmp_pos, ltmp_ind, getPbc() );
560 : }
561 : }
562 :
563 4628 : void MultiColvarBase::setupNonUseSpeciesLinkCells( const unsigned& my_always_active ) {
564 4628 : plumed_assert( !usespecies );
565 5229 : if( nblock==0 || !linkcells.enabled() ) return ;
566 210 : deactivateAllTasks();
567 : std::vector<unsigned> requiredlinkcells;
568 :
569 210 : 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 17115 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
578 11202 : 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 408788 : 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 54 : } else if( nactive_atoms>0 ) {
586 : // Get some parallel info
587 54 : unsigned stride=comm.Get_size();
588 54 : unsigned rank=comm.Get_rank();
589 54 : if( serialCalculation() ) { stride=1; rank=0; }
590 :
591 : unsigned nactive_three=0;
592 199743 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
593 133090 : if( isCurrentlyActive( ablocks[2][i] ) ) nactive_three++;
594 : }
595 :
596 54 : std::vector<Vector> lttmp_pos( nactive_three );
597 54 : std::vector<unsigned> lttmp_ind( nactive_three );
598 :
599 : nactive_three=0;
600 54 : if( allthirdblockintasks ) {
601 199743 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
602 133090 : if( !isCurrentlyActive( ablocks[2][i] ) ) continue;
603 133090 : lttmp_ind[nactive_three]=ablocks[2][i];
604 133090 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
605 66545 : 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 108 : threecells.buildCellLists( lttmp_pos, lttmp_ind, getPbc() );
617 :
618 : // Ensure we only do tasks where atoms are in appropriate link cells
619 54 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
620 54 : std::vector<unsigned> tlinked_atoms( 1+ablocks[2].size() );
621 1845 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
622 1158 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
623 1158 : unsigned natomsper=1; linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
624 1158 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, natomsper, linked_atoms );
625 579 : if( allthirdblockintasks ) {
626 143085 : for(unsigned j=0; j<natomsper; ++j) {
627 355997 : 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 210 : if( !serialCalculation() ) comm.Sum( taskFlags );
639 210 : lockContributors();
640 : }
641 :
642 245988 : void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const {
643 : plumed_dbg_assert( !usespecies && nblock>0 );
644 245988 : if( atoms.size()!=decoder.size() ) atoms.resize( decoder.size() );
645 :
646 245988 : unsigned scode = taskCode;
647 2113404 : for(unsigned i=0; i<decoder.size(); ++i) {
648 540476 : unsigned ind=( scode / decoder[i] );
649 1080952 : atoms[i] = ablocks[i][ind];
650 540476 : scode -= ind*decoder[i];
651 : }
652 245988 : }
653 :
654 421349 : bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ) const {
655 421349 : if( isDensity() ) {
656 19065 : myatoms.setNumberOfAtoms( 1 ); myatoms.setAtom( 0, taskCode ); return true;
657 402284 : } else if( usespecies ) {
658 301396 : std::vector<unsigned> task_atoms(1); task_atoms[0]=taskCode;
659 150698 : unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getPositionOfAtomForLinkCells( taskCode ), linkcells );
660 150698 : return natomsper>1;
661 251586 : } else if( matsums ) {
662 : myatoms.setNumberOfAtoms( getNumberOfAtoms() );
663 53239 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) myatoms.setAtom( i, i );
664 251145 : } else if( allthirdblockintasks ) {
665 39816 : plumed_dbg_assert( ablocks.size()==3 ); std::vector<unsigned> atoms(2); decodeIndexToAtoms( taskCode, atoms );
666 79632 : myatoms.setupAtomsFromLinkCells( atoms, getPositionOfAtomForLinkCells( atoms[0] ), threecells );
667 211329 : } else if( nblock>0 ) {
668 179551 : std::vector<unsigned> atoms( ablocks.size() );
669 359102 : decodeIndexToAtoms( taskCode, atoms ); myatoms.setNumberOfAtoms( ablocks.size() );
670 1174306 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, atoms[i] );
671 : } else {
672 31778 : myatoms.setNumberOfAtoms( ablocks.size() );
673 248996 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, ablocks[i][taskCode] );
674 : }
675 : return true;
676 : }
677 :
678 8765 : void MultiColvarBase::setupActiveTaskSet( std::vector<unsigned>& active_tasks, const std::string& input_label ) {
679 8765 : if( !setup_completed ) {
680 : bool justVolumes=false;
681 1919 : if( usespecies ) {
682 : justVolumes=true;
683 1967 : for(unsigned i=0; i<getNumberOfVessels(); ++i) {
684 1307 : vesselbase::StoreDataVessel* mys=dynamic_cast<vesselbase::StoreDataVessel*>( getPntrToVessel(i) );
685 1307 : if( mys ) continue;
686 799 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(i) );
687 799 : if( !myb ) { justVolumes=false; break; }
688 38 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
689 38 : if( !myv ) { justVolumes=false; break; }
690 : }
691 : }
692 1919 : deactivateAllTasks();
693 1919 : if( justVolumes && mydata ) {
694 88 : if( mydata->getNumberOfDataUsers()==0 ) justVolumes=false;
695 :
696 323 : for(unsigned i=0; i<mydata->getNumberOfDataUsers(); ++i) {
697 49 : MultiColvarBase* myu=dynamic_cast<MultiColvarBase*>( mydata->getDataUser(i) );
698 49 : if( myu ) {
699 98 : myu->setupActiveTaskSet( taskFlags, getLabel() );
700 : } else {
701 0 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
702 : }
703 : }
704 : }
705 1919 : 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 9768257 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
719 : }
720 :
721 : // Now activate all this class
722 1919 : lockContributors();
723 : // Setup the link cells
724 1919 : setupLinkCells();
725 : // Ensures that setup is not performed multiple times during one cycle
726 1919 : setup_completed=true;
727 : }
728 :
729 : // And activate the tasks in input action
730 17530 : if( getLabel()!=input_label ) {
731 : int input_code=-1;
732 134 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
733 122 : if( mybasemulticolvars[i]->getLabel()==input_label ) { input_code=i+1; break; }
734 : }
735 :
736 98 : 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 506063 : 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 : }
746 : }
747 8765 : }
748 :
749 465 : bool MultiColvarBase::filtersUsedAsInput() {
750 : bool inputAreFilters=false;
751 1710 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
752 260 : MultiColvarFilter* myfilt=dynamic_cast<MultiColvarFilter*>( mybasemulticolvars[i] );
753 260 : if( myfilt || mybasemulticolvars[i]->filtersUsedAsInput() ) inputAreFilters=true;
754 : }
755 465 : return inputAreFilters;
756 : }
757 :
758 8716 : void MultiColvarBase::calculate() {
759 : // Recursive function that sets up tasks
760 17432 : setupActiveTaskSet( taskFlags, getLabel() );
761 :
762 : // Check for filters and rerun setup of link cells if there are any
763 8716 : if( mybasemulticolvars.size()>0 && filtersUsedAsInput() ) setupLinkCells();
764 :
765 : // Setup the link cells if we are not using species
766 15360 : 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 4628 : unsigned first_active=std::numeric_limits<unsigned>::max();
770 9658 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
771 9524 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
772 : else {
773 4628 : first_active=i; break;
774 : }
775 : }
776 4628 : setupNonUseSpeciesLinkCells( first_active );
777 : }
778 : // And run all tasks
779 8716 : runAllTasks();
780 8716 : }
781 :
782 212 : void MultiColvarBase::calculateNumericalDerivatives( ActionWithValue* a ) {
783 212 : if( mybasemulticolvars.size()>0 ) plumed_merror("cannot calculate numerical derivatives for this quantity");
784 212 : calculateAtomicNumericalDerivatives( this, 0 );
785 212 : }
786 :
787 2948 : void MultiColvarBase::prepare() {
788 2948 : setup_completed=false; atomsWereRetrieved=false;
789 2948 : }
790 :
791 6526 : void MultiColvarBase::retrieveAtoms() {
792 6526 : if( !atomsWereRetrieved ) { ActionAtomistic::retrieveAtoms(); atomsWereRetrieved=true; }
793 6526 : }
794 :
795 38205 : 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 76410 : 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 39425 : 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 38205 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
810 9573417 : 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 38205 : }
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 202 : 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 38962 : 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 905722 : 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 1811444 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
863 1738684 : 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 646506 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=(mybasemulticolvars[i]->getNumberOfDerivatives() - 9) / 3;
867 1467723 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[katom].second, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
868 978482 : 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 53023 : 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 106046 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
889 54801 : if( usespecies && !normed && iatom==0 ) return mybasedata[mmc]->getTemporyMultiValue(0);
890 :
891 51245 : unsigned oval=0; if( iatom>0 ) oval=1;
892 102490 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(oval);
893 102451 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
894 51206 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
895 78 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
896 : }
897 102490 : mybasedata[mmc]->retrieveDerivatives( atom_lab[katom].second, normed, myder );
898 : return myder;
899 : }
900 :
901 60217341 : 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 120435456 : double weight0=1.0; if( atom_lab[katom].first>0 ) weight0=mybasedata[atom_lab[katom].first-1]->retrieveWeightWithIndex( atom_lab[katom].second );
904 120435456 : 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 62954241 : if( ival<0 ) myatoms.getUnderlyingMultiValue().addTemporyValue( weight0*weighti*val );
907 57480441 : else myatoms.addValue( ival, weight0*weighti*val );
908 :
909 : // Return if we don't need derivatives
910 60217341 : if( doNotCalculateDerivatives() ) return ;
911 : // And virial
912 54045394 : if( ival<0 ) myatoms.addTemporyBoxDerivatives( weight0*weighti*vir );
913 51787247 : else myatoms.addBoxDerivatives( ival, weight0*weighti*vir );
914 :
915 : // Add derivatives of central atom
916 54045394 : if( atom_lab[katom].first>0 ) {
917 774 : addComDerivatives( ival, 0, -weight0*weighti*der, myatoms );
918 2322 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
919 1548 : tmpder[0]=weighti*val; mergeInputDerivatives( ival, 0, 1, 0, tmpder, getInputDerivatives(0, false, myatoms), myatoms );
920 : } else {
921 54044620 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( 0, -der );
922 51786473 : else myatoms.addAtomsDerivatives( ival, 0, -der );
923 : }
924 : // Add derivatives of atom in coordination sphere
925 54045394 : if( atom_lab[jatom].first>0 ) {
926 774 : addComDerivatives( ival, iatom, weight0*weighti*der, myatoms );
927 2322 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
928 1548 : tmpder[0]=weight0*val; mergeInputDerivatives( ival, 0, 1, iatom, tmpder, getInputDerivatives(iatom, false, myatoms), myatoms );
929 : } else {
930 54044620 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
931 51786473 : else myatoms.addAtomsDerivatives( ival, iatom, der );
932 : }
933 : }
934 :
935 1428471 : void MultiColvarBase::addAtomDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
936 1428471 : if( doNotCalculateDerivatives() ) return ;
937 : unsigned jatom=myatoms.getIndex(iatom);
938 :
939 2361790 : if( atom_lab[jatom].first>0 ) {
940 904174 : addComDerivatives( ival, iatom, der, myatoms );
941 : } else {
942 276721 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
943 276721 : else myatoms.addAtomsDerivatives( ival, iatom, der );
944 : }
945 : }
946 :
947 212706 : double MultiColvarBase::calculateWeight( const unsigned& current, const double& weight, AtomValuePack& myvals ) const {
948 212706 : return 1.0;
949 : }
950 :
951 417820 : void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
952 417820 : AtomValuePack myatoms( myvals, this );
953 : // Retrieve the atom list
954 417820 : if( !setupCurrentAtomList( current, myatoms ) ) return;
955 : // Get weight due to dynamic groups
956 417700 : double weight = 1.0;
957 417700 : if( !matsums ) {
958 772150301 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
959 771880062 : 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 244017 : if( allthirdblockintasks && i>1 ) break;
962 244017 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
963 488034 : weight *= mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
964 : }
965 147461 : } else if( usespecies ) {
966 294040 : if( atom_lab[myatoms.getIndex(0)].first>0 ) {
967 18136 : 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 417700 : double multweight = calculateWeight( current, weight, myatoms );
972 835400 : if( weight*multweight<getTolerance() ) {
973 121536 : updateActiveAtoms( myatoms );
974 : return;
975 : }
976 : myatoms.setValue( 0, weight*multweight );
977 : // Deal with derivatives of weights due to dynamic groups
978 386430 : if( !matsums && !doNotCalculateDerivatives() && mybasemulticolvars.size()>0 ) {
979 78832 : MultiValue& outder=myatoms.getUnderlyingMultiValue(); MultiValue myder(0,0);
980 190830 : 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 227121 : 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 8630009 : 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 : }
997 : }
998 : // Retrieve derivative stuff for central atom
999 296164 : if( !doNotCalculateDerivatives() ) {
1000 268667 : if( usespecies && mybasemulticolvars.size()>0 && atom_lab[myatoms.getIndex(0)].first>0 ) {
1001 1605 : unsigned mmc = atom_lab[0].first - 1;
1002 3210 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(0);
1003 3198 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
1004 1593 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
1005 24 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
1006 : }
1007 4815 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(0)].second, false, myder );
1008 1605 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
1009 4815 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[myatoms.getIndex(0)].second, mybasemulticolvars[mmc]->my_tmp_capacks[0] );
1010 : }
1011 : }
1012 : // Compute everything
1013 296164 : double vv=compute( task_index, myatoms ); updateActiveAtoms( myatoms );
1014 : myatoms.setValue( 1, vv );
1015 296164 : return;
1016 : }
1017 :
1018 628082 : void MultiColvarBase::updateActiveAtoms( AtomValuePack& myatoms ) const {
1019 628082 : if( mybasemulticolvars.size()==0 ) myatoms.updateUsingIndices();
1020 141866 : else myatoms.updateDynamicList();
1021 628082 : }
1022 :
1023 2753900 : Vector MultiColvarBase::getCentralAtomPos( const unsigned& taskIndex ) {
1024 2753900 : unsigned curr=getTaskCode( taskIndex );
1025 :
1026 2753900 : if( usespecies || isDensity() ) {
1027 1459911 : 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 17040 : 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 14124079 : for(unsigned i=0; i<ablocks.size(); ++i) {
1039 5170238 : if( use_for_central_atom[i] ) mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(ablocks[i][curr]);
1040 : }
1041 1291859 : return mypos;
1042 : }
1043 : }
1044 :
1045 605385 : void MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsigned& taskIndex, CatomPack& mypack ) {
1046 605385 : unsigned curr=getTaskCode( taskIndex );
1047 :
1048 605385 : if(usespecies) {
1049 464996 : if( mypack.getNumberOfAtomsWithDerivatives()!=1 ) mypack.resize(1);
1050 464996 : mypack.setIndex( 0, basn + curr );
1051 929992 : mypack.setDerivative( 0, Tensor::identity() );
1052 140389 : } else if( nblock>0 ) {
1053 0 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral);
1054 : 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 140389 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral);
1065 : unsigned k=0;
1066 809705 : for(unsigned i=0; i<ablocks.size(); ++i) {
1067 176309 : if( use_for_central_atom[i] ) {
1068 347338 : mypack.setIndex( k, basn + ablocks[i][curr] );
1069 347338 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
1070 173669 : k++;
1071 : }
1072 : }
1073 : }
1074 605385 : }
1075 :
1076 121150993 : Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 ) const {
1077 121150993 : if(usepbc) { return pbcDistance( vec1, vec2 ); }
1078 0 : else { return delta( vec1, vec2 ); }
1079 : }
1080 :
1081 190514 : void MultiColvarBase::applyPbc(std::vector<Vector>& dlist, unsigned int max_index) const {
1082 190514 : if (usepbc) pbcApply(dlist, max_index);
1083 190514 : }
1084 :
1085 1983 : void MultiColvarBase::apply() {
1086 1983 : if( getForcesFromVessels( forcesToApply ) ) setForcesOnAtoms( forcesToApply );
1087 1983 : }
1088 :
1089 : }
1090 5517 : }
|
