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 <algorithm>
27 : #include <string>
28 :
29 : using namespace std;
30 :
31 : namespace PLMD {
32 :
33 : /// Class containing the pointers to the MD data
34 : /// It is templated so that single and double precision versions coexist
35 : /// IT IS STILL UNDOCUMENTED. IT PROBABLY NEEDS A STRONG CLEANUP
36 : template <class T>
37 2788 : class MDAtomsTyped:
38 : public MDAtomsBase
39 : {
40 : T scalep,scalef;
41 : T scaleb,scalev;
42 : T scalec,scalem; // factor to scale charges and masses
43 : int stride;
44 : T *m;
45 : T *c;
46 : T *px; T *py; T *pz;
47 : T *fx; T *fy; T *fz;
48 : T *box;
49 : T *virial;
50 : public:
51 : MDAtomsTyped();
52 : void setm(void*m);
53 : void setc(void*m);
54 : void setBox(void*);
55 : void setp(void*p);
56 : void setVirial(void*);
57 : void setf(void*f);
58 : void setp(void*p,int i);
59 : void setf(void*f,int i);
60 : void setUnits(const Units&,const Units&);
61 8160 : void MD2double(const void*m,double&d)const {
62 8160 : d=double(*(static_cast<const T*>(m)));
63 8160 : }
64 401 : void double2MD(const double&d,void*m)const {
65 401 : *(static_cast<T*>(m))=T(d);
66 401 : }
67 0 : Vector getMDforces(const unsigned index)const {
68 0 : Vector force(fx[stride*index],fy[stride*index],fz[stride*index]);
69 0 : return force/scalef;
70 : }
71 : void getBox(Tensor &)const;
72 : void getPositions(const vector<int>&index,vector<Vector>&positions)const;
73 : void getPositions(const std::set<AtomNumber>&index,const vector<unsigned>&i,vector<Vector>&positions)const;
74 : void getPositions(unsigned j,unsigned k,vector<Vector>&positions)const;
75 : void getLocalPositions(std::vector<Vector>&p)const;
76 : void getMasses(const vector<int>&index,vector<double>&)const;
77 : void getCharges(const vector<int>&index,vector<double>&)const;
78 : void updateVirial(const Tensor&)const;
79 : void updateForces(const vector<int>&index,const vector<Vector>&);
80 : void updateForces(const std::set<AtomNumber>&index,const vector<unsigned>&i,const vector<Vector>&forces);
81 : void rescaleForces(const vector<int>&index,double factor);
82 : unsigned getRealPrecision()const;
83 : };
84 :
85 : template <class T>
86 648 : void MDAtomsTyped<T>::setUnits(const Units& units,const Units& MDUnits) {
87 648 : double lscale=units.getLength()/MDUnits.getLength();
88 648 : double escale=units.getEnergy()/MDUnits.getEnergy();
89 648 : double cscale=units.getCharge()/MDUnits.getCharge();
90 648 : double mscale=units.getMass()/MDUnits.getMass();
91 : // scalep and scaleb are used to convert MD to plumed
92 648 : scalep=1.0/lscale;
93 648 : scaleb=1.0/lscale;
94 : // scalef and scalev are used to convert plumed to MD
95 648 : scalef=escale/lscale;
96 648 : scalev=escale;
97 648 : scalec=1.0/cscale;
98 648 : scalem=1.0/mscale;
99 648 : }
100 :
101 : template <class T>
102 56100 : void MDAtomsTyped<T>::getBox(Tensor&box)const {
103 56100 : 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;
104 3818 : else box.zero();
105 56100 : }
106 :
107 : template <class T>
108 0 : void MDAtomsTyped<T>::getPositions(const vector<int>&index,vector<Vector>&positions)const {
109 : // cannot be parallelized with omp because access to positions is not ordered
110 0 : for(unsigned i=0; i<index.size(); ++i) {
111 0 : positions[index[i]][0]=px[stride*i]*scalep;
112 0 : positions[index[i]][1]=py[stride*i]*scalep;
113 0 : positions[index[i]][2]=pz[stride*i]*scalep;
114 : }
115 0 : }
116 :
117 : template <class T>
118 2010 : void MDAtomsTyped<T>::getPositions(const std::set<AtomNumber>&index,const vector<unsigned>&i, vector<Vector>&positions)const {
119 : // cannot be parallelized with omp because access to positions is not ordered
120 : unsigned k=0;
121 23607 : for(const auto & p : index) {
122 86388 : positions[p.index()][0]=px[stride*i[k]]*scalep;
123 86388 : positions[p.index()][1]=py[stride*i[k]]*scalep;
124 86388 : positions[p.index()][2]=pz[stride*i[k]]*scalep;
125 21597 : k++;
126 : }
127 2010 : }
128 :
129 : template <class T>
130 26595 : void MDAtomsTyped<T>::getPositions(unsigned j,unsigned k,vector<Vector>&positions)const {
131 82055 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(&positions[j],(k-j)))
132 : for(unsigned i=j; i<k; ++i) {
133 4598790 : positions[i][0]=px[stride*i]*scalep;
134 4598790 : positions[i][1]=py[stride*i]*scalep;
135 4598790 : positions[i][2]=pz[stride*i]*scalep;
136 : }
137 26595 : }
138 :
139 :
140 : template <class T>
141 450 : void MDAtomsTyped<T>::getLocalPositions(vector<Vector>&positions)const {
142 952 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(positions))
143 1004 : for(unsigned i=0; i<positions.size(); ++i) {
144 32400 : positions[i][0]=px[stride*i]*scalep;
145 32400 : positions[i][1]=py[stride*i]*scalep;
146 32400 : positions[i][2]=pz[stride*i]*scalep;
147 : }
148 450 : }
149 :
150 :
151 : template <class T>
152 559 : void MDAtomsTyped<T>::getMasses(const vector<int>&index,vector<double>&masses)const {
153 1043360 : if(m) for(unsigned i=0; i<index.size(); ++i) masses[index[i]]=scalem*m[i];
154 0 : else for(unsigned i=0; i<index.size(); ++i) masses[index[i]]=0.0;
155 559 : }
156 :
157 : template <class T>
158 559 : void MDAtomsTyped<T>::getCharges(const vector<int>&index,vector<double>&charges)const {
159 1040798 : if(c) for(unsigned i=0; i<index.size(); ++i) charges[index[i]]=scalec*c[i];
160 3468 : else for(unsigned i=0; i<index.size(); ++i) charges[index[i]]=0.0;
161 559 : }
162 :
163 : template <class T>
164 19076 : void MDAtomsTyped<T>::updateVirial(const Tensor&virial)const {
165 19076 : 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);
166 19076 : }
167 :
168 : template <class T>
169 1926 : void MDAtomsTyped<T>::updateForces(const std::set<AtomNumber>&index,const vector<unsigned>&i,const vector<Vector>&forces) {
170 : unsigned k=0;
171 19617 : for(const auto & p : index) {
172 70764 : fx[stride*i[k]]+=scalef*T(forces[p.index()][0]);
173 70764 : fy[stride*i[k]]+=scalef*T(forces[p.index()][1]);
174 70764 : fz[stride*i[k]]+=scalef*T(forces[p.index()][2]);
175 17691 : k++;
176 : }
177 1926 : }
178 :
179 : template <class T>
180 27252 : void MDAtomsTyped<T>::updateForces(const vector<int>&index,const vector<Vector>&forces) {
181 84030 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(fx,stride*index.size()))
182 59052 : for(unsigned i=0; i<index.size(); ++i) {
183 6935838 : fx[stride*i]+=scalef*T(forces[index[i]][0]);
184 6935838 : fy[stride*i]+=scalef*T(forces[index[i]][1]);
185 6935838 : fz[stride*i]+=scalef*T(forces[index[i]][2]);
186 : }
187 27252 : }
188 :
189 : template <class T>
190 50 : void MDAtomsTyped<T>::rescaleForces(const vector<int>&index,double factor) {
191 50 : if(virial) for(unsigned i=0; i<3; i++)for(unsigned j=0; j<3; j++) virial[3*i+j]*=T(factor);
192 200 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(fx,stride*index.size()))
193 200 : for(unsigned i=0; i<index.size(); ++i) {
194 5400 : fx[stride*i]*=T(factor);
195 5400 : fy[stride*i]*=T(factor);
196 5400 : fz[stride*i]*=T(factor);
197 : }
198 50 : }
199 :
200 : template <class T>
201 635 : unsigned MDAtomsTyped<T>::getRealPrecision()const {
202 635 : return sizeof(T);
203 : }
204 :
205 : template <class T>
206 30656 : void MDAtomsTyped<T>::setp(void*pp) {
207 : T*p=static_cast<T*>(pp);
208 30656 : plumed_assert(stride==0 || stride==3);
209 30656 : px=p;
210 30656 : py=p+1;
211 30656 : pz=p+2;
212 30656 : stride=3;
213 30656 : }
214 :
215 : template <class T>
216 26838 : void MDAtomsTyped<T>::setBox(void*pp) {
217 26838 : box=static_cast<T*>(pp);
218 26838 : }
219 :
220 :
221 : template <class T>
222 30656 : void MDAtomsTyped<T>::setf(void*ff) {
223 : T*f=static_cast<T*>(ff);
224 30656 : plumed_assert(stride==0 || stride==3);
225 30656 : fx=f;
226 30656 : fy=f+1;
227 30656 : fz=f+2;
228 30656 : stride=3;
229 30656 : }
230 :
231 : template <class T>
232 0 : void MDAtomsTyped<T>::setp(void*pp,int i) {
233 : T*p=static_cast<T*>(pp);
234 0 : plumed_assert(stride==0 || stride==1);
235 0 : if(i==0)px=p;
236 0 : if(i==1)py=p;
237 0 : if(i==2)pz=p;
238 0 : stride=1;
239 0 : }
240 :
241 : template <class T>
242 24658 : void MDAtomsTyped<T>::setVirial(void*pp) {
243 24658 : virial=static_cast<T*>(pp);
244 24658 : }
245 :
246 :
247 : template <class T>
248 0 : void MDAtomsTyped<T>::setf(void*ff,int i) {
249 : T*f=static_cast<T*>(ff);
250 0 : plumed_assert(stride==0 || stride==1);
251 0 : if(i==0)fx=f;
252 0 : if(i==1)fy=f;
253 0 : if(i==2)fz=f;
254 0 : stride=1;
255 0 : }
256 :
257 : template <class T>
258 30656 : void MDAtomsTyped<T>::setm(void*m) {
259 30656 : this->m=static_cast<T*>(m);
260 30656 : }
261 :
262 : template <class T>
263 24595 : void MDAtomsTyped<T>::setc(void*c) {
264 24595 : this->c=static_cast<T*>(c);
265 24595 : }
266 :
267 : template <class T>
268 2788 : MDAtomsTyped<T>::MDAtomsTyped():
269 : scalep(1.0),
270 : scalef(1.0),
271 : scaleb(1.0),
272 : scalev(1.0),
273 : scalec(1.0),
274 : scalem(1.0),
275 : stride(0),
276 : m(NULL),
277 : c(NULL),
278 : px(NULL),
279 : py(NULL),
280 : pz(NULL),
281 : fx(NULL),
282 : fy(NULL),
283 : fz(NULL),
284 : box(NULL),
285 2788 : virial(NULL)
286 2788 : {}
287 :
288 2789 : MDAtomsBase* MDAtomsBase::create(unsigned p) {
289 2789 : if(p==sizeof(double)) {
290 2788 : return new MDAtomsTyped<double>;
291 1 : } else if (p==sizeof(float)) {
292 0 : return new MDAtomsTyped<float>;
293 : }
294 : std::string pp;
295 1 : Tools::convert(p,pp);
296 4 : plumed_merror("cannot create an MD interface with sizeof(real)=="+ pp);
297 : return NULL;
298 : }
299 :
300 : }
301 :
|
