Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 2902 : 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 3310 : void MD2double(const void*m,double&d)const {
62 3310 : d=double(*(static_cast<const T*>(m)));
63 3310 : }
64 191 : void double2MD(const double&d,void*m)const {
65 191 : *(static_cast<T*>(m))=T(d);
66 191 : }
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(unsigned j,unsigned k,vector<Vector>&positions)const;
74 : void getLocalPositions(std::vector<Vector>&p)const;
75 : void getMasses(const vector<int>&index,vector<double>&)const;
76 : void getCharges(const vector<int>&index,vector<double>&)const;
77 : void updateVirial(const Tensor&)const;
78 : void updateForces(const vector<int>&index,const vector<Vector>&);
79 : void rescaleForces(const vector<int>&index,double factor);
80 : unsigned getRealPrecision()const;
81 : };
82 :
83 : template <class T>
84 331 : void MDAtomsTyped<T>::setUnits(const Units& units,const Units& MDUnits) {
85 331 : double lscale=units.getLength()/MDUnits.getLength();
86 331 : double escale=units.getEnergy()/MDUnits.getEnergy();
87 331 : double cscale=units.getCharge()/MDUnits.getCharge();
88 331 : double mscale=units.getMass()/MDUnits.getMass();
89 : // scalep and scaleb are used to convert MD to plumed
90 331 : scalep=1.0/lscale;
91 331 : scaleb=1.0/lscale;
92 : // scalef and scalev are used to convert plumed to MD
93 331 : scalef=escale/lscale;
94 331 : scalev=escale;
95 331 : scalec=1.0/cscale;
96 331 : scalem=1.0/mscale;
97 331 : }
98 :
99 : template <class T>
100 40553 : void MDAtomsTyped<T>::getBox(Tensor&box)const {
101 40553 : 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;
102 28 : else box.zero();
103 40553 : }
104 :
105 : template <class T>
106 324 : void MDAtomsTyped<T>::getPositions(const vector<int>&index,vector<Vector>&positions)const {
107 : // cannot be parallelized with omp because access to positions is not ordered
108 14886 : for(unsigned i=0; i<index.size(); ++i) {
109 14562 : positions[index[i]][0]=px[stride*i]*scalep;
110 14562 : positions[index[i]][1]=py[stride*i]*scalep;
111 14562 : positions[index[i]][2]=pz[stride*i]*scalep;
112 : }
113 324 : }
114 :
115 : template <class T>
116 18559 : void MDAtomsTyped<T>::getPositions(unsigned j,unsigned k,vector<Vector>&positions)const {
117 38808 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(&positions[j],(k-j)))
118 20249 : for(unsigned i=j; i<k; ++i) {
119 1506758 : positions[i][0]=px[stride*i]*scalep;
120 1518144 : positions[i][1]=py[stride*i]*scalep;
121 1518695 : positions[i][2]=pz[stride*i]*scalep;
122 : }
123 18559 : }
124 :
125 :
126 : template <class T>
127 0 : void MDAtomsTyped<T>::getLocalPositions(vector<Vector>&positions)const {
128 0 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(positions))
129 0 : for(unsigned i=0; i<positions.size(); ++i) {
130 0 : positions[i][0]=px[stride*i]*scalep;
131 0 : positions[i][1]=py[stride*i]*scalep;
132 0 : positions[i][2]=pz[stride*i]*scalep;
133 : }
134 0 : }
135 :
136 :
137 : template <class T>
138 322 : void MDAtomsTyped<T>::getMasses(const vector<int>&index,vector<double>&masses)const {
139 322 : if(m) for(unsigned i=0; i<index.size(); ++i) masses[index[i]]=scalem*m[i];
140 0 : else for(unsigned i=0; i<index.size(); ++i) masses[index[i]]=0.0;
141 322 : }
142 :
143 : template <class T>
144 322 : void MDAtomsTyped<T>::getCharges(const vector<int>&index,vector<double>&charges)const {
145 322 : if(c) for(unsigned i=0; i<index.size(); ++i) charges[index[i]]=scalec*c[i];
146 36 : else for(unsigned i=0; i<index.size(); ++i) charges[index[i]]=0.0;
147 322 : }
148 :
149 : template <class T>
150 14811 : void MDAtomsTyped<T>::updateVirial(const Tensor&virial)const {
151 14811 : 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);
152 14811 : }
153 :
154 : template <class T>
155 19518 : void MDAtomsTyped<T>::updateForces(const vector<int>&index,const vector<Vector>&forces) {
156 40483 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(fx,stride*index.size()))
157 20965 : for(unsigned i=0; i<index.size(); ++i) {
158 1653020 : fx[stride*i]+=scalef*T(forces[index[i]][0]);
159 1668495 : fy[stride*i]+=scalef*T(forces[index[i]][1]);
160 1661560 : fz[stride*i]+=scalef*T(forces[index[i]][2]);
161 : }
162 19518 : }
163 :
164 : template <class T>
165 50 : void MDAtomsTyped<T>::rescaleForces(const vector<int>&index,double factor) {
166 50 : if(virial) for(unsigned i=0; i<3; i++)for(unsigned j=0; j<3; j++) virial[3*i+j]*=T(factor);
167 150 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(fx,stride*index.size()))
168 100 : for(unsigned i=0; i<index.size(); ++i) {
169 4201 : fx[stride*i]*=T(factor);
170 4201 : fy[stride*i]*=T(factor);
171 4201 : fz[stride*i]*=T(factor);
172 : }
173 50 : }
174 :
175 : template <class T>
176 331 : unsigned MDAtomsTyped<T>::getRealPrecision()const {
177 331 : return sizeof(T);
178 : }
179 :
180 : template <class T>
181 20991 : void MDAtomsTyped<T>::setp(void*pp) {
182 20991 : T*p=static_cast<T*>(pp);
183 20991 : plumed_assert(stride==0 || stride==3);
184 20991 : px=p;
185 20991 : py=p+1;
186 20991 : pz=p+2;
187 20991 : stride=3;
188 20991 : }
189 :
190 : template <class T>
191 20963 : void MDAtomsTyped<T>::setBox(void*pp) {
192 20963 : box=static_cast<T*>(pp);
193 20963 : }
194 :
195 :
196 : template <class T>
197 20991 : void MDAtomsTyped<T>::setf(void*ff) {
198 20991 : T*f=static_cast<T*>(ff);
199 20991 : plumed_assert(stride==0 || stride==3);
200 20991 : fx=f;
201 20991 : fy=f+1;
202 20991 : fz=f+2;
203 20991 : stride=3;
204 20991 : }
205 :
206 : template <class T>
207 0 : void MDAtomsTyped<T>::setp(void*pp,int i) {
208 0 : T*p=static_cast<T*>(pp);
209 0 : plumed_assert(stride==0 || stride==1);
210 0 : if(i==0)px=p;
211 0 : if(i==1)py=p;
212 0 : if(i==2)pz=p;
213 0 : stride=1;
214 0 : }
215 :
216 : template <class T>
217 18841 : void MDAtomsTyped<T>::setVirial(void*pp) {
218 18841 : virial=static_cast<T*>(pp);
219 18841 : }
220 :
221 :
222 : template <class T>
223 0 : void MDAtomsTyped<T>::setf(void*ff,int i) {
224 0 : T*f=static_cast<T*>(ff);
225 0 : plumed_assert(stride==0 || stride==1);
226 0 : if(i==0)fx=f;
227 0 : if(i==1)fy=f;
228 0 : if(i==2)fz=f;
229 0 : stride=1;
230 0 : }
231 :
232 : template <class T>
233 20991 : void MDAtomsTyped<T>::setm(void*m) {
234 20991 : this->m=static_cast<T*>(m);
235 20991 : }
236 :
237 : template <class T>
238 18773 : void MDAtomsTyped<T>::setc(void*c) {
239 18773 : this->c=static_cast<T*>(c);
240 18773 : }
241 :
242 : template <class T>
243 1451 : MDAtomsTyped<T>::MDAtomsTyped():
244 : scalep(1.0),
245 : scalef(1.0),
246 : scaleb(1.0),
247 : scalev(1.0),
248 : scalec(1.0),
249 : scalem(1.0),
250 : stride(0),
251 : m(NULL),
252 : c(NULL),
253 : px(NULL),
254 : py(NULL),
255 : pz(NULL),
256 : fx(NULL),
257 : fy(NULL),
258 : fz(NULL),
259 : box(NULL),
260 1451 : virial(NULL)
261 1451 : {}
262 :
263 1451 : MDAtomsBase* MDAtomsBase::create(unsigned p) {
264 1451 : if(p==sizeof(double)) {
265 1451 : return new MDAtomsTyped<double>;
266 0 : } else if (p==sizeof(float)) {
267 0 : return new MDAtomsTyped<float>;
268 : }
269 0 : std::string pp;
270 0 : Tools::convert(p,pp);
271 0 : plumed_merror("cannot create an MD interface with sizeof(real)=="+ pp);
272 0 : return NULL;
273 : }
274 :
275 2523 : }
276 :
|
