Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 "MDAtoms.h"
23 : #include "tools/Tools.h"
24 : #include "tools/OpenMP.h"
25 : #include "tools/Exception.h"
26 : #include "tools/Units.h"
27 : #include <algorithm>
28 : #include <string>
29 : #include <map>
30 :
31 : using namespace std;
32 :
33 : namespace PLMD {
34 :
35 : /// Class containing the pointers to the MD data
36 : /// It is templated so that single and double precision versions coexist
37 : /// IT IS STILL UNDOCUMENTED. IT PROBABLY NEEDS A STRONG CLEANUP
38 : template <class T>
39 10167 : class MDAtomsTyped:
40 : public MDAtomsBase
41 : {
42 : T scalep,scalef;
43 : T scaleb,scalev;
44 : T scalec,scalem; // factor to scale charges and masses
45 : int stride;
46 : T *m;
47 : T *c;
48 : T *px; T *py; T *pz;
49 : T *fx; T *fy; T *fz;
50 : T *box;
51 : T *virial;
52 : std::map<std::string,T*> extraCV;
53 : std::map<std::string,T*> extraCVForce;
54 : public:
55 : MDAtomsTyped();
56 : void setm(void*m) override;
57 : void setc(void*m) override;
58 : void setBox(void*) override;
59 : void setp(void*p) override;
60 : void setVirial(void*) override;
61 : void setf(void*f) override;
62 : void setp(void*p,int i) override;
63 : void setf(void*f,int i) override;
64 : void setUnits(const Units&,const Units&) override;
65 20 : void setExtraCV(const std::string &name,void*p) override {
66 20 : extraCV[name]=static_cast<T*>(p);
67 20 : }
68 20 : void setExtraCVForce(const std::string &name,void*p) override {
69 20 : extraCVForce[name]=static_cast<T*>(p);
70 20 : }
71 40 : double getExtraCV(const std::string &name) override {
72 40 : return static_cast<double>(*extraCV[name]);
73 : }
74 40 : void updateExtraCVForce(const std::string &name,double f) override {
75 40 : *extraCVForce[name]+=static_cast<T>(f);
76 40 : }
77 12119 : void MD2double(const void*m,double&d)const override {
78 12119 : d=double(*(static_cast<const T*>(m)));
79 12119 : }
80 2303 : void double2MD(const double&d,void*m) const override {
81 2303 : *(static_cast<T*>(m))=T(d);
82 2303 : }
83 0 : Vector getMDforces(const unsigned index) const override {
84 0 : Vector force(fx[stride*index],fy[stride*index],fz[stride*index]);
85 0 : return force/scalef;
86 : }
87 : void getBox(Tensor &) const override;
88 : void getPositions(const vector<int>&index,vector<Vector>&positions) const override;
89 : void getPositions(const std::set<AtomNumber>&index,const vector<unsigned>&i,vector<Vector>&positions) const override;
90 : void getPositions(unsigned j,unsigned k,vector<Vector>&positions) const override;
91 : void getLocalPositions(std::vector<Vector>&p) const override;
92 : void getMasses(const vector<int>&index,vector<double>&) const override;
93 : void getCharges(const vector<int>&index,vector<double>&) const override;
94 : void updateVirial(const Tensor&) const override;
95 : void updateForces(const vector<int>&index,const vector<Vector>&) override;
96 : void updateForces(const std::set<AtomNumber>&index,const vector<unsigned>&i,const vector<Vector>&forces) override;
97 : void rescaleForces(const vector<int>&index,double factor) override;
98 : unsigned getRealPrecision()const override;
99 : };
100 :
101 : template <class T>
102 771 : void MDAtomsTyped<T>::setUnits(const Units& units,const Units& MDUnits) {
103 771 : double lscale=units.getLength()/MDUnits.getLength();
104 771 : double escale=units.getEnergy()/MDUnits.getEnergy();
105 771 : double cscale=units.getCharge()/MDUnits.getCharge();
106 771 : double mscale=units.getMass()/MDUnits.getMass();
107 : // scalep and scaleb are used to convert MD to plumed
108 771 : scalep=1.0/lscale;
109 771 : scaleb=1.0/lscale;
110 : // scalef and scalev are used to convert plumed to MD
111 771 : scalef=escale/lscale;
112 771 : scalev=escale;
113 771 : scalec=1.0/cscale;
114 771 : scalem=1.0/mscale;
115 771 : }
116 :
117 : template <class T>
118 79497 : void MDAtomsTyped<T>::getBox(Tensor&box)const {
119 79497 : if(this->box) for(int i=0; i<3; i++)for(int j=0; j<3; j++) box(i,j)=this->box[3*i+j]*scaleb;
120 3840 : else box.zero();
121 79497 : }
122 :
123 : template <class T>
124 0 : void MDAtomsTyped<T>::getPositions(const vector<int>&index,vector<Vector>&positions)const {
125 : // cannot be parallelized with omp because access to positions is not ordered
126 0 : for(unsigned i=0; i<index.size(); ++i) {
127 0 : positions[index[i]][0]=px[stride*i]*scalep;
128 0 : positions[index[i]][1]=py[stride*i]*scalep;
129 0 : positions[index[i]][2]=pz[stride*i]*scalep;
130 : }
131 0 : }
132 :
133 : template <class T>
134 2177 : void MDAtomsTyped<T>::getPositions(const std::set<AtomNumber>&index,const vector<unsigned>&i, vector<Vector>&positions)const {
135 : // cannot be parallelized with omp because access to positions is not ordered
136 : unsigned k=0;
137 33753 : for(const auto & p : index) {
138 88197 : positions[p.index()][0]=px[stride*i[k]]*scalep;
139 88197 : positions[p.index()][1]=py[stride*i[k]]*scalep;
140 88197 : positions[p.index()][2]=pz[stride*i[k]]*scalep;
141 29399 : k++;
142 : }
143 2177 : }
144 :
145 : template <class T>
146 37595 : void MDAtomsTyped<T>::getPositions(unsigned j,unsigned k,vector<Vector>&positions)const {
147 120382 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(&positions[j],(k-j)))
148 45192 : for(unsigned i=j; i<k; ++i) {
149 6549692 : positions[i][0]=px[stride*i]*scalep;
150 3229392 : positions[i][1]=py[stride*i]*scalep;
151 3157337 : positions[i][2]=pz[stride*i]*scalep;
152 : }
153 37595 : }
154 :
155 :
156 : template <class T>
157 450 : void MDAtomsTyped<T>::getLocalPositions(vector<Vector>&positions)const {
158 952 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(positions))
159 502 : for(unsigned i=0; i<positions.size(); ++i) {
160 31806 : positions[i][0]=px[stride*i]*scalep;
161 15862 : positions[i][1]=py[stride*i]*scalep;
162 15650 : positions[i][2]=pz[stride*i]*scalep;
163 : }
164 450 : }
165 :
166 :
167 : template <class T>
168 675 : void MDAtomsTyped<T>::getMasses(const vector<int>&index,vector<double>&masses)const {
169 1201749 : if(m) for(unsigned i=0; i<index.size(); ++i) masses[index[i]]=scalem*m[i];
170 0 : else for(unsigned i=0; i<index.size(); ++i) masses[index[i]]=0.0;
171 675 : }
172 :
173 : template <class T>
174 675 : void MDAtomsTyped<T>::getCharges(const vector<int>&index,vector<double>&charges)const {
175 1197477 : if(c) for(unsigned i=0; i<index.size(); ++i) charges[index[i]]=scalec*c[i];
176 4272 : else for(unsigned i=0; i<index.size(); ++i) charges[index[i]]=0.0;
177 675 : }
178 :
179 : template <class T>
180 26967 : void MDAtomsTyped<T>::updateVirial(const Tensor&virial)const {
181 26967 : if(this->virial) for(int i=0; i<3; i++)for(int j=0; j<3; j++) this->virial[3*i+j]+=T(virial(i,j)*scalev);
182 26967 : }
183 :
184 : template <class T>
185 2063 : void MDAtomsTyped<T>::updateForces(const std::set<AtomNumber>&index,const vector<unsigned>&i,const vector<Vector>&forces) {
186 : unsigned k=0;
187 27999 : for(const auto & p : index) {
188 95492 : fx[stride*i[k]]+=scalef*T(forces[p.index()][0]);
189 95492 : fy[stride*i[k]]+=scalef*T(forces[p.index()][1]);
190 95492 : fz[stride*i[k]]+=scalef*T(forces[p.index()][2]);
191 23873 : k++;
192 : }
193 2063 : }
194 :
195 : template <class T>
196 38827 : void MDAtomsTyped<T>::updateForces(const vector<int>&index,const vector<Vector>&forces) {
197 123826 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(fx,stride*index.size()))
198 46172 : for(unsigned i=0; i<index.size(); ++i) {
199 10685127 : fx[stride*i]+=scalef*T(forces[index[i]][0]);
200 10731306 : fy[stride*i]+=scalef*T(forces[index[i]][1]);
201 10716738 : fz[stride*i]+=scalef*T(forces[index[i]][2]);
202 : }
203 38827 : }
204 :
205 : template <class T>
206 50 : void MDAtomsTyped<T>::rescaleForces(const vector<int>&index,double factor) {
207 50 : if(virial) for(unsigned i=0; i<3; i++)for(unsigned j=0; j<3; j++) virial[3*i+j]*=T(factor);
208 199 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(fx,stride*index.size()))
209 99 : for(unsigned i=0; i<index.size(); ++i) {
210 3942 : fx[stride*i]*=T(factor);
211 3942 : fy[stride*i]*=T(factor);
212 3942 : fz[stride*i]*=T(factor);
213 : }
214 50 : }
215 :
216 : template <class T>
217 771 : unsigned MDAtomsTyped<T>::getRealPrecision()const {
218 771 : return sizeof(T);
219 : }
220 :
221 : template <class T>
222 42333 : void MDAtomsTyped<T>::setp(void*pp) {
223 : T*p=static_cast<T*>(pp);
224 42333 : plumed_assert(stride==0 || stride==3);
225 42333 : px=p;
226 42333 : py=p+1;
227 42333 : pz=p+2;
228 42333 : stride=3;
229 42333 : }
230 :
231 : template <class T>
232 38493 : void MDAtomsTyped<T>::setBox(void*pp) {
233 38493 : box=static_cast<T*>(pp);
234 38493 : }
235 :
236 :
237 : template <class T>
238 42333 : void MDAtomsTyped<T>::setf(void*ff) {
239 : T*f=static_cast<T*>(ff);
240 42333 : plumed_assert(stride==0 || stride==3);
241 42333 : fx=f;
242 42333 : fy=f+1;
243 42333 : fz=f+2;
244 42333 : stride=3;
245 42333 : }
246 :
247 : template <class T>
248 0 : void MDAtomsTyped<T>::setp(void*pp,int i) {
249 : T*p=static_cast<T*>(pp);
250 0 : plumed_assert(stride==0 || stride==1);
251 0 : if(i==0)px=p;
252 0 : if(i==1)py=p;
253 0 : if(i==2)pz=p;
254 0 : stride=1;
255 0 : }
256 :
257 : template <class T>
258 32835 : void MDAtomsTyped<T>::setVirial(void*pp) {
259 32835 : virial=static_cast<T*>(pp);
260 32835 : }
261 :
262 :
263 : template <class T>
264 0 : void MDAtomsTyped<T>::setf(void*ff,int i) {
265 : T*f=static_cast<T*>(ff);
266 0 : plumed_assert(stride==0 || stride==1);
267 0 : if(i==0)fx=f;
268 0 : if(i==1)fy=f;
269 0 : if(i==2)fz=f;
270 0 : stride=1;
271 0 : }
272 :
273 : template <class T>
274 42333 : void MDAtomsTyped<T>::setm(void*m) {
275 42333 : this->m=static_cast<T*>(m);
276 42333 : }
277 :
278 : template <class T>
279 32740 : void MDAtomsTyped<T>::setc(void*c) {
280 32740 : this->c=static_cast<T*>(c);
281 32740 : }
282 :
283 : template <class T>
284 3389 : MDAtomsTyped<T>::MDAtomsTyped():
285 : scalep(1.0),
286 : scalef(1.0),
287 : scaleb(1.0),
288 : scalev(1.0),
289 : scalec(1.0),
290 : scalem(1.0),
291 : stride(0),
292 : m(NULL),
293 : c(NULL),
294 : px(NULL),
295 : py(NULL),
296 : pz(NULL),
297 : fx(NULL),
298 : fy(NULL),
299 : fz(NULL),
300 : box(NULL),
301 3389 : virial(NULL)
302 3389 : {}
303 :
304 3390 : std::unique_ptr<MDAtomsBase> MDAtomsBase::create(unsigned p) {
305 3390 : if(p==sizeof(double)) {
306 3389 : return std::unique_ptr<MDAtomsTyped<double>>(new MDAtomsTyped<double>);
307 1 : } else if (p==sizeof(float)) {
308 0 : return std::unique_ptr<MDAtomsTyped<float>>(new MDAtomsTyped<float>);
309 : }
310 : std::string pp;
311 1 : Tools::convert(p,pp);
312 4 : plumed_merror("cannot create an MD interface with sizeof(real)=="+ pp);
313 : return NULL;
314 : }
315 :
316 : }
317 :
|
