# | Line 1 | Line 1 | |
---|---|---|
1 | #include "Integrator.hpp" | |
2 | #include "simError.h" | |
3 | < | #include <cmath> |
4 | < | template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) |
5 | < | : T(theInfo, the_ff), fz(NULL), curZPos(NULL), |
6 | < | indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0) |
7 | < | { |
3 | > | #include <math.h> |
4 | ||
5 | + | const double INFINITE_TIME = 10e30; |
6 | + | template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, |
7 | + | ForceFields* the_ff): T(theInfo, the_ff), |
8 | + | fzOut(NULL), |
9 | + | curZconsTime(0), |
10 | + | forcePolicy(NULL), |
11 | + | usingSMD(false), |
12 | + | hasZConsGap(false){ |
13 | //get properties from SimInfo | |
14 | GenericData* data; | |
15 | ZConsParaData* zConsParaData; | |
16 | DoubleData* sampleTime; | |
17 | DoubleData* tolerance; | |
18 | + | DoubleData* gap; |
19 | + | DoubleData* fixtime; |
20 | StringData* policy; | |
21 | StringData* filename; | |
22 | + | IntData* smdFlag; |
23 | double COM[3]; | |
24 | ||
25 | //by default, the direction of constraint is z | |
# | Line 23 | Line 30 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
30 | ||
31 | //estimate the force constant of harmonical potential | |
32 | double Kb = 1.986E-3 ; //in kcal/K | |
26 | – | |
27 | – | double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2; |
28 | – | zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox); |
33 | ||
34 | < | //creat force substraction policy |
34 | > | double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) / |
35 | > | 2; |
36 | > | zForceConst = Kb * info->target_temp / (halfOfLargestBox * halfOfLargestBox); |
37 | > | |
38 | > | //creat force Subtraction policy |
39 | data = info->getProperty(ZCONSFORCEPOLICY_ID); | |
40 | < | if(!data){ |
41 | < | sprintf( painCave.errMsg, |
42 | < | "ZConstraint Warning: User does not set force substraction policy, " |
43 | < | "PolicyByMass is used\n"); |
40 | > | if (!data){ |
41 | > | sprintf(painCave.errMsg, |
42 | > | "ZConstraint Warning: User does not set force Subtraction policy, " |
43 | > | "PolicyByMass is used\n"); |
44 | painCave.isFatal = 0; | |
45 | simError(); | |
46 | ||
47 | < | forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this); |
47 | > | forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
48 | } | |
49 | else{ | |
50 | policy = dynamic_cast<StringData*>(data); | |
51 | < | |
52 | < | if(!policy){ |
53 | < | sprintf( painCave.errMsg, |
54 | < | "ZConstraint Error: Convertion from GenericData to StringData failure, " |
55 | < | "PolicyByMass is used\n"); |
51 | > | |
52 | > | if (!policy){ |
53 | > | sprintf(painCave.errMsg, |
54 | > | "ZConstraint Error: Convertion from GenericData to StringData failure, " |
55 | > | "PolicyByMass is used\n"); |
56 | painCave.isFatal = 0; | |
57 | simError(); | |
58 | ||
59 | < | forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this); |
59 | > | forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
60 | } | |
61 | else{ | |
62 | < | if(policy->getData() == "BYNUMBER") |
63 | < | forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this); |
64 | < | else if(policy->getData() == "BYMASS") |
65 | < | forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this); |
62 | > | if (policy->getData() == "BYNUMBER") |
63 | > | forcePolicy = (ForceSubtractionPolicy *) new PolicyByNumber(this); |
64 | > | else if (policy->getData() == "BYMASS") |
65 | > | forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
66 | else{ | |
67 | < | sprintf( painCave.errMsg, |
68 | < | "ZConstraint Warning: unknown force substraction policy, " |
69 | < | "PolicyByMass is used\n"); |
67 | > | sprintf(painCave.errMsg, |
68 | > | "ZConstraint Warning: unknown force Subtraction policy, " |
69 | > | "PolicyByMass is used\n"); |
70 | painCave.isFatal = 0; | |
71 | simError(); | |
72 | < | forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this); |
73 | < | } |
72 | > | forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
73 | > | } |
74 | } | |
75 | } | |
76 | < | |
77 | < | |
76 | > | |
77 | > | |
78 | //retrieve sample time of z-contraint | |
79 | data = info->getProperty(ZCONSTIME_ID); | |
80 | < | |
81 | < | if(!data) { |
82 | < | |
83 | < | sprintf( painCave.errMsg, |
84 | < | "ZConstraint error: If you use an ZConstraint\n" |
77 | < | " , you must set sample time.\n"); |
80 | > | |
81 | > | if (!data){ |
82 | > | sprintf(painCave.errMsg, |
83 | > | "ZConstraint error: If you use an ZConstraint\n" |
84 | > | " , you must set sample time.\n"); |
85 | painCave.isFatal = 1; | |
86 | < | simError(); |
86 | > | simError(); |
87 | } | |
88 | else{ | |
82 | – | |
89 | sampleTime = dynamic_cast<DoubleData*>(data); | |
84 | – | |
85 | – | if(!sampleTime){ |
90 | ||
91 | < | sprintf( painCave.errMsg, |
92 | < | "ZConstraint error: Can not get property from SimInfo\n"); |
91 | > | if (!sampleTime){ |
92 | > | sprintf(painCave.errMsg, |
93 | > | "ZConstraint error: Can not get property from SimInfo\n"); |
94 | painCave.isFatal = 1; | |
95 | < | simError(); |
91 | < | |
95 | > | simError(); |
96 | } | |
97 | else{ | |
98 | this->zconsTime = sampleTime->getData(); | |
99 | } | |
96 | – | |
100 | } | |
101 | < | |
101 | > | |
102 | //retrieve output filename of z force | |
103 | data = info->getProperty(ZCONSFILENAME_ID); | |
104 | < | if(!data) { |
105 | < | |
106 | < | |
107 | < | sprintf( painCave.errMsg, |
105 | < | "ZConstraint error: If you use an ZConstraint\n" |
106 | < | " , you must set output filename of z-force.\n"); |
104 | > | if (!data){ |
105 | > | sprintf(painCave.errMsg, |
106 | > | "ZConstraint error: If you use an ZConstraint\n" |
107 | > | " , you must set output filename of z-force.\n"); |
108 | painCave.isFatal = 1; | |
109 | < | simError(); |
109 | < | |
109 | > | simError(); |
110 | } | |
111 | else{ | |
112 | – | |
112 | filename = dynamic_cast<StringData*>(data); | |
114 | – | |
115 | – | if(!filename){ |
113 | ||
114 | < | sprintf( painCave.errMsg, |
115 | < | "ZConstraint error: Can not get property from SimInfo\n"); |
114 | > | if (!filename){ |
115 | > | sprintf(painCave.errMsg, |
116 | > | "ZConstraint error: Can not get property from SimInfo\n"); |
117 | painCave.isFatal = 1; | |
118 | < | simError(); |
121 | < | |
118 | > | simError(); |
119 | } | |
120 | else{ | |
121 | this->zconsOutput = filename->getData(); | |
122 | } | |
126 | – | |
123 | } | |
124 | ||
125 | //retrieve tolerance for z-constraint molecuels | |
126 | data = info->getProperty(ZCONSTOL_ID); | |
127 | < | |
128 | < | if(!data) { |
129 | < | |
134 | < | sprintf( painCave.errMsg, |
135 | < | "ZConstraint error: can not get tolerance \n"); |
127 | > | |
128 | > | if (!data){ |
129 | > | sprintf(painCave.errMsg, "ZConstraint error: can not get tolerance \n"); |
130 | painCave.isFatal = 1; | |
131 | < | simError(); |
131 | > | simError(); |
132 | } | |
133 | else{ | |
140 | – | |
134 | tolerance = dynamic_cast<DoubleData*>(data); | |
142 | – | |
143 | – | if(!tolerance){ |
135 | ||
136 | < | sprintf( painCave.errMsg, |
137 | < | "ZConstraint error: Can not get property from SimInfo\n"); |
136 | > | if (!tolerance){ |
137 | > | sprintf(painCave.errMsg, |
138 | > | "ZConstraint error: Can not get property from SimInfo\n"); |
139 | painCave.isFatal = 1; | |
140 | < | simError(); |
149 | < | |
140 | > | simError(); |
141 | } | |
142 | else{ | |
143 | this->zconsTol = tolerance->getData(); | |
144 | } | |
145 | + | } |
146 | ||
147 | + | //quick hack here |
148 | + | data = info->getProperty(ZCONSGAP_ID); |
149 | + | |
150 | + | if (data){ |
151 | + | gap = dynamic_cast<DoubleData*>(data); |
152 | + | |
153 | + | if (!gap){ |
154 | + | sprintf(painCave.errMsg, |
155 | + | "ZConstraint error: Can not get property from SimInfo\n"); |
156 | + | painCave.isFatal = 1; |
157 | + | simError(); |
158 | + | } |
159 | + | else{ |
160 | + | this->hasZConsGap = true; |
161 | + | this->zconsGap = gap->getData(); |
162 | + | } |
163 | } | |
156 | – | |
157 | – | //retrieve index of z-constraint molecules |
158 | – | data = info->getProperty(ZCONSPARADATA_ID); |
159 | – | if(!data) { |
164 | ||
165 | < | sprintf( painCave.errMsg, |
166 | < | "ZConstraint error: If you use an ZConstraint\n" |
167 | < | " , you must set index of z-constraint molecules.\n"); |
165 | > | |
166 | > | |
167 | > | data = info->getProperty(ZCONSFIXTIME_ID); |
168 | > | |
169 | > | if (data){ |
170 | > | fixtime = dynamic_cast<DoubleData*>(data); |
171 | > | if (!fixtime){ |
172 | > | sprintf(painCave.errMsg, |
173 | > | "ZConstraint error: Can not get zconsFixTime from SimInfo\n"); |
174 | > | painCave.isFatal = 1; |
175 | > | simError(); |
176 | > | } |
177 | > | else{ |
178 | > | this->zconsFixTime = fixtime->getData(); |
179 | > | } |
180 | > | } |
181 | > | else if(hasZConsGap){ |
182 | > | sprintf(painCave.errMsg, |
183 | > | "ZConstraint error: must set fixtime if already set zconsGap\n"); |
184 | > | painCave.isFatal = 1; |
185 | > | simError(); |
186 | > | } |
187 | > | |
188 | > | |
189 | > | |
190 | > | data = info->getProperty(ZCONSUSINGSMD_ID); |
191 | > | |
192 | > | if (data){ |
193 | > | smdFlag = dynamic_cast<IntData*>(data); |
194 | > | |
195 | > | if (!smdFlag){ |
196 | > | sprintf(painCave.errMsg, |
197 | > | "ZConstraint error: Can not get property from SimInfo\n"); |
198 | > | painCave.isFatal = 1; |
199 | > | simError(); |
200 | > | } |
201 | > | else{ |
202 | > | this->usingSMD= smdFlag->getData() ? true : false; |
203 | > | } |
204 | > | |
205 | > | } |
206 | > | |
207 | > | |
208 | > | |
209 | > | //retrieve index of z-constraint molecules |
210 | > | data = info->getProperty(ZCONSPARADATA_ID); |
211 | > | if (!data){ |
212 | > | sprintf(painCave.errMsg, |
213 | > | "ZConstraint error: If you use an ZConstraint\n" |
214 | > | " , you must set index of z-constraint molecules.\n"); |
215 | painCave.isFatal = 1; | |
216 | < | simError(); |
216 | > | simError(); |
217 | } | |
218 | else{ | |
168 | – | |
219 | zConsParaData = dynamic_cast<ZConsParaData*>(data); | |
170 | – | |
171 | – | if(!zConsParaData){ |
220 | ||
221 | < | sprintf( painCave.errMsg, |
222 | < | "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n"); |
221 | > | if (!zConsParaData){ |
222 | > | sprintf(painCave.errMsg, |
223 | > | "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n"); |
224 | painCave.isFatal = 1; | |
225 | < | simError(); |
177 | < | |
225 | > | simError(); |
226 | } | |
227 | else{ | |
180 | – | |
228 | parameters = zConsParaData->getData(); | |
229 | ||
230 | //check the range of zconsIndex | |
# | Line 189 | Line 236 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
236 | int totalNumMol; | |
237 | ||
238 | minIndex = (*parameters)[0].zconsIndex; | |
239 | < | if(minIndex < 0){ |
240 | < | sprintf( painCave.errMsg, |
194 | < | "ZConstraint error: index is out of range\n"); |
239 | > | if (minIndex < 0){ |
240 | > | sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); |
241 | painCave.isFatal = 1; | |
242 | < | simError(); |
243 | < | } |
242 | > | simError(); |
243 | > | } |
244 | ||
245 | maxIndex = (*parameters)[parameters->size() - 1].zconsIndex; | |
246 | ||
247 | #ifndef IS_MPI | |
248 | totalNumMol = nMols; | |
249 | #else | |
250 | < | totalNumMol = mpiSim->getTotNmol(); |
250 | > | totalNumMol = mpiSim->getNMolGlobal(); |
251 | #endif | |
252 | < | |
253 | < | if(maxIndex > totalNumMol - 1){ |
254 | < | sprintf( painCave.errMsg, |
209 | < | "ZConstraint error: index is out of range\n"); |
252 | > | |
253 | > | if (maxIndex > totalNumMol - 1){ |
254 | > | sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); |
255 | painCave.isFatal = 1; | |
256 | < | simError(); |
256 | > | simError(); |
257 | } | |
258 | ||
259 | //if user does not specify the zpos for the zconstraint molecule | |
260 | //its initial z coordinate will be used as default | |
261 | < | for(int i = 0; i < parameters->size(); i++){ |
262 | < | |
218 | < | if(!(*parameters)[i].havingZPos){ |
261 | > | for (int i = 0; i < (int) (parameters->size()); i++){ |
262 | > | if (!(*parameters)[i].havingZPos){ |
263 | #ifndef IS_MPI | |
264 | < | for(int j = 0; j < nMols; j++){ |
264 | > | for (int j = 0; j < nMols; j++){ |
265 | if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ | |
266 | molecules[j].getCOM(COM); | |
267 | < | break; |
267 | > | break; |
268 | } | |
269 | } | |
270 | #else | |
271 | < | //query which processor current zconstraint molecule belongs to |
272 | < | int *MolToProcMap; |
271 | > | //query which processor current zconstraint molecule belongs to |
272 | > | int* MolToProcMap; |
273 | int whichNode; | |
274 | < | double initZPos; |
274 | > | |
275 | MolToProcMap = mpiSim->getMolToProcMap(); | |
276 | whichNode = MolToProcMap[(*parameters)[i].zconsIndex]; | |
277 | < | |
277 | > | |
278 | //broadcast the zpos of current z-contraint molecule | |
279 | //the node which contain this | |
280 | < | |
281 | < | if (worldRank == whichNode ){ |
282 | < | |
239 | < | for(int j = 0; j < nMols; j++) |
280 | > | |
281 | > | if (worldRank == whichNode){ |
282 | > | for (int j = 0; j < nMols; j++) |
283 | if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ | |
284 | molecules[j].getCOM(COM); | |
285 | < | break; |
285 | > | break; |
286 | } | |
244 | – | |
287 | } | |
288 | ||
289 | < | MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD); |
289 | > | MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE, whichNode, |
290 | > | MPI_COMM_WORLD); |
291 | #endif | |
292 | < | |
292 | > | |
293 | (*parameters)[i].zPos = COM[whichDirection]; | |
294 | ||
295 | < | sprintf( painCave.errMsg, |
296 | < | "ZConstraint warning: Does not specify zpos for z-constraint molecule " |
297 | < | "initial z coornidate will be used \n"); |
298 | < | painCave.isFatal = 0; |
299 | < | simError(); |
300 | < | |
258 | < | } |
295 | > | sprintf(painCave.errMsg, |
296 | > | "ZConstraint warning: Does not specify zpos for z-constraint molecule " |
297 | > | "initial z coornidate will be used \n"); |
298 | > | painCave.isFatal = 0; |
299 | > | simError(); |
300 | > | } |
301 | } | |
260 | – | |
302 | }//end if (!zConsParaData) | |
303 | + | |
304 | }//end if (!data) | |
305 | < | |
306 | < | // |
305 | > | |
306 | > | // |
307 | #ifdef IS_MPI | |
308 | update(); | |
309 | #else | |
310 | int searchResult; | |
311 | < | |
312 | < | for(int i = 0; i < nMols; i++){ |
271 | < | |
311 | > | |
312 | > | for (int i = 0; i < nMols; i++){ |
313 | searchResult = isZConstraintMol(&molecules[i]); | |
314 | < | |
315 | < | if(searchResult > -1){ |
275 | < | |
314 | > | |
315 | > | if (searchResult > -1){ |
316 | zconsMols.push_back(&molecules[i]); | |
317 | massOfZConsMols.push_back(molecules[i].getTotalMass()); | |
318 | ||
319 | zPos.push_back((*parameters)[searchResult].zPos); | |
280 | – | // cout << "index: "<< (*parameters)[searchResult].zconsIndex |
281 | – | // <<"\tzPos = " << (*parameters)[searchResult].zPos << endl; |
320 | kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); | |
321 | ||
322 | < | molecules[i].getCOM(COM); |
322 | > | if(usingSMD) |
323 | > | cantVel.push_back((*parameters)[searchResult].cantVel); |
324 | > | |
325 | } | |
326 | < | else |
287 | < | { |
288 | < | |
326 | > | else{ |
327 | unconsMols.push_back(&molecules[i]); | |
328 | massOfUnconsMols.push_back(molecules[i].getTotalMass()); | |
291 | – | |
329 | } | |
330 | } | |
331 | ||
332 | < | fz = new double[zconsMols.size()]; |
333 | < | curZPos = new double[zconsMols.size()]; |
334 | < | indexOfZConsMols = new int [zconsMols.size()]; |
332 | > | fz.resize(zconsMols.size()); |
333 | > | curZPos.resize(zconsMols.size()); |
334 | > | indexOfZConsMols.resize(zconsMols.size()); |
335 | ||
299 | – | if(!fz || !curZPos || !indexOfZConsMols){ |
300 | – | sprintf( painCave.errMsg, |
301 | – | "Memory allocation failure in class Zconstraint\n"); |
302 | – | painCave.isFatal = 1; |
303 | – | simError(); |
304 | – | } |
305 | – | |
336 | //determine the states of z-constraint molecules | |
337 | < | for(int i = 0; i < zconsMols.size(); i++){ |
337 | > | for (size_t i = 0; i < zconsMols.size(); i++){ |
338 | indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); | |
339 | ||
340 | zconsMols[i]->getCOM(COM); | |
341 | < | if (fabs(zPos[i] - COM[whichDirection]) < zconsTol) |
341 | > | |
342 | > | if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){ |
343 | states.push_back(zcsFixed); | |
344 | < | else |
344 | > | |
345 | > | if (hasZConsGap) |
346 | > | endFixTime.push_back(info->getTime() + zconsFixTime); |
347 | > | } |
348 | > | else{ |
349 | states.push_back(zcsMoving); | |
350 | + | |
351 | + | if (hasZConsGap) |
352 | + | endFixTime.push_back(INFINITE_TIME); |
353 | + | } |
354 | + | |
355 | + | if(usingSMD) |
356 | + | cantPos.push_back(COM[whichDirection]); |
357 | } | |
358 | < | |
358 | > | |
359 | > | if(usingSMD) |
360 | > | prevCantPos = cantPos; |
361 | #endif | |
362 | ||
363 | + | |
364 | //get total masss of unconstraint molecules | |
365 | double totalMassOfUncons_local; | |
366 | totalMassOfUncons_local = 0; | |
367 | < | |
368 | < | for(int i = 0; i < unconsMols.size(); i++) |
367 | > | |
368 | > | for (size_t i = 0; i < unconsMols.size(); i++) |
369 | totalMassOfUncons_local += unconsMols[i]->getTotalMass(); | |
370 | < | |
370 | > | |
371 | #ifndef IS_MPI | |
372 | totalMassOfUncons = totalMassOfUncons_local; | |
373 | #else | |
374 | < | MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, |
375 | < | MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
374 | > | MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE, |
375 | > | MPI_SUM, MPI_COMM_WORLD); |
376 | #endif | |
377 | ||
333 | – | |
378 | //get total number of unconstrained atoms | |
379 | int nUnconsAtoms_local; | |
380 | nUnconsAtoms_local = 0; | |
381 | < | for(int i = 0; i < unconsMols.size(); i++) |
381 | > | for (int i = 0; i < (int) (unconsMols.size()); i++) |
382 | nUnconsAtoms_local += unconsMols[i]->getNAtoms(); | |
383 | < | |
383 | > | |
384 | #ifndef IS_MPI | |
385 | totNumOfUnconsAtoms = nUnconsAtoms_local; | |
386 | #else | |
387 | < | MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, |
388 | < | MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
387 | > | MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM, |
388 | > | MPI_COMM_WORLD); |
389 | #endif | |
390 | ||
347 | – | // creat zconsWriter |
348 | – | fzOut = new ZConsWriter(zconsOutput.c_str(), parameters); |
349 | – | |
350 | – | if(!fzOut){ |
351 | – | sprintf( painCave.errMsg, |
352 | – | "Memory allocation failure in class Zconstraint\n"); |
353 | – | painCave.isFatal = 1; |
354 | – | simError(); |
355 | – | } |
356 | – | |
391 | forcePolicy->update(); | |
392 | } | |
393 | ||
394 | < | template<typename T> ZConstraint<T>::~ZConstraint() |
361 | < | { |
362 | < | if(fz) |
363 | < | delete[] fz; |
394 | > | template<typename T> ZConstraint<T>::~ZConstraint(){ |
395 | ||
396 | < | if(curZPos) |
366 | < | delete[] curZPos; |
367 | < | |
368 | < | if(indexOfZConsMols) |
369 | < | delete[] indexOfZConsMols; |
370 | < | |
371 | < | if(fzOut) |
396 | > | if (fzOut){ |
397 | delete fzOut; | |
398 | < | |
399 | < | if(forcePolicy) |
398 | > | } |
399 | > | |
400 | > | if (forcePolicy){ |
401 | delete forcePolicy; | |
402 | + | } |
403 | } | |
404 | ||
405 | ||
# | Line 381 | Line 408 | template<typename T> ZConstraint<T>::~ZConstraint() | |
408 | */ | |
409 | ||
410 | #ifdef IS_MPI | |
411 | < | template<typename T> void ZConstraint<T>::update() |
385 | < | { |
411 | > | template<typename T> void ZConstraint<T>::update(){ |
412 | double COM[3]; | |
413 | int index; | |
414 | < | |
414 | > | |
415 | zconsMols.clear(); | |
416 | massOfZConsMols.clear(); | |
417 | zPos.clear(); | |
418 | kz.clear(); | |
419 | < | |
419 | > | cantPos.clear(); |
420 | > | cantVel.clear(); |
421 | > | |
422 | unconsMols.clear(); | |
423 | massOfUnconsMols.clear(); | |
396 | – | |
424 | ||
425 | + | |
426 | //creat zconsMol and unconsMol lists | |
427 | < | for(int i = 0; i < nMols; i++){ |
400 | < | |
427 | > | for (int i = 0; i < nMols; i++){ |
428 | index = isZConstraintMol(&molecules[i]); | |
429 | < | |
430 | < | if(index > -1){ |
404 | < | |
429 | > | |
430 | > | if (index > -1){ |
431 | zconsMols.push_back(&molecules[i]); | |
432 | zPos.push_back((*parameters)[index].zPos); | |
433 | kz.push_back((*parameters)[index].kRatio * zForceConst); | |
434 | + | massOfZConsMols.push_back(molecules[i].getTotalMass()); |
435 | ||
436 | < | massOfZConsMols.push_back(molecules[i].getTotalMass()); |
437 | < | |
438 | < | molecules[i].getCOM(COM); |
436 | > | if(usingSMD) |
437 | > | cantVel.push_back((*parameters)[index].cantVel); |
438 | > | |
439 | } | |
440 | < | else |
414 | < | { |
415 | < | |
440 | > | else{ |
441 | unconsMols.push_back(&molecules[i]); | |
442 | massOfUnconsMols.push_back(molecules[i].getTotalMass()); | |
418 | – | |
443 | } | |
444 | } | |
445 | ||
446 | + | fz.resize(zconsMols.size()); |
447 | + | curZPos.resize(zconsMols.size()); |
448 | + | indexOfZConsMols.resize(zconsMols.size()); |
449 | + | |
450 | + | for (size_t i = 0; i < zconsMols.size(); i++){ |
451 | + | indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); |
452 | + | } |
453 | + | |
454 | //determine the states of z-constraint molecules | |
455 | < | for(int i = 0; i < zconsMols.size(); i++){ |
456 | < | zconsMols[i]->getCOM(COM); |
457 | < | if (fabs(zPos[i] - COM[whichDirection]) < zconsTol) |
455 | > | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
456 | > | |
457 | > | zconsMols[i]->getCOM(COM); |
458 | > | |
459 | > | if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){ |
460 | states.push_back(zcsFixed); | |
461 | < | else |
461 | > | |
462 | > | if (hasZConsGap) |
463 | > | endFixTime.push_back(info->getTime() + zconsFixTime); |
464 | > | } |
465 | > | else{ |
466 | states.push_back(zcsMoving); | |
429 | – | } |
467 | ||
468 | < | |
469 | < | //The reason to declare fz and indexOfZconsMols as pointer to array is |
433 | < | // that we want to make the MPI communication simple |
434 | < | if(fz) |
435 | < | delete[] fz; |
436 | < | |
437 | < | if(curZPos) |
438 | < | delete[] curZPos; |
439 | < | |
440 | < | if(indexOfZConsMols) |
441 | < | delete[] indexOfZConsMols; |
442 | < | |
443 | < | if (zconsMols.size() > 0){ |
444 | < | fz = new double[zconsMols.size()]; |
445 | < | curZPos = new double[zconsMols.size()]; |
446 | < | indexOfZConsMols = new int[zconsMols.size()]; |
447 | < | |
448 | < | if(!fz || !curZPos || !indexOfZConsMols){ |
449 | < | sprintf( painCave.errMsg, |
450 | < | "Memory allocation failure in class Zconstraint\n"); |
451 | < | painCave.isFatal = 1; |
452 | < | simError(); |
468 | > | if (hasZConsGap) |
469 | > | endFixTime.push_back(INFINITE_TIME); |
470 | } | |
454 | – | |
455 | – | for(int i = 0; i < zconsMols.size(); i++){ |
456 | – | indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); |
457 | – | } |
471 | ||
472 | < | } |
473 | < | else{ |
474 | < | fz = NULL; |
475 | < | curZPos = NULL; |
476 | < | indexOfZConsMols = NULL; |
477 | < | } |
478 | < | |
472 | > | if(usingSMD) |
473 | > | cantPos.push_back(COM[whichDirection]); |
474 | > | } |
475 | > | |
476 | > | if(usingSMD) |
477 | > | prevCantPos = cantPos; |
478 | > | |
479 | // | |
480 | forcePolicy->update(); | |
468 | – | |
481 | } | |
482 | ||
483 | #endif | |
# | Line 479 | Line 491 | template<typename T> void ZConstraint<T>::update() | |
491 | * other non-negative values, its index in indexOfAllZConsMols vector | |
492 | */ | |
493 | ||
494 | < | template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol) |
483 | < | { |
494 | > | template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){ |
495 | int index; | |
496 | int low; | |
497 | int high; | |
498 | int mid; | |
499 | ||
500 | index = mol->getGlobalIndex(); | |
501 | < | |
501 | > | |
502 | low = 0; | |
503 | high = parameters->size() - 1; | |
504 | < | |
504 | > | |
505 | //Binary Search (we have sorted the array) | |
506 | < | while(low <= high){ |
507 | < | mid = (low + high) /2; |
506 | > | while (low <= high){ |
507 | > | mid = (low + high) / 2; |
508 | if ((*parameters)[mid].zconsIndex == index) | |
509 | return mid; | |
510 | < | else if ((*parameters)[mid].zconsIndex > index ) |
511 | < | high = mid -1; |
512 | < | else |
513 | < | low = mid + 1; |
510 | > | else if ((*parameters)[mid].zconsIndex > index) |
511 | > | high = mid - 1; |
512 | > | else |
513 | > | low = mid + 1; |
514 | } | |
515 | < | |
515 | > | |
516 | return -1; | |
517 | } | |
518 | ||
519 | template<typename T> void ZConstraint<T>::integrate(){ | |
520 | < | |
520 | > | // creat zconsWriter |
521 | > | fzOut = new ZConsWriter(zconsOutput.c_str(), parameters); |
522 | > | |
523 | > | if (!fzOut){ |
524 | > | sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n"); |
525 | > | painCave.isFatal = 1; |
526 | > | simError(); |
527 | > | } |
528 | > | |
529 | //zero out the velocities of center of mass of unconstrained molecules | |
530 | //and the velocities of center of mass of every single z-constrained molecueles | |
531 | zeroOutVel(); | |
532 | ||
533 | curZconsTime = zconsTime + info->getTime(); | |
515 | – | |
516 | – | T::integrate(); |
534 | ||
535 | + | T::integrate(); |
536 | } | |
519 | – | |
537 | ||
538 | + | |
539 | /** | |
540 | * | |
541 | * | |
# | Line 532 | Line 550 | template<typename T> void ZConstraint<T>::calcForce(in | |
550 | ||
551 | T::calcForce(calcPot, calcStress); | |
552 | ||
535 | – | if (checkZConsState()){ |
536 | – | |
537 | – | #ifdef IS_MPI |
538 | – | if(worldRank == 0){ |
539 | – | #endif |
540 | – | // std::cerr << "\n" |
541 | – | // << "*******************************************\n" |
542 | – | // << " about to call zeroOutVel()\n" |
543 | – | // << "*******************************************\n" |
544 | – | // << "\n"; |
545 | – | #ifdef IS_MPI |
546 | – | } |
547 | – | #endif |
548 | – | zeroOutVel(); |
553 | ||
554 | < | #ifdef IS_MPI |
555 | < | if(worldRank == 0){ |
556 | < | #endif |
557 | < | // std::cerr << "\n" |
558 | < | // << "*******************************************\n" |
559 | < | // << " finished zeroOutVel()\n" |
556 | < | // << "*******************************************\n" |
557 | < | // << "\n"; |
558 | < | #ifdef IS_MPI |
559 | < | } |
560 | < | #endif |
561 | < | |
554 | > | if (hasZConsGap){ |
555 | > | updateZPos(); |
556 | > | } |
557 | > | |
558 | > | if (checkZConsState()){ |
559 | > | zeroOutVel(); |
560 | forcePolicy->update(); | |
561 | } | |
562 | < | |
562 | > | |
563 | zsys = calcZSys(); | |
564 | zSysCOMVel = calcSysCOMVel(); | |
565 | #ifdef IS_MPI | |
566 | < | if(worldRank == 0){ |
566 | > | if (worldRank == 0){ |
567 | #endif | |
568 | < | // cout << "---------------------------------------------------------------------" <<endl; |
569 | < | // cout << "current time: " << info->getTime() << endl; |
570 | < | // cout << "center of mass at z: " << zsys << endl; |
571 | < | // cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl; |
568 | > | //cout << "---------------------------------------------------------------------" <<endl; |
569 | > | //cout << "current time: " << info->getTime() << endl; |
570 | > | //cout << "center of mass at z: " << zsys << endl; |
571 | > | //cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl; |
572 | ||
573 | #ifdef IS_MPI | |
574 | } | |
575 | #endif | |
576 | ||
577 | //do zconstraint force; | |
578 | < | if (haveFixedZMols()) |
578 | > | if (haveFixedZMols()){ |
579 | this->doZconstraintForce(); | |
580 | + | } |
581 | ||
582 | < | //use harmonical poteintial to move the molecules to the specified positions |
583 | < | if (haveMovingZMols()) |
584 | < | this->doHarmonic(); |
582 | > | //use external force to move the molecules to the specified positions |
583 | > | if (haveMovingZMols()){ |
584 | > | if (usingSMD) |
585 | > | this->doHarmonic(cantPos); |
586 | > | else |
587 | > | this->doHarmonic(zPos); |
588 | > | } |
589 | ||
590 | //write out forces and current positions of z-constraint molecules | |
591 | < | if(info->getTime() >= curZconsTime){ |
592 | < | for(int i = 0; i < zconsMols.size(); i++){ |
591 | > | if (info->getTime() >= curZconsTime){ |
592 | > | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
593 | zconsMols[i]->getCOM(COM); | |
594 | < | curZPos[i] = COM[whichDirection]; |
594 | > | curZPos[i] = COM[whichDirection]; |
595 | ||
596 | < | //if the z-constraint molecule is still moving, just record its force |
597 | < | if(states[i] == zcsMoving){ |
598 | < | fz[i] = 0; |
599 | < | Atom** movingZAtoms; |
600 | < | movingZAtoms = zconsMols[i]->getMyAtoms(); |
601 | < | for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
602 | < | movingZAtoms[j]->getFrc(force); |
603 | < | fz[i] += force[whichDirection]; |
596 | > | //if the z-constraint molecule is still moving, just record its force |
597 | > | if (states[i] == zcsMoving){ |
598 | > | fz[i] = 0; |
599 | > | Atom** movingZAtoms; |
600 | > | movingZAtoms = zconsMols[i]->getMyAtoms(); |
601 | > | for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
602 | > | movingZAtoms[j]->getFrc(force); |
603 | > | fz[i] += force[whichDirection]; |
604 | > | } |
605 | } | |
606 | < | } |
607 | < | } |
608 | < | fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos); |
609 | < | curZconsTime += zconsTime; |
606 | > | } |
607 | > | fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0], |
608 | > | &curZPos[0], &zPos[0]); |
609 | > | curZconsTime += zconsTime; |
610 | } | |
611 | ||
612 | zSysCOMVel = calcSysCOMVel(); | |
613 | #ifdef IS_MPI | |
614 | < | if(worldRank == 0){ |
614 | > | if (worldRank == 0){ |
615 | #endif | |
616 | < | // cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl; |
616 | > | //cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl; |
617 | #ifdef IS_MPI | |
618 | } | |
619 | #endif | |
# | Line 619 | Line 623 | template<typename T> void ZConstraint<T>::calcForce(in | |
623 | /** | |
624 | * | |
625 | */ | |
626 | < | |
627 | < | template<typename T> double ZConstraint<T>::calcZSys() |
624 | < | { |
626 | > | |
627 | > | template<typename T> double ZConstraint<T>::calcZSys(){ |
628 | //calculate reference z coordinate for z-constraint molecules | |
629 | double totalMass_local; | |
630 | double totalMass; | |
# | Line 629 | Line 632 | template<typename T> double ZConstraint<T>::calcZSys() | |
632 | double totalMZ; | |
633 | double massOfCurMol; | |
634 | double COM[3]; | |
635 | < | |
635 | > | |
636 | totalMass_local = 0; | |
637 | totalMZ_local = 0; | |
638 | < | |
639 | < | for(int i = 0; i < nMols; i++){ |
638 | > | |
639 | > | for (int i = 0; i < nMols; i++){ |
640 | massOfCurMol = molecules[i].getTotalMass(); | |
641 | molecules[i].getCOM(COM); | |
642 | < | |
642 | > | |
643 | totalMass_local += massOfCurMol; | |
644 | totalMZ_local += massOfCurMol * COM[whichDirection]; | |
642 | – | |
645 | } | |
646 | < | |
647 | < | |
646 | > | |
647 | > | |
648 | #ifdef IS_MPI | |
649 | < | MPI_Allreduce(&totalMass_local, &totalMass, 1, |
650 | < | MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
651 | < | MPI_Allreduce(&totalMZ_local, &totalMZ, 1, |
652 | < | MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
651 | < | #else |
649 | > | MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM, |
650 | > | MPI_COMM_WORLD); |
651 | > | MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
652 | > | #else |
653 | totalMass = totalMass_local; | |
654 | totalMZ = totalMZ_local; | |
655 | < | #endif |
655 | > | #endif |
656 | ||
657 | double zsys; | |
658 | zsys = totalMZ / totalMass; | |
# | Line 662 | Line 663 | template<typename T> double ZConstraint<T>::calcZSys() | |
663 | /** | |
664 | * | |
665 | */ | |
666 | < | template<typename T> void ZConstraint<T>::thermalize( void ){ |
666 | < | |
666 | > | template<typename T> void ZConstraint<T>::thermalize(void){ |
667 | T::thermalize(); | |
668 | zeroOutVel(); | |
669 | } | |
# | Line 673 | Line 673 | template<typename T> void ZConstraint<T>::zeroOutVel() | |
673 | */ | |
674 | ||
675 | template<typename T> void ZConstraint<T>::zeroOutVel(){ | |
676 | – | |
676 | Atom** fixedZAtoms; | |
677 | double COMvel[3]; | |
678 | double vel[3]; | |
679 | double zSysCOMVel; | |
680 | ||
681 | //zero out the velocities of center of mass of fixed z-constrained molecules | |
683 | – | |
684 | – | for(int i = 0; i < zconsMols.size(); i++){ |
682 | ||
683 | < | if (states[i] == zcsFixed){ |
683 | > | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
684 | > | if (states[i] == zcsFixed){ |
685 | > | zconsMols[i]->getCOMvel(COMvel); |
686 | > | //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
687 | ||
688 | – | zconsMols[i]->getCOMvel(COMvel); |
689 | – | //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
690 | – | |
688 | fixedZAtoms = zconsMols[i]->getMyAtoms(); | |
689 | < | |
690 | < | for(int j =0; j < zconsMols[i]->getNAtoms(); j++){ |
689 | > | |
690 | > | for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
691 | fixedZAtoms[j]->getVel(vel); | |
692 | < | vel[whichDirection] -= COMvel[whichDirection]; |
693 | < | fixedZAtoms[j]->setVel(vel); |
692 | > | vel[whichDirection] -= COMvel[whichDirection]; |
693 | > | fixedZAtoms[j]->setVel(vel); |
694 | } | |
695 | ||
696 | < | zconsMols[i]->getCOMvel(COMvel); |
697 | < | //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
696 | > | zconsMols[i]->getCOMvel(COMvel); |
697 | > | //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
698 | } | |
702 | – | |
699 | } | |
700 | ||
701 | < | //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; |
701 | > | //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; |
702 | ||
703 | zSysCOMVel = calcSysCOMVel(); | |
704 | #ifdef IS_MPI | |
705 | < | if(worldRank == 0){ |
705 | > | if (worldRank == 0){ |
706 | #endif | |
707 | < | // cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl; |
707 | > | //cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl; |
708 | #ifdef IS_MPI | |
709 | } | |
710 | #endif | |
711 | < | |
711 | > | |
712 | // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules | |
713 | double MVzOfMovingMols_local; | |
714 | double MVzOfMovingMols; | |
715 | double totalMassOfMovingZMols_local; | |
716 | double totalMassOfMovingZMols; | |
717 | < | |
717 | > | |
718 | MVzOfMovingMols_local = 0; | |
719 | totalMassOfMovingZMols_local = 0; | |
720 | ||
721 | < | for(int i =0; i < unconsMols.size(); i++){ |
721 | > | for (int i = 0; i < (int) (unconsMols.size()); i++){ |
722 | unconsMols[i]->getCOMvel(COMvel); | |
723 | < | MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
723 | > | MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
724 | } | |
725 | ||
726 | < | for(int i = 0; i < zconsMols.size(); i++){ |
726 | > | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
727 | if (states[i] == zcsMoving){ | |
728 | zconsMols[i]->getCOMvel(COMvel); | |
729 | MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; | |
730 | < | totalMassOfMovingZMols_local += massOfZConsMols[i]; |
730 | > | totalMassOfMovingZMols_local += massOfZConsMols[i]; |
731 | } | |
736 | – | |
732 | } | |
733 | ||
734 | #ifndef IS_MPI | |
735 | MVzOfMovingMols = MVzOfMovingMols_local; | |
736 | totalMassOfMovingZMols = totalMassOfMovingZMols_local; | |
737 | #else | |
738 | < | MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
739 | < | MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
738 | > | MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE, |
739 | > | MPI_SUM, MPI_COMM_WORLD); |
740 | > | MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, |
741 | > | MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
742 | #endif | |
743 | ||
744 | double vzOfMovingMols; | |
745 | < | vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols); |
745 | > | vzOfMovingMols = MVzOfMovingMols / |
746 | > | (totalMassOfUncons + totalMassOfMovingZMols); |
747 | ||
748 | //modify the velocites of unconstrained molecules | |
749 | Atom** unconsAtoms; | |
750 | < | for(int i = 0; i < unconsMols.size(); i++){ |
753 | < | |
750 | > | for (int i = 0; i < (int) (unconsMols.size()); i++){ |
751 | unconsAtoms = unconsMols[i]->getMyAtoms(); | |
752 | < | for(int j = 0; j < unconsMols[i]->getNAtoms();j++){ |
752 | > | for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
753 | unconsAtoms[j]->getVel(vel); | |
754 | vel[whichDirection] -= vzOfMovingMols; | |
755 | unconsAtoms[j]->setVel(vel); | |
756 | } | |
760 | – | |
757 | } | |
758 | ||
759 | //modify the velocities of moving z-constrained molecuels | |
760 | Atom** movingZAtoms; | |
761 | < | for(int i = 0; i < zconsMols.size(); i++){ |
762 | < | |
767 | < | if (states[i] ==zcsMoving){ |
768 | < | |
761 | > | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
762 | > | if (states[i] == zcsMoving){ |
763 | movingZAtoms = zconsMols[i]->getMyAtoms(); | |
764 | < | for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
764 | > | for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
765 | movingZAtoms[j]->getVel(vel); | |
766 | vel[whichDirection] -= vzOfMovingMols; | |
767 | < | movingZAtoms[j]->setVel(vel); |
767 | > | movingZAtoms[j]->setVel(vel); |
768 | > | } |
769 | } | |
770 | < | |
776 | < | } |
770 | > | } |
771 | ||
778 | – | } |
772 | ||
780 | – | |
773 | zSysCOMVel = calcSysCOMVel(); | |
774 | #ifdef IS_MPI | |
775 | < | if(worldRank == 0){ |
775 | > | if (worldRank == 0){ |
776 | #endif | |
777 | < | // cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl; |
777 | > | //cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl; |
778 | #ifdef IS_MPI | |
779 | } | |
780 | #endif | |
789 | – | |
781 | } | |
782 | ||
783 | /** | |
# | Line 794 | Line 785 | template<typename T> void ZConstraint<T>::doZconstrain | |
785 | */ | |
786 | ||
787 | template<typename T> void ZConstraint<T>::doZconstraintForce(){ | |
797 | – | |
788 | Atom** zconsAtoms; | |
789 | double totalFZ; | |
790 | double totalFZ_local; | |
801 | – | double COMvel[3]; |
791 | double COM[3]; | |
792 | double force[3]; | |
793 | ||
794 | //constrain the molecules which do not reach the specified positions | |
795 | < | |
795 | > | |
796 | //Zero Out the force of z-contrained molecules | |
797 | totalFZ_local = 0; | |
798 | ||
799 | //calculate the total z-contrained force of fixed z-contrained molecules | |
800 | ||
801 | < | for(int i = 0; i < zconsMols.size(); i++){ |
802 | < | |
801 | > | //cout << "before zero out z-constraint force on fixed z-constraint molecuels " |
802 | > | // << "total force is " << calcTotalForce() << endl; |
803 | > | |
804 | > | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
805 | if (states[i] == zcsFixed){ | |
815 | – | |
806 | zconsMols[i]->getCOM(COM); | |
807 | zconsAtoms = zconsMols[i]->getMyAtoms(); | |
808 | ||
809 | fz[i] = 0; | |
810 | < | for(int j =0; j < zconsMols[i]->getNAtoms(); j++) { |
810 | > | for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
811 | zconsAtoms[j]->getFrc(force); | |
812 | < | fz[i] += force[whichDirection]; |
812 | > | fz[i] += force[whichDirection]; |
813 | } | |
814 | totalFZ_local += fz[i]; | |
815 | ||
816 | //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i] | |
817 | // <<"\tcurrent zpos: " << COM[whichDirection] | |
818 | < | // << "\tcurrent fz: " <<fz[i] << endl; |
829 | < | |
830 | < | |
818 | > | // << "\tcurrent fz: " <<fz[i] << endl; |
819 | } | |
832 | – | |
820 | } | |
821 | ||
822 | //calculate total z-constraint force | |
823 | #ifdef IS_MPI | |
824 | < | MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
824 | > | MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
825 | #else | |
826 | totalFZ = totalFZ_local; | |
827 | #endif | |
828 | ||
829 | < | |
829 | > | |
830 | // apply negative to fixed z-constrained molecues; | |
831 | < | force[0]= 0; |
832 | < | force[1]= 0; |
833 | < | force[2]= 0; |
831 | > | force[0] = 0; |
832 | > | force[1] = 0; |
833 | > | force[2] = 0; |
834 | ||
835 | < | for(int i = 0; i < zconsMols.size(); i++){ |
836 | < | |
850 | < | if (states[i] == zcsFixed){ |
851 | < | |
835 | > | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
836 | > | if (states[i] == zcsFixed){ |
837 | int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms(); | |
838 | zconsAtoms = zconsMols[i]->getMyAtoms(); | |
839 | < | |
840 | < | for(int j =0; j < nAtomOfCurZConsMol; j++) { |
839 | > | |
840 | > | for (int j = 0; j < nAtomOfCurZConsMol; j++){ |
841 | //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; | |
842 | < | force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]); |
842 | > | force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i], |
843 | > | zconsAtoms[j], |
844 | > | fz[i]); |
845 | zconsAtoms[j]->addFrc(force); | |
846 | } | |
860 | – | |
847 | } | |
862 | – | |
848 | } | |
849 | ||
850 | < | // cout << "after zero out z-constraint force on fixed z-constraint molecuels " |
851 | < | // << "total force is " << calcTotalForce() << endl; |
850 | > | //cout << "after zero out z-constraint force on fixed z-constraint molecuels " |
851 | > | // << "total force is " << calcTotalForce() << endl; |
852 | ||
868 | – | force[0]= 0; |
869 | – | force[1]= 0; |
870 | – | force[2]= 0; |
853 | ||
854 | + | force[0] = 0; |
855 | + | force[1] = 0; |
856 | + | force[2] = 0; |
857 | + | |
858 | //modify the forces of unconstrained molecules | |
859 | < | for(int i = 0; i < unconsMols.size(); i++){ |
860 | < | |
861 | < | Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
862 | < | |
863 | < | for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
864 | < | //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); |
865 | < | force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ); |
866 | < | unconsAtoms[j]->addFrc(force); |
867 | < | } |
882 | < | |
859 | > | for (int i = 0; i < (int) (unconsMols.size()); i++){ |
860 | > | Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
861 | > | |
862 | > | for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
863 | > | //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); |
864 | > | force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j], |
865 | > | totalFZ); |
866 | > | unconsAtoms[j]->addFrc(force); |
867 | > | } |
868 | } | |
869 | ||
870 | < | //modify the forces of moving z-constrained molecules |
871 | < | for(int i = 0; i < zconsMols.size(); i++) { |
870 | > | //modify the forces of moving z-constrained molecules |
871 | > | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
872 | if (states[i] == zcsMoving){ | |
888 | – | |
873 | Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); | |
874 | ||
875 | < | for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
875 | > | for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
876 | //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); | |
877 | < | force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ); |
877 | > | force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j], |
878 | > | totalFZ); |
879 | movingZAtoms[j]->addFrc(force); | |
880 | } | |
881 | } | |
882 | } | |
883 | < | |
884 | < | //cout << "after substracting z-constraint force from moving molecuels " |
900 | < | // << "total force is " << calcTotalForce() << endl; |
901 | < | |
883 | > | // cout << "after substracting z-constraint force from moving molecuels " |
884 | > | // << "total force is " << calcTotalForce() << endl; |
885 | } | |
886 | ||
887 | /** | |
# | Line 906 | Line 889 | template<typename T> void ZConstraint<T>::doZconstrain | |
889 | * | |
890 | */ | |
891 | ||
892 | < | template<typename T> void ZConstraint<T>::doHarmonic(){ |
892 | > | template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){ |
893 | double force[3]; | |
894 | double harmonicU; | |
895 | double harmonicF; | |
# | Line 914 | Line 897 | template<typename T> void ZConstraint<T>::doHarmonic() | |
897 | double diff; | |
898 | double totalFZ_local; | |
899 | double totalFZ; | |
900 | < | |
900 | > | |
901 | force[0] = 0; | |
902 | force[1] = 0; | |
903 | force[2] = 0; | |
904 | ||
905 | totalFZ_local = 0; | |
906 | ||
907 | < | for(int i = 0; i < zconsMols.size(); i++) { |
925 | < | |
907 | > | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
908 | if (states[i] == zcsMoving){ | |
909 | zconsMols[i]->getCOM(COM); | |
910 | < | // cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] |
911 | < | // << "\tcurrent zpos: " << COM[whichDirection] << endl; |
910 | > | // cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] |
911 | > | // << "\tcurrent zpos: " << COM[whichDirection] << endl; |
912 | ||
913 | < | diff = COM[whichDirection] -zPos[i]; |
914 | < | |
913 | > | diff = COM[whichDirection] - resPos[i]; |
914 | > | |
915 | harmonicU = 0.5 * kz[i] * diff * diff; | |
916 | info->lrPot += harmonicU; | |
917 | ||
918 | < | harmonicF = - kz[i] * diff; |
918 | > | harmonicF = -kz[i] * diff; |
919 | totalFZ_local += harmonicF; | |
920 | ||
921 | < | //adjust force |
922 | < | |
921 | > | //adjust force |
922 | > | |
923 | Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); | |
924 | ||
925 | < | for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
925 | > | for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
926 | //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms(); | |
927 | < | force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF); |
928 | < | movingZAtoms[j]->addFrc(force); |
929 | < | } |
927 | > | force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], |
928 | > | movingZAtoms[j], |
929 | > | harmonicF); |
930 | > | movingZAtoms[j]->addFrc(force); |
931 | > | } |
932 | } | |
949 | – | |
933 | } | |
934 | ||
935 | #ifndef IS_MPI | |
936 | totalFZ = totalFZ_local; | |
937 | #else | |
938 | < | MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
938 | > | MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
939 | #endif | |
940 | ||
941 | < | force[0]= 0; |
942 | < | force[1]= 0; |
960 | < | force[2]= 0; |
941 | > | //cout << "before substracting harmonic force from moving molecuels " |
942 | > | // << "total force is " << calcTotalForce() << endl; |
943 | ||
944 | + | force[0] = 0; |
945 | + | force[1] = 0; |
946 | + | force[2] = 0; |
947 | + | |
948 | //modify the forces of unconstrained molecules | |
949 | < | for(int i = 0; i < unconsMols.size(); i++){ |
950 | < | |
951 | < | Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
952 | < | |
953 | < | for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
954 | < | //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms; |
955 | < | force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ); |
956 | < | unconsAtoms[j]->addFrc(force); |
957 | < | } |
949 | > | for (int i = 0; i < (int) (unconsMols.size()); i++){ |
950 | > | Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
951 | > | |
952 | > | for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
953 | > | //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms; |
954 | > | force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j], |
955 | > | totalFZ); |
956 | > | unconsAtoms[j]->addFrc(force); |
957 | > | } |
958 | } | |
959 | ||
960 | + | //cout << "after substracting harmonic force from moving molecuels " |
961 | + | // << "total force is " << calcTotalForce() << endl; |
962 | } | |
963 | ||
964 | /** | |
# | Line 980 | Line 968 | template<typename T> bool ZConstraint<T>::checkZConsSt | |
968 | template<typename T> bool ZConstraint<T>::checkZConsState(){ | |
969 | double COM[3]; | |
970 | double diff; | |
971 | < | |
971 | > | |
972 | int changed_local; | |
973 | int changed; | |
974 | < | |
974 | > | |
975 | changed_local = 0; | |
988 | – | |
989 | – | for(int i =0; i < zconsMols.size(); i++){ |
976 | ||
977 | + | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
978 | zconsMols[i]->getCOM(COM); | |
979 | diff = fabs(COM[whichDirection] - zPos[i]); | |
980 | < | if ( diff <= zconsTol && states[i] == zcsMoving){ |
980 | > | if (diff <= zconsTol && states[i] == zcsMoving){ |
981 | states[i] = zcsFixed; | |
982 | < | changed_local = 1; |
982 | > | changed_local = 1; |
983 | > | |
984 | > | if(usingSMD) |
985 | > | prevCantPos = cantPos; |
986 | > | |
987 | > | if (hasZConsGap) |
988 | > | endFixTime[i] = info->getTime() + zconsFixTime; |
989 | } | |
990 | < | else if ( diff > zconsTol && states[i] == zcsFixed){ |
990 | > | else if (diff > zconsTol && states[i] == zcsFixed){ |
991 | states[i] = zcsMoving; | |
992 | < | changed_local = 1; |
992 | > | changed_local = 1; |
993 | > | |
994 | > | if(usingSMD) |
995 | > | cantPos = prevCantPos; |
996 | > | |
997 | > | if (hasZConsGap) |
998 | > | endFixTime[i] = INFINITE_TIME; |
999 | } | |
1001 | – | |
1000 | } | |
1001 | ||
1002 | #ifndef IS_MPI | |
1003 | < | changed =changed_local; |
1003 | > | changed = changed_local; |
1004 | #else | |
1005 | < | MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
1005 | > | MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
1006 | #endif | |
1007 | ||
1008 | return (changed > 0); | |
1011 | – | |
1009 | } | |
1010 | ||
1011 | template<typename T> bool ZConstraint<T>::haveFixedZMols(){ | |
1015 | – | |
1012 | int havingFixed_local; | |
1013 | int havingFixed; | |
1014 | ||
1015 | havingFixed_local = 0; | |
1016 | ||
1017 | < | for(int i = 0; i < zconsMols.size(); i++) |
1017 | > | for (int i = 0; i < (int) (zconsMols.size()); i++) |
1018 | if (states[i] == zcsFixed){ | |
1019 | havingFixed_local = 1; | |
1020 | < | break; |
1020 | > | break; |
1021 | } | |
1022 | ||
1023 | #ifndef IS_MPI | |
1024 | havingFixed = havingFixed_local; | |
1025 | #else | |
1026 | < | MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
1026 | > | MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM, |
1027 | > | MPI_COMM_WORLD); |
1028 | #endif | |
1029 | ||
1030 | return (havingFixed > 0); | |
# | Line 1038 | Line 1035 | template<typename T> bool ZConstraint<T>::haveMovingZM | |
1035 | * | |
1036 | */ | |
1037 | template<typename T> bool ZConstraint<T>::haveMovingZMols(){ | |
1041 | – | |
1038 | int havingMoving_local; | |
1039 | int havingMoving; | |
1040 | ||
1041 | havingMoving_local = 0; | |
1042 | ||
1043 | < | for(int i = 0; i < zconsMols.size(); i++) |
1043 | > | for (int i = 0; i < (int) (zconsMols.size()); i++) |
1044 | if (states[i] == zcsMoving){ | |
1045 | havingMoving_local = 1; | |
1046 | < | break; |
1046 | > | break; |
1047 | } | |
1048 | ||
1049 | #ifndef IS_MPI | |
1050 | havingMoving = havingMoving_local; | |
1051 | #else | |
1052 | < | MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
1052 | > | MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM, |
1053 | > | MPI_COMM_WORLD); |
1054 | #endif | |
1055 | ||
1056 | return (havingMoving > 0); | |
1060 | – | |
1057 | } | |
1058 | ||
1059 | /** | |
1060 | * | |
1061 | */ | |
1062 | ||
1063 | < | template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel() |
1068 | < | { |
1063 | > | template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){ |
1064 | double MVzOfMovingMols_local; | |
1065 | double MVzOfMovingMols; | |
1066 | double totalMassOfMovingZMols_local; | |
1067 | double totalMassOfMovingZMols; | |
1068 | double COMvel[3]; | |
1069 | < | |
1069 | > | |
1070 | MVzOfMovingMols_local = 0; | |
1071 | totalMassOfMovingZMols_local = 0; | |
1072 | ||
1073 | < | for(int i =0; i < unconsMols.size(); i++){ |
1073 | > | for (int i = 0; i < unconsMols.size(); i++){ |
1074 | unconsMols[i]->getCOMvel(COMvel); | |
1075 | < | MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
1075 | > | MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
1076 | } | |
1077 | ||
1078 | < | for(int i = 0; i < zconsMols.size(); i++){ |
1084 | < | |
1078 | > | for (int i = 0; i < zconsMols.size(); i++){ |
1079 | if (states[i] == zcsMoving){ | |
1080 | zconsMols[i]->getCOMvel(COMvel); | |
1081 | MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; | |
1082 | < | totalMassOfMovingZMols_local += massOfZConsMols[i]; |
1082 | > | totalMassOfMovingZMols_local += massOfZConsMols[i]; |
1083 | } | |
1090 | – | |
1084 | } | |
1085 | ||
1086 | #ifndef IS_MPI | |
1087 | MVzOfMovingMols = MVzOfMovingMols_local; | |
1088 | totalMassOfMovingZMols = totalMassOfMovingZMols_local; | |
1089 | #else | |
1090 | < | MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1091 | < | MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1090 | > | MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE, |
1091 | > | MPI_SUM, MPI_COMM_WORLD); |
1092 | > | MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, |
1093 | > | MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
1094 | #endif | |
1095 | ||
1096 | double vzOfMovingMols; | |
1097 | < | vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols); |
1097 | > | vzOfMovingMols = MVzOfMovingMols / |
1098 | > | (totalMassOfUncons + totalMassOfMovingZMols); |
1099 | ||
1100 | return vzOfMovingMols; | |
1101 | } | |
# | Line 1108 | Line 1104 | template<typename T> double ZConstraint<T>::calcMoving | |
1104 | * | |
1105 | */ | |
1106 | ||
1107 | < | template<typename T> double ZConstraint<T>::calcSysCOMVel() |
1112 | < | { |
1107 | > | template<typename T> double ZConstraint<T>::calcSysCOMVel(){ |
1108 | double COMvel[3]; | |
1109 | double tempMVz_local; | |
1110 | double tempMVz; | |
# | Line 1117 | Line 1112 | template<typename T> double ZConstraint<T>::calcSysCOM | |
1112 | double massOfZCons; | |
1113 | ||
1114 | ||
1115 | < | tempMVz_local = 0; |
1115 | > | tempMVz_local = 0; |
1116 | ||
1117 | < | for(int i =0 ; i < nMols; i++){ |
1117 | > | for (int i = 0 ; i < nMols; i++){ |
1118 | molecules[i].getCOMvel(COMvel); | |
1119 | < | tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection]; |
1119 | > | tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection]; |
1120 | } | |
1121 | ||
1122 | massOfZCons_local = 0; | |
1123 | < | |
1124 | < | for(int i = 0; i < massOfZConsMols.size(); i++){ |
1123 | > | |
1124 | > | for (int i = 0; i < (int) (massOfZConsMols.size()); i++){ |
1125 | massOfZCons_local += massOfZConsMols[i]; | |
1126 | } | |
1127 | #ifndef IS_MPI | |
1128 | massOfZCons = massOfZCons_local; | |
1129 | tempMVz = tempMVz_local; | |
1130 | #else | |
1131 | < | MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1132 | < | MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1131 | > | MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM, |
1132 | > | MPI_COMM_WORLD); |
1133 | > | MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
1134 | #endif | |
1135 | ||
1136 | < | return tempMVz /(totalMassOfUncons + massOfZCons); |
1136 | > | return tempMVz / (totalMassOfUncons + massOfZCons); |
1137 | } | |
1138 | ||
1139 | /** | |
# | Line 1145 | Line 1141 | template<typename T> double ZConstraint<T>::calcTotalF | |
1141 | */ | |
1142 | ||
1143 | template<typename T> double ZConstraint<T>::calcTotalForce(){ | |
1148 | – | |
1144 | double force[3]; | |
1145 | double totalForce_local; | |
1146 | double totalForce; | |
1147 | ||
1148 | totalForce_local = 0; | |
1149 | ||
1150 | < | for(int i = 0; i < nAtoms; i++){ |
1150 | > | for (int i = 0; i < nAtoms; i++){ |
1151 | atoms[i]->getFrc(force); | |
1152 | totalForce_local += force[whichDirection]; | |
1153 | } | |
# | Line 1160 | Line 1155 | template<typename T> double ZConstraint<T>::calcTotalF | |
1155 | #ifndef IS_MPI | |
1156 | totalForce = totalForce_local; | |
1157 | #else | |
1158 | < | MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1158 | > | MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM, |
1159 | > | MPI_COMM_WORLD); |
1160 | #endif | |
1161 | ||
1162 | return totalForce; | |
1167 | – | |
1163 | } | |
1164 | ||
1165 | /** | |
# | Line 1175 | Line 1170 | template<typename T> void ZConstraint<T>::PolicyByNumb | |
1170 | //calculate the number of atoms of moving z-constrained molecules | |
1171 | int nMovingZAtoms_local; | |
1172 | int nMovingZAtoms; | |
1173 | < | |
1173 | > | |
1174 | nMovingZAtoms_local = 0; | |
1175 | < | for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++) |
1176 | < | if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)) |
1177 | < | nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms(); |
1178 | < | |
1175 | > | for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++) |
1176 | > | if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){ |
1177 | > | nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms(); |
1178 | > | } |
1179 | > | |
1180 | #ifdef IS_MPI | |
1181 | < | MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
1181 | > | MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, |
1182 | > | MPI_COMM_WORLD); |
1183 | #else | |
1184 | nMovingZAtoms = nMovingZAtoms_local; | |
1185 | #endif | |
1186 | totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms; | |
1190 | – | |
1191 | – | #ifdef IS_MPI |
1192 | – | if(worldRank == 0){ |
1193 | – | #endif |
1194 | – | // std::cerr << "\n" |
1195 | – | // << "*******************************************\n" |
1196 | – | // << " fiished Policy by numbr()\n" |
1197 | – | // << "*******************************************\n" |
1198 | – | // << "\n"; |
1199 | – | #ifdef IS_MPI |
1200 | – | } |
1201 | – | #endif |
1187 | } | |
1188 | ||
1189 | < | template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){ |
1189 | > | template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, |
1190 | > | Atom* atom, |
1191 | > | double totalForce){ |
1192 | return totalForce / mol->getNAtoms(); | |
1193 | } | |
1194 | ||
1195 | < | template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, double totalForce){ |
1195 | > | template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, |
1196 | > | double totalForce){ |
1197 | return totalForce / totNumOfMovingAtoms; | |
1198 | } | |
1199 | ||
1200 | < | template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){ |
1201 | < | return totalForce / mol->getNAtoms(); |
1200 | > | template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, |
1201 | > | Atom* atom, |
1202 | > | double totalForce){ |
1203 | > | return totalForce / mol->getNAtoms(); |
1204 | } | |
1205 | ||
1206 | < | template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, double totalForce){ |
1206 | > | template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, |
1207 | > | double totalForce){ |
1208 | return totalForce / zconsIntegrator->totNumOfUnconsAtoms; | |
1209 | } | |
1210 | ||
# | Line 1225 | Line 1216 | template<typename T> void ZConstraint<T>::PolicyByMass | |
1216 | //calculate the number of atoms of moving z-constrained molecules | |
1217 | double massOfMovingZAtoms_local; | |
1218 | double massOfMovingZAtoms; | |
1219 | < | |
1219 | > | |
1220 | massOfMovingZAtoms_local = 0; | |
1221 | < | for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++) |
1222 | < | if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)) |
1223 | < | massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass(); |
1224 | < | |
1221 | > | for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++) |
1222 | > | if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){ |
1223 | > | massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass(); |
1224 | > | } |
1225 | > | |
1226 | #ifdef IS_MPI | |
1227 | < | MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1227 | > | MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE, |
1228 | > | MPI_SUM, MPI_COMM_WORLD); |
1229 | #else | |
1230 | massOfMovingZAtoms = massOfMovingZAtoms_local; | |
1231 | #endif | |
1232 | < | totMassOfMovingAtoms = massOfMovingZAtoms_local + zconsIntegrator->totalMassOfUncons; |
1232 | > | totMassOfMovingAtoms = massOfMovingZAtoms + |
1233 | > | zconsIntegrator->totalMassOfUncons; |
1234 | } | |
1235 | ||
1236 | < | template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){ |
1236 | > | template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, |
1237 | > | Atom* atom, |
1238 | > | double totalForce){ |
1239 | return totalForce * atom->getMass() / mol->getTotalMass(); | |
1240 | } | |
1241 | ||
1242 | < | template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols( Atom* atom, double totalForce){ |
1243 | < | return totalForce * atom->getMass() / totMassOfMovingAtoms; |
1242 | > | template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom, |
1243 | > | double totalForce){ |
1244 | > | return totalForce * atom->getMass() / totMassOfMovingAtoms; |
1245 | } | |
1246 | ||
1247 | < | template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){ |
1247 | > | template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, |
1248 | > | Atom* atom, |
1249 | > | double totalForce){ |
1250 | return totalForce * atom->getMass() / mol->getTotalMass(); | |
1251 | } | |
1252 | ||
1253 | < | template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, double totalForce){ |
1254 | < | return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons; |
1253 | > | template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, |
1254 | > | double totalForce){ |
1255 | > | return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons; |
1256 | } | |
1257 | ||
1258 | + | template<typename T> void ZConstraint<T>::updateZPos(){ |
1259 | + | double curTime; |
1260 | + | double COM[3]; |
1261 | + | |
1262 | + | curTime = info->getTime(); |
1263 | + | |
1264 | + | for (size_t i = 0; i < zconsMols.size(); i++){ |
1265 | + | |
1266 | + | if (states[i] == zcsFixed && curTime >= endFixTime[i]){ |
1267 | + | zPos[i] += zconsGap; |
1268 | + | |
1269 | + | if (usingSMD){ |
1270 | + | zconsMols[i]->getCOM(COM); |
1271 | + | cantPos[i] = COM[whichDirection]; |
1272 | + | } |
1273 | + | |
1274 | + | } |
1275 | + | |
1276 | + | } |
1277 | + | |
1278 | + | } |
1279 | + | |
1280 | + | template<typename T> void ZConstraint<T>::updateCantPos(){ |
1281 | + | double curTime; |
1282 | + | double dt; |
1283 | + | |
1284 | + | curTime = info->getTime(); |
1285 | + | dt = info->dt; |
1286 | + | |
1287 | + | for (size_t i = 0; i < zconsMols.size(); i++){ |
1288 | + | if (states[i] == zcsMoving){ |
1289 | + | cantPos[i] += cantVel[i] * dt; |
1290 | + | } |
1291 | + | } |
1292 | + | |
1293 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |