Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2018 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 "CatomPack.h"
32 : #include "CatomPack.h"
33 : #include <vector>
34 : #include <string>
35 : #include <limits>
36 :
37 : using namespace std;
38 :
39 : namespace PLMD {
40 : namespace multicolvar {
41 :
42 320 : void MultiColvarBase::registerKeywords( Keywords& keys ) {
43 320 : Action::registerKeywords( keys );
44 320 : ActionWithValue::registerKeywords( keys );
45 320 : ActionAtomistic::registerKeywords( keys );
46 320 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
47 320 : ActionWithVessel::registerKeywords( keys );
48 : keys.add("hidden","NL_STRIDE","the frequency with which the neighbor list should be updated. Between neighbour list update steps all quantities "
49 320 : "that contributed less than TOL at the previous neighbor list update step are ignored.");
50 : 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 "
51 : "regular collective variables. The label is used to reference the full set of quantities calculated by "
52 : "the action. This is usual when using \\ref multicolvarfunction. Generally when doing this the previously calculated "
53 : "multicolvar will be referenced using the DATA keyword rather than ARG.\n\n"
54 : "This Action can be used to calculate the following scalar quantities directly. These quantities are calculated by "
55 : "employing the keywords listed below. "
56 : "These quantities can then be referenced elsewhere in the input file by using this Action's label "
57 : "followed by a dot and the name of the quantity. Some amongst them can be calculated multiple times "
58 : "with different parameters. In this case the quantities calculated can be referenced elsewhere in the "
59 : "input by using the name of the quantity followed by a numerical identifier "
60 : "e.g. <em>label</em>.lessthan-1, <em>label</em>.lessthan-2 etc. When doing this and, for clarity we have "
61 : "made the label of the components customizable. As such by using the LABEL keyword in the description of the keyword "
62 320 : "input you can customize the component name");
63 : keys.reserve("atoms-3","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
64 : "one coordination number for each of the atoms specified. Each of these coordination numbers specifies how many of the "
65 : "other specified atoms are within a certain cutoff of the central atom. You can specify the atoms here as another multicolvar "
66 : "action or using a MultiColvarFilter or ActionVolume action. When you do so the quantity is calculated for those atoms specified "
67 : "in the previous multicolvar. This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
68 320 : "coordination number more than four for example");
69 : 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 "
70 : "one coordination number for each of the atoms specified in SPECIESA. Each of these cooordination numbers specifies how many "
71 : "of the atoms specifies using SPECIESB is within the specified cutoff. As with the species keyword the input can also be specified "
72 320 : "using the label of another multicolvar");
73 : 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 "
74 320 : "the documentation for that keyword");
75 320 : }
76 :
77 264 : MultiColvarBase::MultiColvarBase(const ActionOptions&ao):
78 : Action(ao),
79 : ActionAtomistic(ao),
80 : ActionWithValue(ao),
81 : ActionWithVessel(ao),
82 : usepbc(false),
83 : allthirdblockintasks(false),
84 : uselinkforthree(false),
85 : linkcells(comm),
86 : threecells(comm),
87 : setup_completed(false),
88 : atomsWereRetrieved(false),
89 : matsums(false),
90 : usespecies(false),
91 264 : nblock(0)
92 : {
93 264 : if( keywords.exists("NOPBC") ) {
94 264 : bool nopbc=!usepbc; parseFlag("NOPBC",nopbc);
95 264 : usepbc=!nopbc;
96 : }
97 264 : if( keywords.exists("SPECIESA") ) { matsums=usespecies=true; }
98 264 : }
99 :
100 350 : bool MultiColvarBase::parseMultiColvarAtomList(const std::string& key, const int& num, std::vector<AtomNumber>& t) {
101 350 : std::vector<std::string> mlabs;
102 350 : if( num<0 ) parseVector(key,mlabs);
103 0 : else parseNumberedVector(key,num,mlabs);
104 :
105 350 : if( mlabs.size()==0 ) return false;
106 :
107 440 : std::string mname; unsigned found_mcolv=mlabs.size();
108 303 : for(unsigned i=0; i<mlabs.size(); ++i) {
109 222 : MultiColvarBase* mycolv = plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlabs[i]);
110 222 : if(!mycolv) { found_mcolv=i; break; }
111 : // Check all base multicolvars are of same type
112 83 : if( i==0 ) {
113 81 : mname = mycolv->getName();
114 81 : if( mycolv->isPeriodic() ) error("multicolvar functions don't work with this multicolvar");
115 : } else {
116 2 : if( mname!=mycolv->getName() ) error("All input multicolvars must be of same type");
117 : }
118 : // And track which variable stores each colvar
119 83 : for(unsigned j=0; j<mycolv->getFullNumberOfTasks(); ++j) atom_lab.push_back( std::pair<unsigned,unsigned>( mybasemulticolvars.size()+1, j ) );
120 : // And store the multicolvar base
121 83 : mybasemulticolvars.push_back( mycolv );
122 : // And create a basedata stash
123 83 : mybasedata.push_back( mybasemulticolvars[mybasemulticolvars.size()-1]->buildDataStashes( this ) );
124 : // Check if weight has derivatives
125 83 : if( mybasemulticolvars[ mybasemulticolvars.size()-1 ]->weightHasDerivatives ) weightHasDerivatives=true;
126 83 : plumed_assert( mybasemulticolvars.size()==mybasedata.size() );
127 : }
128 : // Have we conventional atoms to read in
129 220 : if( found_mcolv==0 ) {
130 139 : std::vector<AtomNumber> tt; ActionAtomistic::interpretAtomList( mlabs, tt );
131 139 : for(unsigned i=0; i<tt.size(); ++i) { atom_lab.push_back( std::pair<unsigned,unsigned>( 0, t.size() + i ) ); }
132 139 : log.printf(" keyword %s takes atoms : ", key.c_str() );
133 139 : for(unsigned i=0; i<tt.size(); ++i) { t.push_back( tt[i] ); log.printf("%d ",tt[i].serial() ); }
134 139 : log.printf("\n");
135 81 : } else if( found_mcolv==mlabs.size() ) {
136 81 : if( checkNumericalDerivatives() ) error("cannot use NUMERICAL_DERIVATIVES keyword with dynamic groups as input");
137 81 : log.printf(" keyword %s takes dynamic groups of atoms constructed from multicolvars labelled : ", key.c_str() );
138 81 : for(unsigned i=0; i<mlabs.size(); ++i) log.printf("%s ",mlabs[i].c_str() );
139 81 : log.printf("\n");
140 0 : } else if( found_mcolv<mlabs.size() ) {
141 0 : error("cannot mix multicolvar input and atom input in one line");
142 : }
143 570 : return true;
144 : }
145 :
146 68 : void MultiColvarBase::readTwoGroups( const std::string& key0, const std::string& key1, const std::string& key2, std::vector<AtomNumber>& all_atoms ) {
147 68 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) ); ablocks.resize( 2 );
148 :
149 68 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
150 20 : nblock=atom_lab.size(); for(unsigned i=0; i<2; ++i) ablocks[i].resize(nblock);
151 20 : resizeBookeepingArray( nblock, nblock );
152 20 : for(unsigned i=0; i<nblock; ++i) ablocks[0][i]=ablocks[1][i]=i;
153 5430 : for(unsigned i=1; i<nblock; ++i) {
154 4216245 : for(unsigned j=0; j<i; ++j) {
155 4210835 : bookeeping(i,j).first=getFullNumberOfTasks();
156 4210835 : addTaskToList( i*nblock + j );
157 4210835 : bookeeping(i,j).second=getFullNumberOfTasks();
158 : }
159 : }
160 : } else {
161 48 : parseMultiColvarAtomList(key1,-1,all_atoms);
162 48 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=i;
163 48 : parseMultiColvarAtomList(key2,-1,all_atoms);
164 48 : 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;
165 :
166 48 : if( ablocks[0].size()>ablocks[1].size() ) nblock = ablocks[0].size();
167 48 : else nblock=ablocks[1].size();
168 :
169 48 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
170 65 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
171 1126 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
172 1109 : bookeeping(i,j).first=getFullNumberOfTasks();
173 1109 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 ) {
174 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
175 0 : atom_lab[ablocks[0][i]].second!=atom_lab[ablocks[1][j]].second ) addTaskToList( i*nblock + j );
176 1109 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] ) addTaskToList( i*nblock + j );
177 1109 : bookeeping(i,j).second=getFullNumberOfTasks();
178 : }
179 : }
180 : }
181 68 : }
182 :
183 4 : void MultiColvarBase::readGroupKeywords( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& key3,
184 : const bool& no_third_dim_accum, const bool& symmetric, std::vector<AtomNumber>& all_atoms ) {
185 4 : plumed_dbg_assert( keywords.exists(key0) && keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) ); ablocks.resize( 3 );
186 :
187 4 : if( parseMultiColvarAtomList(key0,-1,all_atoms) ) {
188 2 : if( no_third_dim_accum ) {
189 2 : nblock=atom_lab.size(); ablocks[0].resize(nblock); ablocks[1].resize( nblock );
190 2 : for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=ablocks[1][i]=i;
191 2 : resizeBookeepingArray( nblock, nblock );
192 2 : if( symmetric ) {
193 0 : for(unsigned i=1; i<nblock; ++i) {
194 0 : for(unsigned j=0; j<i; ++j) {
195 0 : bookeeping(i,j).first=getFullNumberOfTasks();
196 0 : addTaskToList( i*nblock + j );
197 0 : bookeeping(i,j).second=getFullNumberOfTasks();
198 : }
199 : }
200 : } else {
201 69 : for(unsigned i=0; i<nblock; ++i) {
202 4172 : for(unsigned j=0; j<nblock; ++j) {
203 4105 : if( i==j ) continue ;
204 4038 : bookeeping(i,j).first=getFullNumberOfTasks();
205 4038 : addTaskToList( i*nblock + j );
206 4038 : bookeeping(i,j).second=getFullNumberOfTasks();
207 : }
208 : }
209 : }
210 2 : if( !parseMultiColvarAtomList(key3,-1,all_atoms) ) error("missing required keyword " + key3 + " in input");
211 2 : ablocks[2].resize( atom_lab.size() - ablocks[0].size() );
212 2 : for(unsigned i=0; i<ablocks[2].size(); ++i) ablocks[2][i]=ablocks[0].size() + i;
213 : } else {
214 0 : nblock=atom_lab.size(); for(unsigned i=0; i<3; ++i) ablocks[i].resize(nblock);
215 0 : resizeBookeepingArray( nblock, nblock );
216 0 : for(unsigned i=0; i<nblock; ++i) { ablocks[0][i]=i; ablocks[1][i]=i; ablocks[2][i]=i; }
217 0 : if( symmetric ) {
218 0 : for(unsigned i=2; i<nblock; ++i) {
219 0 : for(unsigned j=1; j<i; ++j) {
220 0 : bookeeping(i,j).first=getFullNumberOfTasks();
221 0 : for(unsigned k=0; k<j; ++k) addTaskToList( i*nblock*nblock + j*nblock + k );
222 0 : bookeeping(i,j).second=getFullNumberOfTasks();
223 : }
224 : }
225 : } else {
226 0 : for(unsigned i=0; i<nblock; ++i) {
227 0 : for(unsigned j=0; j<nblock; ++j) {
228 0 : if( i==j ) continue;
229 0 : bookeeping(i,j).first=getFullNumberOfTasks();
230 0 : for(unsigned k=0; k<nblock; ++k) {
231 0 : if( i!=k && j!=k ) addTaskToList( i*nblock*nblock + j*nblock + k );
232 : }
233 0 : bookeeping(i,j).first=getFullNumberOfTasks();
234 : }
235 : }
236 : }
237 : }
238 : } else {
239 2 : readThreeGroups( key1, key2, key3, true, no_third_dim_accum, all_atoms );
240 : }
241 :
242 4 : }
243 :
244 11 : void MultiColvarBase::readThreeGroups( const std::string& key1, const std::string& key2, const std::string& key3,
245 : const bool& allow2, const bool& no_third_dim_accum, std::vector<AtomNumber>& all_atoms ) {
246 11 : plumed_dbg_assert( keywords.exists(key1) && keywords.exists(key2) && keywords.exists(key3) ); ablocks.resize( 3 );
247 :
248 11 : bool readkey1=parseMultiColvarAtomList(key1,-1,all_atoms);
249 11 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<ablocks[0].size(); ++i) ablocks[0][i]=i;
250 11 : bool readkey2=parseMultiColvarAtomList(key2,-1,all_atoms);
251 22 : if( !readkey1 && !readkey2 ) return ;
252 11 : 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;
253 :
254 11 : resizeBookeepingArray( ablocks[0].size(), ablocks[1].size() );
255 11 : bool readkey3=parseMultiColvarAtomList(key3,-1,all_atoms);
256 11 : if( !readkey3 && !allow2 ) {
257 0 : error("missing atom specification " + key3);
258 11 : } else if( !readkey3 ) {
259 2 : if( ablocks[1].size()>ablocks[0].size() ) nblock=ablocks[1].size();
260 0 : else nblock=ablocks[0].size();
261 :
262 2 : ablocks[2].resize( ablocks[1].size() );
263 2 : for(unsigned i=0; i<ablocks[1].size(); ++i) ablocks[2][i]=ablocks[0].size() + i;
264 4 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
265 198 : for(unsigned j=1; j<ablocks[1].size(); ++j) {
266 196 : bookeeping(i,j).first=getFullNumberOfTasks();
267 9898 : for(unsigned k=0; k<j; ++k) {
268 9702 : if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
269 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
270 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
271 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
272 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 &&
273 0 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) addTaskToList( nblock*nblock*i + nblock*j + k );
274 29106 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
275 19404 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
276 19404 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) addTaskToList( nblock*nblock*i + nblock*j + k );
277 : }
278 196 : bookeeping(i,j).second=getFullNumberOfTasks();
279 : }
280 : }
281 : } else {
282 9 : ablocks[2].resize( atom_lab.size() - ablocks[1].size() - ablocks[0].size() );
283 9 : for(unsigned i=0; i<ablocks[2].size(); ++i) ablocks[2][i] = ablocks[0].size() + ablocks[1].size() + i;
284 :
285 9 : if( ablocks[1].size()>ablocks[0].size() ) nblock=ablocks[1].size();
286 9 : else nblock=ablocks[0].size();
287 9 : if( ablocks[2].size()>nblock ) nblock=ablocks[2].size();
288 :
289 9 : unsigned kcount; if( no_third_dim_accum ) kcount=1; else kcount=ablocks[2].size();
290 :
291 38 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
292 238 : for(unsigned j=0; j<ablocks[1].size(); ++j) {
293 209 : bookeeping(i,j).first=getFullNumberOfTasks();
294 418 : for(unsigned k=0; k<kcount; ++k) {
295 209 : if( no_third_dim_accum ) addTaskToList( nblock*i + j );
296 0 : else if( atom_lab[ablocks[0][i]].first>0 && atom_lab[ablocks[1][j]].first>0 && atom_lab[ablocks[2][k]].first>0 ) {
297 0 : if( mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel() &&
298 0 : mybasemulticolvars[atom_lab[ablocks[0][i]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
299 0 : mybasemulticolvars[atom_lab[ablocks[1][j]].first-1]->getLabel()!=mybasemulticolvars[atom_lab[ablocks[2][k]].first-1]->getLabel() &&
300 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 &&
301 0 : atom_lab[ablocks[1][j]].second!=atom_lab[ablocks[2][k]].second ) addTaskToList( nblock*nblock*i + nblock*j + k );
302 0 : } else if( all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[1][j]].second] &&
303 0 : all_atoms[atom_lab[ablocks[0][i]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] &&
304 0 : all_atoms[atom_lab[ablocks[1][j]].second]!=all_atoms[atom_lab[ablocks[2][k]].second] ) addTaskToList( nblock*nblock*i + nblock*j + k );
305 : }
306 209 : bookeeping(i,j).second=getFullNumberOfTasks();
307 : }
308 : }
309 : }
310 : }
311 :
312 12264184 : void MultiColvarBase::addTaskToList( const unsigned& taskCode ) {
313 12264184 : plumed_assert( getNumberOfVessels()==0 );
314 12264184 : ActionWithVessel::addTaskToList( taskCode );
315 12264184 : }
316 :
317 81 : void MultiColvarBase::resizeBookeepingArray( const unsigned& num1, const unsigned& num2 ) {
318 81 : bookeeping.resize( num1, num2 );
319 5626 : for(unsigned i=0; i<num1; ++i) {
320 5545 : for(unsigned j=0; j<num2; ++j) { bookeeping(i,j).first=0; bookeeping(i,j).second=0; }
321 : }
322 81 : }
323 :
324 221 : void MultiColvarBase::setupMultiColvarBase( const std::vector<AtomNumber>& atoms ) {
325 221 : if( !matsums && atom_lab.size()==0 ) error("No atoms have been read in");
326 221 : std::vector<AtomNumber> all_atoms;
327 : // Setup decoder array
328 221 : if( !usespecies && nblock>0 ) {
329 :
330 46 : ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
331 46 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
332 46 : if( ablocks.size()==3 ) {
333 13 : allthirdblockintasks=uselinkforthree=true;
334 111 : for(unsigned i=0; i<bookeeping.nrows(); ++i) {
335 4610 : for(unsigned j=0; j<bookeeping.ncols(); ++j) {
336 4512 : unsigned ntper = bookeeping(i,j).second - bookeeping(i,j).first;
337 4512 : if( i==j && ntper==0 ) {
338 69 : continue;
339 4443 : } else if( ntper == 1 && allthirdblockintasks ) {
340 4249 : allthirdblockintasks=true;
341 194 : } else if( ntper != ablocks[2].size() ) {
342 194 : allthirdblockintasks=uselinkforthree=false;
343 : } else {
344 0 : allthirdblockintasks=false;
345 : }
346 : }
347 : }
348 : }
349 :
350 46 : if( allthirdblockintasks ) {
351 11 : decoder.resize(2); plumed_assert( ablocks.size()==3 );
352 : // Check if number of atoms is too large
353 11 : if( pow( double(nblock), 2.0 )>std::numeric_limits<unsigned>::max() ) error("number of atoms in groups is too big for PLUMED to handle");
354 : } else {
355 35 : decoder.resize( ablocks.size() );
356 : // Check if number of atoms is too large
357 35 : 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");
358 : }
359 46 : unsigned code=1; for(unsigned i=0; i<decoder.size(); ++i) { decoder[decoder.size()-1-i]=code; code *= nblock; }
360 175 : } else if( !usespecies ) {
361 76 : ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
362 76 : numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
363 99 : } else if( keywords.exists("SPECIESA") ) {
364 70 : plumed_assert( atom_lab.size()==0 && all_atoms.size()==0 );
365 70 : ablocks.resize( 1 ); bool readspecies=parseMultiColvarAtomList("SPECIES", -1, all_atoms);
366 70 : if( readspecies ) {
367 62 : ablocks[0].resize( atom_lab.size() ); for(unsigned i=0; i<atom_lab.size(); ++i) { addTaskToList(i); ablocks[0][i]=i; }
368 : } else {
369 8 : if( !parseMultiColvarAtomList("SPECIESA", -1, all_atoms) ) error("missing SPECIES/SPECIESA keyword");
370 8 : unsigned nat1=atom_lab.size();
371 8 : if( !parseMultiColvarAtomList("SPECIESB", -1, all_atoms) ) error("missing SPECIESB keyword");
372 8 : unsigned nat2=atom_lab.size() - nat1;
373 :
374 8 : for(unsigned i=0; i<nat1; ++i) addTaskToList(i);
375 8 : ablocks[0].resize( nat2 );
376 718 : for(unsigned i=0; i<nat2; ++i) {
377 710 : bool found=false; unsigned inum;
378 6241 : for(unsigned j=0; j<nat1; ++j) {
379 5597 : if( atom_lab[nat1+i].first>0 && atom_lab[j].first>0 ) {
380 8160 : if( atom_lab[nat1+i].first==atom_lab[j].first &&
381 2720 : mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
382 0 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) { found=true; inum=j; break; }
383 2877 : } else if( atom_lab[nat1+i].first>0 ) {
384 0 : if( mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
385 0 : all_atoms[atom_lab[j].second] ) { found=true; inum=nat1 + i; break; }
386 2877 : } else if( atom_lab[j].first>0 ) {
387 0 : if( all_atoms[atom_lab[nat1+i].second]==
388 0 : mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ) { found=true; inum=nat1+i; break; }
389 2877 : } else if( all_atoms[atom_lab[nat1+i].second]==all_atoms[atom_lab[j].second] ) { found=true; inum=j; break; }
390 : }
391 : // This prevents mistakes being made in colvar setup
392 710 : if( found ) { ablocks[0][i]=inum; }
393 644 : else { ablocks[0][i]=nat1 + i; }
394 : }
395 : }
396 : }
397 221 : if( mybasemulticolvars.size()>0 ) { for(unsigned i=0; i<mybasedata.size(); ++i) mybasedata[i]->resizeTemporyMultiValues(2); }
398 :
399 : // Copy lists of atoms involved from base multicolvars
400 442 : std::vector<AtomNumber> tmp_atoms;
401 304 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
402 83 : tmp_atoms=mybasemulticolvars[i]->getAbsoluteIndexes();
403 83 : for(unsigned j=0; j<tmp_atoms.size(); ++j) all_atoms.push_back( tmp_atoms[j] );
404 : }
405 : // Copy atom lists from input
406 221 : if( !usespecies ) {
407 122 : for(unsigned i=0; i<atoms.size(); ++i) all_atoms.push_back( atoms[i] );
408 : }
409 :
410 : // Now make sure we get all the atom positions
411 221 : ActionAtomistic::requestAtoms( all_atoms );
412 : // And setup dependencies
413 221 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) addDependency( mybasemulticolvars[i] );
414 :
415 : // Setup underlying ActionWithVessel
416 442 : readVesselKeywords();
417 221 : }
418 :
419 19 : void MultiColvarBase::setAtomsForCentralAtom( const std::vector<bool>& catom_ind ) {
420 19 : unsigned nat=0; plumed_assert( catom_ind.size()==ablocks.size() );
421 88 : for(unsigned i=0; i<catom_ind.size(); ++i) {
422 69 : use_for_central_atom[i]=catom_ind[i];
423 69 : if( use_for_central_atom[i] ) nat++;
424 : }
425 19 : plumed_dbg_assert( nat>0 ); ncentral=nat;
426 19 : numberForCentralAtom = 1.0 / static_cast<double>( nat );
427 19 : }
428 :
429 339 : void MultiColvarBase::turnOnDerivatives() {
430 339 : ActionWithValue::turnOnDerivatives();
431 339 : needsDerivatives();
432 339 : forcesToApply.resize( getNumberOfDerivatives() );
433 339 : }
434 :
435 146 : void MultiColvarBase::setLinkCellCutoff( const double& lcut, double tcut ) {
436 146 : plumed_assert( usespecies || ablocks.size()<4 );
437 146 : if( tcut<0 ) tcut=lcut;
438 146 : linkcells.setCutoff( lcut );
439 146 : threecells.setCutoff( tcut );
440 146 : }
441 :
442 4 : double MultiColvarBase::getLinkCellCutoff() const {
443 4 : return linkcells.getCutoff();
444 : }
445 :
446 1509 : void MultiColvarBase::setupLinkCells() {
447 3018 : if( (!usespecies && nblock==0) || !linkcells.enabled() ) return ;
448 : // Retrieve any atoms that haven't already been retrieved
449 900 : for(std::vector<MultiColvarBase*>::iterator p=mybasemulticolvars.begin(); p!=mybasemulticolvars.end(); ++p) {
450 168 : (*p)->retrieveAtoms();
451 : }
452 732 : retrieveAtoms();
453 :
454 : unsigned iblock;
455 732 : if( usespecies ) {
456 556 : iblock=0;
457 176 : } else if( ablocks.size()<4 ) {
458 176 : iblock=1;
459 : } else {
460 0 : plumed_error();
461 : }
462 :
463 : // Count number of currently active atoms
464 732 : nactive_atoms=0;
465 197122 : for(unsigned i=0; i<ablocks[iblock].size(); ++i) {
466 196390 : if( isCurrentlyActive( ablocks[iblock][i] ) ) nactive_atoms++;
467 : }
468 :
469 732 : if( nactive_atoms>0 ) {
470 732 : std::vector<Vector> ltmp_pos( nactive_atoms );
471 1464 : std::vector<unsigned> ltmp_ind( nactive_atoms );
472 :
473 732 : nactive_atoms=0;
474 732 : if( usespecies ) {
475 177372 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
476 176816 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
477 176720 : ltmp_ind[nactive_atoms]=ablocks[0][i];
478 176720 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ltmp_ind[nactive_atoms] );
479 176720 : nactive_atoms++;
480 : }
481 : } else {
482 19750 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
483 19574 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
484 14654 : ltmp_ind[nactive_atoms]=i;
485 14654 : ltmp_pos[nactive_atoms]=getPositionOfAtomForLinkCells( ablocks[1][i] );
486 14654 : nactive_atoms++;
487 : }
488 : }
489 :
490 : // Build the lists for the link cells
491 1464 : linkcells.buildCellLists( ltmp_pos, ltmp_ind, getPbc() );
492 : }
493 : }
494 :
495 5214 : void MultiColvarBase::setupNonUseSpeciesLinkCells( const unsigned& my_always_active ) {
496 5214 : plumed_assert( !usespecies );
497 10428 : if( nblock==0 || !linkcells.enabled() ) return ;
498 664 : deactivateAllTasks();
499 :
500 664 : if( !uselinkforthree && nactive_atoms>0 ) {
501 : // Get some parallel info
502 126 : unsigned stride=comm.Get_size();
503 126 : unsigned rank=comm.Get_rank();
504 126 : if( serialCalculation() ) { stride=1; rank=0; }
505 :
506 : // Ensure we only do tasks where atoms are in appropriate link cells
507 126 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
508 5577 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
509 5451 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
510 2472 : unsigned natomsper=1; linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
511 2472 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), natomsper, linked_atoms );
512 204805 : for(unsigned j=0; j<natomsper; ++j) {
513 202333 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) taskFlags[k]=1;
514 : }
515 126 : }
516 538 : } else if( nactive_atoms>0 ) {
517 : // Get some parallel info
518 538 : unsigned stride=comm.Get_size();
519 538 : unsigned rank=comm.Get_rank();
520 538 : if( serialCalculation() ) { stride=1; rank=0; }
521 :
522 538 : unsigned nactive_three=0;
523 21538 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
524 21000 : if( isCurrentlyActive( ablocks[2][i] ) ) nactive_three++;
525 : }
526 :
527 538 : std::vector<Vector> lttmp_pos( nactive_three );
528 1076 : std::vector<unsigned> lttmp_ind( nactive_three );
529 :
530 538 : nactive_three=0;
531 538 : if( allthirdblockintasks ) {
532 21538 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
533 21000 : if( !isCurrentlyActive( ablocks[2][i] ) ) continue;
534 21000 : lttmp_ind[nactive_three]=ablocks[2][i];
535 21000 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
536 21000 : nactive_three++;
537 : }
538 : } else {
539 0 : for(unsigned i=0; i<ablocks[2].size(); ++i) {
540 0 : if( !isCurrentlyActive( ablocks[2][i] ) ) continue;
541 0 : lttmp_ind[nactive_three]=i;
542 0 : lttmp_pos[nactive_three]=getPositionOfAtomForLinkCells( ablocks[2][i] );
543 0 : nactive_three++;
544 : }
545 : }
546 : // Build the list of the link cells
547 538 : threecells.buildCellLists( lttmp_pos, lttmp_ind, getPbc() );
548 :
549 : // Ensure we only do tasks where atoms are in appropriate link cells
550 1076 : std::vector<unsigned> linked_atoms( 1+ablocks[1].size() );
551 1076 : std::vector<unsigned> tlinked_atoms( 1+ablocks[2].size() );
552 5692 : for(unsigned i=rank; i<ablocks[0].size(); i+=stride) {
553 5154 : if( !isCurrentlyActive( ablocks[0][i] ) ) continue;
554 5154 : unsigned natomsper=1; linked_atoms[0]=my_always_active; // Note we always check atom 0 because it is simpler than changing LinkCells.cpp
555 5154 : linkcells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), natomsper, linked_atoms );
556 5154 : if( allthirdblockintasks ) {
557 59796 : for(unsigned j=0; j<natomsper; ++j) {
558 54642 : for(unsigned k=bookeeping(i,linked_atoms[j]).first; k<bookeeping(i,linked_atoms[j]).second; ++k) taskFlags[k]=1;
559 : }
560 : } else {
561 0 : unsigned ntatomsper=1; tlinked_atoms[0]=lttmp_ind[0];
562 0 : threecells.retrieveNeighboringAtoms( getPositionOfAtomForLinkCells( ablocks[0][i] ), ntatomsper, tlinked_atoms );
563 0 : for(unsigned j=0; j<natomsper; ++j) {
564 0 : for(unsigned k=0; k<ntatomsper; ++k) taskFlags[bookeeping(i,linked_atoms[j]).first+tlinked_atoms[k]]=1;
565 : }
566 : }
567 538 : }
568 : }
569 664 : if( !serialCalculation() ) comm.Sum( taskFlags );
570 664 : lockContributors();
571 : }
572 :
573 269144 : void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const {
574 : plumed_dbg_assert( !usespecies && nblock>0 );
575 269144 : if( atoms.size()!=decoder.size() ) atoms.resize( decoder.size() );
576 :
577 269132 : unsigned scode = taskCode;
578 855761 : for(unsigned i=0; i<decoder.size(); ++i) {
579 586673 : unsigned ind=( scode / decoder[i] );
580 586684 : atoms[i] = ablocks[i][ind];
581 586589 : scode -= ind*decoder[i];
582 : }
583 269123 : }
584 :
585 423461 : bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ) const {
586 423461 : if( isDensity() ) {
587 17146 : myatoms.setNumberOfAtoms( 1 ); myatoms.setAtom( 0, taskCode ); return true;
588 406258 : } else if( usespecies ) {
589 136815 : std::vector<unsigned> task_atoms(1); task_atoms[0]=taskCode;
590 136891 : unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getLinkCellPosition(task_atoms), linkcells );
591 136884 : return natomsper>1;
592 269443 : } else if( matsums ) {
593 303 : myatoms.setNumberOfAtoms( getNumberOfAtoms() );
594 303 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) myatoms.setAtom( i, i );
595 269140 : } else if( allthirdblockintasks ) {
596 54567 : plumed_dbg_assert( ablocks.size()==3 ); std::vector<unsigned> atoms(2); decodeIndexToAtoms( taskCode, atoms );
597 54563 : myatoms.setupAtomsFromLinkCells( atoms, getLinkCellPosition(atoms), threecells );
598 214573 : } else if( nblock>0 ) {
599 179087 : std::vector<unsigned> atoms( ablocks.size() );
600 179111 : decodeIndexToAtoms( taskCode, atoms ); myatoms.setNumberOfAtoms( ablocks.size() );
601 179101 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, atoms[i] );
602 : } else {
603 35486 : myatoms.setNumberOfAtoms( ablocks.size() );
604 35516 : for(unsigned i=0; i<ablocks.size(); ++i) myatoms.setAtom( i, ablocks[i][taskCode] );
605 : }
606 269509 : return true;
607 : }
608 :
609 8781 : void MultiColvarBase::setupActiveTaskSet( std::vector<unsigned>& active_tasks, const std::string& input_label ) {
610 8781 : if( !setup_completed ) {
611 1501 : bool justVolumes=false;
612 1501 : if( usespecies ) {
613 555 : justVolumes=true;
614 817 : for(unsigned i=0; i<getNumberOfVessels(); ++i) {
615 555 : vesselbase::StoreDataVessel* mys=dynamic_cast<vesselbase::StoreDataVessel*>( getPntrToVessel(i) );
616 555 : if( mys ) continue;
617 304 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(i) );
618 304 : if( !myb ) { justVolumes=false; break; }
619 17 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
620 17 : if( !myv ) { justVolumes=false; break; }
621 : }
622 : }
623 1501 : deactivateAllTasks();
624 1501 : if( justVolumes && mydata ) {
625 251 : if( mydata->getNumberOfDataUsers()==0 ) justVolumes=false;
626 :
627 269 : for(unsigned i=0; i<mydata->getNumberOfDataUsers(); ++i) {
628 18 : MultiColvarBase* myu=dynamic_cast<MultiColvarBase*>( mydata->getDataUser(i) );
629 18 : if( myu ) {
630 18 : myu->setupActiveTaskSet( taskFlags, getLabel() );
631 : } else {
632 0 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
633 : }
634 : }
635 : }
636 1501 : if( justVolumes ) {
637 56 : for(unsigned j=0; j<getNumberOfVessels(); ++j) {
638 28 : vesselbase::BridgeVessel* myb=dynamic_cast<vesselbase::BridgeVessel*>( getPntrToVessel(j) );
639 28 : if( !myb ) continue ;
640 11 : ActionVolume* myv=dynamic_cast<ActionVolume*>( myb->getOutputAction() );
641 11 : if( !myv ) continue ;
642 11 : myv->retrieveAtoms(); myv->setupRegions();
643 :
644 61515 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
645 61504 : if( myv->inVolumeOfInterest(i) ) taskFlags[i]=1;
646 : }
647 : }
648 : } else {
649 1473 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
650 : }
651 :
652 : // Now activate all this class
653 1501 : lockContributors();
654 : // Setup the link cells
655 1501 : setupLinkCells();
656 : // Ensures that setup is not performed multiple times during one cycle
657 1501 : setup_completed=true;
658 : }
659 :
660 : // And activate the tasks in input action
661 8781 : if( getLabel()!=input_label ) {
662 18 : int input_code=-1;
663 20 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
664 20 : if( mybasemulticolvars[i]->getLabel()==input_label ) { input_code=i+1; break; }
665 : }
666 :
667 18 : MultiValue my_tvals( getNumberOfQuantities(), getNumberOfDerivatives() );
668 18 : AtomValuePack mytmp_atoms( my_tvals, this );
669 61548 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
670 61530 : if( !taskIsCurrentlyActive(i) ) continue;
671 3123 : setupCurrentAtomList( getTaskCode(i), mytmp_atoms );
672 238948 : for(unsigned j=0; j<mytmp_atoms.getNumberOfAtoms(); ++j) {
673 235825 : unsigned itask=mytmp_atoms.getIndex(j);
674 235825 : if( atom_lab[itask].first==input_code ) active_tasks[ atom_lab[itask].second ]=1;
675 : }
676 18 : }
677 : }
678 8781 : }
679 :
680 326 : bool MultiColvarBase::filtersUsedAsInput() {
681 326 : bool inputAreFilters=false;
682 508 : for(unsigned i=0; i<mybasemulticolvars.size(); ++i) {
683 182 : MultiColvarFilter* myfilt=dynamic_cast<MultiColvarFilter*>( mybasemulticolvars[i] );
684 182 : if( myfilt || mybasemulticolvars[i]->filtersUsedAsInput() ) inputAreFilters=true;
685 : }
686 326 : return inputAreFilters;
687 : }
688 :
689 8763 : void MultiColvarBase::calculate() {
690 : // Recursive function that sets up tasks
691 8763 : setupActiveTaskSet( taskFlags, getLabel() );
692 :
693 : // Check for filters and rerun setup of link cells if there are any
694 8763 : if( mybasemulticolvars.size()>0 && filtersUsedAsInput() ) setupLinkCells();
695 :
696 : // Setup the link cells if we are not using species
697 8763 : if( !usespecies && ablocks.size()>1 ) {
698 : // This loop finds the first active atom, which is always checked because
699 : // of a peculiarity in linkcells
700 5214 : unsigned first_active=std::numeric_limits<unsigned>::max();
701 5348 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
702 5348 : if( !isCurrentlyActive( ablocks[1][i] ) ) continue;
703 : else {
704 5214 : first_active=i; break;
705 : }
706 : }
707 5214 : setupNonUseSpeciesLinkCells( first_active );
708 : }
709 : // And run all tasks
710 8763 : runAllTasks();
711 8763 : }
712 :
713 217 : void MultiColvarBase::calculateNumericalDerivatives( ActionWithValue* a ) {
714 217 : if( mybasemulticolvars.size()>0 ) plumed_merror("cannot calculate numerical derivatives for this quantity");
715 217 : calculateAtomicNumericalDerivatives( this, 0 );
716 217 : }
717 :
718 2471 : void MultiColvarBase::prepare() {
719 2471 : setup_completed=false; atomsWereRetrieved=false;
720 2471 : }
721 :
722 5600 : void MultiColvarBase::retrieveAtoms() {
723 5600 : if( !atomsWereRetrieved ) { ActionAtomistic::retrieveAtoms(); atomsWereRetrieved=true; }
724 5600 : }
725 :
726 36973 : void MultiColvarBase::mergeInputDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
727 : const unsigned& jatom, const std::vector<double>& der,
728 : MultiValue& myder, AtomValuePack& myatoms ) const {
729 36973 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
730 : plumed_dbg_assert( ival<myatoms.getUnderlyingMultiValue().getNumberOfValues() );
731 : plumed_dbg_assert( start<myder.getNumberOfValues() && end<=myder.getNumberOfValues() );
732 : plumed_dbg_assert( der.size()==myder.getNumberOfValues() && jatom<myatoms.getNumberOfAtoms() );
733 : // Convert input atom to local index
734 36973 : unsigned katom = myatoms.getIndex( jatom ); plumed_dbg_assert( katom<atom_lab.size() ); plumed_dbg_assert( atom_lab[katom].first>0 );
735 : // Find base colvar
736 36973 : unsigned mmc=atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
737 : // Get start of indices for this atom
738 36973 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
739 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
740 :
741 36973 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
742 4782403 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
743 4745430 : unsigned jder=myder.getActiveIndex(j);
744 4745430 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
745 4423473 : unsigned kder=basen+jder;
746 110152908 : for(unsigned icomp=start; icomp<end; ++icomp) {
747 105729435 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
748 : }
749 : } else {
750 321957 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
751 7507008 : for(unsigned icomp=start; icomp<end; ++icomp) {
752 7185051 : myvals.addDerivative( ival, kder, der[icomp]*myder.getDerivative( icomp, jder ) );
753 : }
754 : }
755 : }
756 36973 : }
757 :
758 93393 : void MultiColvarBase::addComDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
759 : plumed_dbg_assert( ival<static_cast<int>(myatoms.getUnderlyingMultiValue().getNumberOfValues()) && iatom<myatoms.getNumberOfAtoms() );
760 : // Convert input atom to local index
761 93393 : unsigned katom = myatoms.getIndex( iatom ); plumed_dbg_assert( atom_lab[katom].first>0 );
762 : // Find base colvar
763 93393 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
764 : // Get start of indices for this atom
765 93393 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=(mybasemulticolvars[i]->getNumberOfDerivatives() - 9) / 3;
766 93393 : multicolvar::CatomPack atom0=mybasemulticolvars[mmc]->getCentralAtomPack( basen, atom_lab[katom].second );
767 93393 : myatoms.addComDerivatives( ival, der, atom0 );
768 93393 : }
769 :
770 1114367 : void MultiColvarBase::getInputData( const unsigned& ind, const bool& normed,
771 : const multicolvar::AtomValuePack& myatoms,
772 : std::vector<double>& orient ) const {
773 : // Converint input atom to local index
774 1114367 : unsigned katom = myatoms.getIndex(ind); plumed_dbg_assert( atom_lab[katom].first>0 );
775 : // Find base colvar
776 1114367 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
777 : // Check if orient is the correct size
778 1114367 : if( orient.size()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ) orient.resize( mybasemulticolvars[mmc]->getNumberOfQuantities() );
779 : // Retrieve the value
780 1114367 : mybasedata[mmc]->retrieveValueWithIndex( atom_lab[katom].second, normed, orient );
781 1114367 : }
782 :
783 44762 : MultiValue& MultiColvarBase::getInputDerivatives( const unsigned& iatom, const bool& normed, const multicolvar::AtomValuePack& myatoms ) const {
784 : // Converint input atom to local index
785 44762 : unsigned katom = myatoms.getIndex(iatom); plumed_dbg_assert( atom_lab[katom].first>0 );
786 : // Find base colvar
787 44762 : unsigned mmc = atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
788 44762 : if( usespecies && !normed && iatom==0 ) return mybasedata[mmc]->getTemporyMultiValue(0);
789 :
790 43860 : unsigned oval=0; if( iatom>0 ) oval=1;
791 43860 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(oval);
792 87694 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
793 43834 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
794 26 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
795 : }
796 43860 : mybasedata[mmc]->retrieveDerivatives( atom_lab[katom].second, normed, myder );
797 43860 : return myder;
798 : }
799 :
800 58692014 : void MultiColvarBase::accumulateSymmetryFunction( const int& ival, const unsigned& iatom, const double& val, const Vector& der, const Tensor& vir, multicolvar::AtomValuePack& myatoms ) const {
801 58692014 : plumed_dbg_assert( usespecies ); unsigned katom=myatoms.getIndex(0), jatom=myatoms.getIndex(iatom);
802 58702075 : double weight0=1.0; if( atom_lab[katom].first>0 ) weight0=mybasedata[atom_lab[katom].first-1]->retrieveWeightWithIndex( atom_lab[katom].second );
803 58703001 : double weighti=1.0; if( atom_lab[jatom].first>0 ) weighti=mybasedata[atom_lab[jatom].first-1]->retrieveWeightWithIndex( atom_lab[jatom].second );
804 : // Accumulate the value
805 58704478 : if( ival<0 ) myatoms.getUnderlyingMultiValue().addTemporyValue( weight0*weighti*val );
806 56023828 : else myatoms.addValue( ival, weight0*weighti*val );
807 :
808 : // Return if we don't need derivatives
809 117407317 : if( doNotCalculateDerivatives() ) return ;
810 : // And virial
811 52142426 : if( ival<0 ) myatoms.addTemporyBoxDerivatives( weight0*weighti*vir );
812 49956430 : else myatoms.addBoxDerivatives( ival, weight0*weighti*vir );
813 :
814 : // Add derivatives of central atom
815 52141304 : if( atom_lab[katom].first>0 ) {
816 774 : addComDerivatives( ival, 0, -weight0*weighti*der, myatoms );
817 774 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
818 774 : tmpder[0]=weighti*val; mergeInputDerivatives( ival, 0, 1, 0, tmpder, getInputDerivatives(0, false, myatoms), myatoms );
819 : } else {
820 52140648 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( 0, -der );
821 49954657 : else myatoms.addAtomsDerivatives( ival, 0, -der );
822 : }
823 : // Add derivatives of atom in coordination sphere
824 52140708 : if( atom_lab[jatom].first>0 ) {
825 774 : addComDerivatives( ival, iatom, weight0*weighti*der, myatoms );
826 774 : std::vector<double> tmpder( mybasemulticolvars[atom_lab[katom].first - 1]->getNumberOfQuantities(), 0. );
827 774 : tmpder[0]=weight0*val; mergeInputDerivatives( ival, 0, 1, iatom, tmpder, getInputDerivatives(iatom, false, myatoms), myatoms );
828 : } else {
829 52139904 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
830 49953914 : else myatoms.addAtomsDerivatives( ival, iatom, der );
831 : }
832 : }
833 :
834 2113722 : void MultiColvarBase::addAtomDerivatives( const int& ival, const unsigned& iatom, const Vector& der, multicolvar::AtomValuePack& myatoms ) const {
835 4229133 : if( doNotCalculateDerivatives() ) return ;
836 1843473 : unsigned jatom=myatoms.getIndex(iatom);
837 :
838 1842111 : if( atom_lab[jatom].first>0 ) {
839 91845 : addComDerivatives( ival, iatom, der, myatoms );
840 : } else {
841 1750591 : if( ival<0 ) myatoms.addTemporyAtomsDerivatives( iatom, der );
842 1750591 : else myatoms.addAtomsDerivatives( ival, iatom, der );
843 : }
844 : }
845 :
846 246733 : double MultiColvarBase::calculateWeight( const unsigned& current, const double& weight, AtomValuePack& myvals ) const {
847 246733 : return 1.0;
848 : }
849 :
850 420442 : void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
851 420442 : AtomValuePack myatoms( myvals, this );
852 : // Retrieve the atom list
853 420464 : if( !setupCurrentAtomList( current, myatoms ) ) return;
854 : // Get weight due to dynamic groups
855 420274 : double weight = 1.0;
856 420274 : if( !matsums ) {
857 1958267 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
858 1672232 : if( atom_lab[myatoms.getIndex(i)].first==0 ) continue;
859 : // Only need to do first two atoms for thigns like TopologyMatrix, HbondMatrix, Bridge and so on
860 242089 : if( allthirdblockintasks && i>1 ) break;
861 241089 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
862 242089 : weight *= mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
863 : }
864 133991 : } else if( usespecies ) {
865 133674 : if( atom_lab[myatoms.getIndex(0)].first>0 ) {
866 8104 : if( mybasedata[atom_lab[myatoms.getIndex(0)].first-1]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(0)].second )<epsilon ) weight=0.;
867 : }
868 : }
869 : // Do a quick check on the size of this contribution
870 420334 : double multweight = calculateWeight( current, weight, myatoms );
871 420154 : if( weight*multweight<getTolerance() ) {
872 99739 : updateActiveAtoms( myatoms );
873 99739 : return;
874 : }
875 320456 : myatoms.setValue( 0, weight*multweight );
876 : // Deal with derivatives of weights due to dynamic groups
877 320587 : if( !matsums && !doNotCalculateDerivatives() && mybasemulticolvars.size()>0 ) {
878 37536 : MultiValue& outder=myatoms.getUnderlyingMultiValue(); MultiValue myder(0,0);
879 111315 : for(unsigned i=0; i<myatoms.getNumberOfAtoms(); ++i) {
880 : // Neglect any atoms without differentiable weights
881 73779 : if( atom_lab[myatoms.getIndex(i)].first==0 ) continue;
882 :
883 : // Retrieve derivatives
884 73779 : unsigned mmc = atom_lab[myatoms.getIndex(i)].first - 1;
885 73779 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() || myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
886 37536 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
887 : }
888 73779 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(i)].second, false, myder );
889 :
890 : // Retrieve the prefactor (product of all other weights)
891 73779 : double prefactor = multweight*weight / mybasedata[mmc]->retrieveWeightWithIndex( atom_lab[myatoms.getIndex(i)].second );
892 : // And accumulate the derivatives
893 73779 : for(unsigned j=0; j<myder.getNumberActive(); ++j) { unsigned jder=myder.getActiveIndex(j); outder.addDerivative( 0, jder, prefactor*myder.getDerivative(0,jder) ); }
894 73779 : myder.clearAll();
895 37536 : }
896 : }
897 : // Compute everything
898 320572 : double vv=doCalculation( task_index, myatoms );
899 320536 : myatoms.setValue( 1, vv );
900 320549 : return;
901 : }
902 :
903 224917 : double MultiColvarBase::doCalculation( const unsigned& taskIndex, AtomValuePack& myatoms ) const {
904 224917 : if( usespecies && mybasemulticolvars.size()>0 && atom_lab[myatoms.getIndex(0)].first>0 ) {
905 8026 : unsigned mmc = atom_lab[0].first - 1;
906 8026 : MultiValue& myder=mybasedata[mmc]->getTemporyMultiValue(0);
907 16038 : if( myder.getNumberOfValues()!=mybasemulticolvars[mmc]->getNumberOfQuantities() ||
908 8012 : myder.getNumberOfDerivatives()!=mybasemulticolvars[mmc]->getNumberOfDerivatives() ) {
909 14 : myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
910 : }
911 8026 : mybasedata[mmc]->retrieveDerivatives( atom_lab[myatoms.getIndex(0)].second, false, myder );
912 : }
913 224939 : double val=compute( taskIndex, myatoms ); updateActiveAtoms( myatoms );
914 225152 : return val;
915 : }
916 :
917 520607 : void MultiColvarBase::updateActiveAtoms( AtomValuePack& myatoms ) const {
918 520607 : if( mybasemulticolvars.size()==0 ) myatoms.updateUsingIndices();
919 137920 : else myatoms.updateDynamicList();
920 520763 : }
921 :
922 2201754 : Vector MultiColvarBase::getCentralAtomPos( const unsigned& taskIndex ) {
923 2201754 : unsigned curr=getTaskCode( taskIndex );
924 :
925 2201796 : if( usespecies || isDensity() ) {
926 908927 : return getPositionOfAtomForLinkCells(curr);
927 1292878 : } else if( nblock>0 ) {
928 : // double factor=1.0/static_cast<double>( ablocks.size() );
929 2100 : Vector mypos; mypos.zero();
930 2100 : std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
931 6300 : for(unsigned i=0; i<ablocks.size(); ++i) {
932 4200 : if( use_for_central_atom[i] ) mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(atoms[i]);
933 : }
934 2100 : return mypos;
935 : } else {
936 1290778 : Vector mypos; mypos.zero();
937 5134501 : for(unsigned i=0; i<ablocks.size(); ++i) {
938 3843719 : if( use_for_central_atom[i] ) mypos+=numberForCentralAtom*getPositionOfAtomForLinkCells(ablocks[i][curr]);
939 : }
940 1290779 : return mypos;
941 : }
942 : }
943 :
944 216011 : CatomPack MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsigned& taskIndex ) {
945 216011 : unsigned curr=getTaskCode( taskIndex );
946 :
947 216035 : CatomPack mypack;
948 215949 : if(usespecies) {
949 76673 : mypack.resize(1);
950 76674 : mypack.setIndex( 0, basn + curr );
951 76673 : mypack.setDerivative( 0, Tensor::identity() );
952 139276 : } else if( nblock>0 ) {
953 0 : mypack.resize(ncentral); unsigned k=0;
954 0 : std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
955 0 : for(unsigned i=0; i<ablocks.size(); ++i) {
956 0 : if( use_for_central_atom[i] ) {
957 0 : mypack.setIndex( k, basn + atoms[i] );
958 0 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
959 0 : k++;
960 : }
961 0 : }
962 : } else {
963 139276 : mypack.resize(ncentral); unsigned k=0;
964 311808 : for(unsigned i=0; i<ablocks.size(); ++i) {
965 172541 : if( use_for_central_atom[i] ) {
966 169927 : mypack.setIndex( k, basn + ablocks[i][curr] );
967 169996 : mypack.setDerivative( k, numberForCentralAtom*Tensor::identity() );
968 169919 : k++;
969 : }
970 : }
971 : }
972 215993 : return mypack;
973 : }
974 :
975 2127473 : Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 ) const {
976 2127473 : if(usepbc) { return pbcDistance( vec1, vec2 ); }
977 0 : else { return delta( vec1, vec2 ); }
978 : }
979 :
980 191420 : void MultiColvarBase::applyPbc(std::vector<Vector>& dlist, unsigned int max_index) const {
981 191420 : if (usepbc) pbcApply(dlist, max_index);
982 191417 : }
983 :
984 1536 : void MultiColvarBase::apply() {
985 1536 : if( getForcesFromVessels( forcesToApply ) ) setForcesOnAtoms( forcesToApply );
986 1536 : }
987 :
988 : }
989 2523 : }
|
