Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "MultiColvarBase.h"
23 : #include "ActionVolume.h"
24 : #include "MultiColvarFilter.h"
25 : #include "vesselbase/Vessel.h"
26 : #include "vesselbase/BridgeVessel.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/ActionSet.h"
29 : #include "tools/Pbc.h"
30 : #include "AtomValuePack.h"
31 : #include <vector>
32 : #include <string>
33 : #include <limits>
34 :
35 : namespace PLMD {
36 : namespace multicolvar {
37 :
38 425 : void MultiColvarBase::registerKeywords( Keywords& keys ) {
39 425 : Action::registerKeywords( keys );
40 425 : ActionWithValue::registerKeywords( keys );
41 425 : ActionAtomistic::registerKeywords( keys );
42 850 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
43 425 : ActionWithVessel::registerKeywords( keys );
44 850 : keys.add("hidden","NL_STRIDE","the frequency with which the neighbor list should be updated. Between neighbour list update steps all quantities "
45 : "that contributed less than TOL at the previous neighbor list update step are ignored.");
46 425 : keys.setComponentsIntroduction("When the label of this action is used as the input for a second you are not referring to a scalar quantity as you are in "
47 : "regular collective variables. The label is used to reference the full set of quantities calculated by "
48 : "the action. This is usual when using \\ref multicolvarfunction. Generally when doing this the previously calculated "
49 : "multicolvar will be referenced using the DATA keyword rather than ARG.\n\n"
50 : "This Action can be used to calculate the following scalar quantities directly. These quantities are calculated by "
51 : "employing the keywords listed below. "
52 : "These quantities can then be referenced elsewhere in the input file by using this Action's label "
53 : "followed by a dot and the name of the quantity. Some of them can be calculated multiple times "
54 : "with different parameters. In this case the quantities calculated can be referenced elsewhere in the "
55 : "input by using the name of the quantity followed by a numerical identifier "
56 : "e.g. <em>label</em>.lessthan-1, <em>label</em>.lessthan-2 etc. When doing this and, for clarity we have "
57 : "made it so that the user can set a particular label for each of the components. As such by using the LABEL keyword in the description of the keyword "
58 : "input you can customize the component name");
59 850 : keys.reserve("atoms-3","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
60 : "one coordination number for each of the atoms specified. Each of these coordination numbers specifies how many of the "
61 : "other specified atoms are within a certain cutoff of the central atom. You can specify the atoms here as another multicolvar "
62 : "action or using a MultiColvarFilter or ActionVolume action. When you do so the quantity is calculated for those atoms specified "
63 : "in the previous multicolvar. This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
64 : "coordination number more than four for example");
65 850 : keys.reserve("atoms-4","SPECIESA","this keyword is used for colvars such as the coordination number. In that context it species that plumed should calculate "
66 : "one coordination number for each of the atoms specified in SPECIESA. Each of these coordination numbers specifies how many "
67 : "of the atoms specifies using SPECIESB is within the specified cutoff. As with the species keyword the input can also be specified "
68 : "using the label of another multicolvar");
69 850 : keys.reserve("atoms-4","SPECIESB","this keyword is used for colvars such as the coordination number. It must appear with SPECIESA. For a full explanation see "
70 : "the documentation for that keyword");
71 1275 : keys.add("hidden","ALL_INPUT_SAME_TYPE","remove this keyword to remove certain checks in the input on the sanity of your input file. See code for details");
72 425 : }
73 :
74 350 : MultiColvarBase::MultiColvarBase(const ActionOptions&ao):
75 : Action(ao),
76 : ActionAtomistic(ao),
77 : ActionWithValue(ao),
78 : ActionWithVessel(ao),
79 350 : usepbc(false),
80 350 : allthirdblockintasks(false),
81 350 : uselinkforthree(false),
82 350 : linkcells(comm),
83 350 : threecells(comm),
84 350 : setup_completed(false),
85 350 : atomsWereRetrieved(false),
86 350 : matsums(false),
87 350 : usespecies(false),
88 700 : nblock(0)
89 : {
90 700 : if( keywords.exists("NOPBC") ) {
91 350 : bool nopbc=!usepbc; parseFlag("NOPBC",nopbc);
92 350 : usepbc=!nopbc;
93 : }
94 700 : if( keywords.exists("SPECIESA") ) { matsums=usespecies=true; }
95 350 : }
96 :
97 48 : void MultiColvarBase::readAtomsLikeKeyword( const std::string & key, const int& natoms, std::vector<AtomNumber>& all_atoms ) {
98 48 : plumed_assert( !usespecies );
99 48 : if( all_atoms.size()>0 ) return;
100 :
101 : std::vector<AtomNumber> t;
102 48 : for(int i=1;; ++i ) {
103 970 : parseAtomList(key, i, t );
104 970 : if( t.empty() ) break;
105 :
106 922 : log.printf(" Colvar %d is calculated from atoms : ", i);
107 3450 : for(unsigned j=0; j<t.size(); ++j) log.printf("%d ",t[j].serial() );
108 922 : log.printf("\n");
109 :
110 922 : if( i==1 && natoms<0 ) { ablocks.resize(t.size()); }
111 916 : else if( i==1 ) ablocks.resize(natoms);
112 922 : if( t.size()!=ablocks.size() ) {
113 0 : std::string ss; Tools::convert(i,ss);
114 0 : error(key + ss + " keyword has the wrong number of atoms");
115 : }
116 3450 : for(unsigned j=0; j<ablocks.size(); ++j) {
117 2528 : ablocks[j].push_back( ablocks.size()*(i-1)+j ); all_atoms.push_back( t[j] );
118 2528 : atom_lab.push_back( std::pair<unsigned,unsigned>( 0, ablocks.size()*(i-1)+j ) );
119 : }
120 922 : t.resize(0);
121 922 : }
122 48 : if( all_atoms.size()>0 ) {
123 48 : nblock=0;
124 970 : for(unsigned i=0; i<ablocks[0].size(); ++i) addTaskToList( i );
125 : }
126 : }
127 :
128 468 : bool MultiColvarBase::parseMultiColvarAtomList(const std::string& key, const int& num, std::vector<AtomNumber>& t) {
129 : std::vector<std::string> mlabs;
130 468 : if( num<0 ) parseVector(key,mlabs);
131 0 : else parseNumberedVector(key,num,mlabs);
132 :
133 468 : if( mlabs.size()==0 ) return false;
134 :
135 321 : std::string mname; unsigned found_mcolv=mlabs.size();
136 448 : for(unsigned i=0; i<mlabs.size(); ++i) {
137 326 : MultiColvarBase* mycolv = plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlabs[i]);
138 326 : if(!mycolv) { found_mcolv=i; break; }
139 : // Check all base multicolvars are of same type
140 127 : if( i==0 ) {
141 122 : mname = mycolv->getName();
142 244 : if( keywords.exists("ALL_INPUT_SAME_TYPE") && mycolv->isPeriodic() ) error("multicolvar functions don't work with this multicolvar");
143 : } else {
144 10 : if( keywords.exists("ALL_INPUT_SAME_TYPE") && mname!=mycolv->getName() ) error("All input multicolvars must be of same type");
145 : }
146 : // And track which variable stores each colvar
147 4458438 : for(unsigned j=0; j<mycolv->getFullNumberOfTasks(); ++j) atom_lab.push_back( std::pair<unsigned,unsigned>( mybasemulticolvars.size()+1, j ) );
148 : // And store the multicolvar base
149 127 : mybasemulticolvars.push_back( mycolv );
150 : // And create a basedata stash
151 127 : mybasedata.push_back( mybasemulticolvars[mybasemulticolvars.size()-1]->buildDataStashes( this ) );
152 : // Check if weight has derivatives
153 127 : if( mybasemulticolvars[ mybasemulticolvars.size()-1 ]->weightHasDerivatives ) weightHasDerivatives=true;
154 127 : plumed_assert( mybasemulticolvars.size()==mybasedata.size() );
155 : }
156 : // Have we conventional atoms to read in
157 321 : if( found_mcolv==0 ) {
158 199 : std::vector<AtomNumber> tt; ActionAtomistic::interpretAtomList( mlabs, tt );
159 89460 : for(unsigned i=0; i<tt.size(); ++i) { atom_lab.push_back( std::pair<unsigned,unsigned>( 0, t.size() + i ) ); }
160 199 : log.printf(" keyword %s takes atoms : ", key.c_str() );
161 89460 : for(unsigned i=0; i<tt.size(); ++i) { t.push_back( tt[i] ); log.printf("%d ",tt[i].serial() ); }
162 199 : log.printf("\n");
163 122 : } else if( found_mcolv==mlabs.size() ) {
164 122 : if( checkNumericalDerivatives() ) error("cannot use NUMERICAL_DERIVATIVES keyword with dynamic groups as input");
165 122 : log.printf(" keyword %s takes dynamic groups of atoms constructed from multicolvars labelled : ", key.c_str() );
166 249 : for(unsigned i=0; i<mlabs.size(); ++i) log.printf("%s ",mlabs[i].c_str() );
167 122 : log.printf("\n");
168 0 : } else if( found_mcolv<mlabs.size() ) {
169 0 : error("cannot mix multicolvar input and atom input in one line");
170 : }
171 : return true;
172 468 : }
173 :
174 76 : void MultiColvarBase::readTwoGroups( const std::string& key0, const std::string& key1, const std::string& key2, std::vector<AtomNumber>& all_atoms ) {
175 76 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) ); ablocks.resize( 2 );
176 :
177 76 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
178 78 : nblock=atom_lab.size(); for(unsigned i=0; i<2; ++i) ablocks[i].resize(nblock);
179 26 : resizeBookeepingArray( nblock, nblock );
180 5486 : for(unsigned i=0; i<nblock; ++i) ablocks[0][i]=ablocks[1][i]=i;
181 5460 : for(unsigned i=1; i<nblock; ++i) {
182 4216329 : for(unsigned j=0; j<i; ++j) {
183 4210895 : bookeeping(i,j).first=getFullNumberOfTasks();
184 4210895 : addTaskToList( i*nblock + j );
185 4210895 : bookeeping(i,j).second=getFullNumberOfTasks();
186 : }
187 : }
188 : } else {
189 50 : parseMultiColvarAtomList(key1,-1,all_atoms);
190 67 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=i;
191 50 : parseMultiColvarAtomList(key2,-1,all_atoms);
192 1119 : 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;
193 :
194 50 : if( ablocks[0].size()>ablocks[1].size() ) nblock = ablocks[0].size();
195 50 : else nblock=ablocks[1].size();
196 :
197 50 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
198 67 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
199 1126 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
200 1109 : bookeeping(i,j).first=getFullNumberOfTasks();
201 1109 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 ) {
202 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
203 0 : atom_lab[ablocks[0][i]].second!=atom_lab[ablocks[1][j]].second ) addTaskToList( i*nblock + j );
204 1109 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] ) addTaskToList( i*nblock + j );
205 1109 : bookeeping(i,j).second=getFullNumberOfTasks();
206 : }
207 : }
208 : }
209 76 : }
210 :
211 12 : void MultiColvarBase::readGroupKeywords( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& key3,
212 : const bool& no_third_dim_accum, const bool& symmetric, std::vector<AtomNumber>& all_atoms ) {
213 12 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) ); ablocks.resize( 3 );
214 :
215 12 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
216 10 : if( no_third_dim_accum ) {
217 10 : nblock=atom_lab.size(); ablocks[0].resize(nblock); ablocks[1].resize( nblock );
218 1481 : for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=ablocks[1][i]=i;
219 10 : resizeBookeepingArray( nblock, nblock );
220 10 : if( symmetric ) {
221 : // This ensures that later parts of the code don't switch off allthirdblockintasks
222 1343 : for(unsigned i=0; i<nblock; ++i) { bookeeping(i,i).first=0; bookeeping(i,i).second=1; }
223 1337 : for(unsigned i=1; i<nblock; ++i) {
224 222447 : for(unsigned j=0; j<i; ++j) {
225 221116 : bookeeping(j,i).first=bookeeping(i,j).first=getFullNumberOfTasks();
226 221116 : addTaskToList( i*nblock + j );
227 221116 : bookeeping(j,i).second=bookeeping(i,j).second=getFullNumberOfTasks();
228 : }
229 : }
230 : } else {
231 138 : for(unsigned i=0; i<nblock; ++i) {
232 8344 : for(unsigned j=0; j<nblock; ++j) {
233 8210 : if( i==j ) continue ;
234 8076 : bookeeping(i,j).first=getFullNumberOfTasks();
235 8076 : addTaskToList( i*nblock + j );
236 8076 : bookeeping(i,j).second=getFullNumberOfTasks();
237 : }
238 : }
239 : }
240 10 : if( !parseMultiColvarAtomList(key3,-1,all_atoms) ) error("missing required keyword " + key3 + " in input");
241 10 : ablocks[2].resize( atom_lab.size() - ablocks[0].size() );
242 49845 : for(unsigned i=0; i<ablocks[2].size(); ++i) ablocks[2][i]=ablocks[0].size() + i;
243 : } else {
244 0 : nblock=atom_lab.size(); for(unsigned i=0; i<3; ++i) ablocks[i].resize(nblock);
245 0 : resizeBookeepingArray( nblock, nblock );
246 0 : for(unsigned i=0; i<nblock; ++i) { ablocks[0][i]=i; ablocks[1][i]=i; ablocks[2][i]=i; }
247 0 : if( symmetric ) {
248 0 : for(unsigned i=2; i<nblock; ++i) {
249 0 : for(unsigned j=1; j<i; ++j) {
250 0 : bookeeping(i,j).first=getFullNumberOfTasks();
251 0 : for(unsigned k=0; k<j; ++k) addTaskToList( i*nblock*nblock + j*nblock + k );
252 0 : bookeeping(i,j).second=getFullNumberOfTasks();
253 : }
254 : }
255 : } else {
256 0 : for(unsigned i=0; i<nblock; ++i) {
257 0 : for(unsigned j=0; j<nblock; ++j) {
258 0 : if( i==j ) continue;
259 0 : bookeeping(i,j).first=getFullNumberOfTasks();
260 0 : for(unsigned k=0; k<nblock; ++k) {
261 0 : if( i!=k && j!=k ) addTaskToList( i*nblock*nblock + j*nblock + k );
262 : }
263 0 : bookeeping(i,j).first=getFullNumberOfTasks();
264 : }
265 : }
266 : }
267 : }
268 : } else {
269 2 : readThreeGroups( key1, key2, key3, true, no_third_dim_accum, all_atoms );
270 : }
271 :
272 12 : }
273 :
274 10 : void MultiColvarBase::readThreeGroups( const std::string& key1, const std::string& key2, const std::string& key3,
275 : const bool& allow2, const bool& no_third_dim_accum, std::vector<AtomNumber>& all_atoms ) {
276 10 : plumed_dbg_assert( keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) ); ablocks.resize( 3 );
277 :
278 10 : bool readkey1=parseMultiColvarAtomList(key1,-1,all_atoms);
279 31 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=i;
280 10 : bool readkey2=parseMultiColvarAtomList(key2,-1,all_atoms);
281 10 : if( !readkey1 && !readkey2 ) return ;
282 225 : 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;
283 :
284 10 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
285 10 : bool readkey3=parseMultiColvarAtomList(key3,-1,all_atoms);
286 10 : if( !readkey3 && !allow2 ) {
287 0 : error("missing atom specification " + key3);
288 10 : } else if( !readkey3 ) {
289 2 : if( ablocks[1].size()>ablocks[0].size() ) nblock=ablocks[1].size();
290 0 : else nblock=ablocks[0].size();
291 :
292 2 : ablocks[2].resize( ablocks[1].size() );
293 200 : for(unsigned i=0; i<ablocks[1].size(); ++i) ablocks[2][i]=ablocks[0].size() + i;
294 4 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
295 198 : for(unsigned j=1; j<ablocks[1].size(); ++j) {
296 196 : bookeeping(i,j).first=getFullNumberOfTasks();
297 9898 : for(unsigned k=0; k<j; ++k) {
298 9702 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
299 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
300 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
301 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
302 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 &&
303 0 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) addTaskToList( nblock*nblock*i + nblock*j + k );
304 9702 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
305 9702 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
306 9702 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) addTaskToList( nblock*nblock*i + nblock*j + k );
307 : }
308 196 : bookeeping(i,j).second=getFullNumberOfTasks();
309 : }
310 : }
311 : } else {
312 8 : ablocks[2].resize( atom_lab.size() - ablocks[1].size() - ablocks[0].size() );
313 3182 : for(unsigned i=0; i<ablocks[2].size(); ++i) ablocks[2][i] = ablocks[0].size() + ablocks[1].size() + i;
314 :
315 8 : if( ablocks[1].size()>ablocks[0].size() ) nblock=ablocks[1].size();
316 8 : else nblock=ablocks[0].size();
317 8 : if( ablocks[2].size()>nblock ) nblock=ablocks[2].size();
318 :
319 8 : unsigned kcount; if( no_third_dim_accum ) kcount=1; else kcount=ablocks[2].size();
320 :
321 27 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
322 128 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
323 109 : bookeeping(i,j).first=getFullNumberOfTasks();
324 218 : for(unsigned k=0; k<kcount; ++k) {
325 109 : if( no_third_dim_accum ) addTaskToList( nblock*i + j );
326 0 : else if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
327 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
328 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
329 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
330 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 &&
331 0 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) addTaskToList( nblock*nblock*i + nblock*j + k );
332 0 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
333 0 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
334 0 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) addTaskToList( nblock*nblock*i + nblock*j + k );
335 : }
336 109 : bookeeping(i,j).second=getFullNumberOfTasks();
337 : }
338 : }
339 : }
340 : }
341 :
342 4 : void MultiColvarBase::buildSets() {
343 : std::vector<AtomNumber> fake_atoms;
344 8 : if( !parseMultiColvarAtomList("DATA",-1,fake_atoms) ) error("missing DATA keyword");
345 4 : if( fake_atoms.size()>0 ) error("no atoms should appear in the specification for this object. Input should be other multicolvars");
346 :
347 4 : nblock = mybasemulticolvars[0]->getFullNumberOfTasks();
348 11 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
349 7 : if( mybasemulticolvars[i]->getFullNumberOfTasks()!=nblock ) {
350 0 : error("mismatch between numbers of tasks in various base multicolvars");
351 : }
352 : }
353 4 : ablocks.resize( mybasemulticolvars.size() ); usespecies=false;
354 11 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
355 7 : ablocks[i].resize( nblock );
356 27 : for(unsigned j=0; j<nblock; ++j) ablocks[i][j]=i*nblock+j;
357 : }
358 15 : for(unsigned i=0; i<nblock; ++i) {
359 11 : if( mybasemulticolvars.size()<4 ) {
360 11 : unsigned cvcode=0, tmpc=1;
361 31 : for(unsigned j=0; j<ablocks.size(); ++j) { cvcode += i*tmpc; tmpc *= nblock; }
362 11 : addTaskToList( cvcode );
363 : } else {
364 0 : addTaskToList( i );
365 : }
366 : }
367 4 : mybasedata[0]->resizeTemporyMultiValues( mybasemulticolvars.size() ); setupMultiColvarBase( fake_atoms );
368 4 : }
369 :
370 12512606 : void MultiColvarBase::addTaskToList( const unsigned& taskCode ) {
371 12512606 : plumed_assert( getNumberOfVessels()==0 );
372 12512606 : ActionWithVessel::addTaskToList( taskCode );
373 12512606 : }
374 :
375 96 : void MultiColvarBase::resizeBookeepingArray( const unsigned& num1, const unsigned& num2 ) {
376 96 : bookeeping.resize( num1, num2 );
377 7065 : for(unsigned i=0; i<num1; ++i) {
378 8887414 : for(unsigned j=0; j<num2; ++j) { bookeeping(i,j).first=0; bookeeping(i,j).second=0; }
379 : }
380 96 : }
381 :
382 308 : void MultiColvarBase::setupMultiColvarBase( const std::vector<AtomNumber>& atoms ) {
383 308 : if( !matsums && atom_lab.size()==0 ) error("No atoms have been read in");
384 : std::vector<AtomNumber> all_atoms;
385 : // Setup decoder array
386 308 : if( !usespecies && nblock>0 ) {
387 :
388 63 : ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
389 63 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
390 63 : if( ablocks.size()==3 ) {
391 20 : allthirdblockintasks=uselinkforthree=true;
392 1512 : for(unsigned i=0; i<bookeeping.nrows(); ++i) {
393 453578 : for(unsigned j=0; j<bookeeping.ncols(); ++j) {
394 452086 : unsigned ntper = bookeeping(i,j).second - bookeeping(i,j).first;
395 452086 : if( i==j && ntper==0 ) {
396 136 : continue;
397 451950 : } else if( ntper == 1 && allthirdblockintasks ) {
398 451756 : allthirdblockintasks=true;
399 194 : } else if( ntper != ablocks[2].size() ) {
400 194 : allthirdblockintasks=uselinkforthree=false;
401 : } else {
402 0 : allthirdblockintasks=false;
403 : }
404 : }
405 : }
406 : }
407 :
408 63 : if( allthirdblockintasks ) {
409 18 : decoder.resize(2); plumed_assert( ablocks.size()==3 );
410 : // Check if number of atoms is too large
411 18 : if( std::pow( double(nblock), 2.0 )>std::numeric_limits<unsigned>::max() ) error("number of atoms in groups is too big for PLUMED to handle");
412 : } else {
413 45 : decoder.resize( ablocks.size() );
414 : // Check if number of atoms is too large
415 45 : if( std::pow( double(nblock), double(ablocks.size()) )>std::numeric_limits<unsigned>::max() ) error("number of atoms in groups is too big for PLUMED to handle");
416 : }
417 190 : unsigned code=1; for(unsigned i=0; i<decoder.size(); ++i) { decoder[decoder.size()-1-i]=code; code *= nblock; }
418 245 : } else if( !usespecies ) {
419 87 : ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
420 87 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
421 316 : } else if( keywords.exists("SPECIESA") ) {
422 100 : plumed_assert( atom_lab.size()==0 && all_atoms.size()==0 );
423 100 : ablocks.resize( 1 ); bool readspecies=parseMultiColvarAtomList("SPECIES", -1, all_atoms);
424 100 : if( readspecies ) {
425 41493 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<atom_lab.size(); ++i) { addTaskToList(i); ablocks[0][i]=i; }
426 : } else {
427 38 : if( !parseMultiColvarAtomList("SPECIESA", -1, all_atoms) ) error("missing SPECIES/SPECIESA keyword");
428 19 : unsigned nat1=atom_lab.size();
429 38 : if( !parseMultiColvarAtomList("SPECIESB", -1, all_atoms) ) error("missing SPECIESB keyword");
430 19 : unsigned nat2=atom_lab.size() - nat1;
431 :
432 662 : for(unsigned i=0; i<nat1; ++i) addTaskToList(i);
433 19 : ablocks[0].resize( nat2 );
434 3437 : for(unsigned i=0; i<nat2; ++i) {
435 : bool found=false; unsigned inum=0;
436 265449 : for(unsigned j=0; j<nat1; ++j) {
437 262113 : if( atom_lab[nat1+i].first>0 && atom_lab[j].first>0 ) {
438 252720 : if( atom_lab[nat1+i].first==atom_lab[j].first &&
439 0 : mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
440 252720 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) { found=true; inum=j; break; }
441 9393 : } else if( atom_lab[nat1+i].first>0 ) {
442 0 : if( mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
443 0 : all_atoms[atom_lab[j].second] ) { found=true; inum=nat1 + i; break; }
444 9393 : } else if( atom_lab[j].first>0 ) {
445 0 : if( all_atoms[atom_lab[nat1+i].second]==
446 0 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) { found=true; inum=nat1+i; break; }
447 9393 : } else if( all_atoms[atom_lab[nat1+i].second]==all_atoms[atom_lab[j].second] ) { found=true; inum=j; break; }
448 : }
449 : // This prevents mistakes being made in colvar setup
450 82 : if( found ) { ablocks[0][i]=inum; }
451 3336 : else { ablocks[0][i]=nat1 + i; }
452 : }
453 : }
454 : }
455 308 : if( mybasemulticolvars.size()>0 ) {
456 246 : for(unsigned i=0; i<mybasedata.size(); ++i) {
457 127 : mybasedata[i]->resizeTemporyMultiValues(2); mybasemulticolvars[i]->my_tmp_capacks.resize(2);
458 : }
459 : }
460 :
461 : // Copy lists of atoms involved from base multicolvars
462 : std::vector<AtomNumber> tmp_atoms;
463 435 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
464 127 : tmp_atoms=mybasemulticolvars[i]->getAbsoluteIndexes();
465 416316 : for(unsigned j=0; j<tmp_atoms.size(); ++j) all_atoms.push_back( tmp_atoms[j] );
466 : }
467 : // Copy atom lists from input
468 59575 : for(unsigned i=0; i<atoms.size(); ++i) all_atoms.push_back( atoms[i] );
469 :
470 : // Now make sure we get all the atom positions
471 308 : ActionAtomistic::requestAtoms( all_atoms );
472 : // And setup dependencies
473 435 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) addDependency( mybasemulticolvars[i] );
474 :
475 : // Setup underlying ActionWithVessel
476 308 : readVesselKeywords();
477 308 : }
478 :
479 19 : void MultiColvarBase::setAtomsForCentralAtom( const std::vector<bool>& catom_ind ) {
480 19 : unsigned nat=0; plumed_assert( catom_ind.size()==ablocks.size() );
481 89 : for(unsigned i=0; i<catom_ind.size(); ++i) {
482 70 : use_for_central_atom[i]=catom_ind[i];
483 70 : if( use_for_central_atom[i] ) nat++;
484 : }
485 19 : plumed_dbg_assert( nat>0 ); ncentral=nat;
486 19 : numberForCentralAtom = 1.0 / static_cast<double>( nat );
487 19 : }
488 :
489 429 : void MultiColvarBase::turnOnDerivatives() {
490 429 : ActionWithValue::turnOnDerivatives();
491 429 : needsDerivatives();
492 429 : forcesToApply.resize( getNumberOfDerivatives() );
493 429 : }
494 :
495 191 : void MultiColvarBase::setLinkCellCutoff( const double& lcut, double tcut ) {
496 191 : plumed_assert( usespecies || ablocks.size()<4 );
497 191 : if( tcut<0 ) tcut=lcut;
498 :
499 191 : if( !linkcells.enabled() ) {
500 191 : linkcells.setCutoff( lcut );
501 191 : threecells.setCutoff( tcut );
502 : } else {
503 0 : if( lcut>linkcells.getCutoff() ) linkcells.setCutoff( lcut );
504 0 : if( tcut>threecells.getCutoff() ) threecells.setCutoff( tcut );
505 : }
506 191 : }
507 :
508 8 : double MultiColvarBase::getLinkCellCutoff() const {
509 8 : return linkcells.getCutoff();
510 : }
511 :
512 1938 : void MultiColvarBase::setupLinkCells() {
513 1938 : if( (!usespecies && nblock==0) || !linkcells.enabled() ) return ;
514 : // Retrieve any atoms that haven't already been retrieved
515 1318 : for(std::vector<MultiColvarBase*>::iterator p=mybasemulticolvars.begin(); p!=mybasemulticolvars.end(); ++p) {
516 203 : (*p)->retrieveAtoms();
517 : }
518 1115 : retrieveAtoms();
519 :
520 : unsigned iblock;
521 1115 : if( usespecies ) {
522 : iblock=0;
523 217 : } else if( ablocks.size()<4 ) {
524 : iblock=1;
525 : } else {
526 0 : plumed_error();
527 : }
528 :
529 : // Count number of currently active atoms
530 1115 : nactive_atoms=0;
531 387665 : for(unsigned i=0; i<ablocks[iblock].size(); ++i) {
532 386550 : if( isCurrentlyActive( ablocks[iblock][i] ) ) nactive_atoms++;
533 : }
534 :
535 1115 : if( nactive_atoms>0 ) {
536 1115 : std::vector<Vector> ltmp_pos( nactive_atoms );
537 1115 : std::vector<unsigned> ltmp_ind( nactive_atoms );
538 :
539 1115 : nactive_atoms=0;
540 1115 : if( usespecies ) {
541 366350 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
542 365452 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
543 365356 : ltmp_ind[nactive_atoms]=ablocks[0][i];
544 365356 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ltmp_ind[nactive_atoms] );
545 365356 : nactive_atoms++;
546 : }
547 : } else {
548 21315 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
549 21098 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
550 16178 : ltmp_ind[nactive_atoms]=i;
551 16178 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ablocks[1][i] );
552 16178 : nactive_atoms++;
553 : }
554 : }
555 :
556 : // Build the lists for the link cells
557 1115 : linkcells.buildCellLists( ltmp_pos, ltmp_ind, getPbc() );
558 : }
559 : }
560 :
561 4628 : void MultiColvarBase::setupNonUseSpeciesLinkCells( const unsigned& my_always_active ) {
562 4628 : plumed_assert( !usespecies );
563 4628 : if( nblock==0 || !linkcells.enabled() ) return ;
564 210 : deactivateAllTasks();
565 : std::vector<unsigned> requiredlinkcells;
566 :
567 210 : if( !uselinkforthree && nactive_atoms>0 ) {
568 : // Get some parallel info
569 156 : unsigned stride=comm.Get_size();
570 156 : unsigned rank=comm.Get_rank();
571 156 : if( serialCalculation() ) { stride=1; rank=0; }
572 :
573 : // Ensure we only do tasks where atoms are in appropriate link cells
574 156 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
575 5757 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
576 5601 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
577 2622 : unsigned natomsper=1; linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
578 2622 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, natomsper, linked_atoms );
579 205705 : for(unsigned j=0; j<natomsper; ++j) {
580 353118 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) taskFlags[k]=1;
581 : }
582 : }
583 210 : } else if( nactive_atoms>0 ) {
584 : // Get some parallel info
585 54 : unsigned stride=comm.Get_size();
586 54 : unsigned rank=comm.Get_rank();
587 54 : if( serialCalculation() ) { stride=1; rank=0; }
588 :
589 : unsigned nactive_three=0;
590 66599 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
591 66545 : if( isCurrentlyActive( ablocks[2][i] ) ) nactive_three++;
592 : }
593 :
594 54 : std::vector<Vector> lttmp_pos( nactive_three );
595 54 : std::vector<unsigned> lttmp_ind( nactive_three );
596 :
597 : nactive_three=0;
598 54 : if( allthirdblockintasks ) {
599 66599 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
600 66545 : if( !isCurrentlyActive( ablocks[2][i] ) ) continue;
601 66545 : lttmp_ind[nactive_three]=ablocks[2][i];
602 66545 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
603 66545 : nactive_three++;
604 : }
605 : } else {
606 0 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
607 0 : if( !isCurrentlyActive( ablocks[2][i] ) ) continue;
608 0 : lttmp_ind[nactive_three]=i;
609 0 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
610 0 : nactive_three++;
611 : }
612 : }
613 : // Build the list of the link cells
614 54 : threecells.buildCellLists( lttmp_pos, lttmp_ind, getPbc() );
615 :
616 : // Ensure we only do tasks where atoms are in appropriate link cells
617 54 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
618 54 : std::vector<unsigned> tlinked_atoms( 1+ablocks[2].size() );
619 633 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
620 579 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
621 579 : unsigned natomsper=1; linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
622 579 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, natomsper, linked_atoms );
623 579 : if( allthirdblockintasks ) {
624 71832 : for(unsigned j=0; j<natomsper; ++j) {
625 142372 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) taskFlags[k]=1;
626 : }
627 : } else {
628 0 : unsigned ntatomsper=1; tlinked_atoms[0]=lttmp_ind[0];
629 0 : threecells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), requiredlinkcells, ntatomsper, tlinked_atoms );
630 0 : for(unsigned j=0; j<natomsper; ++j) {
631 0 : for(unsigned k=0; k<ntatomsper; ++k) taskFlags[bookeeping(i,linked_atoms[j]).first+tlinked_atoms[k]]=1;
632 : }
633 : }
634 : }
635 : }
636 210 : if( !serialCalculation() ) comm.Sum( taskFlags );
637 210 : lockContributors();
638 : }
639 :
640 245988 : void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const {
641 : plumed_dbg_assert( !usespecies && nblock>0 );
642 245988 : if( atoms.size()!=decoder.size() ) atoms.resize( decoder.size() );
643 :
644 245988 : unsigned scode = taskCode;
645 786464 : for(unsigned i=0; i<decoder.size(); ++i) {
646 540476 : unsigned ind=( scode / decoder[i] );
647 540476 : atoms[i] = ablocks[i][ind];
648 540476 : scode -= ind*decoder[i];
649 : }
650 245988 : }
651 :
652 451407 : bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ) const {
653 451407 : if( isDensity() ) {
654 19070 : myatoms.setNumberOfAtoms( 1 ); myatoms.setAtom( 0, taskCode ); return true;
655 432337 : } else if( usespecies ) {
656 180751 : std::vector<unsigned> task_atoms(1); task_atoms[0]=taskCode;
657 180751 : unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getPositionOfAtomForLinkCells( taskCode ), linkcells );
658 180751 : return natomsper>1;
659 251586 : } else if( matsums ) {
660 : myatoms.setNumberOfAtoms( getNumberOfAtoms() );
661 52798 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) myatoms.setAtom( i, i );
662 251145 : } else if( allthirdblockintasks ) {
663 39816 : plumed_dbg_assert( ablocks.size()==3 ); std::vector<unsigned> atoms(2); decodeIndexToAtoms( taskCode, atoms );
664 39816 : myatoms.setupAtomsFromLinkCells( atoms, getPositionOfAtomForLinkCells( atoms[0] ), threecells );
665 211329 : } else if( nblock>0 ) {
666 179551 : std::vector<unsigned> atoms( ablocks.size() );
667 179551 : decodeIndexToAtoms( taskCode, atoms ); myatoms.setNumberOfAtoms( ablocks.size() );
668 587153 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, atoms[i] );
669 : } else {
670 31778 : myatoms.setNumberOfAtoms( ablocks.size() );
671 124498 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, ablocks[i][taskCode] );
672 : }
673 : return true;
674 : }
675 :
676 9073 : void MultiColvarBase::setupActiveTaskSet( std::vector<unsigned>& active_tasks, const std::string& input_label ) {
677 9073 : if( !setup_completed ) {
678 : bool justVolumes=false;
679 1930 : if( usespecies ) {
680 : justVolumes=true;
681 1437 : for(unsigned i=0; i<getNumberOfVessels(); ++i) {
682 1317 : vesselbase::StoreDataVessel* mys=dynamic_cast<vesselbase::StoreDataVessel*>( getPntrToVessel(i) );
683 1317 : if( mys ) continue;
684 809 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(i) );
685 809 : if( !myb ) { justVolumes=false; break; }
686 38 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
687 38 : if( !myv ) { justVolumes=false; break; }
688 : }
689 : }
690 1930 : deactivateAllTasks();
691 1930 : if( justVolumes && mydata ) {
692 88 : if( mydata->getNumberOfDataUsers()==0 ) justVolumes=false;
693 :
694 137 : for(unsigned i=0; i<mydata->getNumberOfDataUsers(); ++i) {
695 49 : MultiColvarBase* myu=dynamic_cast<MultiColvarBase*>( mydata->getDataUser(i) );
696 49 : if( myu ) {
697 49 : myu->setupActiveTaskSet( taskFlags, getLabel() );
698 : } else {
699 0 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
700 : }
701 : }
702 : }
703 1930 : if( justVolumes ) {
704 160 : for(unsigned j=0; j<getNumberOfVessels(); ++j) {
705 80 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(j) );
706 80 : if( !myb ) continue ;
707 32 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
708 32 : if( !myv ) continue ;
709 32 : myv->retrieveAtoms(); myv->setupRegions();
710 :
711 148161 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
712 148129 : if( myv->inVolumeOfInterest(i) ) taskFlags[i]=1;
713 : }
714 : }
715 : } else {
716 4886605 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
717 : }
718 :
719 : // Now activate all this class
720 1930 : lockContributors();
721 : // Setup the link cells
722 1930 : setupLinkCells();
723 : // Ensures that setup is not performed multiple times during one cycle
724 1930 : setup_completed=true;
725 : }
726 :
727 : // And activate the tasks in input action
728 9073 : if( getLabel()!=input_label ) {
729 : int input_code=-1;
730 61 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
731 61 : if( mybasemulticolvars[i]->getLabel()==input_label ) { input_code=i+1; break; }
732 : }
733 :
734 49 : MultiValue my_tvals( getNumberOfQuantities(), getNumberOfDerivatives() );
735 49 : AtomValuePack mytmp_atoms( my_tvals, this );
736 148234 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
737 148185 : if( !taskIsCurrentlyActive(i) ) continue;
738 3529 : setupCurrentAtomList( getTaskCode(i), mytmp_atoms );
739 254796 : for(unsigned j=0; j<mytmp_atoms.getNumberOfAtoms(); ++j) {
740 : unsigned itask=mytmp_atoms.getIndex(j);
741 251267 : if( atom_lab[itask].first==input_code ) active_tasks[ atom_lab[itask].second ]=1;
742 : }
743 : }
744 49 : }
745 9073 : }
746 :
747 467 : bool MultiColvarBase::filtersUsedAsInput() {
748 : bool inputAreFilters=false;
749 728 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
750 261 : MultiColvarFilter* myfilt=dynamic_cast<MultiColvarFilter*>( mybasemulticolvars[i] );
751 261 : if( myfilt || mybasemulticolvars[i]->filtersUsedAsInput() ) inputAreFilters=true;
752 : }
753 467 : return inputAreFilters;
754 : }
755 :
756 9024 : void MultiColvarBase::calculate() {
757 : // Recursive function that sets up tasks
758 9024 : setupActiveTaskSet( taskFlags, getLabel() );
759 :
760 : // Check for filters and rerun setup of link cells if there are any
761 9024 : if( mybasemulticolvars.size()>0 && filtersUsedAsInput() ) setupLinkCells();
762 :
763 : // Setup the link cells if we are not using species
764 9024 : if( !usespecies && ablocks.size()>1 ) {
765 : // This loop finds the first active atom, which is always checked because
766 : // of a peculiarity in linkcells
767 4628 : unsigned first_active=std::numeric_limits<unsigned>::max();
768 4762 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
769 4762 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
770 : else {
771 4628 : first_active=i; break;
772 : }
773 : }
774 4628 : setupNonUseSpeciesLinkCells( first_active );
775 : }
776 : // And run all tasks
777 9024 : runAllTasks();
778 9024 : }
779 :
780 213 : void MultiColvarBase::calculateNumericalDerivatives( ActionWithValue* a ) {
781 213 : if( mybasemulticolvars.size()>0 ) plumed_merror("cannot calculate numerical derivatives for this quantity");
782 213 : calculateAtomicNumericalDerivatives( this, 0 );
783 213 : }
784 :
785 2959 : void MultiColvarBase::prepare() {
786 2959 : setup_completed=false; atomsWereRetrieved=false;
787 2959 : }
788 :
789 6548 : void MultiColvarBase::retrieveAtoms() {
790 6548 : if( !atomsWereRetrieved ) { ActionAtomistic::retrieveAtoms(); atomsWereRetrieved=true; }
791 6548 : }
792 :
793 38245 : void MultiColvarBase::mergeInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
794 : const unsigned& jatom, const std::vector<double>& der,
795 : MultiValue& myder, AtomValuePack& myatoms ) const {
796 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
797 : plumed_dbg_assert( ival<myatoms.getUnderlyingMultiValue().getNumberOfValues() );
798 : plumed_dbg_assert( start<myder.getNumberOfValues() && end<=myder.getNumberOfValues() );
799 : plumed_dbg_assert( der.size()==myder.getNumberOfValues() && jatom<myatoms.getNumberOfAtoms() );
800 : // Convert input atom to local index
801 : unsigned katom = myatoms.getIndex( jatom ); plumed_dbg_assert( katom<atom_lab.size() ); plumed_dbg_assert( atom_lab[katom].first>0 );
802 : // Find base colvar
803 38245 : unsigned mmc=atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
804 : // Get start of indices for this atom
805 38855 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
806 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
807 38245 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
808 4805851 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
809 : unsigned jder=myder.getActiveIndex(j);
810 4767606 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
811 4434561 : unsigned kder=basen+jder;
812 110197260 : for(unsigned icomp=start; icomp<end; ++icomp) {
813 105762699 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
814 : }
815 : } else {
816 333045 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
817 7551360 : for(unsigned icomp=start; icomp<end; ++icomp) {
818 7218315 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
819 : }
820 : }
821 : }
822 38245 : }
823 :
824 106 : void MultiColvarBase::splitInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
825 : const unsigned& jatom, const std::vector<double>& der,
826 : MultiValue& myder, AtomValuePack& myatoms ) const {
827 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
828 : plumed_dbg_assert( ival<myder.getNumberOfValues() );
829 : plumed_dbg_assert( start<myvals.getNumberOfValues() && end<=myvals.getNumberOfValues() );
830 : plumed_dbg_assert( der.size()==myatoms.getUnderlyingMultiValue().getNumberOfValues() && jatom<myatoms.getNumberOfAtoms() );
831 : // Convert input atom to local index
832 : unsigned katom = myatoms.getIndex( jatom ); plumed_dbg_assert( katom<atom_lab.size() ); plumed_dbg_assert( atom_lab[katom].first>0 );
833 : // Find base colvar
834 106 : unsigned mmc=atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
835 : // Get start of indices for this atom
836 154 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
837 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
838 106 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
839 19534 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
840 : unsigned jder=myder.getActiveIndex(j);
841 19428 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
842 18474 : unsigned kder=basen+jder; plumed_assert( kder<myvals.getNumberOfDerivatives() );
843 39942 : for(unsigned icomp=start; icomp<end; ++icomp) {
844 21468 : myvals.addDerivative( icomp, kder, der[icomp]*myder.getDerivative( ival, jder ) );
845 : }
846 : } else {
847 954 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
848 3942 : for(unsigned icomp=start; icomp<end; ++icomp) {
849 2988 : myvals.addDerivative( icomp, kder, der[icomp]*myder.getDerivative( ival, jder ) );
850 : }
851 : }
852 : }
853 106 : }
854 :
855 905762 : void MultiColvarBase::addComDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
856 : plumed_dbg_assert( ival<static_cast<int>(myatoms.getUnderlyingMultiValue().getNumberOfValues()) && iatom<myatoms.getNumberOfAtoms() );
857 : // Convert input atom to local index
858 : unsigned katom = myatoms.getIndex( iatom ); plumed_dbg_assert( atom_lab[katom].first>0 );
859 : // Find base colvar
860 905762 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
861 905762 : if( usespecies && iatom==0 ) { myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[0] ); return; }
862 :
863 : // Get start of indices for this atom
864 646526 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=(mybasemulticolvars[i]->getNumberOfDerivatives() - 9) / 3;
865 489261 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[katom].second, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
866 489261 : myatoms.addComDerivatives( ival, der, mybasemulticolvars[mmc]->my_tmp_capacks[1] );
867 : }
868 :
869 1122658 : void MultiColvarBase::getInputData( const unsigned& ind, const bool& normed,
870 : const multicolvar::AtomValuePack& myatoms,
871 : std::vector<double>& orient ) const {
872 : // Converint input atom to local index
873 : unsigned katom = myatoms.getIndex(ind); plumed_dbg_assert( atom_lab[katom].first>0 );
874 : // Find base colvar
875 1122658 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
876 : // Check if orient is the correct size
877 1122658 : if( orient.size()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ) orient.resize( mybasemulticolvars[mmc]->getNumberOfQuantities() );
878 : // Retrieve the value
879 1122658 : mybasedata[mmc]->retrieveValueWithIndex( atom_lab[katom].second, normed, orient );
880 1122658 : }
881 :
882 53063 : MultiValue& MultiColvarBase::getInputDerivatives( const unsigned& iatom, const bool& normed, const multicolvar::AtomValuePack& myatoms ) const {
883 : // Converint input atom to local index
884 : unsigned katom = myatoms.getIndex(iatom); plumed_dbg_assert( atom_lab[katom].first>0 );
885 : // Find base colvar
886 53063 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
887 53063 : if( usespecies && !normed && iatom==0 ) return mybasedata[mmc]->getTemporyMultiValue(0);
888 :
889 51265 : unsigned oval=0; if( iatom>0 ) oval=1;
890 51265 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(oval);
891 102491 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
892 51226 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
893 39 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
894 : }
895 51265 : mybasedata[mmc]->retrieveDerivatives( atom_lab[katom].second, normed, myder );
896 : return myder;
897 : }
898 :
899 64359697 : void MultiColvarBase::accumulateSymmetryFunction( const int& ival, const unsigned& iatom, const double& val, const Vector& der, const Tensor& vir, multicolvar::AtomValuePack& myatoms ) const {
900 : plumed_dbg_assert( usespecies ); unsigned katom=myatoms.getIndex(0), jatom=myatoms.getIndex(iatom);
901 64359697 : double weight0=1.0; if( atom_lab[katom].first>0 ) weight0=mybasedata[atom_lab[katom].first-1]->retrieveWeightWithIndex( atom_lab[katom].second );
902 64359697 : double weighti=1.0; if( atom_lab[jatom].first>0 ) weighti=mybasedata[atom_lab[jatom].first-1]->retrieveWeightWithIndex( atom_lab[jatom].second );
903 : // Accumulate the value
904 64359697 : if( ival<0 ) myatoms.getUnderlyingMultiValue().addTemporyValue( weight0*weighti*val );
905 61622797 : else myatoms.addValue( ival, weight0*weighti*val );
906 :
907 : // Return if we don't need derivatives
908 64359697 : if( doNotCalculateDerivatives() ) return ;
909 : // And virial
910 57992318 : if( ival<0 ) myatoms.addTemporyBoxDerivatives( weight0*weighti*vir );
911 55734171 : else myatoms.addBoxDerivatives( ival, weight0*weighti*vir );
912 :
913 : // Add derivatives of central atom
914 57992318 : if( atom_lab[katom].first>0 ) {
915 794 : addComDerivatives( ival, 0, -weight0*weighti*der, myatoms );
916 794 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
917 794 : tmpder[0]=weighti*val; mergeInputDerivatives( ival, 0, 1, 0, tmpder, getInputDerivatives(0, false, myatoms), myatoms );
918 : } else {
919 57991524 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( 0, -der );
920 55733377 : else myatoms.addAtomsDerivatives( ival, 0, -der );
921 : }
922 : // Add derivatives of atom in coordination sphere
923 57992318 : if( atom_lab[jatom].first>0 ) {
924 794 : addComDerivatives( ival, iatom, weight0*weighti*der, myatoms );
925 794 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
926 794 : tmpder[0]=weight0*val; mergeInputDerivatives( ival, 0, 1, iatom, tmpder, getInputDerivatives(iatom, false, myatoms), myatoms );
927 : } else {
928 57991524 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
929 55733377 : else myatoms.addAtomsDerivatives( ival, iatom, der );
930 : }
931 : }
932 :
933 1428471 : void MultiColvarBase::addAtomDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
934 1428471 : if( doNotCalculateDerivatives() ) return ;
935 : unsigned jatom=myatoms.getIndex(iatom);
936 :
937 1180895 : if( atom_lab[jatom].first>0 ) {
938 904174 : addComDerivatives( ival, iatom, der, myatoms );
939 : } else {
940 276721 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
941 276721 : else myatoms.addAtomsDerivatives( ival, iatom, der );
942 : }
943 : }
944 :
945 242764 : double MultiColvarBase::calculateWeight( const unsigned& current, const double& weight, AtomValuePack& myvals ) const {
946 242764 : return 1.0;
947 : }
948 :
949 447878 : void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
950 447878 : AtomValuePack myatoms( myvals, this );
951 : // Retrieve the atom list
952 447878 : if( !setupCurrentAtomList( current, myatoms ) ) return;
953 : // Get weight due to dynamic groups
954 447758 : double weight = 1.0;
955 447758 : if( !matsums ) {
956 386210280 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
957 385940036 : if( atom_lab[myatoms.getIndex(i)].first==0 ) continue;
958 : // Only need to do first two atoms for things like TopologyMatrix, HbondMatrix, Bridge and so on
959 244017 : if( allthirdblockintasks && i>1 ) break;
960 244017 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
961 244017 : weight *= mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
962 : }
963 177514 : } else if( usespecies ) {
964 177073 : if( atom_lab[myatoms.getIndex(0)].first>0 ) {
965 9073 : if( mybasedata[atom_lab[myatoms.getIndex(0)].first-1]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(0)].second )<epsilon ) weight=0.;
966 : }
967 : }
968 : // Do a quick check on the size of this contribution
969 447758 : double multweight = calculateWeight( current, weight, myatoms );
970 447758 : if( weight*multweight<getTolerance() ) {
971 121536 : updateActiveAtoms( myatoms );
972 : return;
973 : }
974 : myatoms.setValue( 0, weight*multweight );
975 : // Deal with derivatives of weights due to dynamic groups
976 326222 : if( !matsums && !doNotCalculateDerivatives() && mybasemulticolvars.size()>0 ) {
977 39416 : MultiValue& outder=myatoms.getUnderlyingMultiValue(); MultiValue myder(0,0);
978 115123 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
979 : // Neglect any atoms without differentiable weights
980 75707 : if( atom_lab[myatoms.getIndex(i)].first==0 ) continue;
981 :
982 : // Retrieve derivatives
983 75707 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
984 75707 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() || myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
985 39416 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
986 : }
987 75707 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(i)].second, false, myder );
988 :
989 : // Retrieve the prefactor (product of all other weights)
990 75707 : double prefactor = multweight*weight / mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
991 : // And accumulate the derivatives
992 2927141 : for(unsigned j=0; j<myder.getNumberActive(); ++j) { unsigned jder=myder.getActiveIndex(j); outder.addDerivative( 0, jder, prefactor*myder.getDerivative(0,jder) ); }
993 75707 : myder.clearAll();
994 : }
995 39416 : }
996 : // Retrieve derivative stuff for central atom
997 326222 : if( !doNotCalculateDerivatives() ) {
998 206915 : if( usespecies && mybasemulticolvars.size()>0 && atom_lab[myatoms.getIndex(0)].first>0 ) {
999 1610 : unsigned mmc = atom_lab[0].first - 1;
1000 1610 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(0);
1001 3208 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
1002 1598 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
1003 12 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
1004 : }
1005 1610 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(0)].second, false, myder );
1006 1610 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
1007 1610 : mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[myatoms.getIndex(0)].second, mybasemulticolvars[mmc]->my_tmp_capacks[0] );
1008 : }
1009 : }
1010 : // Compute everything
1011 326222 : double vv=compute( task_index, myatoms ); updateActiveAtoms( myatoms );
1012 : myatoms.setValue( 1, vv );
1013 326222 : return;
1014 : }
1015 :
1016 658140 : void MultiColvarBase::updateActiveAtoms( AtomValuePack& myatoms ) const {
1017 658140 : if( mybasemulticolvars.size()==0 ) myatoms.updateUsingIndices();
1018 141871 : else myatoms.updateDynamicList();
1019 658140 : }
1020 :
1021 2753935 : Vector MultiColvarBase::getCentralAtomPos( const unsigned& taskIndex ) {
1022 2753935 : unsigned curr=getTaskCode( taskIndex );
1023 :
1024 2753935 : if( usespecies || isDensity() ) {
1025 1459946 : return getPositionOfAtomForLinkCells(curr);
1026 1293989 : } else if( nblock>0 ) {
1027 : // double factor=1.0/static_cast<double>( ablocks.size() );
1028 2130 : Vector mypos; mypos.zero();
1029 2130 : std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
1030 6390 : for(unsigned i=0; i<ablocks.size(); ++i) {
1031 4260 : if( use_for_central_atom[i] ) mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(atoms[i]);
1032 : }
1033 2130 : return mypos;
1034 : } else {
1035 1291859 : Vector mypos; mypos.zero();
1036 5138646 : for(unsigned i=0; i<ablocks.size(); ++i) {
1037 3846787 : if( use_for_central_atom[i] ) mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(ablocks[i][curr]);
1038 : }
1039 1291859 : return mypos;
1040 : }
1041 : }
1042 :
1043 605410 : void MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsigned& taskIndex, CatomPack& mypack ) {
1044 605410 : unsigned curr=getTaskCode( taskIndex );
1045 :
1046 605410 : if(usespecies) {
1047 464996 : if( mypack.getNumberOfAtomsWithDerivatives()!=1 ) mypack.resize(1);
1048 464996 : mypack.setIndex( 0, basn + curr );
1049 464996 : mypack.setDerivative( 0, Tensor::identity() );
1050 140414 : } else if( nblock>0 ) {
1051 0 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral);
1052 0 : unsigned k=0;
1053 0 : std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
1054 0 : for(unsigned i=0; i<ablocks.size(); ++i) {
1055 0 : if( use_for_central_atom[i] ) {
1056 0 : mypack.setIndex( k, basn + atoms[i] );
1057 0 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
1058 0 : k++;
1059 : }
1060 : }
1061 : } else {
1062 140414 : if( mypack.getNumberOfAtomsWithDerivatives()!=ncentral ) mypack.resize(ncentral);
1063 140414 : unsigned k=0;
1064 316748 : for(unsigned i=0; i<ablocks.size(); ++i) {
1065 176334 : if( use_for_central_atom[i] ) {
1066 173694 : mypack.setIndex( k, basn + ablocks[i][curr] );
1067 173694 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
1068 173694 : k++;
1069 : }
1070 : }
1071 : }
1072 605410 : }
1073 :
1074 121150993 : Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 ) const {
1075 121150993 : if(usepbc) { return pbcDistance( vec1, vec2 ); }
1076 0 : else { return delta( vec1, vec2 ); }
1077 : }
1078 :
1079 220567 : void MultiColvarBase::applyPbc(std::vector<Vector>& dlist, unsigned int max_index) const {
1080 220567 : if (usepbc) pbcApply(dlist, max_index);
1081 220567 : }
1082 :
1083 1994 : void MultiColvarBase::apply() {
1084 1994 : if( getForcesFromVessels( forcesToApply ) ) setForcesOnAtoms( forcesToApply );
1085 1994 : }
1086 :
1087 : }
1088 : }
|
