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 : namespace PLMD {
27 :
28 : template<typename T>
29 1075805 : static void getPointer(const TypesafePtr & p, const std::vector<unsigned>& shape, const unsigned& start, const unsigned& stride, T*&pp ) {
30 1075805 : if( shape.size()==1 && stride==1 ) {
31 71969 : pp=p.get<T*>( {shape[0]} );
32 1003836 : } else if( shape.size()==1 ) {
33 762201 : auto p_=p.get<T*>( {shape[0],stride} );
34 762183 : pp = p_+start;
35 241635 : } else if( shape.size()==2 ) {
36 241635 : pp=p.get<T*>( {shape[0],shape[1]} );
37 : }
38 1075787 : }
39 :
40 : template <class T>
41 : class DataPassingObjectTyped : public DataPassingObject {
42 : private:
43 : /// A pointer to the value
44 : TypesafePtr v;
45 : /// A pointer to the force
46 : TypesafePtr f;
47 : public:
48 : /// This convers a number from the MD code into a double
49 : double MD2double(const TypesafePtr &) const override ;
50 : /// This is used when you want to save the passed object to a double variable in PLUMED rather than the pointer
51 : /// this can be used even when you don't pass a pointer from the MD code
52 : void saveValueAsDouble( const TypesafePtr & val ) override;
53 : /// Set the pointer to the value
54 : void setValuePointer( const TypesafePtr & p, const std::vector<unsigned>& shape, const bool& isconst ) override;
55 : /// Set the pointer to the force
56 : void setForcePointer( const TypesafePtr & p, const std::vector<unsigned>& shape ) override;
57 : /// This gets the data in the pointer and passes it to the output value
58 : void share_data( std::vector<double>& values ) const override ;
59 : /// Share the data and put it in the value from sequential data
60 : void share_data( const unsigned& j, const unsigned& k, Value* value ) override;
61 : /// Share the data and put it in the value from a scattered data
62 : void share_data( const std::vector<AtomNumber>&index, const std::vector<unsigned>& i, Value* value ) override;
63 : /// Pass the force from the value to the output value
64 : void add_force( Value* vv ) override;
65 : void add_force( const std::vector<int>& index, Value* value ) override;
66 : void add_force( const std::vector<AtomNumber>& index, const std::vector<unsigned>& i, Value* value ) override;
67 : /// Rescale the force on the output value
68 : void rescale_force( const unsigned& n, const double& factor, Value* value ) override;
69 : /// This transfers everything to the output
70 : void setData( Value* value ) override;
71 : };
72 :
73 10990 : std::unique_ptr<DataPassingObject> DataPassingObject::create(unsigned n) {
74 10990 : if(n==sizeof(double)) {
75 10973 : return std::make_unique<DataPassingObjectTyped<double>>();
76 17 : } else if(n==sizeof(float)) {
77 17 : return std::make_unique<DataPassingObjectTyped<float>>();
78 : }
79 : std::string pp;
80 0 : Tools::convert(n,pp);
81 0 : plumed_merror("cannot create an MD interface with sizeof(real)=="+ pp);
82 : return NULL;
83 : }
84 :
85 : template <class T>
86 0 : double DataPassingObjectTyped<T>::MD2double(const TypesafePtr & m) const {
87 0 : double d=double(m.template get<T>());
88 0 : return d;
89 : }
90 :
91 : template <class T>
92 5576 : void DataPassingObjectTyped<T>::saveValueAsDouble( const TypesafePtr & val ) {
93 5576 : hasbackup=true;
94 5576 : bvalue=double(val.template get<T>());
95 5576 : }
96 :
97 : template <>
98 5 : void DataPassingObjectTyped<float>::saveValueAsDouble( const TypesafePtr & val ) {
99 5 : hasbackup=true;
100 5 : bvalue=double(val.template get<float>());
101 : // The following is to avoid extra digits in case the MD code uses floats
102 : // e.g.: float f=0.002 when converted to double becomes 0.002000000094995
103 : // To avoid this, we keep only up to 6 significant digits after first one.
104 : // This is done only for quantities that the user writes in base ten, namely the timestep,
105 : // because the resulting grid is coarser than the float itself and would discard
106 : // information from quantities that the MD code computes, such as the energy.
107 : // std::to_chars would express this better, as it writes the shortest decimal that reads
108 : // back as the same float, but floating point <charconv> requires GCC 11 or libc++ 14
109 5 : if( roundToBaseTen && bvalue != 0.0 ) {
110 3 : double magnitude=std::pow(10,std::floor(std::log10(std::fabs(bvalue))));
111 3 : bvalue = std::round(bvalue/magnitude*1e6)/1e6*magnitude;
112 : }
113 5 : }
114 :
115 : template <class T>
116 814524 : void DataPassingObjectTyped<T>::setValuePointer( const TypesafePtr & val, const std::vector<unsigned>& shape, const bool& isconst ) {
117 814524 : if( shape.size()==0 ) {
118 631382 : if( isconst ) {
119 631254 : val.get<const T*>();
120 : } else {
121 128 : val.get<T*>(); // just check type and discard pointer
122 : }
123 183142 : } else if( shape.size()==1 ) {
124 46410 : if( isconst )
125 46390 : val.get<const T*>({shape[0]});
126 : else
127 20 : val.get<T*>({shape[0]}); // just check type and discard pointer
128 136732 : } else if( shape.size()==2 ) {
129 136732 : if( isconst )
130 136729 : val.get<const T*>({shape[0],shape[1]});
131 : else
132 3 : val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
133 : }
134 814523 : v=val.copy();
135 814523 : }
136 :
137 : template <class T>
138 555645 : void DataPassingObjectTyped<T>::setForcePointer( const TypesafePtr & val, const std::vector<unsigned>& shape ) {
139 555645 : if( shape.size()==0 ) {
140 422652 : val.get<T*>(); // just check type and discard pointer
141 132993 : } else if( shape.size()==1 )
142 0 : val.get<T*>({shape[0]}); // just check type and discard pointer
143 132993 : else if( shape.size()==2 )
144 132993 : val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
145 555629 : f=val.copy();
146 555627 : }
147 :
148 : template <class T>
149 23914 : void DataPassingObjectTyped<T>::setData( Value* value ) {
150 23914 : if( value->getRank()==0 ) {
151 703 : *v.template get<T*>() = static_cast<T>(value->get()) / unit;
152 703 : return;
153 : }
154 : T* pp;
155 23211 : getPointer( v, value->getShape(), start, stride, pp );
156 23211 : unsigned nvals=value->getNumberOfValues();
157 888006 : for(unsigned i=0; i<nvals; ++i) {
158 864795 : pp[i] = T( value->get(i) );
159 : }
160 : }
161 :
162 : template <class T>
163 536288 : void DataPassingObjectTyped<T>::share_data( const unsigned& j, const unsigned& k, Value* value ) {
164 536288 : if( value->getRank()==0 ) {
165 8852 : if( hasbackup ) {
166 8780 : value->set( unit*bvalue );
167 : } else {
168 72 : value->set( unit*double(v.template get<T>()) );
169 : }
170 8852 : return;
171 : }
172 527436 : std::vector<unsigned> s(value->getShape());
173 527436 : if( s.size()==1 ) {
174 392741 : s[0]=k-j;
175 : }
176 : const T* pp;
177 527436 : getPointer( v, s, start, stride, pp );
178 527420 : std::vector<double> & d=value->data;
179 527420 : #pragma omp parallel for num_threads(value->getGoodNumThreads(j,k))
180 : for(unsigned i=j; i<k; ++i) {
181 : d[i]=unit*pp[i*stride];
182 : }
183 : }
184 :
185 : template <class T>
186 1350 : void DataPassingObjectTyped<T>::share_data( std::vector<double>& values ) const {
187 1350 : std::vector<unsigned> maxel(1,values.size());
188 : const T* pp;
189 1350 : getPointer( v, maxel, start, stride, pp );
190 1350 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(values))
191 : for(unsigned i=0; i<values.size(); ++i) {
192 : values[i]=unit*pp[i*stride];
193 : }
194 1350 : }
195 :
196 : template <class T>
197 71881 : void DataPassingObjectTyped<T>::share_data( const std::vector<AtomNumber>&index, const std::vector<unsigned>& i, Value* value ) {
198 : plumed_dbg_assert( value->getRank()==1 );
199 71881 : std::vector<unsigned> maxel(1,index.size());
200 : #ifndef NDEBUG
201 : // bounds are only checked in debug mode since they require this extra step that is potentially expensive
202 : maxel[0]=(i.size()>0?*std::max_element(i.begin(),i.end())+1:0);
203 : #else
204 71881 : maxel[0]=0;
205 : #endif
206 : const T* pp;
207 71881 : getPointer( v, maxel, start, stride, pp );
208 : // cannot be parallelized with omp because access to data is not ordered
209 : unsigned k=0;
210 2513334 : for(const auto & p : index) {
211 2441453 : value->data[p.index()]=unit*pp[i[k]*stride];
212 2441453 : k++;
213 : }
214 71881 : }
215 :
216 : template <class T>
217 106968 : void DataPassingObjectTyped<T>::add_force( Value* value ) {
218 106968 : if( value->getRank()==0 ) {
219 40 : *f.template get<T*>() += funit*static_cast<T>(value->getForce(0));
220 40 : return;
221 : }
222 : T* pp;
223 106928 : getPointer( f, value->getShape(), start, stride, pp );
224 106928 : unsigned nvals=value->getNumberOfValues();
225 106928 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,nvals))
226 : for(unsigned i=0; i<nvals; ++i) {
227 : pp[i*stride] += funit*T(value->getForce(i));
228 : }
229 : }
230 :
231 : template <class T>
232 281234 : void DataPassingObjectTyped<T>::add_force( const std::vector<int>& index, Value* value ) {
233 281234 : plumed_assert( value->getRank()==1 );
234 281234 : std::vector<unsigned> s(1,index.size());
235 : T* pp;
236 281234 : getPointer( f, s, start, stride, pp );
237 281232 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,index.size()))
238 : for(unsigned i=0; i<index.size(); ++i) {
239 : pp[i*stride] += funit*T(value->getForce(index[i]));
240 : }
241 281232 : }
242 :
243 : template <class T>
244 63615 : void DataPassingObjectTyped<T>::add_force( const std::vector<AtomNumber>& index, const std::vector<unsigned>& i, Value* value ) {
245 : plumed_dbg_assert( value->getRank()==1 );
246 63615 : std::vector<unsigned> maxel(1,index.size());
247 : #ifndef NDEBUG
248 : // bounds are only checked in debug mode since they require this extra step that is potentially expensive
249 : maxel[0]=(i.size()>0?*std::max_element(i.begin(),i.end())+1:0);
250 : #else
251 63615 : maxel[0]=0;
252 : #endif
253 : T* pp;
254 63615 : getPointer( f, maxel, start, stride, pp );
255 : unsigned k=0;
256 1861299 : for(const auto & p : index) {
257 1797684 : pp[stride*i[k]] += funit*T(value->getForce(p.index()));
258 1797684 : k++;
259 : }
260 63615 : }
261 :
262 : template <class T>
263 150 : void DataPassingObjectTyped<T>::rescale_force( const unsigned& n, const double& factor, Value* value ) {
264 150 : plumed_assert( value->getRank()>0 );
265 150 : std::vector<unsigned> s( value->getShape() );
266 150 : if( s.size()==1 ) {
267 150 : s[0] = n;
268 : }
269 : T* pp;
270 150 : getPointer( f, s, start, stride, pp );
271 : // The MD code might not pass this buffer on every rank (e.g. the virial)
272 150 : if( !pp ) {
273 : return;
274 : }
275 150 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,n))
276 : for(unsigned i=0; i<n; ++i) {
277 : pp[i*stride] *= T(factor);
278 : }
279 : }
280 :
281 : }
|
