Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "DataPassingObject.h"
23 : #include "tools/OpenMP.h"
24 : #include "tools/Tools.h"
25 :
26 : #include <algorithm>
27 :
28 : namespace PLMD {
29 :
30 : #ifdef _OPENMP
31 : /// Number of threads used to copy n elements between the arrays of the MD code and PLUMED.
32 : /// These loops do very little work per element, so for the system sizes that are typical in MD
33 : /// starting the threads (and waking them up, as PLUMED is called once per step) costs more than
34 : /// the copy itself: with 4 threads, copying the positions of 7000 atoms took 5 to 10 times longer
35 : /// than with one thread. Threads are thus only used when each one has at least this many elements.
36 : /// It is called outside the pragmas (and only with OpenMP), so that the call is visible also to static analysis.
37 1046095 : static unsigned getCopyNumThreads( const std::size_t& n ) {
38 : constexpr std::size_t minElementsPerThread=32768;
39 1046095 : const std::size_t m=n/minElementsPerThread;
40 1046095 : if( m<2 ) {
41 : return 1;
42 : }
43 10 : return std::min<std::size_t>( m, OpenMP::getNumThreads() );
44 : }
45 : #endif
46 :
47 : template<typename T>
48 1283721 : static void getPointer(const TypesafePtr & p, const std::vector<std::size_t>& shape, const unsigned& start, const unsigned& stride, T*&pp ) {
49 1283721 : if( shape.size()==1 && stride==1 ) {
50 82348 : pp=p.get<T*>( {shape[0]} );
51 1201373 : } else if( shape.size()==1 ) {
52 912108 : auto p_=p.get<T*>( {shape[0],stride} );
53 912087 : pp = p_+start;
54 289265 : } else if( shape.size()==2 ) {
55 : // start and stride are only used for vectors that are interleaved in the same array (e.g. positions)
56 289265 : plumed_massert( start==0 && stride==1, "start and stride cannot be used when passing a matrix");
57 289265 : pp=p.get<T*>( {shape[0],shape[1]} );
58 : } else {
59 0 : plumed_merror("data with rank larger than two cannot be passed between PLUMED and the MD code");
60 : }
61 1283700 : }
62 :
63 : template <class T>
64 : class DataPassingObjectTyped : public DataPassingObject {
65 : private:
66 : /// A pointer to the value
67 : TypesafePtr v;
68 : /// A pointer to the force
69 : TypesafePtr f;
70 : public:
71 : /// This convers a number from the MD code into a double
72 : double MD2double(const TypesafePtr &) const override ;
73 : /// This is used when you want to save the passed object to a double variable in PLUMED rather than the pointer
74 : /// this can be used even when you don't pass a pointer from the MD code
75 : void saveValueAsDouble( const TypesafePtr & val ) override;
76 : /// Set the pointer to the value
77 : void setValuePointer( const TypesafePtr & p, const std::vector<std::size_t>& shape, const bool& isconst ) override;
78 : /// Set the pointer to the force
79 : void setForcePointer( const TypesafePtr & p, const std::vector<std::size_t>& shape ) override;
80 : /// This gets the data in the pointer and passes it to the output value
81 : void share_data( std::vector<double>& values ) const override ;
82 : /// Share the data and put it in the value from sequential data
83 : void share_data( const unsigned& j, const unsigned& k, Value* value ) override;
84 : /// Share the data and put it in the value from a scattered data
85 : void share_data( const std::vector<AtomNumber>&index, const std::vector<unsigned>& i, const std::size_t& nlocal, Value* value ) override;
86 : /// Check if the MD code passed an array to store the forces on this quantity
87 546415 : bool hasForcePointer() const override {
88 546415 : return static_cast<bool>(f);
89 : }
90 : /// Pass the force from the value to the output value
91 : void add_force( Value* vv ) override;
92 : void add_force( const std::vector<int>& index, Value* value ) override;
93 : void add_force( const std::vector<AtomNumber>& index, const std::vector<unsigned>& i, const std::size_t& nlocal, Value* value ) override;
94 : /// Rescale the force on the output value
95 : void rescale_force( const unsigned& n, const double& factor, Value* value ) override;
96 : /// This transfers everything to the output
97 : void setData( Value* value ) override;
98 : };
99 :
100 15586 : std::unique_ptr<DataPassingObject> DataPassingObject::create(unsigned n) {
101 15586 : if(n==sizeof(double)) {
102 15562 : return std::make_unique<DataPassingObjectTyped<double>>();
103 24 : } else if(n==sizeof(float)) {
104 24 : return std::make_unique<DataPassingObjectTyped<float>>();
105 : }
106 : std::string pp;
107 0 : Tools::convert(n,pp);
108 0 : plumed_merror("cannot create an MD interface with sizeof(real)=="+ pp);
109 : return NULL;
110 : }
111 :
112 : template <class T>
113 0 : double DataPassingObjectTyped<T>::MD2double(const TypesafePtr & m) const {
114 0 : double d=double(m.template get<T>());
115 0 : return d;
116 : }
117 :
118 : template <class T>
119 6406 : void DataPassingObjectTyped<T>::saveValueAsDouble( const TypesafePtr & val ) {
120 6406 : hasbackup=true;
121 6406 : bvalue=double(val.template get<T>());
122 6406 : }
123 :
124 : template <>
125 6 : void DataPassingObjectTyped<float>::saveValueAsDouble( const TypesafePtr & val ) {
126 6 : hasbackup=true;
127 6 : bvalue=double(val.template get<float>());
128 : // The following is to avoid extra digits in case the MD code uses floats
129 : // e.g.: float f=0.002 when converted to double becomes 0.002000000094995
130 : // To avoid this, we keep only up to 6 significant digits after first one.
131 : // This is done only for quantities that the user writes in base ten, namely the timestep,
132 : // because the resulting grid is coarser than the float itself and would discard
133 : // information from quantities that the MD code computes, such as the energy.
134 : // std::to_chars would express this better, as it writes the shortest decimal that reads
135 : // back as the same float, but floating point <charconv> requires GCC 11 or libc++ 14
136 6 : if( roundToBaseTen && bvalue != 0.0 ) {
137 4 : double magnitude=std::pow(10,std::floor(std::log10(std::fabs(bvalue))));
138 4 : bvalue = std::round(bvalue/magnitude*1e6)/1e6*magnitude;
139 : }
140 6 : }
141 :
142 : template <class T>
143 969796 : void DataPassingObjectTyped<T>::setValuePointer( const TypesafePtr & val, const std::vector<std::size_t>& shape, const bool& isconst ) {
144 969796 : if( shape.size()==0 ) {
145 754840 : if( isconst ) {
146 754709 : val.get<const T*>();
147 : } else {
148 131 : val.get<T*>(); // just check type and discard pointer
149 : }
150 214956 : } else if( shape.size()==1 ) {
151 52965 : if( isconst )
152 52942 : val.get<const T*>({shape[0]});
153 : else
154 23 : val.get<T*>({shape[0]}); // just check type and discard pointer
155 161991 : } else if( shape.size()==2 ) {
156 161991 : if( isconst )
157 161988 : val.get<const T*>({shape[0],shape[1]});
158 : else
159 3 : val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
160 : } else {
161 0 : plumed_merror("data with rank larger than two cannot be passed between PLUMED and the MD code");
162 : }
163 969795 : v=val.copy();
164 969795 : }
165 :
166 : template <class T>
167 657478 : void DataPassingObjectTyped<T>::setForcePointer( const TypesafePtr & val, const std::vector<std::size_t>& shape ) {
168 657478 : if( shape.size()==0 ) {
169 499326 : val.get<T*>(); // just check type and discard pointer
170 158152 : } else if( shape.size()==1 )
171 0 : val.get<T*>({shape[0]}); // just check type and discard pointer
172 158152 : else if( shape.size()==2 )
173 158152 : val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
174 : else {
175 0 : plumed_merror("data with rank larger than two cannot be passed between PLUMED and the MD code");
176 : }
177 657462 : f=val.copy();
178 657460 : }
179 :
180 : template <class T>
181 27190 : void DataPassingObjectTyped<T>::setData( Value* value ) {
182 27190 : if( value->getRank()==0 ) {
183 703 : *v.template get<T*>() = static_cast<T>(value->get()) / unit;
184 703 : return;
185 : }
186 26487 : T* pp=nullptr;
187 26487 : getPointer( v, value->getShape(), start, stride, pp );
188 26487 : unsigned nvals=value->getNumberOfValues();
189 1019046 : for(unsigned i=0; i<nvals; ++i) {
190 992559 : pp[i] = T( value->get(i) );
191 : }
192 : }
193 :
194 : template <class T>
195 609437 : void DataPassingObjectTyped<T>::share_data( const unsigned& j, const unsigned& k, Value* value ) {
196 609437 : if( value->getRank()==0 ) {
197 11215 : if( hasbackup ) {
198 10925 : value->set( unit*bvalue );
199 : } else {
200 290 : value->set( unit*double(v.template get<T>()) );
201 : }
202 11215 : return;
203 : }
204 598222 : std::vector<std::size_t> s(value->getShape());
205 598222 : if( s.size()==1 ) {
206 : // elements j to k-1 are read, so the array passed by the MD code must contain at least k elements
207 438300 : s[0]=k;
208 : }
209 598222 : const T* pp=nullptr;
210 598222 : getPointer( v, s, start, stride, pp );
211 598206 : std::vector<double> & d=value->data;
212 : #ifdef _OPENMP
213 598206 : const unsigned nt=getCopyNumThreads(k-j);
214 : #endif
215 598206 : #pragma omp parallel for num_threads(nt)
216 : for(unsigned i=j; i<k; ++i) {
217 : d[i]=unit*pp[i*stride];
218 : }
219 598222 : }
220 :
221 : template <class T>
222 1422 : void DataPassingObjectTyped<T>::share_data( std::vector<double>& values ) const {
223 1422 : std::vector<std::size_t> maxel(1,values.size());
224 1422 : const T* pp=nullptr;
225 1422 : getPointer( v, maxel, start, stride, pp );
226 : #ifdef _OPENMP
227 1422 : const unsigned nt=getCopyNumThreads(values.size());
228 : #endif
229 1422 : #pragma omp parallel for num_threads(nt)
230 : for(unsigned i=0; i<values.size(); ++i) {
231 : values[i]=unit*pp[i*stride];
232 : }
233 1422 : }
234 :
235 : template <class T>
236 110015 : void DataPassingObjectTyped<T>::share_data( const std::vector<AtomNumber>&index, const std::vector<unsigned>& i, const std::size_t& nlocal, Value* value ) {
237 : plumed_dbg_assert( value->getRank()==1 );
238 : plumed_dbg_assert( index.size()==i.size() );
239 : plumed_dbg_assert( i.size()==0 || *std::max_element(i.begin(),i.end())<nlocal );
240 : // the size of the array passed by the MD code is checked against the number of local atoms
241 110015 : std::vector<std::size_t> maxel(1,nlocal);
242 110015 : const T* pp=nullptr;
243 110015 : getPointer( v, maxel, start, stride, pp );
244 : // cannot be parallelized with omp because access to data is not ordered
245 : unsigned k=0;
246 3244681 : for(const auto & p : index) {
247 3134669 : value->data[p.index()]=unit*pp[i[k]*stride];
248 3134669 : k++;
249 : }
250 110015 : }
251 :
252 : template <class T>
253 129371 : void DataPassingObjectTyped<T>::add_force( Value* value ) {
254 129371 : if( value->getRank()==0 ) {
255 40 : *f.template get<T*>() += funit*static_cast<T>(value->getForce(0));
256 40 : return;
257 : }
258 129331 : T* pp=nullptr;
259 129331 : getPointer( f, value->getShape(), start, stride, pp );
260 129331 : unsigned nvals=value->getNumberOfValues();
261 : #ifdef _OPENMP
262 129331 : const unsigned nt=getCopyNumThreads(nvals);
263 : #endif
264 129331 : #pragma omp parallel for num_threads(nt)
265 : for(unsigned i=0; i<nvals; ++i) {
266 : pp[i*stride] += funit*T(value->getForce(i));
267 : }
268 : }
269 :
270 : template <class T>
271 316988 : void DataPassingObjectTyped<T>::add_force( const std::vector<int>& index, Value* value ) {
272 316988 : plumed_assert( value->getRank()==1 );
273 316988 : std::vector<std::size_t> s(1,index.size());
274 316988 : T* pp=nullptr;
275 316988 : getPointer( f, s, start, stride, pp );
276 : #ifdef _OPENMP
277 316986 : const unsigned nt=getCopyNumThreads(index.size());
278 : #endif
279 316986 : #pragma omp parallel for num_threads(nt)
280 : for(unsigned i=0; i<index.size(); ++i) {
281 : pp[i*stride] += funit*T(value->getForce(index[i]));
282 : }
283 316988 : }
284 :
285 : template <class T>
286 101106 : void DataPassingObjectTyped<T>::add_force( const std::vector<AtomNumber>& index, const std::vector<unsigned>& i, const std::size_t& nlocal, Value* value ) {
287 : plumed_dbg_assert( value->getRank()==1 );
288 : plumed_dbg_assert( index.size()==i.size() );
289 : plumed_dbg_assert( i.size()==0 || *std::max_element(i.begin(),i.end())<nlocal );
290 : // the size of the array passed by the MD code is checked against the number of local atoms
291 101106 : std::vector<std::size_t> maxel(1,nlocal);
292 101106 : T* pp=nullptr;
293 101106 : getPointer( f, maxel, start, stride, pp );
294 : unsigned k=0;
295 2081682 : for(const auto & p : index) {
296 1980576 : pp[stride*i[k]] += funit*T(value->getForce(p.index()));
297 1980576 : k++;
298 : }
299 101106 : }
300 :
301 : template <class T>
302 150 : void DataPassingObjectTyped<T>::rescale_force( const unsigned& n, const double& factor, Value* value ) {
303 150 : plumed_assert( value->getRank()>0 );
304 150 : std::vector<std::size_t> s( value->getShape() );
305 150 : if( s.size()==1 ) {
306 150 : s[0] = n;
307 : }
308 150 : T* pp=nullptr;
309 150 : getPointer( f, s, start, stride, pp );
310 : // The MD code might not pass this buffer on every rank (e.g. the virial)
311 150 : if( !pp ) {
312 : return;
313 : }
314 : #ifdef _OPENMP
315 150 : const unsigned nt=getCopyNumThreads(n);
316 : #endif
317 150 : #pragma omp parallel for num_threads(nt)
318 : for(unsigned i=0; i<n; ++i) {
319 : pp[i*stride] *= T(factor);
320 : }
321 150 : }
322 :
323 : }
|
