Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 : #ifndef __PLUMED_colvar_MultiColvarTemplate_h
23 : #define __PLUMED_colvar_MultiColvarTemplate_h
24 :
25 : #include "core/ActionWithVector.h"
26 : #include "core/ParallelTaskManager.h"
27 : #include "tools/ColvarOutput.h"
28 : #include "ColvarInput.h"
29 :
30 : namespace PLMD {
31 :
32 : class Colvar;
33 :
34 : namespace colvar {
35 :
36 : struct MultiColvarInput {
37 : bool usepbc;
38 : unsigned mode;
39 : unsigned nindices_per_task;
40 : //this is so simple that it better to use the Rule of 0
41 : //and the openacc helpers
42 : void toACCDevice()const {
43 : #pragma acc enter data copyin(this[0:1], usepbc, mode, nindices_per_task)
44 : }
45 : void removeFromACCDevice() const {
46 : #pragma acc exit data delete(nindices_per_task, mode, usepbc, this[0:1])
47 : }
48 : };
49 :
50 : template <class CV, typename myPTM=defaultPTM>
51 : class MultiColvarTemplate : public ActionWithVector {
52 : public:
53 : using input_type = MultiColvarInput;
54 : using mytype=MultiColvarTemplate<CV,myPTM>;
55 : using PTM =typename myPTM::template PTM<mytype>;
56 : using precision = typename cvprecision<CV>::type;
57 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
58 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
59 : constexpr static size_t virialSize = 9;
60 : private:
61 : /// The parallel task manager
62 : PTM taskmanager;
63 : /// An index that decides what we are calculating
64 : unsigned mode;
65 : /// Are we using pbc to calculate the CVs
66 : bool usepbc;
67 : /// Do we reassemble the molecule
68 : bool wholemolecules;
69 : /// Variable that is turned true if charges are needed to compute the underlying CVs
70 : /// Add the command:
71 : /// keys.add("hidden","NEEDS_CHARGES","Inform the parallel task manager that we need the charges to compute this CV");
72 : /// in the registerKeywords method for your colvar to turn this to true
73 : bool needs_charges;
74 : /// The number of atoms per task
75 : unsigned natoms_per_task;
76 : public:
77 : static void registerKeywords(Keywords&);
78 : explicit MultiColvarTemplate(const ActionOptions&);
79 : unsigned getNumberOfDerivatives() override ;
80 : void addValueWithDerivatives( const std::vector<std::size_t>& shape=std::vector<std::size_t>() ) override ;
81 : void addComponentWithDerivatives( const std::string& name, const std::vector<std::size_t>& shape=std::vector<std::size_t>() ) override ;
82 : void getAtomicInputData( std::vector<double>& inputdata ) const override ;
83 : void getAtomicInputData( std::vector<float>& inputdata ) const override ;
84 : void calculate() override;
85 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
86 : static void performTask( std::size_t task_index,
87 : const MultiColvarInput& actiondata,
88 : ParallelActionsInput& input,
89 : ParallelActionsOutput& output );
90 : static int getNumberOfValuesPerTask( std::size_t task_index, const MultiColvarInput& actiondata );
91 : static void getForceIndices( std::size_t task_index,
92 : std::size_t colno,
93 : std::size_t ntotal_force,
94 : const MultiColvarInput& actiondata,
95 : const ParallelActionsInput& input,
96 : ForceIndexHolder force_indices );
97 : };
98 :
99 : template <class CV, typename myPTM>
100 681 : void MultiColvarTemplate<CV,myPTM>::registerKeywords(Keywords& keys ) {
101 681 : CV::registerKeywords( keys );
102 681 : PTM::registerKeywords( keys );
103 1362 : keys.addInputKeyword("optional","MASK","vector","the label for a sparse vector that should be used to determine which elements of the vector should be computed");
104 681 : unsigned nkeys = keys.size();
105 8768 : for(unsigned i=0; i<nkeys; ++i) {
106 16174 : if( keys.style( keys.getKeyword(i), "atoms" ) ) {
107 1850 : keys.reset_style( keys.getKeyword(i), "numbered" );
108 : }
109 : }
110 1362 : if( keys.outputComponentExists(".#!value") ) {
111 2160 : keys.setValueDescription("vector","the " + keys.getDisplayName() + " for each set of specified atoms");
112 : }
113 681 : }
114 :
115 : template <class CV, typename myPTM>
116 321 : MultiColvarTemplate<CV,myPTM>::MultiColvarTemplate(const ActionOptions&ao):
117 : Action(ao),
118 : ActionWithVector(ao),
119 321 : taskmanager(this),
120 321 : mode(0),
121 321 : usepbc(true),
122 321 : wholemolecules(false),
123 321 : needs_charges(keywords.exists("NEEDS_CHARGES")) {
124 321 : std::vector<AtomNumber> all_atoms;
125 321 : if( getName()=="POSITION_VECTOR" || getName()=="MASS_VECTOR" || getName()=="CHARGE_VECTOR" ) {
126 96 : parseAtomList( "ATOMS", all_atoms );
127 : }
128 321 : if( all_atoms.size()>0 ) {
129 46 : natoms_per_task=1;
130 : } else {
131 275 : std::vector<AtomNumber> t;
132 232191 : for(int i=1;; ++i ) {
133 232191 : CV::parseAtomList( i, t, this );
134 232191 : if( t.empty() ) {
135 : break;
136 : }
137 :
138 231916 : if( i==1 ) {
139 275 : natoms_per_task=t.size();
140 : }
141 231916 : if( t.size()!=natoms_per_task ) {
142 : std::string ss;
143 0 : Tools::convert(i,ss);
144 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
145 : }
146 698316 : for(unsigned j=0; j<natoms_per_task; ++j) {
147 466400 : all_atoms.push_back( t[j] );
148 : }
149 231916 : t.resize(0);
150 : }
151 275 : }
152 321 : if( all_atoms.size()==0 ) {
153 0 : error("No atoms have been specified");
154 : }
155 321 : requestAtoms(all_atoms);
156 321 : if( keywords.exists("NOPBC") ) {
157 321 : bool nopbc=!usepbc;
158 321 : parseFlag("NOPBC",nopbc);
159 321 : usepbc=!nopbc;
160 : }
161 321 : if( keywords.exists("WHOLEMOLECULES") ) {
162 48 : parseFlag("WHOLEMOLECULES",wholemolecules);
163 48 : if( wholemolecules ) {
164 1 : usepbc=false;
165 : }
166 : }
167 321 : if( usepbc ) {
168 278 : log.printf(" using periodic boundary conditions\n");
169 : } else {
170 43 : log.printf(" without periodic boundary conditions\n");
171 : }
172 :
173 : // Setup the values
174 321 : mode = CV::getModeAndSetupValues( this );
175 : // This sets up an array in the parallel task manager to hold all the indices
176 : // Sets up the index list in the task manager
177 321 : taskmanager.setupParallelTaskManager( 3*natoms_per_task + virialSize, virialSize );
178 321 : taskmanager.setActionInput( MultiColvarInput{ usepbc, mode, natoms_per_task });
179 321 : }
180 :
181 : template <class CV, typename myPTM>
182 2560 : unsigned MultiColvarTemplate<CV,myPTM>::getNumberOfDerivatives() {
183 2560 : return 3*getNumberOfAtoms()+9;
184 : }
185 :
186 : template <class CV, typename myPTM>
187 31974 : void MultiColvarTemplate<CV,myPTM>::calculate() {
188 31974 : if( wholemolecules ) {
189 1 : makeWhole();
190 : }
191 31974 : taskmanager.runAllTasks();
192 31974 : }
193 :
194 : template <class CV, typename myPTM>
195 24168 : void MultiColvarTemplate<CV,myPTM>::applyNonZeroRankForces( std::vector<double>& outforces ) {
196 24168 : taskmanager.applyForces( outforces );
197 24168 : }
198 :
199 : template <class CV, typename myPTM>
200 184 : void MultiColvarTemplate<CV,myPTM>::addValueWithDerivatives( const std::vector<std::size_t>& shape ) {
201 184 : std::vector<std::size_t> s(1);
202 184 : s[0]=getNumberOfAtoms() / natoms_per_task;
203 184 : addValue( s );
204 184 : }
205 :
206 : template <class T, typename myPTM>
207 425 : void MultiColvarTemplate<T,myPTM>::addComponentWithDerivatives( const std::string& compName, const std::vector<std::size_t>& shape ) {
208 425 : std::vector<std::size_t> s(1);
209 425 : s[0]=getNumberOfAtoms() / natoms_per_task;
210 425 : addComponent( compName, s );
211 425 : }
212 :
213 : template <class CV, typename myPTM>
214 31974 : void MultiColvarTemplate<CV,myPTM>::getAtomicInputData( std::vector<double>& inputdata ) const {
215 31974 : std::size_t ntasks = getConstPntrToComponent(0)->getNumberOfStoredValues();
216 31974 : if( inputdata.size()!=5*natoms_per_task*ntasks ) {
217 300 : inputdata.resize( 5*natoms_per_task*ntasks );
218 : }
219 :
220 : std::size_t k=0;
221 1727128 : for(unsigned i=0; i<ntasks; ++i) {
222 4953888 : for(unsigned j=0; j<natoms_per_task; ++j) {
223 3258734 : Vector mypos( getPosition( natoms_per_task*i + j ) );
224 3258734 : inputdata[k] = mypos[0];
225 3258734 : k++;
226 3258734 : inputdata[k] = mypos[1];
227 3258734 : k++;
228 3258734 : inputdata[k] = mypos[2];
229 3258734 : k++;
230 : }
231 4953888 : for(unsigned j=0; j<natoms_per_task; ++j) {
232 3258734 : inputdata[k] = getMass( natoms_per_task*i + j );
233 3258734 : k++;
234 : }
235 1695154 : if( !chargesWereSet && !needs_charges ) {
236 15 : k += natoms_per_task;
237 15 : continue ;
238 : }
239 4953843 : for(unsigned j=0; j<natoms_per_task; ++j) {
240 3258704 : inputdata[k] = getCharge( natoms_per_task*i + j );
241 3258704 : k++;
242 : }
243 : }
244 31974 : }
245 :
246 : template <class CV, typename myPTM>
247 0 : void MultiColvarTemplate<CV,myPTM>::getAtomicInputData( std::vector<float>& inputdata ) const {
248 0 : std::size_t ntasks = getConstPntrToComponent(0)->getNumberOfStoredValues();
249 0 : if( inputdata.size()!=5*natoms_per_task*ntasks ) {
250 0 : inputdata.resize( 5*natoms_per_task*ntasks );
251 : }
252 :
253 : std::size_t k=0;
254 0 : for(unsigned i=0; i<ntasks; ++i) {
255 0 : for(unsigned j=0; j<natoms_per_task; ++j) {
256 0 : Vector mypos( getPosition( natoms_per_task*i + j ) );
257 0 : inputdata[k] = mypos[0];
258 0 : k++;
259 0 : inputdata[k] = mypos[1];
260 0 : k++;
261 0 : inputdata[k] = mypos[2];
262 0 : k++;
263 : }
264 0 : for(unsigned j=0; j<natoms_per_task; ++j) {
265 0 : inputdata[k] = getMass( natoms_per_task*i + j );
266 0 : k++;
267 : }
268 0 : if( !chargesWereSet && !needs_charges ) {
269 0 : k += natoms_per_task;
270 0 : continue ;
271 : }
272 0 : for(unsigned j=0; j<natoms_per_task; ++j) {
273 0 : inputdata[k] = getCharge( natoms_per_task*i + j );
274 0 : k++;
275 : }
276 : }
277 0 : }
278 :
279 : template <class CV, typename myPTM>
280 2724769 : void MultiColvarTemplate<CV,myPTM>::performTask( std::size_t task_index,
281 : const MultiColvarInput& actiondata,
282 : ParallelActionsInput& input,
283 : ParallelActionsOutput& output ) {
284 2724769 : std::size_t pos_start = 5*actiondata.nindices_per_task*task_index;
285 : using V3=VectorT<precision>;
286 2724769 : if( actiondata.usepbc ) {
287 2288376 : if( actiondata.nindices_per_task==1 ) {
288 : //this may be changed to input.pbc.apply() en mass or only on this one
289 49576 : V3 fpos=input.pbc->distance(V3(0.0,0.0,0.0),
290 49576 : V3(input.inputdata[pos_start],
291 24788 : input.inputdata[pos_start+1],
292 24788 : input.inputdata[pos_start+2]) );
293 24788 : input.inputdata[pos_start] =fpos[0];
294 24788 : input.inputdata[pos_start+1]=fpos[1];
295 24788 : input.inputdata[pos_start+2]=fpos[2];
296 : } else {
297 : //make whole?
298 : std::size_t apos_start = pos_start;
299 : //if accidentaly nindices_per_task is 0, this will work by looping on all possible unsigned integers!!!!
300 4541226 : for(unsigned j=0; j<actiondata.nindices_per_task-1; ++j) {
301 2277638 : V3 first(input.inputdata[apos_start],
302 2277638 : input.inputdata[apos_start+1],
303 2277638 : input.inputdata[apos_start+2]);
304 2277638 : V3 second(input.inputdata[apos_start+3],
305 2277638 : input.inputdata[apos_start+4],
306 2277638 : input.inputdata[apos_start+5]);
307 : //calling the pbc here gives problems
308 2277638 : second=first+input.pbc->distance(first,second);
309 2277638 : input.inputdata[apos_start+3]=second[0];
310 2277638 : input.inputdata[apos_start+4]=second[1];
311 2277638 : input.inputdata[apos_start+5]=second[2];
312 : apos_start += 3;
313 : }
314 : }
315 436393 : } else if( actiondata.nindices_per_task==1 ) {
316 : //isn't this equivalent to x = x-0?
317 : //why this is needed?
318 : V3 fpos=delta(V3(0.0,0.0,0.0),
319 162779 : V3(input.inputdata[pos_start],
320 162779 : input.inputdata[pos_start+1],
321 162779 : input.inputdata[pos_start+2]));
322 : input.inputdata[pos_start]=fpos[0];
323 162779 : input.inputdata[pos_start+1]=fpos[1];
324 162779 : input.inputdata[pos_start+2]=fpos[2];
325 : }
326 2724769 : const size_t mass_start = pos_start + 3*actiondata.nindices_per_task;
327 2724769 : const size_t charge_start = mass_start + actiondata.nindices_per_task;
328 2724769 : const size_t local_ndev = 3*actiondata.nindices_per_task+virialSize;
329 :
330 2724769 : ColvarOutput cvout { output.values,
331 : local_ndev,
332 : output.derivatives.data() };
333 5449538 : CV::calculateCV( ColvarInput{actiondata.mode,
334 : actiondata.nindices_per_task,
335 2724769 : input.inputdata+pos_start,
336 2724769 : input.inputdata+mass_start,
337 2724769 : input.inputdata+charge_start,
338 2724769 : *input.pbc},
339 : cvout );
340 2724769 : }
341 :
342 : template <class CV, typename myPTM>
343 1242940 : int MultiColvarTemplate<CV,myPTM>::getNumberOfValuesPerTask( std::size_t task_index,
344 : const MultiColvarInput& actiondata ) {
345 1242940 : return 1;
346 : }
347 :
348 : template <class CV, typename myPTM>
349 1242940 : void MultiColvarTemplate<CV,myPTM>::getForceIndices( std::size_t task_index,
350 : std::size_t colno,
351 : std::size_t ntotal_force,
352 : const MultiColvarInput& actiondata,
353 : const ParallelActionsInput& input,
354 : ForceIndexHolder force_indices ) {
355 2653789 : for(unsigned i=0; i<input.ncomponents; ++i) {
356 : std::size_t m=0;
357 1410849 : std::size_t base = 3*task_index*actiondata.nindices_per_task;
358 4083385 : for(unsigned j=0; j<actiondata.nindices_per_task; ++j) {
359 2672536 : force_indices.indices[i][m] = base + m;
360 2672536 : ++m;
361 2672536 : force_indices.indices[i][m] = base + m;
362 2672536 : ++m;
363 2672536 : force_indices.indices[i][m] = base + m;
364 2672536 : ++m;
365 : }
366 1410849 : force_indices.threadsafe_derivatives_end[i] = 3*actiondata.nindices_per_task;
367 1410849 : force_indices.indices[i][m+0] = ntotal_force - 9;
368 1410849 : force_indices.indices[i][m+1] = ntotal_force - 8;
369 1410849 : force_indices.indices[i][m+2] = ntotal_force - 7;
370 1410849 : force_indices.indices[i][m+3] = ntotal_force - 6;
371 1410849 : force_indices.indices[i][m+4] = ntotal_force - 5;
372 1410849 : force_indices.indices[i][m+5] = ntotal_force - 4;
373 1410849 : force_indices.indices[i][m+6] = ntotal_force - 3;
374 1410849 : force_indices.indices[i][m+7] = ntotal_force - 2;
375 1410849 : force_indices.indices[i][m+8] = ntotal_force - 1;
376 1410849 : force_indices.tot_indices[i] = 3*actiondata.nindices_per_task + virialSize;
377 : }
378 1242940 : }
379 :
380 : } // namespace colvar
381 : } // namespace PLMD
382 : #endif
|
