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 T>
51 : class MultiColvarTemplate : public ActionWithVector {
52 : public:
53 : using input_type = MultiColvarInput;
54 : using PTM = ParallelTaskManager<MultiColvarTemplate<T>>;
55 : constexpr static size_t virialSize = 9;
56 : private:
57 : /// The parallel task manager
58 : PTM taskmanager;
59 : /// An index that decides what we are calculating
60 : unsigned mode;
61 : /// Are we using pbc to calculate the CVs
62 : bool usepbc;
63 : /// Do we reassemble the molecule
64 : bool wholemolecules;
65 : /// The number of atoms per task
66 : unsigned natoms_per_task;
67 : public:
68 : static void registerKeywords(Keywords&);
69 : explicit MultiColvarTemplate(const ActionOptions&);
70 : unsigned getNumberOfDerivatives() override ;
71 : void addValueWithDerivatives( const std::vector<std::size_t>& shape=std::vector<std::size_t>() ) override ;
72 : void addComponentWithDerivatives( const std::string& name, const std::vector<std::size_t>& shape=std::vector<std::size_t>() ) override ;
73 : void getInputData( std::vector<double>& inputdata ) const override ;
74 : void calculate() override;
75 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
76 : static void performTask( std::size_t task_index,
77 : const MultiColvarInput& actiondata,
78 : ParallelActionsInput& input,
79 : ParallelActionsOutput& output );
80 : static int getNumberOfValuesPerTask( std::size_t task_index, const MultiColvarInput& actiondata );
81 : static void getForceIndices( std::size_t task_index,
82 : std::size_t colno,
83 : std::size_t ntotal_force,
84 : const MultiColvarInput& actiondata,
85 : const ParallelActionsInput& input,
86 : ForceIndexHolder force_indices );
87 : };
88 :
89 : template <class T>
90 472 : void MultiColvarTemplate<T>::registerKeywords(Keywords& keys ) {
91 472 : T::registerKeywords( keys );
92 472 : PTM::registerKeywords( keys );
93 944 : 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");
94 472 : unsigned nkeys = keys.size();
95 5774 : for(unsigned i=0; i<nkeys; ++i) {
96 10604 : if( keys.style( keys.getKeyword(i), "atoms" ) ) {
97 1296 : keys.reset_style( keys.getKeyword(i), "numbered" );
98 : }
99 : }
100 944 : if( keys.outputComponentExists(".#!value") ) {
101 708 : keys.setValueDescription("vector","the " + keys.getDisplayName() + " for each set of specified atoms");
102 : }
103 472 : }
104 :
105 : template <class T>
106 226 : MultiColvarTemplate<T>::MultiColvarTemplate(const ActionOptions&ao):
107 : Action(ao),
108 : ActionWithVector(ao),
109 226 : taskmanager(this),
110 226 : mode(0),
111 226 : usepbc(true),
112 226 : wholemolecules(false) {
113 : std::vector<AtomNumber> all_atoms;
114 594 : if( getName()=="POSITION_VECTOR" || getName()=="MASS_VECTOR" || getName()=="CHARGE_VECTOR" ) {
115 84 : parseAtomList( "ATOMS", all_atoms );
116 : }
117 226 : if( all_atoms.size()>0 ) {
118 40 : natoms_per_task=1;
119 : } else {
120 : std::vector<AtomNumber> t;
121 192624 : for(int i=1;; ++i ) {
122 192624 : T::parseAtomList( i, t, this );
123 192624 : if( t.empty() ) {
124 : break;
125 : }
126 :
127 192438 : if( i==1 ) {
128 186 : natoms_per_task=t.size();
129 : }
130 192438 : if( t.size()!=natoms_per_task ) {
131 : std::string ss;
132 0 : Tools::convert(i,ss);
133 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
134 : }
135 579625 : for(unsigned j=0; j<natoms_per_task; ++j) {
136 387187 : all_atoms.push_back( t[j] );
137 : }
138 192438 : t.resize(0);
139 : }
140 : }
141 226 : if( all_atoms.size()==0 ) {
142 0 : error("No atoms have been specified");
143 : }
144 226 : requestAtoms(all_atoms);
145 226 : if( keywords.exists("NOPBC") ) {
146 226 : bool nopbc=!usepbc;
147 226 : parseFlag("NOPBC",nopbc);
148 226 : usepbc=!nopbc;
149 : }
150 226 : if( keywords.exists("WHOLEMOLECULES") ) {
151 42 : parseFlag("WHOLEMOLECULES",wholemolecules);
152 42 : if( wholemolecules ) {
153 1 : usepbc=false;
154 : }
155 : }
156 226 : if( usepbc ) {
157 191 : log.printf(" using periodic boundary conditions\n");
158 : } else {
159 35 : log.printf(" without periodic boundary conditions\n");
160 : }
161 :
162 : // Setup the values
163 226 : mode = T::getModeAndSetupValues( this );
164 : // This sets up an array in the parallel task manager to hold all the indices
165 : // Sets up the index list in the task manager
166 226 : taskmanager.setupParallelTaskManager( 3*natoms_per_task + virialSize, virialSize );
167 226 : taskmanager.setActionInput( MultiColvarInput{ usepbc, mode, natoms_per_task });
168 226 : }
169 :
170 : template <class T>
171 1946 : unsigned MultiColvarTemplate<T>::getNumberOfDerivatives() {
172 1946 : return 3*getNumberOfAtoms()+9;
173 : }
174 :
175 : template <class T>
176 25806 : void MultiColvarTemplate<T>::calculate() {
177 25806 : if( wholemolecules ) {
178 1 : makeWhole();
179 : }
180 25806 : taskmanager.runAllTasks();
181 25806 : }
182 :
183 : template <class T>
184 21985 : void MultiColvarTemplate<T>::applyNonZeroRankForces( std::vector<double>& outforces ) {
185 21985 : taskmanager.applyForces( outforces );
186 21985 : }
187 :
188 : template <class T>
189 119 : void MultiColvarTemplate<T>::addValueWithDerivatives( const std::vector<std::size_t>& shape ) {
190 119 : std::vector<std::size_t> s(1);
191 119 : s[0]=getNumberOfAtoms() / natoms_per_task;
192 119 : addValue( s );
193 119 : }
194 :
195 : template <class T>
196 334 : void MultiColvarTemplate<T>::addComponentWithDerivatives( const std::string& compName, const std::vector<std::size_t>& shape ) {
197 334 : std::vector<std::size_t> s(1);
198 334 : s[0]=getNumberOfAtoms() / natoms_per_task;
199 334 : addComponent( compName, s );
200 334 : }
201 :
202 : template <class T>
203 25806 : void MultiColvarTemplate<T>::getInputData( std::vector<double>& inputdata ) const {
204 25806 : std::size_t ntasks = getConstPntrToComponent(0)->getNumberOfStoredValues();
205 25806 : if( inputdata.size()!=5*natoms_per_task*ntasks ) {
206 219 : inputdata.resize( 5*natoms_per_task*ntasks );
207 : }
208 :
209 : std::size_t k=0;
210 1439574 : for(unsigned i=0; i<ntasks; ++i) {
211 4136461 : for(unsigned j=0; j<natoms_per_task; ++j) {
212 2722693 : Vector mypos( getPosition( natoms_per_task*i + j ) );
213 2722693 : inputdata[k] = mypos[0];
214 2722693 : k++;
215 2722693 : inputdata[k] = mypos[1];
216 2722693 : k++;
217 2722693 : inputdata[k] = mypos[2];
218 2722693 : k++;
219 : }
220 4136461 : for(unsigned j=0; j<natoms_per_task; ++j) {
221 2722693 : inputdata[k] = getMass( natoms_per_task*i + j );
222 2722693 : k++;
223 : }
224 4136461 : for(unsigned j=0; j<natoms_per_task; ++j) {
225 2722693 : inputdata[k] = getCharge( natoms_per_task*i + j );
226 2722693 : k++;
227 : }
228 : }
229 25806 : }
230 :
231 : template <class T>
232 2230892 : void MultiColvarTemplate<T>::performTask( std::size_t task_index,
233 : const MultiColvarInput& actiondata,
234 : ParallelActionsInput& input,
235 : ParallelActionsOutput& output ) {
236 2230892 : std::size_t pos_start = 5*actiondata.nindices_per_task*task_index;
237 2230892 : if( actiondata.usepbc ) {
238 1823243 : if( actiondata.nindices_per_task==1 ) {
239 : //this may be changed to input.pbc.apply() en mass or only on this one
240 49576 : Vector fpos=input.pbc->distance(Vector(0.0,0.0,0.0),
241 49576 : Vector(input.inputdata[pos_start],
242 24788 : input.inputdata[pos_start+1],
243 24788 : input.inputdata[pos_start+2]) );
244 24788 : input.inputdata[pos_start] =fpos[0];
245 24788 : input.inputdata[pos_start+1]=fpos[1];
246 24788 : input.inputdata[pos_start+2]=fpos[2];
247 : } else {
248 : //make whole?
249 : std::size_t apos_start = pos_start;
250 : //if accidentaly nindices_per_task is 0, this will work by looping on all possible unsigned integers!!!!
251 3608846 : for(unsigned j=0; j<actiondata.nindices_per_task-1; ++j) {
252 1810391 : Vector first(input.inputdata[apos_start],
253 1810391 : input.inputdata[apos_start+1],
254 1810391 : input.inputdata[apos_start+2]);
255 1810391 : Vector second(input.inputdata[apos_start+3],
256 1810391 : input.inputdata[apos_start+4],
257 1810391 : input.inputdata[apos_start+5]);
258 : //calling the pbc here gives problems
259 1810391 : second=first+input.pbc->distance(first,second);
260 1810391 : input.inputdata[apos_start+3]=second[0];
261 1810391 : input.inputdata[apos_start+4]=second[1];
262 1810391 : input.inputdata[apos_start+5]=second[2];
263 : apos_start += 3;
264 : }
265 : }
266 407649 : } else if( actiondata.nindices_per_task==1 ) {
267 : //isn't this equivalent to x = x-0?
268 : //why this is needed?
269 : Vector fpos=delta(Vector(0.0,0.0,0.0),
270 134291 : Vector(input.inputdata[pos_start],
271 134291 : input.inputdata[pos_start+1],
272 134291 : input.inputdata[pos_start+2]));
273 : input.inputdata[pos_start]=fpos[0];
274 134291 : input.inputdata[pos_start+1]=fpos[1];
275 134291 : input.inputdata[pos_start+2]=fpos[2];
276 : }
277 2230892 : const size_t mass_start = pos_start + 3*actiondata.nindices_per_task;
278 2230892 : const size_t charge_start = mass_start + actiondata.nindices_per_task;
279 2230892 : const size_t local_ndev = 3*actiondata.nindices_per_task+virialSize;
280 :
281 2230892 : ColvarOutput cvout { output.values,
282 : local_ndev,
283 : output.derivatives.data() };
284 4461784 : T::calculateCV( ColvarInput{actiondata.mode,
285 : actiondata.nindices_per_task,
286 2230892 : input.inputdata+pos_start,
287 2230892 : input.inputdata+mass_start,
288 2230892 : input.inputdata+charge_start,
289 2230892 : *input.pbc},
290 : cvout );
291 2230892 : }
292 :
293 : template <class T>
294 1030449 : int MultiColvarTemplate<T>::getNumberOfValuesPerTask( std::size_t task_index,
295 : const MultiColvarInput& actiondata ) {
296 1030449 : return 1;
297 : }
298 :
299 : template <class T>
300 1030449 : void MultiColvarTemplate<T>::getForceIndices( std::size_t task_index,
301 : std::size_t colno,
302 : std::size_t ntotal_force,
303 : const MultiColvarInput& actiondata,
304 : const ParallelActionsInput& input,
305 : ForceIndexHolder force_indices ) {
306 2228186 : for(unsigned i=0; i<input.ncomponents; ++i) {
307 : std::size_t m=0;
308 1197737 : std::size_t base = 3*task_index*actiondata.nindices_per_task;
309 3442683 : for(unsigned j=0; j<actiondata.nindices_per_task; ++j) {
310 2244946 : force_indices.indices[i][m] = base + m;
311 2244946 : ++m;
312 2244946 : force_indices.indices[i][m] = base + m;
313 2244946 : ++m;
314 2244946 : force_indices.indices[i][m] = base + m;
315 2244946 : ++m;
316 : }
317 1197737 : force_indices.threadsafe_derivatives_end[i] = 3*actiondata.nindices_per_task;
318 1197737 : force_indices.indices[i][m+0] = ntotal_force - 9;
319 1197737 : force_indices.indices[i][m+1] = ntotal_force - 8;
320 1197737 : force_indices.indices[i][m+2] = ntotal_force - 7;
321 1197737 : force_indices.indices[i][m+3] = ntotal_force - 6;
322 1197737 : force_indices.indices[i][m+4] = ntotal_force - 5;
323 1197737 : force_indices.indices[i][m+5] = ntotal_force - 4;
324 1197737 : force_indices.indices[i][m+6] = ntotal_force - 3;
325 1197737 : force_indices.indices[i][m+7] = ntotal_force - 2;
326 1197737 : force_indices.indices[i][m+8] = ntotal_force - 1;
327 1197737 : force_indices.tot_indices[i] = 3*actiondata.nindices_per_task + virialSize;
328 : }
329 1030449 : }
330 :
331 : } // namespace colvar
332 : } // namespace PLMD
333 : #endif
|
