Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2020 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #ifndef __PLUMED_secondarystructure_SecondaryStructureBase_h
23 : #define __PLUMED_secondarystructure_SecondaryStructureBase_h
24 :
25 : #include "core/ActionWithVector.h"
26 : #include "core/ActionShortcut.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/ActionSet.h"
29 : #include "core/GenericMolInfo.h"
30 : #include "core/ParallelTaskManager.h"
31 : #include "tools/ColvarOutput.h"
32 : #include "tools/TypeUtils.h"
33 : #include <vector>
34 : #include <type_traits>
35 :
36 : namespace PLMD {
37 : namespace secondarystructure {
38 :
39 : namespace helpers {
40 : //this makes the variable needstype only needed to be declared explicitly
41 : //needsType activate the reading of "TYPE" during the parsing
42 : template <typename T, typename=void>
43 : constexpr bool needsType=false;
44 : //specialization to intercept needsType
45 : template <typename T>
46 : constexpr bool needsType<T,std::void_t<decltype(T::needsType)>> =T::needsType;
47 :
48 : //this makes the variable needstype only needed to be declared explicitly
49 : //needsBondLength activate the reading of "TYPE" during the parsing
50 : template <typename T, typename=void>
51 : constexpr bool needsBondLength=false;
52 : //specialization to intercept needsBondLength
53 : template <typename T>
54 : constexpr bool needsBondLength<T,std::void_t<decltype(T::needsBondLength)>> =T::needsBondLength;
55 : }
56 :
57 : /// Base action for calculating things like AlphRMSD, AntibetaRMSD, etc
58 : template <class T, typename myPTM=defaultPTM>
59 : class SecondaryStructureBase: public ActionWithVector {
60 : public:
61 : using input_type = T;
62 : using mytype=SecondaryStructureBase<T,myPTM>;
63 : using PTM =typename myPTM::template PTM<mytype>;
64 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
65 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
66 : using precision = cvprecision_t<T>;
67 : static constexpr size_t virialSize=9;
68 : static constexpr unsigned customGatherStep=3;
69 : static constexpr unsigned customGatherStopBefore=virialSize;
70 : /// The largest number of atoms in a segment (the size of the position buffer in performTask).
71 : /// ALPHARMSD, ANTIBETARMSD and PARABETARMSD all use segments of exactly this size
72 : static constexpr unsigned maxAtomsInSegment=30;
73 :
74 : private:
75 : PTM taskmanager;
76 : static constexpr bool needsType = helpers::needsType<T>;
77 : static constexpr bool needsBondLength = helpers::needsBondLength<T>;
78 : public:
79 : static void registerKeywords( Keywords& keys );
80 : static void readBackboneAtoms( ActionShortcut* action, PlumedMain& plumed, const std::string& backnames, std::vector<unsigned>& chain_lengths, std::vector<std::string>& all_atoms );
81 : static bool readShortcutWords( std::string& ltmap, ActionShortcut* action );
82 : explicit SecondaryStructureBase(const ActionOptions&);
83 : unsigned getNumberOfDerivatives() override ;
84 : void calculate() override;
85 : void getAtomicInputData( std::vector<double>& inputdata ) const override ;
86 : void getAtomicInputData( std::vector<float>& inputdata ) const override ;
87 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
88 : static void performTask( unsigned task_index, const T& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output );
89 : static int getNumberOfValuesPerTask( std::size_t task_index, const T& actiondata );
90 : static void getForceIndices( std::size_t task_index, std::size_t colno, std::size_t ntotal_force, const T& actiondata, const ParallelActionsInput& input, ForceIndexHolder force_indices );
91 : };
92 :
93 : //this is used to get some of the static non template dependant functionalities for Secondary structure
94 : using SecondaryStructureInput=SecondaryStructureBase<Empty>;
95 :
96 : template <class T,typename myPTM>
97 198 : unsigned SecondaryStructureBase<T,myPTM>::getNumberOfDerivatives() {
98 198 : return 3*getNumberOfAtoms()+virialSize;
99 : }
100 :
101 : template <class T,typename myPTM>
102 60 : bool SecondaryStructureBase<T,myPTM>::readShortcutWords( std::string& ltmap, ActionShortcut* action ) {
103 120 : action->parse("LESS_THAN",ltmap);
104 60 : if( ltmap.length()==0 ) {
105 : std::string nn, mm, d_0, r_0;
106 28 : action->parse("R_0",r_0);
107 14 : if( r_0.length()==0 ) {
108 : // The default is 0.08 nm. The distances are in the units that were set with UNITS, so it has to be converted
109 4 : Tools::convert( 0.08/action->plumed.getUnits().getLength(), r_0 );
110 : }
111 14 : action->parse("NN",nn);
112 14 : action->parse("D_0",d_0);
113 14 : action->parse("MM",mm);
114 70 : ltmap = "RATIONAL R_0=" + r_0 + " D_0=" + d_0 + " NN=" + nn + " MM=" + mm;
115 : return false;
116 : }
117 : return true;
118 : }
119 :
120 : template <class T,typename myPTM>
121 497 : void SecondaryStructureBase<T,myPTM>::registerKeywords( Keywords& keys ) {
122 497 : ActionWithVector::registerKeywords( keys );
123 221 : T::registerKeywords( keys );
124 497 : PTM::registerKeywords( keys );
125 497 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions");
126 994 : keys.addInputKeyword("optional","MASK","vector","a vector which is used to determine which elements of the secondary structure variable should be computed");
127 497 : keys.add("atoms","ATOMS","this is the full list of atoms that we are investigating");
128 497 : keys.add("numbered","SEGMENT","this is the lists of atoms in the segment that are being considered");
129 : if constexpr ( needsBondLength /*keys.getDisplayName()=="SECONDARY_STRUCTURE_DRMSD"*/ ) {
130 148 : keys.add("compulsory","BONDLENGTH","0.17","the length to use for bonds");
131 : }
132 497 : keys.add("numbered","STRUCTURE","the reference structure");
133 : if constexpr ( needsType /*keys.getDisplayName()=="SECONDARY_STRUCTURE_RMSD"*/ ) {
134 73 : keys.add("compulsory","TYPE","OPTIMAL","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE. "
135 : "For more details on the OPTIMAL and SIMPLE methods see RMSD.");
136 : } else if constexpr ( !needsBondLength /*keys.getDisplayName()!="SECONDARY_STRUCTURE_DRMSD"*/ ) {
137 276 : keys.add("compulsory","TYPE","DRMSD","the manner in which RMSD alignment is performed. Should be OPTIMAL, SIMPLE or DRMSD. "
138 : "For more details on the OPTIMAL and SIMPLE methods see RMSD. For more details on the "
139 : "DRMSD method see DRMSD.");
140 : }
141 497 : keys.addFlag("VERBOSE",false,"write a more detailed output");
142 994 : keys.reset_style("VERBOSE","hidden");
143 : if constexpr ( !needsBondLength && !needsType ) {
144 276 : keys.add("residues","RESIDUES","this command is used to specify the set of residues that could conceivably form part of the secondary structure. "
145 : "It is possible to use residues numbers as the various chains and residues should have been identified else using an instance of the "
146 : "MOLINFO action. If you wish to use all the residues from all the chains in your system you can do so by "
147 : "specifying all. Alternatively, if you wish to use a subset of the residues you can specify the particular residues "
148 : "you are interested in as a list of numbers. Please be aware that to form secondary structure elements your chain "
149 : "must contain at least N residues, where N is dependent on the particular secondary structure you are interested in. "
150 : "As such if you define portions of the chain with fewer than N residues the code will crash.");
151 276 : keys.add("optional","LESS_THAN","calculate the number of a residue segments that are within a certain target distance of this secondary structure type. "
152 : "This quantity is calculated using \\f$\\sum_i \\sigma(s_i)\\f$, where \\f$\\sigma(s)\\f$ is a switching function.");
153 276 : keys.add("optional","R_0","The r_0 parameter of the switching function.");
154 276 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
155 276 : keys.add("compulsory","NN","8","The n parameter of the switching function");
156 276 : keys.add("compulsory","MM","12","The m parameter of the switching function");
157 : }
158 497 : keys.addFlag("ALIGN_STRANDS",false,"ensure that the two halves of a beta sheet are not broken by the periodic boundaries before doing alignment");
159 994 : keys.addOutputComponent("struct","default","vector","the vectors containing the rmsd distances between the residues and each of the reference structures");
160 994 : keys.addOutputComponent("lessthan","default","scalar","the number blocks of residues that have an RMSD from the secondary structure that is less than the threshold");
161 497 : keys.needsAction("SECONDARY_STRUCTURE_RMSD");
162 497 : keys.needsAction("SECONDARY_STRUCTURE_DRMSD");
163 497 : keys.needsAction("LESS_THAN");
164 497 : keys.needsAction("SUM");
165 497 : keys.addDOI("10.1021/ct900202f");
166 497 : }
167 :
168 : template <class T,typename myPTM>
169 60 : void SecondaryStructureBase<T,myPTM>::readBackboneAtoms( ActionShortcut* action, PlumedMain& plumed, const std::string& moltype, std::vector<unsigned>& chain_lengths, std::vector<std::string>& all_atoms ) {
170 60 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(action);
171 60 : if( ! moldat ) {
172 0 : action->error("Unable to find MOLINFO in input");
173 : }
174 :
175 60 : std::vector<std::string> resstrings;
176 120 : action->parseVector( "RESIDUES", resstrings );
177 60 : if(resstrings.size()==0) {
178 0 : action->error("residues are not defined, check the keyword RESIDUES");
179 60 : } else if( Tools::caseInSensStringCompare(resstrings[0], "all") ) {
180 : resstrings[0]="all";
181 56 : action->log.printf(" examining all possible secondary structure combinations\n");
182 : } else {
183 4 : action->log.printf(" examining secondary structure in residue positions : %s ",resstrings[0].c_str() );
184 6 : for(unsigned i=1; i<resstrings.size(); ++i) {
185 2 : action->log.printf(", %s",resstrings[i].c_str() );
186 : }
187 4 : action->log.printf("\n");
188 : }
189 60 : std::vector< std::vector<AtomNumber> > backatoms;
190 60 : moldat->getBackbone( resstrings, moltype, backatoms );
191 :
192 60 : chain_lengths.resize( backatoms.size() );
193 758 : for(unsigned i=0; i<backatoms.size(); ++i) {
194 698 : chain_lengths[i]=backatoms[i].size();
195 28908 : for(unsigned j=0; j<backatoms[i].size(); ++j) {
196 : std::string bat_str;
197 28210 : Tools::convert( backatoms[i][j].serial(), bat_str );
198 28210 : all_atoms.push_back( bat_str );
199 : }
200 : }
201 60 : if( all_atoms.empty() ) {
202 0 : action->error("no backbone atoms were found for the residues specified with RESIDUES");
203 : }
204 60 : }
205 :
206 : template <class T,typename myPTM>
207 76 : SecondaryStructureBase<T,myPTM>::SecondaryStructureBase(const ActionOptions&ao):
208 : Action(ao),
209 : ActionWithVector(ao),
210 76 : taskmanager(this) {
211 76 : if( plumed.usingNaturalUnits() ) {
212 0 : error("cannot use this collective variable when using natural units");
213 : }
214 :
215 76 : input_type myinput;
216 76 : parseFlag("NOPBC",myinput.nopbc);
217 84 : std::string alignType="";
218 : if constexpr ( needsType ) {
219 40 : parse("TYPE",alignType);
220 : }
221 76 : log.printf(" distances from secondary structure elements are calculated using %s algorithm\n", alignType.c_str() );
222 76 : log<<" Bibliography "<<plumed.cite("Pietrucci and Laio, J. Chem. Theory Comput. 5, 2197 (2009)");
223 76 : log<<"\n";
224 :
225 76 : parseFlag("ALIGN_STRANDS",myinput.align_strands);
226 76 : bool verbose_output=false;
227 87 : parseFlag("VERBOSE",verbose_output);
228 76 : log.printf(" ensuring atoms 7 and 22 in each residue are not separated by pbc before doing alignment\n");
229 :
230 : // Read in the atoms
231 76 : std::vector<AtomNumber> all_atoms;
232 152 : parseAtomList("ATOMS",all_atoms);
233 76 : if( all_atoms.size()==0 ) {
234 0 : error("no atoms were specified -- use ATOMS");
235 : }
236 76 : requestAtoms( all_atoms );
237 :
238 76 : std::vector<std::vector<unsigned> > colvar_atoms;
239 199232 : for(unsigned i=1;; ++i) {
240 199308 : std::vector<unsigned> newatoms;
241 398616 : if( !parseNumberedVector("SEGMENT",i,newatoms) ) {
242 : break;
243 : }
244 : // The segments are copied into a rectangular matrix and the positions of their atoms into a
245 : // buffer of fixed size, so everything is checked here rather than in performTask
246 199238 : if( newatoms.size()<2 || newatoms.size()>maxAtomsInSegment ) {
247 12 : error("SEGMENT" + std::to_string(i) + " contains " + std::to_string(newatoms.size()) + " atoms. A segment must contain between 2 and "
248 : + std::to_string(maxAtomsInSegment) + " atoms");
249 : }
250 199236 : if( !colvar_atoms.empty() && newatoms.size()!=colvar_atoms[0].size() ) {
251 10 : error("all the SEGMENT keywords must contain the same number of atoms but SEGMENT1 contains "
252 : + std::to_string(colvar_atoms[0].size()) + " atoms and SEGMENT" + std::to_string(i) + " contains " + std::to_string(newatoms.size()) );
253 : }
254 6175664 : for(unsigned ii=0; ii<newatoms.size(); ++ii) {
255 5976432 : if( newatoms[ii]>=all_atoms.size() ) {
256 18 : error("SEGMENT" + std::to_string(i) + " contains index " + std::to_string(newatoms[ii]) + " but only " + std::to_string(all_atoms.size())
257 : + " atoms were specified with ATOMS (indices start from zero)");
258 : }
259 : }
260 199232 : if( verbose_output ) {
261 0 : log.printf(" Secondary structure segment %u contains atoms : ", static_cast<unsigned>(colvar_atoms.size()+1));
262 0 : for(unsigned ii=0; ii<newatoms.size(); ++ii) {
263 0 : log.printf("%d ",all_atoms[newatoms[ii]].serial() );
264 : }
265 0 : log.printf("\n");
266 : }
267 199232 : colvar_atoms.push_back( newatoms );
268 : }
269 70 : if( colvar_atoms.size()==0 ) {
270 0 : error("did not find any SEGMENT keywords in input");
271 : }
272 70 : myinput.colvar_atoms.resize( colvar_atoms.size(), colvar_atoms[0].size() );
273 199300 : for(unsigned i=0; i<colvar_atoms.size(); ++i) {
274 6175650 : for(unsigned j=0; j<colvar_atoms[i].size(); ++j) {
275 5976420 : myinput.colvar_atoms[i][j] = colvar_atoms[i][j];
276 : }
277 : }
278 :
279 34 : double bondlength=0.0;
280 : if constexpr ( needsBondLength ) {
281 34 : parse("BONDLENGTH",bondlength);
282 34 : bondlength=bondlength/getUnits().getLength();
283 : }
284 : // Read in the reference structure
285 144 : for(unsigned ii=1;; ++ii) {
286 165 : std::vector<double> cstruct;
287 330 : if( !parseNumberedVector("STRUCTURE",ii,cstruct) ) {
288 : break ;
289 : }
290 98 : if( cstruct.size()%3!=0 || cstruct.size()/3!=colvar_atoms[0].size() ) {
291 12 : error("STRUCTURE" + std::to_string(ii) + " contains " + std::to_string(cstruct.size()) + " numbers. It must contain the three coordinates of each of the "
292 : + std::to_string(colvar_atoms[0].size()) + " atoms in a segment");
293 : }
294 99 : std::vector<Vector> structure( cstruct.size()/3 );
295 2736 : for(unsigned i=0; i<structure.size(); ++i) {
296 10560 : for(unsigned j=0; j<3; ++j) {
297 7920 : structure[i][j] = 0.1*cstruct[3*i+j]/getUnits().getLength();
298 : }
299 : }
300 96 : myinput.setReferenceStructure( alignType, bondlength, structure );
301 : }
302 :
303 : // And create values to hold everything
304 67 : if( myinput.nstructures==0 ) {
305 1 : error("no reference structure was found in input -- use STRUCTURE");
306 : }
307 66 : if( myinput.align_strands && colvar_atoms[0].size()!=maxAtomsInSegment ) {
308 4 : error("ALIGN_STRANDS assumes that the two strands are the first and second halves of a segment containing " + std::to_string(maxAtomsInSegment)
309 : + " atoms, but the segments contain " + std::to_string(colvar_atoms[0].size()) + " atoms");
310 : }
311 65 : std::vector<std::size_t> shape(1);
312 65 : shape[0]=colvar_atoms.size();
313 65 : if( myinput.nstructures==1 ) {
314 37 : addValue( shape );
315 37 : setNotPeriodic();
316 : } else {
317 : std::string num;
318 84 : for(unsigned i=0; i<myinput.nstructures; ++i) {
319 56 : Tools::convert( i+1, num );
320 56 : addComponent( "struct-" + num, shape );
321 112 : componentIsNotPeriodic( "struct-" + num );
322 : }
323 : }
324 158 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
325 93 : getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
326 : }
327 65 : myinput.nindices_per_task = colvar_atoms[0].size();
328 65 : taskmanager.setupParallelTaskManager( 3*colvar_atoms[0].size() + virialSize, 3*getNumberOfAtoms() + 9 );
329 65 : taskmanager.setActionInput( myinput );
330 188 : }
331 :
332 : template <class T,typename myPTM>
333 534 : void SecondaryStructureBase<T,myPTM>::calculate() {
334 534 : taskmanager.runAllTasks();
335 534 : }
336 :
337 : template <class T,typename myPTM>
338 534 : void SecondaryStructureBase<T,myPTM>::getAtomicInputData( std::vector<double>& inputdata ) const {
339 534 : if( inputdata.size()!=3*getNumberOfAtoms() ) {
340 63 : inputdata.resize( 3*getNumberOfAtoms() );
341 : }
342 :
343 : std::size_t k=0;
344 179646 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
345 179112 : Vector mypos( getPosition(i) );
346 179112 : inputdata[k] = mypos[0];
347 179112 : k++;
348 179112 : inputdata[k] = mypos[1];
349 179112 : k++;
350 179112 : inputdata[k] = mypos[2];
351 179112 : k++;
352 : }
353 534 : }
354 :
355 : template <class T,typename myPTM>
356 0 : void SecondaryStructureBase<T,myPTM>::getAtomicInputData( std::vector<float>& inputdata ) const {
357 0 : if( inputdata.size()!=3*getNumberOfAtoms() ) {
358 0 : inputdata.resize( 3*getNumberOfAtoms() );
359 : }
360 :
361 : std::size_t k=0;
362 0 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
363 0 : Vector mypos( getPosition(i) );
364 0 : inputdata[k] = mypos[0];
365 0 : k++;
366 0 : inputdata[k] = mypos[1];
367 0 : k++;
368 0 : inputdata[k] = mypos[2];
369 0 : k++;
370 : }
371 0 : }
372 :
373 : template <class T,typename myPTM>
374 188954 : void SecondaryStructureBase<T,myPTM>::performTask( unsigned task_index, const T& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output ) {
375 : // Empty is only used to reach the static shortcut helpers (see the
376 : // SecondaryStructureInput alias); this body is never called for it. MSVC
377 : // instantiates member bodies eagerly, hence the guard.
378 : if constexpr ( !std::is_same_v<T,Empty> ) {
379 : // The constructor checks that natoms<=maxAtomsInSegment
380 188954 : std::array<Vector,maxAtomsInSegment> pos;
381 :
382 5851886 : for(unsigned i=0; i<actiondata.natoms; ++i) {
383 5662932 : const unsigned atno = actiondata.colvar_atoms(task_index,i);
384 5662932 : pos[i][0] = input.inputdata[3*atno+0];
385 5662932 : pos[i][1] = input.inputdata[3*atno+1];
386 5662932 : pos[i][2] = input.inputdata[3*atno+2];
387 : }
388 188954 : if( !actiondata.nopbc ) {
389 : // Reconstruct the segment by making each pair of consecutive atoms whole. The two strands of a sheet
390 : // (atoms 0-14 and atoms 15-29) are made whole separately, because the last atom of the first strand
391 : // and the first atom of the second strand are not bonded
392 188954 : const unsigned nstrand = actiondata.align_strands ? actiondata.natoms/2 : actiondata.natoms;
393 565422 : for(unsigned start=0; start<actiondata.natoms; start+=nstrand) {
394 5662932 : for(unsigned i=start; i+1<start+nstrand; ++i) {
395 5286464 : const Vector & first (pos[i]);
396 5286464 : Vector & second (pos[i+1]);
397 5286464 : second=first+input.pbc->distance(first,second);
398 : }
399 : }
400 : // This aligns the two strands if this is required: atom 22 is moved to the periodic image that is
401 : // closest to atom 7 and the rest of the second strand follows
402 188954 : if( actiondata.align_strands ) {
403 187514 : Vector distance=input.pbc->distance( pos[6],pos[21] );
404 : distance += pos[6] - pos[21];
405 3000224 : for(unsigned i=actiondata.natoms/2; i<actiondata.natoms; ++i) {
406 : pos[i]+=distance;
407 : }
408 : }
409 : }
410 :
411 : // Create a holder for the derivatives. The number of derivatives per component has to be the one that
412 : // was given to setupParallelTaskManager (3*natoms+virialSize), and the view must contain just the atoms of
413 : // the segment: the size of the buffer is only an upper limit.
414 188954 : const unsigned rs = actiondata.nstructures;
415 188954 : ColvarOutput rmsd_output( output.values,
416 188954 : 3*actiondata.natoms+virialSize,
417 : output.derivatives.data() );
418 : // And now calculate the DRMSD
419 476202 : for(unsigned i=0; i<rs; ++i) {
420 287248 : T::calculateDistance( i, input.noderiv, actiondata, View{pos.data(),actiondata.natoms}, rmsd_output );
421 : }
422 : }
423 188954 : }
424 :
425 : template <class T,typename myPTM>
426 456 : void SecondaryStructureBase<T,myPTM>::applyNonZeroRankForces( std::vector<double>& outforces ) {
427 456 : taskmanager.applyForces( outforces );
428 456 : }
429 :
430 : template <class T,typename myPTM>
431 92638 : int SecondaryStructureBase<T,myPTM>::getNumberOfValuesPerTask( std::size_t task_index,
432 : const T& actiondata ) {
433 92638 : return 1;
434 : }
435 :
436 : template <class T,typename myPTM>
437 92638 : void SecondaryStructureBase<T,myPTM>::getForceIndices( std::size_t task_index,
438 : std::size_t /* colno */,
439 : std::size_t ntotal_force,
440 : const T& actiondata,
441 : const ParallelActionsInput& input,
442 : ForceIndexHolder force_indices ) {
443 : // Guarded for the same reason as performTask().
444 : if constexpr ( !std::is_same_v<T,Empty> ) {
445 232668 : for(unsigned i=0; i<input.ncomponents; ++i) {
446 : std::size_t m = 0;
447 4337150 : for(unsigned j=0; j<actiondata.nindices_per_task; ++j) {
448 4197120 : std::size_t base = 3*actiondata.colvar_atoms[task_index][j];
449 4197120 : force_indices.indices[i][m] = base + 0;
450 4197120 : ++m;
451 4197120 : force_indices.indices[i][m] = base + 1;
452 4197120 : ++m;
453 4197120 : force_indices.indices[i][m] = base + 2;
454 4197120 : ++m;
455 : }
456 1400300 : for(unsigned n=ntotal_force-virialSize; n<ntotal_force; ++n) {
457 1260270 : force_indices.indices[i][m] = n;
458 1260270 : ++m;
459 : }
460 140030 : force_indices.threadsafe_derivatives_end[i] = 0;
461 140030 : force_indices.tot_indices[i] = m;
462 : }
463 : }
464 92638 : }
465 :
466 : }
467 : }
468 :
469 : #endif
|
