1 |
< |
#include "integrators/Integrator.hpp" |
2 |
< |
#include "utils/simError.h" |
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 |
< |
DoubleGenericData* sampleTime; |
17 |
< |
DoubleGenericData* tolerance; |
18 |
< |
DoubleGenericData* gap; |
19 |
< |
DoubleGenericData* fixtime; |
20 |
< |
StringGenericData* policy; |
21 |
< |
StringGenericData* filename; |
22 |
< |
IntGenericData* smdFlag; |
23 |
< |
double COM[3]; |
24 |
< |
|
25 |
< |
//by default, the direction of constraint is z |
26 |
< |
// 0 --> x |
27 |
< |
// 1 --> y |
28 |
< |
// 2 --> z |
29 |
< |
whichDirection = 2; |
30 |
< |
|
31 |
< |
//estimate the force constant of harmonical potential |
32 |
< |
double Kb = 1.986E-3 ; //in kcal/K |
33 |
< |
|
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->getPropertyByName(ZCONSFORCEPOLICY_ID); |
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 = (ForceSubtractionPolicy *) new PolicyByMass(this); |
48 |
< |
} |
49 |
< |
else{ |
50 |
< |
policy = dynamic_cast<StringGenericData*>(data); |
51 |
< |
|
52 |
< |
if (!policy){ |
53 |
< |
sprintf(painCave.errMsg, |
54 |
< |
"ZConstraint Error: Convertion from GenericData to StringGenericData failure, " |
55 |
< |
"PolicyByMass is used\n"); |
56 |
< |
painCave.isFatal = 0; |
57 |
< |
simError(); |
58 |
< |
|
59 |
< |
forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
60 |
< |
} |
61 |
< |
else{ |
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 Subtraction policy, " |
69 |
< |
"PolicyByMass is used\n"); |
70 |
< |
painCave.isFatal = 0; |
71 |
< |
simError(); |
72 |
< |
forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
73 |
< |
} |
74 |
< |
} |
75 |
< |
} |
76 |
< |
|
77 |
< |
|
78 |
< |
//retrieve sample time of z-contraint |
79 |
< |
data = info->getPropertyByName(ZCONSTIME_ID); |
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(); |
87 |
< |
} |
88 |
< |
else{ |
89 |
< |
sampleTime = dynamic_cast<DoubleGenericData*>(data); |
90 |
< |
|
91 |
< |
if (!sampleTime){ |
92 |
< |
sprintf(painCave.errMsg, |
93 |
< |
"ZConstraint error: Can not get property from SimInfo\n"); |
94 |
< |
painCave.isFatal = 1; |
95 |
< |
simError(); |
96 |
< |
} |
97 |
< |
else{ |
98 |
< |
this->zconsTime = sampleTime->getData(); |
99 |
< |
} |
100 |
< |
} |
101 |
< |
|
102 |
< |
//retrieve output filename of z force |
103 |
< |
data = info->getPropertyByName(ZCONSFILENAME_ID); |
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(); |
110 |
< |
} |
111 |
< |
else{ |
112 |
< |
filename = dynamic_cast<StringGenericData*>(data); |
113 |
< |
|
114 |
< |
if (!filename){ |
115 |
< |
sprintf(painCave.errMsg, |
116 |
< |
"ZConstraint error: Can not get property from SimInfo\n"); |
117 |
< |
painCave.isFatal = 1; |
118 |
< |
simError(); |
119 |
< |
} |
120 |
< |
else{ |
121 |
< |
this->zconsOutput = filename->getData(); |
122 |
< |
} |
123 |
< |
} |
124 |
< |
|
125 |
< |
//retrieve tolerance for z-constraint molecuels |
126 |
< |
data = info->getPropertyByName(ZCONSTOL_ID); |
127 |
< |
|
128 |
< |
if (!data){ |
129 |
< |
sprintf(painCave.errMsg, "ZConstraint error: can not get tolerance \n"); |
130 |
< |
painCave.isFatal = 1; |
131 |
< |
simError(); |
132 |
< |
} |
133 |
< |
else{ |
134 |
< |
tolerance = dynamic_cast<DoubleGenericData*>(data); |
135 |
< |
|
136 |
< |
if (!tolerance){ |
137 |
< |
sprintf(painCave.errMsg, |
138 |
< |
"ZConstraint error: Can not get property from SimInfo\n"); |
139 |
< |
painCave.isFatal = 1; |
140 |
< |
simError(); |
141 |
< |
} |
142 |
< |
else{ |
143 |
< |
this->zconsTol = tolerance->getData(); |
144 |
< |
} |
145 |
< |
} |
146 |
< |
|
147 |
< |
//quick hack here |
148 |
< |
data = info->getPropertyByName(ZCONSGAP_ID); |
149 |
< |
|
150 |
< |
if (data){ |
151 |
< |
gap = dynamic_cast<DoubleGenericData*>(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 |
< |
} |
164 |
< |
|
165 |
< |
|
166 |
< |
|
167 |
< |
data = info->getPropertyByName(ZCONSFIXTIME_ID); |
168 |
< |
|
169 |
< |
if (data){ |
170 |
< |
fixtime = dynamic_cast<DoubleGenericData*>(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->getPropertyByName(ZCONSUSINGSMD_ID); |
191 |
< |
|
192 |
< |
if (data){ |
193 |
< |
smdFlag = dynamic_cast<IntGenericData*>(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->getPropertyByName(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(); |
217 |
< |
} |
218 |
< |
else{ |
219 |
< |
zConsParaData = dynamic_cast<ZConsParaData*>(data); |
220 |
< |
|
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(); |
226 |
< |
} |
227 |
< |
else{ |
228 |
< |
parameters = zConsParaData->getData(); |
229 |
< |
|
230 |
< |
//check the range of zconsIndex |
231 |
< |
//and the minimum value of index is the first one (we already sorted the data) |
232 |
< |
//the maximum value of index is the last one |
233 |
< |
|
234 |
< |
int maxIndex; |
235 |
< |
int minIndex; |
236 |
< |
int totalNumMol; |
237 |
< |
|
238 |
< |
minIndex = (*parameters)[0].zconsIndex; |
239 |
< |
if (minIndex < 0){ |
240 |
< |
sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); |
241 |
< |
painCave.isFatal = 1; |
242 |
< |
simError(); |
243 |
< |
} |
244 |
< |
|
245 |
< |
maxIndex = (*parameters)[parameters->size() - 1].zconsIndex; |
246 |
< |
|
247 |
< |
#ifndef IS_MPI |
248 |
< |
totalNumMol = nMols; |
249 |
< |
#else |
250 |
< |
totalNumMol = mpiSim->getNMolGlobal(); |
251 |
< |
#endif |
252 |
< |
|
253 |
< |
if (maxIndex > totalNumMol - 1){ |
254 |
< |
sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); |
255 |
< |
painCave.isFatal = 1; |
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 < (int) (parameters->size()); i++){ |
262 |
< |
if (!(*parameters)[i].havingZPos){ |
263 |
< |
#ifndef IS_MPI |
264 |
< |
for (int j = 0; j < nMols; j++){ |
265 |
< |
if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ |
266 |
< |
molecules[j].getCOM(COM); |
267 |
< |
break; |
268 |
< |
} |
269 |
< |
} |
270 |
< |
#else |
271 |
< |
//query which processor current zconstraint molecule belongs to |
272 |
< |
int* MolToProcMap; |
273 |
< |
int whichNode; |
274 |
< |
|
275 |
< |
MolToProcMap = mpiSim->getMolToProcMap(); |
276 |
< |
whichNode = MolToProcMap[(*parameters)[i].zconsIndex]; |
277 |
< |
|
278 |
< |
//broadcast the zpos of current z-contraint molecule |
279 |
< |
//the node which contain this |
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; |
286 |
< |
} |
287 |
< |
} |
288 |
< |
|
289 |
< |
MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE, whichNode, |
290 |
< |
MPI_COMM_WORLD); |
291 |
< |
#endif |
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 |
< |
} |
301 |
< |
} |
302 |
< |
}//end if (!zConsParaData) |
303 |
< |
|
304 |
< |
}//end if (!data) |
305 |
< |
|
306 |
< |
// |
307 |
< |
#ifdef IS_MPI |
308 |
< |
update(); |
309 |
< |
#else |
310 |
< |
int searchResult; |
311 |
< |
|
312 |
< |
for (int i = 0; i < nMols; i++){ |
313 |
< |
searchResult = isZConstraintMol(&molecules[i]); |
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); |
320 |
< |
kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); |
321 |
< |
|
322 |
< |
if(usingSMD) |
323 |
< |
cantVel.push_back((*parameters)[searchResult].cantVel); |
324 |
< |
|
325 |
< |
} |
326 |
< |
else{ |
327 |
< |
unconsMols.push_back(&molecules[i]); |
328 |
< |
massOfUnconsMols.push_back(molecules[i].getTotalMass()); |
329 |
< |
} |
330 |
< |
} |
331 |
< |
|
332 |
< |
fz.resize(zconsMols.size()); |
333 |
< |
curZPos.resize(zconsMols.size()); |
334 |
< |
indexOfZConsMols.resize(zconsMols.size()); |
335 |
< |
|
336 |
< |
//determine the states of z-constraint molecules |
337 |
< |
for (size_t i = 0; i < zconsMols.size(); i++){ |
338 |
< |
indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); |
339 |
< |
|
340 |
< |
zconsMols[i]->getCOM(COM); |
341 |
< |
|
342 |
< |
if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){ |
343 |
< |
states.push_back(zcsFixed); |
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 |
< |
|
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 (size_t i = 0; i < unconsMols.size(); i++) |
369 |
< |
totalMassOfUncons_local += unconsMols[i]->getTotalMass(); |
370 |
< |
|
371 |
< |
#ifndef IS_MPI |
372 |
< |
totalMassOfUncons = totalMassOfUncons_local; |
373 |
< |
#else |
374 |
< |
MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE, |
375 |
< |
MPI_SUM, MPI_COMM_WORLD); |
376 |
< |
#endif |
377 |
< |
|
378 |
< |
//get total number of unconstrained atoms |
379 |
< |
int nUnconsAtoms_local; |
380 |
< |
nUnconsAtoms_local = 0; |
381 |
< |
for (int i = 0; i < (int) (unconsMols.size()); i++) |
382 |
< |
nUnconsAtoms_local += unconsMols[i]->getNAtoms(); |
383 |
< |
|
384 |
< |
#ifndef IS_MPI |
385 |
< |
totNumOfUnconsAtoms = nUnconsAtoms_local; |
386 |
< |
#else |
387 |
< |
MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM, |
388 |
< |
MPI_COMM_WORLD); |
389 |
< |
#endif |
390 |
< |
|
391 |
< |
forcePolicy->update(); |
392 |
< |
} |
393 |
< |
|
394 |
< |
template<typename T> ZConstraint<T>::~ZConstraint(){ |
395 |
< |
|
396 |
< |
if (fzOut){ |
397 |
< |
delete fzOut; |
398 |
< |
} |
399 |
< |
|
400 |
< |
if (forcePolicy){ |
401 |
< |
delete forcePolicy; |
402 |
< |
} |
403 |
< |
} |
404 |
< |
|
405 |
< |
|
406 |
< |
/** |
407 |
< |
* |
408 |
< |
*/ |
409 |
< |
|
410 |
< |
#ifdef IS_MPI |
411 |
< |
template<typename T> void ZConstraint<T>::update(){ |
412 |
< |
double COM[3]; |
413 |
< |
int index; |
414 |
< |
|
415 |
< |
zconsMols.clear(); |
416 |
< |
massOfZConsMols.clear(); |
417 |
< |
zPos.clear(); |
418 |
< |
kz.clear(); |
419 |
< |
cantPos.clear(); |
420 |
< |
cantVel.clear(); |
421 |
< |
|
422 |
< |
unconsMols.clear(); |
423 |
< |
massOfUnconsMols.clear(); |
424 |
< |
|
425 |
< |
|
426 |
< |
//creat zconsMol and unconsMol lists |
427 |
< |
for (int i = 0; i < nMols; i++){ |
428 |
< |
index = isZConstraintMol(&molecules[i]); |
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 |
< |
if(usingSMD) |
437 |
< |
cantVel.push_back((*parameters)[index].cantVel); |
438 |
< |
|
439 |
< |
} |
440 |
< |
else{ |
441 |
< |
unconsMols.push_back(&molecules[i]); |
442 |
< |
massOfUnconsMols.push_back(molecules[i].getTotalMass()); |
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 < (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 |
< |
|
462 |
< |
if (hasZConsGap) |
463 |
< |
endFixTime.push_back(info->getTime() + zconsFixTime); |
464 |
< |
} |
465 |
< |
else{ |
466 |
< |
states.push_back(zcsMoving); |
467 |
< |
|
468 |
< |
if (hasZConsGap) |
469 |
< |
endFixTime.push_back(INFINITE_TIME); |
470 |
< |
} |
471 |
< |
|
472 |
< |
if(usingSMD) |
473 |
< |
cantPos.push_back(COM[whichDirection]); |
474 |
< |
} |
475 |
< |
|
476 |
< |
if(usingSMD) |
477 |
< |
prevCantPos = cantPos; |
478 |
< |
|
479 |
< |
forcePolicy->update(); |
480 |
< |
} |
481 |
< |
|
482 |
< |
#endif |
483 |
< |
|
484 |
< |
/** |
485 |
< |
* Function Name: isZConstraintMol |
486 |
< |
* Parameter |
487 |
< |
* Molecule* mol |
488 |
< |
* Return value: |
489 |
< |
* -1, if the molecule is not z-constraint molecule, |
490 |
< |
* other non-negative values, its index in indexOfAllZConsMols vector |
491 |
< |
*/ |
492 |
< |
|
493 |
< |
template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){ |
494 |
< |
int index; |
495 |
< |
int low; |
496 |
< |
int high; |
497 |
< |
int mid; |
498 |
< |
|
499 |
< |
index = mol->getGlobalIndex(); |
500 |
< |
|
501 |
< |
low = 0; |
502 |
< |
high = parameters->size() - 1; |
503 |
< |
|
504 |
< |
//Binary Search (we have sorted the array) |
505 |
< |
while (low <= high){ |
506 |
< |
mid = (low + high) / 2; |
507 |
< |
if ((*parameters)[mid].zconsIndex == index) |
508 |
< |
return mid; |
509 |
< |
else if ((*parameters)[mid].zconsIndex > index) |
510 |
< |
high = mid - 1; |
511 |
< |
else |
512 |
< |
low = mid + 1; |
513 |
< |
} |
514 |
< |
|
515 |
< |
return -1; |
516 |
< |
} |
517 |
< |
|
518 |
< |
template<typename T> void ZConstraint<T>::integrate(){ |
519 |
< |
// creat zconsWriter |
520 |
< |
fzOut = new ZConsWriter(zconsOutput.c_str(), parameters); |
521 |
< |
|
522 |
< |
if (!fzOut){ |
523 |
< |
sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n"); |
524 |
< |
painCave.isFatal = 1; |
525 |
< |
simError(); |
526 |
< |
} |
527 |
< |
|
528 |
< |
//zero out the velocities of center of mass of unconstrained molecules |
529 |
< |
//and the velocities of center of mass of every single z-constrained molecueles |
530 |
< |
zeroOutVel(); |
531 |
< |
|
532 |
< |
curZconsTime = zconsTime + info->getTime(); |
533 |
< |
|
534 |
< |
T::integrate(); |
535 |
< |
} |
536 |
< |
|
537 |
< |
template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){ |
538 |
< |
double zsys; |
539 |
< |
double COM[3]; |
540 |
< |
double force[3]; |
541 |
< |
double zSysCOMVel; |
542 |
< |
|
543 |
< |
T::calcForce(calcPot, calcStress); |
544 |
< |
|
545 |
< |
|
546 |
< |
if (hasZConsGap){ |
547 |
< |
updateZPos(); |
548 |
< |
} |
549 |
< |
|
550 |
< |
if (checkZConsState()){ |
551 |
< |
zeroOutVel(); |
552 |
< |
forcePolicy->update(); |
553 |
< |
} |
554 |
< |
|
555 |
< |
zsys = calcZSys(); |
556 |
< |
zSysCOMVel = calcSysCOMVel(); |
557 |
< |
#ifdef IS_MPI |
558 |
< |
if (worldRank == 0){ |
559 |
< |
#endif |
560 |
< |
|
561 |
< |
#ifdef IS_MPI |
562 |
< |
} |
563 |
< |
#endif |
564 |
< |
|
565 |
< |
//do zconstraint force; |
566 |
< |
if (haveFixedZMols()){ |
567 |
< |
this->doZconstraintForce(); |
568 |
< |
} |
569 |
< |
|
570 |
< |
//use external force to move the molecules to the specified positions |
571 |
< |
if (haveMovingZMols()){ |
572 |
< |
if (usingSMD) |
573 |
< |
this->doHarmonic(cantPos); |
574 |
< |
else |
575 |
< |
this->doHarmonic(zPos); |
576 |
< |
} |
577 |
< |
|
578 |
< |
//write out forces and current positions of z-constraint molecules |
579 |
< |
if (info->getTime() >= curZconsTime){ |
580 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
581 |
< |
zconsMols[i]->getCOM(COM); |
582 |
< |
curZPos[i] = COM[whichDirection]; |
583 |
< |
|
584 |
< |
//if the z-constraint molecule is still moving, just record its force |
585 |
< |
if (states[i] == zcsMoving){ |
586 |
< |
fz[i] = 0; |
587 |
< |
Atom** movingZAtoms; |
588 |
< |
movingZAtoms = zconsMols[i]->getMyAtoms(); |
589 |
< |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
590 |
< |
movingZAtoms[j]->getFrc(force); |
591 |
< |
fz[i] += force[whichDirection]; |
592 |
< |
} |
593 |
< |
} |
594 |
< |
} |
595 |
< |
fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0], |
596 |
< |
&curZPos[0], &zPos[0]); |
597 |
< |
curZconsTime += zconsTime; |
598 |
< |
} |
599 |
< |
|
600 |
< |
zSysCOMVel = calcSysCOMVel(); |
601 |
< |
#ifdef IS_MPI |
602 |
< |
if (worldRank == 0){ |
603 |
< |
#endif |
604 |
< |
#ifdef IS_MPI |
605 |
< |
} |
606 |
< |
#endif |
607 |
< |
} |
608 |
< |
|
609 |
< |
|
610 |
< |
template<typename T> double ZConstraint<T>::calcZSys(){ |
611 |
< |
//calculate reference z coordinate for z-constraint molecules |
612 |
< |
double totalMass_local; |
613 |
< |
double totalMass; |
614 |
< |
double totalMZ_local; |
615 |
< |
double totalMZ; |
616 |
< |
double massOfCurMol; |
617 |
< |
double COM[3]; |
618 |
< |
|
619 |
< |
totalMass_local = 0; |
620 |
< |
totalMZ_local = 0; |
621 |
< |
|
622 |
< |
for (int i = 0; i < nMols; i++){ |
623 |
< |
massOfCurMol = molecules[i].getTotalMass(); |
624 |
< |
molecules[i].getCOM(COM); |
625 |
< |
|
626 |
< |
totalMass_local += massOfCurMol; |
627 |
< |
totalMZ_local += massOfCurMol * COM[whichDirection]; |
628 |
< |
} |
629 |
< |
|
630 |
< |
|
631 |
< |
#ifdef IS_MPI |
632 |
< |
MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM, |
633 |
< |
MPI_COMM_WORLD); |
634 |
< |
MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
635 |
< |
#else |
636 |
< |
totalMass = totalMass_local; |
637 |
< |
totalMZ = totalMZ_local; |
638 |
< |
#endif |
639 |
< |
|
640 |
< |
double zsys; |
641 |
< |
zsys = totalMZ / totalMass; |
642 |
< |
|
643 |
< |
return zsys; |
644 |
< |
} |
645 |
< |
|
646 |
< |
template<typename T> void ZConstraint<T>::thermalize(void){ |
647 |
< |
T::thermalize(); |
648 |
< |
zeroOutVel(); |
649 |
< |
} |
650 |
< |
|
651 |
< |
template<typename T> void ZConstraint<T>::zeroOutVel(){ |
652 |
< |
Atom** fixedZAtoms; |
653 |
< |
double COMvel[3]; |
654 |
< |
double vel[3]; |
655 |
< |
double zSysCOMVel; |
656 |
< |
|
657 |
< |
//zero out the velocities of center of mass of fixed z-constrained molecules |
658 |
< |
|
659 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
660 |
< |
if (states[i] == zcsFixed){ |
661 |
< |
zconsMols[i]->getCOMvel(COMvel); |
662 |
< |
//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
663 |
< |
|
664 |
< |
fixedZAtoms = zconsMols[i]->getMyAtoms(); |
665 |
< |
|
666 |
< |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
667 |
< |
vel = fixedZAtoms[j]->getVel(); |
668 |
< |
vel[whichDirection] -= COMvel[whichDirection]; |
669 |
< |
fixedZAtoms[j]->setVel(vel); |
670 |
< |
} |
671 |
< |
|
672 |
< |
zconsMols[i]->getCOMvel(COMvel); |
673 |
< |
} |
674 |
< |
} |
675 |
< |
|
676 |
< |
zSysCOMVel = calcSysCOMVel(); |
677 |
< |
#ifdef IS_MPI |
678 |
< |
if (worldRank == 0){ |
679 |
< |
#endif |
680 |
< |
#ifdef IS_MPI |
681 |
< |
} |
682 |
< |
#endif |
683 |
< |
|
684 |
< |
// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules |
685 |
< |
double MVzOfMovingMols_local; |
686 |
< |
double MVzOfMovingMols; |
687 |
< |
double totalMassOfMovingZMols_local; |
688 |
< |
double totalMassOfMovingZMols; |
689 |
< |
|
690 |
< |
MVzOfMovingMols_local = 0; |
691 |
< |
totalMassOfMovingZMols_local = 0; |
692 |
< |
|
693 |
< |
for (int i = 0; i < (int) (unconsMols.size()); i++){ |
694 |
< |
unconsMols[i]->getCOMvel(COMvel); |
695 |
< |
MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
696 |
< |
} |
697 |
< |
|
698 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
699 |
< |
if (states[i] == zcsMoving){ |
700 |
< |
zconsMols[i]->getCOMvel(COMvel); |
701 |
< |
MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; |
702 |
< |
totalMassOfMovingZMols_local += massOfZConsMols[i]; |
703 |
< |
} |
704 |
< |
} |
705 |
< |
|
706 |
< |
#ifndef IS_MPI |
707 |
< |
MVzOfMovingMols = MVzOfMovingMols_local; |
708 |
< |
totalMassOfMovingZMols = totalMassOfMovingZMols_local; |
709 |
< |
#else |
710 |
< |
MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE, |
711 |
< |
MPI_SUM, MPI_COMM_WORLD); |
712 |
< |
MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, |
713 |
< |
MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
714 |
< |
#endif |
715 |
< |
|
716 |
< |
double vzOfMovingMols; |
717 |
< |
vzOfMovingMols = MVzOfMovingMols / |
718 |
< |
(totalMassOfUncons + totalMassOfMovingZMols); |
719 |
< |
|
720 |
< |
//modify the velocites of unconstrained molecules |
721 |
< |
Atom** unconsAtoms; |
722 |
< |
for (int i = 0; i < (int) (unconsMols.size()); i++){ |
723 |
< |
unconsAtoms = unconsMols[i]->getMyAtoms(); |
724 |
< |
for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
725 |
< |
vel = unconsAtoms[j]->getVel(); |
726 |
< |
vel[whichDirection] -= vzOfMovingMols; |
727 |
< |
unconsAtoms[j]->setVel(vel); |
728 |
< |
} |
729 |
< |
} |
730 |
< |
|
731 |
< |
//modify the velocities of moving z-constrained molecuels |
732 |
< |
Atom** movingZAtoms; |
733 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
734 |
< |
if (states[i] == zcsMoving){ |
735 |
< |
movingZAtoms = zconsMols[i]->getMyAtoms(); |
736 |
< |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
737 |
< |
vel = movingZAtoms[j]->getVel(); |
738 |
< |
vel[whichDirection] -= vzOfMovingMols; |
739 |
< |
movingZAtoms[j]->setVel(vel); |
740 |
< |
} |
741 |
< |
} |
742 |
< |
} |
743 |
< |
|
744 |
< |
|
745 |
< |
zSysCOMVel = calcSysCOMVel(); |
746 |
< |
#ifdef IS_MPI |
747 |
< |
if (worldRank == 0){ |
748 |
< |
#endif |
749 |
< |
#ifdef IS_MPI |
750 |
< |
} |
751 |
< |
#endif |
752 |
< |
} |
753 |
< |
|
754 |
< |
|
755 |
< |
template<typename T> void ZConstraint<T>::doZconstraintForce(){ |
756 |
< |
Atom** zconsAtoms; |
757 |
< |
double totalFZ; |
758 |
< |
double totalFZ_local; |
759 |
< |
double COM[3]; |
760 |
< |
double force[3]; |
761 |
< |
|
762 |
< |
//constrain the molecules which do not reach the specified positions |
763 |
< |
|
764 |
< |
//Zero Out the force of z-contrained molecules |
765 |
< |
totalFZ_local = 0; |
766 |
< |
|
767 |
< |
//calculate the total z-contrained force of fixed z-contrained molecules |
768 |
< |
|
769 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
770 |
< |
if (states[i] == zcsFixed){ |
771 |
< |
zconsMols[i]->getCOM(COM); |
772 |
< |
zconsAtoms = zconsMols[i]->getMyAtoms(); |
773 |
< |
|
774 |
< |
fz[i] = 0; |
775 |
< |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
776 |
< |
zconsAtoms[j]->getFrc(force); |
777 |
< |
fz[i] += force[whichDirection]; |
778 |
< |
} |
779 |
< |
totalFZ_local += fz[i]; |
780 |
< |
|
781 |
< |
} |
782 |
< |
} |
783 |
< |
|
784 |
< |
//calculate total z-constraint force |
785 |
< |
#ifdef IS_MPI |
786 |
< |
MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
787 |
< |
#else |
788 |
< |
totalFZ = totalFZ_local; |
789 |
< |
#endif |
790 |
< |
|
791 |
< |
|
792 |
< |
// apply negative to fixed z-constrained molecues; |
793 |
< |
force[0] = 0; |
794 |
< |
force[1] = 0; |
795 |
< |
force[2] = 0; |
796 |
< |
|
797 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
798 |
< |
if (states[i] == zcsFixed){ |
799 |
< |
int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms(); |
800 |
< |
zconsAtoms = zconsMols[i]->getMyAtoms(); |
801 |
< |
|
802 |
< |
for (int j = 0; j < nAtomOfCurZConsMol; j++){ |
803 |
< |
//force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; |
804 |
< |
force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i], |
805 |
< |
zconsAtoms[j], |
806 |
< |
fz[i]); |
807 |
< |
zconsAtoms[j]->addFrc(force); |
808 |
< |
} |
809 |
< |
} |
810 |
< |
} |
811 |
< |
|
812 |
< |
force[0] = 0; |
813 |
< |
force[1] = 0; |
814 |
< |
force[2] = 0; |
815 |
< |
|
816 |
< |
//modify the forces of unconstrained molecules |
817 |
< |
for (int i = 0; i < (int) (unconsMols.size()); i++){ |
818 |
< |
Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
819 |
< |
|
820 |
< |
for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
821 |
< |
//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); |
822 |
< |
force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j], |
823 |
< |
totalFZ); |
824 |
< |
unconsAtoms[j]->addFrc(force); |
825 |
< |
} |
826 |
< |
} |
827 |
< |
|
828 |
< |
//modify the forces of moving z-constrained molecules |
829 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
830 |
< |
if (states[i] == zcsMoving){ |
831 |
< |
Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
832 |
< |
|
833 |
< |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
834 |
< |
//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); |
835 |
< |
force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j], |
836 |
< |
totalFZ); |
837 |
< |
movingZAtoms[j]->addFrc(force); |
838 |
< |
} |
839 |
< |
} |
840 |
< |
} |
841 |
< |
} |
842 |
< |
|
843 |
< |
|
844 |
< |
template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){ |
845 |
< |
double force[3]; |
846 |
< |
double harmonicU; |
847 |
< |
double harmonicF; |
848 |
< |
double COM[3]; |
849 |
< |
double diff; |
850 |
< |
double totalFZ_local; |
851 |
< |
double totalFZ; |
852 |
< |
|
853 |
< |
force[0] = 0; |
854 |
< |
force[1] = 0; |
855 |
< |
force[2] = 0; |
856 |
< |
|
857 |
< |
totalFZ_local = 0; |
858 |
< |
|
859 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
860 |
< |
if (states[i] == zcsMoving){ |
861 |
< |
zconsMols[i]->getCOM(COM); |
862 |
< |
|
863 |
< |
diff = COM[whichDirection] - resPos[i]; |
864 |
< |
|
865 |
< |
harmonicU = 0.5 * kz[i] * diff * diff; |
866 |
< |
info->lrPot += harmonicU; |
867 |
< |
|
868 |
< |
harmonicF = -kz[i] * diff; |
869 |
< |
totalFZ_local += harmonicF; |
870 |
< |
|
871 |
< |
//adjust force |
872 |
< |
|
873 |
< |
Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
874 |
< |
|
875 |
< |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
876 |
< |
force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], |
877 |
< |
movingZAtoms[j], |
878 |
< |
harmonicF); |
879 |
< |
movingZAtoms[j]->addFrc(force); |
880 |
< |
} |
881 |
< |
} |
882 |
< |
} |
883 |
< |
|
884 |
< |
#ifndef IS_MPI |
885 |
< |
totalFZ = totalFZ_local; |
886 |
< |
#else |
887 |
< |
MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
888 |
< |
#endif |
889 |
< |
|
890 |
< |
force[0] = 0; |
891 |
< |
force[1] = 0; |
892 |
< |
force[2] = 0; |
893 |
< |
|
894 |
< |
//modify the forces of unconstrained molecules |
895 |
< |
for (int i = 0; i < (int) (unconsMols.size()); i++){ |
896 |
< |
Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
897 |
< |
|
898 |
< |
for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
899 |
< |
//force[whichDirection] = - totalFZ /totNumOfUnconsAtoms; |
900 |
< |
force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j], |
901 |
< |
totalFZ); |
902 |
< |
unconsAtoms[j]->addFrc(force); |
903 |
< |
} |
904 |
< |
} |
905 |
< |
|
906 |
< |
} |
907 |
< |
|
908 |
< |
template<typename T> bool ZConstraint<T>::checkZConsState(){ |
909 |
< |
double COM[3]; |
910 |
< |
double diff; |
911 |
< |
|
912 |
< |
int changed_local; |
913 |
< |
int changed; |
914 |
< |
|
915 |
< |
changed_local = 0; |
916 |
< |
|
917 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
918 |
< |
zconsMols[i]->getCOM(COM); |
919 |
< |
diff = fabs(COM[whichDirection] - zPos[i]); |
920 |
< |
if (diff <= zconsTol && states[i] == zcsMoving){ |
921 |
< |
states[i] = zcsFixed; |
922 |
< |
changed_local = 1; |
923 |
< |
|
924 |
< |
if(usingSMD) |
925 |
< |
prevCantPos = cantPos; |
926 |
< |
|
927 |
< |
if (hasZConsGap) |
928 |
< |
endFixTime[i] = info->getTime() + zconsFixTime; |
929 |
< |
} |
930 |
< |
else if (diff > zconsTol && states[i] == zcsFixed){ |
931 |
< |
states[i] = zcsMoving; |
932 |
< |
changed_local = 1; |
933 |
< |
|
934 |
< |
if(usingSMD) |
935 |
< |
cantPos = prevCantPos; |
936 |
< |
|
937 |
< |
if (hasZConsGap) |
938 |
< |
endFixTime[i] = INFINITE_TIME; |
939 |
< |
} |
940 |
< |
} |
941 |
< |
|
942 |
< |
#ifndef IS_MPI |
943 |
< |
changed = changed_local; |
944 |
< |
#else |
945 |
< |
MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
946 |
< |
#endif |
947 |
< |
|
948 |
< |
return (changed > 0); |
949 |
< |
} |
950 |
< |
|
951 |
< |
template<typename T> bool ZConstraint<T>::haveFixedZMols(){ |
952 |
< |
int havingFixed_local; |
953 |
< |
int havingFixed; |
954 |
< |
|
955 |
< |
havingFixed_local = 0; |
956 |
< |
|
957 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++) |
958 |
< |
if (states[i] == zcsFixed){ |
959 |
< |
havingFixed_local = 1; |
960 |
< |
break; |
961 |
< |
} |
962 |
< |
|
963 |
< |
#ifndef IS_MPI |
964 |
< |
havingFixed = havingFixed_local; |
965 |
< |
#else |
966 |
< |
MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM, |
967 |
< |
MPI_COMM_WORLD); |
968 |
< |
#endif |
969 |
< |
|
970 |
< |
return (havingFixed > 0); |
971 |
< |
} |
972 |
< |
|
973 |
< |
|
974 |
< |
template<typename T> bool ZConstraint<T>::haveMovingZMols(){ |
975 |
< |
int havingMoving_local; |
976 |
< |
int havingMoving; |
977 |
< |
|
978 |
< |
havingMoving_local = 0; |
979 |
< |
|
980 |
< |
for (int i = 0; i < (int) (zconsMols.size()); i++) |
981 |
< |
if (states[i] == zcsMoving){ |
982 |
< |
havingMoving_local = 1; |
983 |
< |
break; |
984 |
< |
} |
985 |
< |
|
986 |
< |
#ifndef IS_MPI |
987 |
< |
havingMoving = havingMoving_local; |
988 |
< |
#else |
989 |
< |
MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM, |
990 |
< |
MPI_COMM_WORLD); |
991 |
< |
#endif |
992 |
< |
|
993 |
< |
return (havingMoving > 0); |
994 |
< |
} |
995 |
< |
|
996 |
< |
|
997 |
< |
template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){ |
998 |
< |
double MVzOfMovingMols_local; |
999 |
< |
double MVzOfMovingMols; |
1000 |
< |
double totalMassOfMovingZMols_local; |
1001 |
< |
double totalMassOfMovingZMols; |
1002 |
< |
double COMvel[3]; |
1003 |
< |
|
1004 |
< |
MVzOfMovingMols_local = 0; |
1005 |
< |
totalMassOfMovingZMols_local = 0; |
1006 |
< |
|
1007 |
< |
for (int i = 0; i < unconsMols.size(); i++){ |
1008 |
< |
unconsMols[i]->getCOMvel(COMvel); |
1009 |
< |
MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
1010 |
< |
} |
1011 |
< |
|
1012 |
< |
for (int i = 0; i < zconsMols.size(); i++){ |
1013 |
< |
if (states[i] == zcsMoving){ |
1014 |
< |
zconsMols[i]->getCOMvel(COMvel); |
1015 |
< |
MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; |
1016 |
< |
totalMassOfMovingZMols_local += massOfZConsMols[i]; |
1017 |
< |
} |
1018 |
< |
} |
1019 |
< |
|
1020 |
< |
#ifndef IS_MPI |
1021 |
< |
MVzOfMovingMols = MVzOfMovingMols_local; |
1022 |
< |
totalMassOfMovingZMols = totalMassOfMovingZMols_local; |
1023 |
< |
#else |
1024 |
< |
MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE, |
1025 |
< |
MPI_SUM, MPI_COMM_WORLD); |
1026 |
< |
MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, |
1027 |
< |
MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
1028 |
< |
#endif |
1029 |
< |
|
1030 |
< |
double vzOfMovingMols; |
1031 |
< |
vzOfMovingMols = MVzOfMovingMols / |
1032 |
< |
(totalMassOfUncons + totalMassOfMovingZMols); |
1033 |
< |
|
1034 |
< |
return vzOfMovingMols; |
1035 |
< |
} |
1036 |
< |
|
1037 |
< |
template<typename T> double ZConstraint<T>::calcSysCOMVel(){ |
1038 |
< |
double COMvel[3]; |
1039 |
< |
double tempMVz_local; |
1040 |
< |
double tempMVz; |
1041 |
< |
double massOfZCons_local; |
1042 |
< |
double massOfZCons; |
1043 |
< |
|
1044 |
< |
|
1045 |
< |
tempMVz_local = 0; |
1046 |
< |
|
1047 |
< |
for (int i = 0 ; i < nMols; i++){ |
1048 |
< |
molecules[i].getCOMvel(COMvel); |
1049 |
< |
tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection]; |
1050 |
< |
} |
1051 |
< |
|
1052 |
< |
massOfZCons_local = 0; |
1053 |
< |
|
1054 |
< |
for (int i = 0; i < (int) (massOfZConsMols.size()); i++){ |
1055 |
< |
massOfZCons_local += massOfZConsMols[i]; |
1056 |
< |
} |
1057 |
< |
#ifndef IS_MPI |
1058 |
< |
massOfZCons = massOfZCons_local; |
1059 |
< |
tempMVz = tempMVz_local; |
1060 |
< |
#else |
1061 |
< |
MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM, |
1062 |
< |
MPI_COMM_WORLD); |
1063 |
< |
MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
1064 |
< |
#endif |
1065 |
< |
|
1066 |
< |
return tempMVz / (totalMassOfUncons + massOfZCons); |
1067 |
< |
} |
1068 |
< |
|
1069 |
< |
template<typename T> double ZConstraint<T>::calcTotalForce(){ |
1070 |
< |
double force[3]; |
1071 |
< |
double totalForce_local; |
1072 |
< |
double totalForce; |
1073 |
< |
|
1074 |
< |
totalForce_local = 0; |
1075 |
< |
|
1076 |
< |
for (int i = 0; i < nAtoms; i++){ |
1077 |
< |
atoms[i]->getFrc(force); |
1078 |
< |
totalForce_local += force[whichDirection]; |
1079 |
< |
} |
1080 |
< |
|
1081 |
< |
#ifndef IS_MPI |
1082 |
< |
totalForce = totalForce_local; |
1083 |
< |
#else |
1084 |
< |
MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM, |
1085 |
< |
MPI_COMM_WORLD); |
1086 |
< |
#endif |
1087 |
< |
|
1088 |
< |
return totalForce; |
1089 |
< |
} |
1090 |
< |
|
1091 |
< |
template<typename T> void ZConstraint<T>::PolicyByNumber::update(){ |
1092 |
< |
//calculate the number of atoms of moving z-constrained molecules |
1093 |
< |
int nMovingZAtoms_local; |
1094 |
< |
int nMovingZAtoms; |
1095 |
< |
|
1096 |
< |
nMovingZAtoms_local = 0; |
1097 |
< |
for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++) |
1098 |
< |
if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){ |
1099 |
< |
nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms(); |
1100 |
< |
} |
1101 |
< |
|
1102 |
< |
#ifdef IS_MPI |
1103 |
< |
MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, |
1104 |
< |
MPI_COMM_WORLD); |
1105 |
< |
#else |
1106 |
< |
nMovingZAtoms = nMovingZAtoms_local; |
1107 |
< |
#endif |
1108 |
< |
totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms; |
1109 |
< |
} |
1110 |
< |
|
1111 |
< |
template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, |
1112 |
< |
Atom* atom, |
1113 |
< |
double totalForce){ |
1114 |
< |
return totalForce / mol->getNAtoms(); |
1115 |
< |
} |
1116 |
< |
|
1117 |
< |
template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, |
1118 |
< |
double totalForce){ |
1119 |
< |
return totalForce / totNumOfMovingAtoms; |
1120 |
< |
} |
1121 |
< |
|
1122 |
< |
template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, |
1123 |
< |
Atom* atom, |
1124 |
< |
double totalForce){ |
1125 |
< |
return totalForce / mol->getNAtoms(); |
1126 |
< |
} |
1127 |
< |
|
1128 |
< |
template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, |
1129 |
< |
double totalForce){ |
1130 |
< |
return totalForce / zconsIntegrator->totNumOfUnconsAtoms; |
1131 |
< |
} |
1132 |
< |
|
1133 |
< |
|
1134 |
< |
template<typename T> void ZConstraint<T>::PolicyByMass::update(){ |
1135 |
< |
//calculate the number of atoms of moving z-constrained molecules |
1136 |
< |
double massOfMovingZAtoms_local; |
1137 |
< |
double massOfMovingZAtoms; |
1138 |
< |
|
1139 |
< |
massOfMovingZAtoms_local = 0; |
1140 |
< |
for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++) |
1141 |
< |
if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){ |
1142 |
< |
massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass(); |
1143 |
< |
} |
1144 |
< |
|
1145 |
< |
#ifdef IS_MPI |
1146 |
< |
MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE, |
1147 |
< |
MPI_SUM, MPI_COMM_WORLD); |
1148 |
< |
#else |
1149 |
< |
massOfMovingZAtoms = massOfMovingZAtoms_local; |
1150 |
< |
#endif |
1151 |
< |
totMassOfMovingAtoms = massOfMovingZAtoms + |
1152 |
< |
zconsIntegrator->totalMassOfUncons; |
1153 |
< |
} |
1154 |
< |
|
1155 |
< |
template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, |
1156 |
< |
Atom* atom, |
1157 |
< |
double totalForce){ |
1158 |
< |
return totalForce * atom->getMass() / mol->getTotalMass(); |
1159 |
< |
} |
1160 |
< |
|
1161 |
< |
template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom, |
1162 |
< |
double totalForce){ |
1163 |
< |
return totalForce * atom->getMass() / totMassOfMovingAtoms; |
1164 |
< |
} |
1165 |
< |
|
1166 |
< |
template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, |
1167 |
< |
Atom* atom, |
1168 |
< |
double totalForce){ |
1169 |
< |
return totalForce * atom->getMass() / mol->getTotalMass(); |
1170 |
< |
} |
1171 |
< |
|
1172 |
< |
template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, |
1173 |
< |
double totalForce){ |
1174 |
< |
return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons; |
1175 |
< |
} |
1176 |
< |
|
1177 |
< |
template<typename T> void ZConstraint<T>::updateZPos(){ |
1178 |
< |
double curTime; |
1179 |
< |
double COM[3]; |
1180 |
< |
|
1181 |
< |
curTime = info->getTime(); |
1182 |
< |
|
1183 |
< |
for (size_t i = 0; i < zconsMols.size(); i++){ |
1184 |
< |
|
1185 |
< |
if (states[i] == zcsFixed && curTime >= endFixTime[i]){ |
1186 |
< |
zPos[i] += zconsGap; |
1187 |
< |
|
1188 |
< |
if (usingSMD){ |
1189 |
< |
zconsMols[i]->getCOM(COM); |
1190 |
< |
cantPos[i] = COM[whichDirection]; |
1191 |
< |
} |
1192 |
< |
|
1193 |
< |
} |
1194 |
< |
|
1195 |
< |
} |
1196 |
< |
|
1197 |
< |
} |
1198 |
< |
|
1199 |
< |
template<typename T> void ZConstraint<T>::updateCantPos(){ |
1200 |
< |
double curTime; |
1201 |
< |
double dt; |
1202 |
< |
|
1203 |
< |
curTime = info->getTime(); |
1204 |
< |
dt = info->dt; |
1205 |
< |
|
1206 |
< |
for (size_t i = 0; i < zconsMols.size(); i++){ |
1207 |
< |
if (states[i] == zcsMoving){ |
1208 |
< |
cantPos[i] += cantVel[i] * dt; |
1209 |
< |
} |
1210 |
< |
} |
1211 |
< |
|
1212 |
< |
} |
1 |
> |
#include "integrators/Integrator.hpp" |
2 |
> |
|
3 |
> |
#include "utils/simError.h" |
4 |
> |
|
5 |
> |
#include <math.h> |
6 |
> |
|
7 |
> |
|
8 |
> |
|
9 |
> |
const double INFINITE_TIME = 10e30; |
10 |
> |
|
11 |
> |
template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, |
12 |
> |
|
13 |
> |
ForceFields* the_ff): T(theInfo, the_ff), |
14 |
> |
|
15 |
> |
fzOut(NULL), |
16 |
> |
|
17 |
> |
curZconsTime(0), |
18 |
> |
|
19 |
> |
forcePolicy(NULL), |
20 |
> |
|
21 |
> |
usingSMD(false), |
22 |
> |
|
23 |
> |
hasZConsGap(false){ |
24 |
> |
|
25 |
> |
//get properties from SimInfo |
26 |
> |
|
27 |
> |
GenericData* data; |
28 |
> |
|
29 |
> |
ZConsParaData* zConsParaData; |
30 |
> |
|
31 |
> |
DoubleGenericData* sampleTime; |
32 |
> |
|
33 |
> |
DoubleGenericData* tolerance; |
34 |
> |
|
35 |
> |
DoubleGenericData* gap; |
36 |
> |
|
37 |
> |
DoubleGenericData* fixtime; |
38 |
> |
|
39 |
> |
StringGenericData* policy; |
40 |
> |
|
41 |
> |
StringGenericData* filename; |
42 |
> |
|
43 |
> |
IntGenericData* smdFlag; |
44 |
> |
|
45 |
> |
double COM[3]; |
46 |
> |
|
47 |
> |
|
48 |
> |
|
49 |
> |
//by default, the direction of constraint is z |
50 |
> |
|
51 |
> |
// 0 --> x |
52 |
> |
|
53 |
> |
// 1 --> y |
54 |
> |
|
55 |
> |
// 2 --> z |
56 |
> |
|
57 |
> |
whichDirection = 2; |
58 |
> |
|
59 |
> |
|
60 |
> |
|
61 |
> |
//estimate the force constant of harmonical potential |
62 |
> |
|
63 |
> |
double Kb = 1.986E-3 ; //in kcal/K |
64 |
> |
|
65 |
> |
|
66 |
> |
|
67 |
> |
double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) / |
68 |
> |
|
69 |
> |
2; |
70 |
> |
|
71 |
> |
zForceConst = Kb * info->target_temp / (halfOfLargestBox * halfOfLargestBox); |
72 |
> |
|
73 |
> |
|
74 |
> |
|
75 |
> |
//creat force Subtraction policy |
76 |
> |
|
77 |
> |
data = info->getPropertyByName(ZCONSFORCEPOLICY_ID); |
78 |
> |
|
79 |
> |
if (!data){ |
80 |
> |
|
81 |
> |
sprintf(painCave.errMsg, |
82 |
> |
|
83 |
> |
"ZConstraint Warning: User does not set force Subtraction policy, " |
84 |
> |
|
85 |
> |
"PolicyByMass is used\n"); |
86 |
> |
|
87 |
> |
painCave.isFatal = 0; |
88 |
> |
|
89 |
> |
simError(); |
90 |
> |
|
91 |
> |
|
92 |
> |
|
93 |
> |
forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
94 |
> |
|
95 |
> |
} |
96 |
> |
|
97 |
> |
else{ |
98 |
> |
|
99 |
> |
policy = dynamic_cast<StringGenericData*>(data); |
100 |
> |
|
101 |
> |
|
102 |
> |
|
103 |
> |
if (!policy){ |
104 |
> |
|
105 |
> |
sprintf(painCave.errMsg, |
106 |
> |
|
107 |
> |
"ZConstraint Error: Convertion from GenericData to StringGenericData failure, " |
108 |
> |
|
109 |
> |
"PolicyByMass is used\n"); |
110 |
> |
|
111 |
> |
painCave.isFatal = 0; |
112 |
> |
|
113 |
> |
simError(); |
114 |
> |
|
115 |
> |
|
116 |
> |
|
117 |
> |
forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
118 |
> |
|
119 |
> |
} |
120 |
> |
|
121 |
> |
else{ |
122 |
> |
|
123 |
> |
if (policy->getData() == "BYNUMBER") |
124 |
> |
|
125 |
> |
forcePolicy = (ForceSubtractionPolicy *) new PolicyByNumber(this); |
126 |
> |
|
127 |
> |
else if (policy->getData() == "BYMASS") |
128 |
> |
|
129 |
> |
forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
130 |
> |
|
131 |
> |
else{ |
132 |
> |
|
133 |
> |
sprintf(painCave.errMsg, |
134 |
> |
|
135 |
> |
"ZConstraint Warning: unknown force Subtraction policy, " |
136 |
> |
|
137 |
> |
"PolicyByMass is used\n"); |
138 |
> |
|
139 |
> |
painCave.isFatal = 0; |
140 |
> |
|
141 |
> |
simError(); |
142 |
> |
|
143 |
> |
forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); |
144 |
> |
|
145 |
> |
} |
146 |
> |
|
147 |
> |
} |
148 |
> |
|
149 |
> |
} |
150 |
> |
|
151 |
> |
|
152 |
> |
|
153 |
> |
|
154 |
> |
|
155 |
> |
//retrieve sample time of z-contraint |
156 |
> |
|
157 |
> |
data = info->getPropertyByName(ZCONSTIME_ID); |
158 |
> |
|
159 |
> |
|
160 |
> |
|
161 |
> |
if (!data){ |
162 |
> |
|
163 |
> |
sprintf(painCave.errMsg, |
164 |
> |
|
165 |
> |
"ZConstraint error: If you use an ZConstraint\n" |
166 |
> |
|
167 |
> |
" , you must set sample time.\n"); |
168 |
> |
|
169 |
> |
painCave.isFatal = 1; |
170 |
> |
|
171 |
> |
simError(); |
172 |
> |
|
173 |
> |
} |
174 |
> |
|
175 |
> |
else{ |
176 |
> |
|
177 |
> |
sampleTime = dynamic_cast<DoubleGenericData*>(data); |
178 |
> |
|
179 |
> |
|
180 |
> |
|
181 |
> |
if (!sampleTime){ |
182 |
> |
|
183 |
> |
sprintf(painCave.errMsg, |
184 |
> |
|
185 |
> |
"ZConstraint error: Can not get property from SimInfo\n"); |
186 |
> |
|
187 |
> |
painCave.isFatal = 1; |
188 |
> |
|
189 |
> |
simError(); |
190 |
> |
|
191 |
> |
} |
192 |
> |
|
193 |
> |
else{ |
194 |
> |
|
195 |
> |
this->zconsTime = sampleTime->getData(); |
196 |
> |
|
197 |
> |
} |
198 |
> |
|
199 |
> |
} |
200 |
> |
|
201 |
> |
|
202 |
> |
|
203 |
> |
//retrieve output filename of z force |
204 |
> |
|
205 |
> |
data = info->getPropertyByName(ZCONSFILENAME_ID); |
206 |
> |
|
207 |
> |
if (!data){ |
208 |
> |
|
209 |
> |
sprintf(painCave.errMsg, |
210 |
> |
|
211 |
> |
"ZConstraint error: If you use an ZConstraint\n" |
212 |
> |
|
213 |
> |
" , you must set output filename of z-force.\n"); |
214 |
> |
|
215 |
> |
painCave.isFatal = 1; |
216 |
> |
|
217 |
> |
simError(); |
218 |
> |
|
219 |
> |
} |
220 |
> |
|
221 |
> |
else{ |
222 |
> |
|
223 |
> |
filename = dynamic_cast<StringGenericData*>(data); |
224 |
> |
|
225 |
> |
|
226 |
> |
|
227 |
> |
if (!filename){ |
228 |
> |
|
229 |
> |
sprintf(painCave.errMsg, |
230 |
> |
|
231 |
> |
"ZConstraint error: Can not get property from SimInfo\n"); |
232 |
> |
|
233 |
> |
painCave.isFatal = 1; |
234 |
> |
|
235 |
> |
simError(); |
236 |
> |
|
237 |
> |
} |
238 |
> |
|
239 |
> |
else{ |
240 |
> |
|
241 |
> |
this->zconsOutput = filename->getData(); |
242 |
> |
|
243 |
> |
} |
244 |
> |
|
245 |
> |
} |
246 |
> |
|
247 |
> |
|
248 |
> |
|
249 |
> |
//retrieve tolerance for z-constraint molecuels |
250 |
> |
|
251 |
> |
data = info->getPropertyByName(ZCONSTOL_ID); |
252 |
> |
|
253 |
> |
|
254 |
> |
|
255 |
> |
if (!data){ |
256 |
> |
|
257 |
> |
sprintf(painCave.errMsg, "ZConstraint error: can not get tolerance \n"); |
258 |
> |
|
259 |
> |
painCave.isFatal = 1; |
260 |
> |
|
261 |
> |
simError(); |
262 |
> |
|
263 |
> |
} |
264 |
> |
|
265 |
> |
else{ |
266 |
> |
|
267 |
> |
tolerance = dynamic_cast<DoubleGenericData*>(data); |
268 |
> |
|
269 |
> |
|
270 |
> |
|
271 |
> |
if (!tolerance){ |
272 |
> |
|
273 |
> |
sprintf(painCave.errMsg, |
274 |
> |
|
275 |
> |
"ZConstraint error: Can not get property from SimInfo\n"); |
276 |
> |
|
277 |
> |
painCave.isFatal = 1; |
278 |
> |
|
279 |
> |
simError(); |
280 |
> |
|
281 |
> |
} |
282 |
> |
|
283 |
> |
else{ |
284 |
> |
|
285 |
> |
this->zconsTol = tolerance->getData(); |
286 |
> |
|
287 |
> |
} |
288 |
> |
|
289 |
> |
} |
290 |
> |
|
291 |
> |
|
292 |
> |
|
293 |
> |
//quick hack here |
294 |
> |
|
295 |
> |
data = info->getPropertyByName(ZCONSGAP_ID); |
296 |
> |
|
297 |
> |
|
298 |
> |
|
299 |
> |
if (data){ |
300 |
> |
|
301 |
> |
gap = dynamic_cast<DoubleGenericData*>(data); |
302 |
> |
|
303 |
> |
|
304 |
> |
|
305 |
> |
if (!gap){ |
306 |
> |
|
307 |
> |
sprintf(painCave.errMsg, |
308 |
> |
|
309 |
> |
"ZConstraint error: Can not get property from SimInfo\n"); |
310 |
> |
|
311 |
> |
painCave.isFatal = 1; |
312 |
> |
|
313 |
> |
simError(); |
314 |
> |
|
315 |
> |
} |
316 |
> |
|
317 |
> |
else{ |
318 |
> |
|
319 |
> |
this->hasZConsGap = true; |
320 |
> |
|
321 |
> |
this->zconsGap = gap->getData(); |
322 |
> |
|
323 |
> |
} |
324 |
> |
|
325 |
> |
} |
326 |
> |
|
327 |
> |
|
328 |
> |
|
329 |
> |
|
330 |
> |
|
331 |
> |
|
332 |
> |
|
333 |
> |
data = info->getPropertyByName(ZCONSFIXTIME_ID); |
334 |
> |
|
335 |
> |
|
336 |
> |
|
337 |
> |
if (data){ |
338 |
> |
|
339 |
> |
fixtime = dynamic_cast<DoubleGenericData*>(data); |
340 |
> |
|
341 |
> |
if (!fixtime){ |
342 |
> |
|
343 |
> |
sprintf(painCave.errMsg, |
344 |
> |
|
345 |
> |
"ZConstraint error: Can not get zconsFixTime from SimInfo\n"); |
346 |
> |
|
347 |
> |
painCave.isFatal = 1; |
348 |
> |
|
349 |
> |
simError(); |
350 |
> |
|
351 |
> |
} |
352 |
> |
|
353 |
> |
else{ |
354 |
> |
|
355 |
> |
this->zconsFixTime = fixtime->getData(); |
356 |
> |
|
357 |
> |
} |
358 |
> |
|
359 |
> |
} |
360 |
> |
|
361 |
> |
else if(hasZConsGap){ |
362 |
> |
|
363 |
> |
sprintf(painCave.errMsg, |
364 |
> |
|
365 |
> |
"ZConstraint error: must set fixtime if already set zconsGap\n"); |
366 |
> |
|
367 |
> |
painCave.isFatal = 1; |
368 |
> |
|
369 |
> |
simError(); |
370 |
> |
|
371 |
> |
} |
372 |
> |
|
373 |
> |
|
374 |
> |
|
375 |
> |
|
376 |
> |
|
377 |
> |
|
378 |
> |
|
379 |
> |
data = info->getPropertyByName(ZCONSUSINGSMD_ID); |
380 |
> |
|
381 |
> |
|
382 |
> |
|
383 |
> |
if (data){ |
384 |
> |
|
385 |
> |
smdFlag = dynamic_cast<IntGenericData*>(data); |
386 |
> |
|
387 |
> |
|
388 |
> |
|
389 |
> |
if (!smdFlag){ |
390 |
> |
|
391 |
> |
sprintf(painCave.errMsg, |
392 |
> |
|
393 |
> |
"ZConstraint error: Can not get property from SimInfo\n"); |
394 |
> |
|
395 |
> |
painCave.isFatal = 1; |
396 |
> |
|
397 |
> |
simError(); |
398 |
> |
|
399 |
> |
} |
400 |
> |
|
401 |
> |
else{ |
402 |
> |
|
403 |
> |
this->usingSMD= smdFlag->getData() ? true : false; |
404 |
> |
|
405 |
> |
} |
406 |
> |
|
407 |
> |
|
408 |
> |
|
409 |
> |
} |
410 |
> |
|
411 |
> |
|
412 |
> |
|
413 |
> |
|
414 |
> |
|
415 |
> |
|
416 |
> |
|
417 |
> |
//retrieve index of z-constraint molecules |
418 |
> |
|
419 |
> |
data = info->getPropertyByName(ZCONSPARADATA_ID); |
420 |
> |
|
421 |
> |
if (!data){ |
422 |
> |
|
423 |
> |
sprintf(painCave.errMsg, |
424 |
> |
|
425 |
> |
"ZConstraint error: If you use an ZConstraint\n" |
426 |
> |
|
427 |
> |
" , you must set index of z-constraint molecules.\n"); |
428 |
> |
|
429 |
> |
painCave.isFatal = 1; |
430 |
> |
|
431 |
> |
simError(); |
432 |
> |
|
433 |
> |
} |
434 |
> |
|
435 |
> |
else{ |
436 |
> |
|
437 |
> |
zConsParaData = dynamic_cast<ZConsParaData*>(data); |
438 |
> |
|
439 |
> |
|
440 |
> |
|
441 |
> |
if (!zConsParaData){ |
442 |
> |
|
443 |
> |
sprintf(painCave.errMsg, |
444 |
> |
|
445 |
> |
"ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n"); |
446 |
> |
|
447 |
> |
painCave.isFatal = 1; |
448 |
> |
|
449 |
> |
simError(); |
450 |
> |
|
451 |
> |
} |
452 |
> |
|
453 |
> |
else{ |
454 |
> |
|
455 |
> |
parameters = zConsParaData->getData(); |
456 |
> |
|
457 |
> |
|
458 |
> |
|
459 |
> |
//check the range of zconsIndex |
460 |
> |
|
461 |
> |
//and the minimum value of index is the first one (we already sorted the data) |
462 |
> |
|
463 |
> |
//the maximum value of index is the last one |
464 |
> |
|
465 |
> |
|
466 |
> |
|
467 |
> |
int maxIndex; |
468 |
> |
|
469 |
> |
int minIndex; |
470 |
> |
|
471 |
> |
int totalNumMol; |
472 |
> |
|
473 |
> |
|
474 |
> |
|
475 |
> |
minIndex = (*parameters)[0].zconsIndex; |
476 |
> |
|
477 |
> |
if (minIndex < 0){ |
478 |
> |
|
479 |
> |
sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); |
480 |
> |
|
481 |
> |
painCave.isFatal = 1; |
482 |
> |
|
483 |
> |
simError(); |
484 |
> |
|
485 |
> |
} |
486 |
> |
|
487 |
> |
|
488 |
> |
|
489 |
> |
maxIndex = (*parameters)[parameters->size() - 1].zconsIndex; |
490 |
> |
|
491 |
> |
|
492 |
> |
|
493 |
> |
#ifndef IS_MPI |
494 |
> |
|
495 |
> |
totalNumMol = nMols; |
496 |
> |
|
497 |
> |
#else |
498 |
> |
|
499 |
> |
totalNumMol = mpiSim->getNMolGlobal(); |
500 |
> |
|
501 |
> |
#endif |
502 |
> |
|
503 |
> |
|
504 |
> |
|
505 |
> |
if (maxIndex > totalNumMol - 1){ |
506 |
> |
|
507 |
> |
sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); |
508 |
> |
|
509 |
> |
painCave.isFatal = 1; |
510 |
> |
|
511 |
> |
simError(); |
512 |
> |
|
513 |
> |
} |
514 |
> |
|
515 |
> |
|
516 |
> |
|
517 |
> |
//if user does not specify the zpos for the zconstraint molecule |
518 |
> |
|
519 |
> |
//its initial z coordinate will be used as default |
520 |
> |
|
521 |
> |
for (int i = 0; i < (int) (parameters->size()); i++){ |
522 |
> |
|
523 |
> |
if (!(*parameters)[i].havingZPos){ |
524 |
> |
|
525 |
> |
#ifndef IS_MPI |
526 |
> |
|
527 |
> |
for (int j = 0; j < nMols; j++){ |
528 |
> |
|
529 |
> |
if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ |
530 |
> |
|
531 |
> |
molecules[j].getCOM(COM); |
532 |
> |
|
533 |
> |
break; |
534 |
> |
|
535 |
> |
} |
536 |
> |
|
537 |
> |
} |
538 |
> |
|
539 |
> |
#else |
540 |
> |
|
541 |
> |
//query which processor current zconstraint molecule belongs to |
542 |
> |
|
543 |
> |
int* MolToProcMap; |
544 |
> |
|
545 |
> |
int whichNode; |
546 |
> |
|
547 |
> |
|
548 |
> |
|
549 |
> |
MolToProcMap = mpiSim->getMolToProcMap(); |
550 |
> |
|
551 |
> |
whichNode = MolToProcMap[(*parameters)[i].zconsIndex]; |
552 |
> |
|
553 |
> |
|
554 |
> |
|
555 |
> |
//broadcast the zpos of current z-contraint molecule |
556 |
> |
|
557 |
> |
//the node which contain this |
558 |
> |
|
559 |
> |
|
560 |
> |
|
561 |
> |
if (worldRank == whichNode){ |
562 |
> |
|
563 |
> |
for (int j = 0; j < nMols; j++) |
564 |
> |
|
565 |
> |
if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ |
566 |
> |
|
567 |
> |
molecules[j].getCOM(COM); |
568 |
> |
|
569 |
> |
break; |
570 |
> |
|
571 |
> |
} |
572 |
> |
|
573 |
> |
} |
574 |
> |
|
575 |
> |
|
576 |
> |
|
577 |
> |
MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE, whichNode, |
578 |
> |
|
579 |
> |
MPI_COMM_WORLD); |
580 |
> |
|
581 |
> |
#endif |
582 |
> |
|
583 |
> |
|
584 |
> |
|
585 |
> |
(*parameters)[i].zPos = COM[whichDirection]; |
586 |
> |
|
587 |
> |
|
588 |
> |
|
589 |
> |
sprintf(painCave.errMsg, |
590 |
> |
|
591 |
> |
"ZConstraint warning: Does not specify zpos for z-constraint molecule " |
592 |
> |
|
593 |
> |
"initial z coornidate will be used \n"); |
594 |
> |
|
595 |
> |
painCave.isFatal = 0; |
596 |
> |
|
597 |
> |
simError(); |
598 |
> |
|
599 |
> |
} |
600 |
> |
|
601 |
> |
} |
602 |
> |
|
603 |
> |
}//end if (!zConsParaData) |
604 |
> |
|
605 |
> |
|
606 |
> |
|
607 |
> |
}//end if (!data) |
608 |
> |
|
609 |
> |
|
610 |
> |
|
611 |
> |
// |
612 |
> |
|
613 |
> |
#ifdef IS_MPI |
614 |
> |
|
615 |
> |
update(); |
616 |
> |
|
617 |
> |
#else |
618 |
> |
|
619 |
> |
int searchResult; |
620 |
> |
|
621 |
> |
|
622 |
> |
|
623 |
> |
for (int i = 0; i < nMols; i++){ |
624 |
> |
|
625 |
> |
searchResult = isZConstraintMol(&molecules[i]); |
626 |
> |
|
627 |
> |
|
628 |
> |
|
629 |
> |
if (searchResult > -1){ |
630 |
> |
|
631 |
> |
zconsMols.push_back(&molecules[i]); |
632 |
> |
|
633 |
> |
massOfZConsMols.push_back(molecules[i].getTotalMass()); |
634 |
> |
|
635 |
> |
|
636 |
> |
|
637 |
> |
zPos.push_back((*parameters)[searchResult].zPos); |
638 |
> |
|
639 |
> |
kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); |
640 |
> |
|
641 |
> |
|
642 |
> |
|
643 |
> |
if(usingSMD) |
644 |
> |
|
645 |
> |
cantVel.push_back((*parameters)[searchResult].cantVel); |
646 |
> |
|
647 |
> |
|
648 |
> |
|
649 |
> |
} |
650 |
> |
|
651 |
> |
else{ |
652 |
> |
|
653 |
> |
unconsMols.push_back(&molecules[i]); |
654 |
> |
|
655 |
> |
massOfUnconsMols.push_back(molecules[i].getTotalMass()); |
656 |
> |
|
657 |
> |
} |
658 |
> |
|
659 |
> |
} |
660 |
> |
|
661 |
> |
|
662 |
> |
|
663 |
> |
fz.resize(zconsMols.size()); |
664 |
> |
|
665 |
> |
curZPos.resize(zconsMols.size()); |
666 |
> |
|
667 |
> |
indexOfZConsMols.resize(zconsMols.size()); |
668 |
> |
|
669 |
> |
|
670 |
> |
|
671 |
> |
//determine the states of z-constraint molecules |
672 |
> |
|
673 |
> |
for (size_t i = 0; i < zconsMols.size(); i++){ |
674 |
> |
|
675 |
> |
indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); |
676 |
> |
|
677 |
> |
|
678 |
> |
|
679 |
> |
zconsMols[i]->getCOM(COM); |
680 |
> |
|
681 |
> |
|
682 |
> |
|
683 |
> |
if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){ |
684 |
> |
|
685 |
> |
states.push_back(zcsFixed); |
686 |
> |
|
687 |
> |
|
688 |
> |
|
689 |
> |
if (hasZConsGap) |
690 |
> |
|
691 |
> |
endFixTime.push_back(info->getTime() + zconsFixTime); |
692 |
> |
|
693 |
> |
} |
694 |
> |
|
695 |
> |
else{ |
696 |
> |
|
697 |
> |
states.push_back(zcsMoving); |
698 |
> |
|
699 |
> |
|
700 |
> |
|
701 |
> |
if (hasZConsGap) |
702 |
> |
|
703 |
> |
endFixTime.push_back(INFINITE_TIME); |
704 |
> |
|
705 |
> |
} |
706 |
> |
|
707 |
> |
|
708 |
> |
|
709 |
> |
if(usingSMD) |
710 |
> |
|
711 |
> |
cantPos.push_back(COM[whichDirection]); |
712 |
> |
|
713 |
> |
} |
714 |
> |
|
715 |
> |
|
716 |
> |
|
717 |
> |
if(usingSMD) |
718 |
> |
|
719 |
> |
prevCantPos = cantPos; |
720 |
> |
|
721 |
> |
#endif |
722 |
> |
|
723 |
> |
|
724 |
> |
|
725 |
> |
|
726 |
> |
|
727 |
> |
//get total masss of unconstraint molecules |
728 |
> |
|
729 |
> |
double totalMassOfUncons_local; |
730 |
> |
|
731 |
> |
totalMassOfUncons_local = 0; |
732 |
> |
|
733 |
> |
|
734 |
> |
|
735 |
> |
for (size_t i = 0; i < unconsMols.size(); i++) |
736 |
> |
|
737 |
> |
totalMassOfUncons_local += unconsMols[i]->getTotalMass(); |
738 |
> |
|
739 |
> |
|
740 |
> |
|
741 |
> |
#ifndef IS_MPI |
742 |
> |
|
743 |
> |
totalMassOfUncons = totalMassOfUncons_local; |
744 |
> |
|
745 |
> |
#else |
746 |
> |
|
747 |
> |
MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE, |
748 |
> |
|
749 |
> |
MPI_SUM, MPI_COMM_WORLD); |
750 |
> |
|
751 |
> |
#endif |
752 |
> |
|
753 |
> |
|
754 |
> |
|
755 |
> |
//get total number of unconstrained atoms |
756 |
> |
|
757 |
> |
int nUnconsAtoms_local; |
758 |
> |
|
759 |
> |
nUnconsAtoms_local = 0; |
760 |
> |
|
761 |
> |
for (int i = 0; i < (int) (unconsMols.size()); i++) |
762 |
> |
|
763 |
> |
nUnconsAtoms_local += unconsMols[i]->getNAtoms(); |
764 |
> |
|
765 |
> |
|
766 |
> |
|
767 |
> |
#ifndef IS_MPI |
768 |
> |
|
769 |
> |
totNumOfUnconsAtoms = nUnconsAtoms_local; |
770 |
> |
|
771 |
> |
#else |
772 |
> |
|
773 |
> |
MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM, |
774 |
> |
|
775 |
> |
MPI_COMM_WORLD); |
776 |
> |
|
777 |
> |
#endif |
778 |
> |
|
779 |
> |
|
780 |
> |
|
781 |
> |
forcePolicy->update(); |
782 |
> |
|
783 |
> |
} |
784 |
> |
|
785 |
> |
|
786 |
> |
|
787 |
> |
template<typename T> ZConstraint<T>::~ZConstraint(){ |
788 |
> |
|
789 |
> |
|
790 |
> |
|
791 |
> |
if (fzOut){ |
792 |
> |
|
793 |
> |
delete fzOut; |
794 |
> |
|
795 |
> |
} |
796 |
> |
|
797 |
> |
|
798 |
> |
|
799 |
> |
if (forcePolicy){ |
800 |
> |
|
801 |
> |
delete forcePolicy; |
802 |
> |
|
803 |
> |
} |
804 |
> |
|
805 |
> |
} |
806 |
> |
|
807 |
> |
|
808 |
> |
|
809 |
> |
|
810 |
> |
|
811 |
> |
/** |
812 |
> |
|
813 |
> |
* |
814 |
> |
|
815 |
> |
*/ |
816 |
> |
|
817 |
> |
|
818 |
> |
|
819 |
> |
#ifdef IS_MPI |
820 |
> |
|
821 |
> |
template<typename T> void ZConstraint<T>::update(){ |
822 |
> |
|
823 |
> |
double COM[3]; |
824 |
> |
|
825 |
> |
int index; |
826 |
> |
|
827 |
> |
|
828 |
> |
|
829 |
> |
zconsMols.clear(); |
830 |
> |
|
831 |
> |
massOfZConsMols.clear(); |
832 |
> |
|
833 |
> |
zPos.clear(); |
834 |
> |
|
835 |
> |
kz.clear(); |
836 |
> |
|
837 |
> |
cantPos.clear(); |
838 |
> |
|
839 |
> |
cantVel.clear(); |
840 |
> |
|
841 |
> |
|
842 |
> |
|
843 |
> |
unconsMols.clear(); |
844 |
> |
|
845 |
> |
massOfUnconsMols.clear(); |
846 |
> |
|
847 |
> |
|
848 |
> |
|
849 |
> |
|
850 |
> |
|
851 |
> |
//creat zconsMol and unconsMol lists |
852 |
> |
|
853 |
> |
for (int i = 0; i < nMols; i++){ |
854 |
> |
|
855 |
> |
index = isZConstraintMol(&molecules[i]); |
856 |
> |
|
857 |
> |
|
858 |
> |
|
859 |
> |
if (index > -1){ |
860 |
> |
|
861 |
> |
zconsMols.push_back(&molecules[i]); |
862 |
> |
|
863 |
> |
zPos.push_back((*parameters)[index].zPos); |
864 |
> |
|
865 |
> |
kz.push_back((*parameters)[index].kRatio * zForceConst); |
866 |
> |
|
867 |
> |
massOfZConsMols.push_back(molecules[i].getTotalMass()); |
868 |
> |
|
869 |
> |
|
870 |
> |
|
871 |
> |
if(usingSMD) |
872 |
> |
|
873 |
> |
cantVel.push_back((*parameters)[index].cantVel); |
874 |
> |
|
875 |
> |
|
876 |
> |
|
877 |
> |
} |
878 |
> |
|
879 |
> |
else{ |
880 |
> |
|
881 |
> |
unconsMols.push_back(&molecules[i]); |
882 |
> |
|
883 |
> |
massOfUnconsMols.push_back(molecules[i].getTotalMass()); |
884 |
> |
|
885 |
> |
} |
886 |
> |
|
887 |
> |
} |
888 |
> |
|
889 |
> |
|
890 |
> |
|
891 |
> |
fz.resize(zconsMols.size()); |
892 |
> |
|
893 |
> |
curZPos.resize(zconsMols.size()); |
894 |
> |
|
895 |
> |
indexOfZConsMols.resize(zconsMols.size()); |
896 |
> |
|
897 |
> |
|
898 |
> |
|
899 |
> |
for (size_t i = 0; i < zconsMols.size(); i++){ |
900 |
> |
|
901 |
> |
indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); |
902 |
> |
|
903 |
> |
} |
904 |
> |
|
905 |
> |
|
906 |
> |
|
907 |
> |
//determine the states of z-constraint molecules |
908 |
> |
|
909 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
910 |
> |
|
911 |
> |
|
912 |
> |
|
913 |
> |
zconsMols[i]->getCOM(COM); |
914 |
> |
|
915 |
> |
|
916 |
> |
|
917 |
> |
if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){ |
918 |
> |
|
919 |
> |
states.push_back(zcsFixed); |
920 |
> |
|
921 |
> |
|
922 |
> |
|
923 |
> |
if (hasZConsGap) |
924 |
> |
|
925 |
> |
endFixTime.push_back(info->getTime() + zconsFixTime); |
926 |
> |
|
927 |
> |
} |
928 |
> |
|
929 |
> |
else{ |
930 |
> |
|
931 |
> |
states.push_back(zcsMoving); |
932 |
> |
|
933 |
> |
|
934 |
> |
|
935 |
> |
if (hasZConsGap) |
936 |
> |
|
937 |
> |
endFixTime.push_back(INFINITE_TIME); |
938 |
> |
|
939 |
> |
} |
940 |
> |
|
941 |
> |
|
942 |
> |
|
943 |
> |
if(usingSMD) |
944 |
> |
|
945 |
> |
cantPos.push_back(COM[whichDirection]); |
946 |
> |
|
947 |
> |
} |
948 |
> |
|
949 |
> |
|
950 |
> |
|
951 |
> |
if(usingSMD) |
952 |
> |
|
953 |
> |
prevCantPos = cantPos; |
954 |
> |
|
955 |
> |
|
956 |
> |
|
957 |
> |
forcePolicy->update(); |
958 |
> |
|
959 |
> |
} |
960 |
> |
|
961 |
> |
|
962 |
> |
|
963 |
> |
#endif |
964 |
> |
|
965 |
> |
|
966 |
> |
|
967 |
> |
/** |
968 |
> |
|
969 |
> |
* Function Name: isZConstraintMol |
970 |
> |
|
971 |
> |
* Parameter |
972 |
> |
|
973 |
> |
* Molecule* mol |
974 |
> |
|
975 |
> |
* Return value: |
976 |
> |
|
977 |
> |
* -1, if the molecule is not z-constraint molecule, |
978 |
> |
|
979 |
> |
* other non-negative values, its index in indexOfAllZConsMols vector |
980 |
> |
|
981 |
> |
*/ |
982 |
> |
|
983 |
> |
|
984 |
> |
|
985 |
> |
template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){ |
986 |
> |
|
987 |
> |
int index; |
988 |
> |
|
989 |
> |
int low; |
990 |
> |
|
991 |
> |
int high; |
992 |
> |
|
993 |
> |
int mid; |
994 |
> |
|
995 |
> |
|
996 |
> |
|
997 |
> |
index = mol->getGlobalIndex(); |
998 |
> |
|
999 |
> |
|
1000 |
> |
|
1001 |
> |
low = 0; |
1002 |
> |
|
1003 |
> |
high = parameters->size() - 1; |
1004 |
> |
|
1005 |
> |
|
1006 |
> |
|
1007 |
> |
//Binary Search (we have sorted the array) |
1008 |
> |
|
1009 |
> |
while (low <= high){ |
1010 |
> |
|
1011 |
> |
mid = (low + high) / 2; |
1012 |
> |
|
1013 |
> |
if ((*parameters)[mid].zconsIndex == index) |
1014 |
> |
|
1015 |
> |
return mid; |
1016 |
> |
|
1017 |
> |
else if ((*parameters)[mid].zconsIndex > index) |
1018 |
> |
|
1019 |
> |
high = mid - 1; |
1020 |
> |
|
1021 |
> |
else |
1022 |
> |
|
1023 |
> |
low = mid + 1; |
1024 |
> |
|
1025 |
> |
} |
1026 |
> |
|
1027 |
> |
|
1028 |
> |
|
1029 |
> |
return -1; |
1030 |
> |
|
1031 |
> |
} |
1032 |
> |
|
1033 |
> |
|
1034 |
> |
|
1035 |
> |
template<typename T> void ZConstraint<T>::integrate(){ |
1036 |
> |
|
1037 |
> |
// creat zconsWriter |
1038 |
> |
|
1039 |
> |
fzOut = new ZConsWriter(zconsOutput.c_str(), parameters); |
1040 |
> |
|
1041 |
> |
|
1042 |
> |
|
1043 |
> |
if (!fzOut){ |
1044 |
> |
|
1045 |
> |
sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n"); |
1046 |
> |
|
1047 |
> |
painCave.isFatal = 1; |
1048 |
> |
|
1049 |
> |
simError(); |
1050 |
> |
|
1051 |
> |
} |
1052 |
> |
|
1053 |
> |
|
1054 |
> |
|
1055 |
> |
//zero out the velocities of center of mass of unconstrained molecules |
1056 |
> |
|
1057 |
> |
//and the velocities of center of mass of every single z-constrained molecueles |
1058 |
> |
|
1059 |
> |
zeroOutVel(); |
1060 |
> |
|
1061 |
> |
|
1062 |
> |
|
1063 |
> |
curZconsTime = zconsTime + info->getTime(); |
1064 |
> |
|
1065 |
> |
|
1066 |
> |
|
1067 |
> |
T::integrate(); |
1068 |
> |
|
1069 |
> |
} |
1070 |
> |
|
1071 |
> |
|
1072 |
> |
|
1073 |
> |
template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){ |
1074 |
> |
|
1075 |
> |
double zsys; |
1076 |
> |
|
1077 |
> |
double COM[3]; |
1078 |
> |
|
1079 |
> |
double force[3]; |
1080 |
> |
|
1081 |
> |
double zSysCOMVel; |
1082 |
> |
|
1083 |
> |
|
1084 |
> |
|
1085 |
> |
T::calcForce(calcPot, calcStress); |
1086 |
> |
|
1087 |
> |
|
1088 |
> |
|
1089 |
> |
|
1090 |
> |
|
1091 |
> |
if (hasZConsGap){ |
1092 |
> |
|
1093 |
> |
updateZPos(); |
1094 |
> |
|
1095 |
> |
} |
1096 |
> |
|
1097 |
> |
|
1098 |
> |
|
1099 |
> |
if (checkZConsState()){ |
1100 |
> |
|
1101 |
> |
zeroOutVel(); |
1102 |
> |
|
1103 |
> |
forcePolicy->update(); |
1104 |
> |
|
1105 |
> |
} |
1106 |
> |
|
1107 |
> |
|
1108 |
> |
|
1109 |
> |
zsys = calcZSys(); |
1110 |
> |
|
1111 |
> |
zSysCOMVel = calcSysCOMVel(); |
1112 |
> |
|
1113 |
> |
#ifdef IS_MPI |
1114 |
> |
|
1115 |
> |
if (worldRank == 0){ |
1116 |
> |
|
1117 |
> |
#endif |
1118 |
> |
|
1119 |
> |
|
1120 |
> |
|
1121 |
> |
#ifdef IS_MPI |
1122 |
> |
|
1123 |
> |
} |
1124 |
> |
|
1125 |
> |
#endif |
1126 |
> |
|
1127 |
> |
|
1128 |
> |
|
1129 |
> |
//do zconstraint force; |
1130 |
> |
|
1131 |
> |
if (haveFixedZMols()){ |
1132 |
> |
|
1133 |
> |
this->doZconstraintForce(); |
1134 |
> |
|
1135 |
> |
} |
1136 |
> |
|
1137 |
> |
|
1138 |
> |
|
1139 |
> |
//use external force to move the molecules to the specified positions |
1140 |
> |
|
1141 |
> |
if (haveMovingZMols()){ |
1142 |
> |
|
1143 |
> |
if (usingSMD) |
1144 |
> |
|
1145 |
> |
this->doHarmonic(cantPos); |
1146 |
> |
|
1147 |
> |
else |
1148 |
> |
|
1149 |
> |
this->doHarmonic(zPos); |
1150 |
> |
|
1151 |
> |
} |
1152 |
> |
|
1153 |
> |
|
1154 |
> |
|
1155 |
> |
//write out forces and current positions of z-constraint molecules |
1156 |
> |
|
1157 |
> |
if (info->getTime() >= curZconsTime){ |
1158 |
> |
|
1159 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
1160 |
> |
|
1161 |
> |
zconsMols[i]->getCOM(COM); |
1162 |
> |
|
1163 |
> |
curZPos[i] = COM[whichDirection]; |
1164 |
> |
|
1165 |
> |
|
1166 |
> |
|
1167 |
> |
//if the z-constraint molecule is still moving, just record its force |
1168 |
> |
|
1169 |
> |
if (states[i] == zcsMoving){ |
1170 |
> |
|
1171 |
> |
fz[i] = 0; |
1172 |
> |
|
1173 |
> |
Atom** movingZAtoms; |
1174 |
> |
|
1175 |
> |
movingZAtoms = zconsMols[i]->getMyAtoms(); |
1176 |
> |
|
1177 |
> |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
1178 |
> |
|
1179 |
> |
movingZAtoms[j]->getFrc(force); |
1180 |
> |
|
1181 |
> |
fz[i] += force[whichDirection]; |
1182 |
> |
|
1183 |
> |
} |
1184 |
> |
|
1185 |
> |
} |
1186 |
> |
|
1187 |
> |
} |
1188 |
> |
|
1189 |
> |
fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0], |
1190 |
> |
|
1191 |
> |
&curZPos[0], &zPos[0]); |
1192 |
> |
|
1193 |
> |
curZconsTime += zconsTime; |
1194 |
> |
|
1195 |
> |
} |
1196 |
> |
|
1197 |
> |
|
1198 |
> |
|
1199 |
> |
zSysCOMVel = calcSysCOMVel(); |
1200 |
> |
|
1201 |
> |
#ifdef IS_MPI |
1202 |
> |
|
1203 |
> |
if (worldRank == 0){ |
1204 |
> |
|
1205 |
> |
#endif |
1206 |
> |
|
1207 |
> |
#ifdef IS_MPI |
1208 |
> |
|
1209 |
> |
} |
1210 |
> |
|
1211 |
> |
#endif |
1212 |
> |
|
1213 |
> |
} |
1214 |
> |
|
1215 |
> |
|
1216 |
> |
|
1217 |
> |
|
1218 |
> |
|
1219 |
> |
template<typename T> double ZConstraint<T>::calcZSys(){ |
1220 |
> |
|
1221 |
> |
//calculate reference z coordinate for z-constraint molecules |
1222 |
> |
|
1223 |
> |
double totalMass_local; |
1224 |
> |
|
1225 |
> |
double totalMass; |
1226 |
> |
|
1227 |
> |
double totalMZ_local; |
1228 |
> |
|
1229 |
> |
double totalMZ; |
1230 |
> |
|
1231 |
> |
double massOfCurMol; |
1232 |
> |
|
1233 |
> |
double COM[3]; |
1234 |
> |
|
1235 |
> |
|
1236 |
> |
|
1237 |
> |
totalMass_local = 0; |
1238 |
> |
|
1239 |
> |
totalMZ_local = 0; |
1240 |
> |
|
1241 |
> |
|
1242 |
> |
|
1243 |
> |
for (int i = 0; i < nMols; i++){ |
1244 |
> |
|
1245 |
> |
massOfCurMol = molecules[i].getTotalMass(); |
1246 |
> |
|
1247 |
> |
molecules[i].getCOM(COM); |
1248 |
> |
|
1249 |
> |
|
1250 |
> |
|
1251 |
> |
totalMass_local += massOfCurMol; |
1252 |
> |
|
1253 |
> |
totalMZ_local += massOfCurMol * COM[whichDirection]; |
1254 |
> |
|
1255 |
> |
} |
1256 |
> |
|
1257 |
> |
|
1258 |
> |
|
1259 |
> |
|
1260 |
> |
|
1261 |
> |
#ifdef IS_MPI |
1262 |
> |
|
1263 |
> |
MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM, |
1264 |
> |
|
1265 |
> |
MPI_COMM_WORLD); |
1266 |
> |
|
1267 |
> |
MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
1268 |
> |
|
1269 |
> |
#else |
1270 |
> |
|
1271 |
> |
totalMass = totalMass_local; |
1272 |
> |
|
1273 |
> |
totalMZ = totalMZ_local; |
1274 |
> |
|
1275 |
> |
#endif |
1276 |
> |
|
1277 |
> |
|
1278 |
> |
|
1279 |
> |
double zsys; |
1280 |
> |
|
1281 |
> |
zsys = totalMZ / totalMass; |
1282 |
> |
|
1283 |
> |
|
1284 |
> |
|
1285 |
> |
return zsys; |
1286 |
> |
|
1287 |
> |
} |
1288 |
> |
|
1289 |
> |
|
1290 |
> |
|
1291 |
> |
template<typename T> void ZConstraint<T>::thermalize(void){ |
1292 |
> |
|
1293 |
> |
T::thermalize(); |
1294 |
> |
|
1295 |
> |
zeroOutVel(); |
1296 |
> |
|
1297 |
> |
} |
1298 |
> |
|
1299 |
> |
|
1300 |
> |
|
1301 |
> |
template<typename T> void ZConstraint<T>::zeroOutVel(){ |
1302 |
> |
|
1303 |
> |
Atom** fixedZAtoms; |
1304 |
> |
|
1305 |
> |
double COMvel[3]; |
1306 |
> |
|
1307 |
> |
double vel[3]; |
1308 |
> |
|
1309 |
> |
double zSysCOMVel; |
1310 |
> |
|
1311 |
> |
|
1312 |
> |
|
1313 |
> |
//zero out the velocities of center of mass of fixed z-constrained molecules |
1314 |
> |
|
1315 |
> |
|
1316 |
> |
|
1317 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
1318 |
> |
|
1319 |
> |
if (states[i] == zcsFixed){ |
1320 |
> |
|
1321 |
> |
zconsMols[i]->getCOMvel(COMvel); |
1322 |
> |
|
1323 |
> |
//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << std::endl; |
1324 |
> |
|
1325 |
> |
|
1326 |
> |
|
1327 |
> |
fixedZAtoms = zconsMols[i]->getMyAtoms(); |
1328 |
> |
|
1329 |
> |
|
1330 |
> |
|
1331 |
> |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
1332 |
> |
|
1333 |
> |
vel = fixedZAtoms[j]->getVel(); |
1334 |
> |
|
1335 |
> |
vel[whichDirection] -= COMvel[whichDirection]; |
1336 |
> |
|
1337 |
> |
fixedZAtoms[j]->setVel(vel); |
1338 |
> |
|
1339 |
> |
} |
1340 |
> |
|
1341 |
> |
|
1342 |
> |
|
1343 |
> |
zconsMols[i]->getCOMvel(COMvel); |
1344 |
> |
|
1345 |
> |
} |
1346 |
> |
|
1347 |
> |
} |
1348 |
> |
|
1349 |
> |
|
1350 |
> |
|
1351 |
> |
zSysCOMVel = calcSysCOMVel(); |
1352 |
> |
|
1353 |
> |
#ifdef IS_MPI |
1354 |
> |
|
1355 |
> |
if (worldRank == 0){ |
1356 |
> |
|
1357 |
> |
#endif |
1358 |
> |
|
1359 |
> |
#ifdef IS_MPI |
1360 |
> |
|
1361 |
> |
} |
1362 |
> |
|
1363 |
> |
#endif |
1364 |
> |
|
1365 |
> |
|
1366 |
> |
|
1367 |
> |
// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules |
1368 |
> |
|
1369 |
> |
double MVzOfMovingMols_local; |
1370 |
> |
|
1371 |
> |
double MVzOfMovingMols; |
1372 |
> |
|
1373 |
> |
double totalMassOfMovingZMols_local; |
1374 |
> |
|
1375 |
> |
double totalMassOfMovingZMols; |
1376 |
> |
|
1377 |
> |
|
1378 |
> |
|
1379 |
> |
MVzOfMovingMols_local = 0; |
1380 |
> |
|
1381 |
> |
totalMassOfMovingZMols_local = 0; |
1382 |
> |
|
1383 |
> |
|
1384 |
> |
|
1385 |
> |
for (int i = 0; i < (int) (unconsMols.size()); i++){ |
1386 |
> |
|
1387 |
> |
unconsMols[i]->getCOMvel(COMvel); |
1388 |
> |
|
1389 |
> |
MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
1390 |
> |
|
1391 |
> |
} |
1392 |
> |
|
1393 |
> |
|
1394 |
> |
|
1395 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
1396 |
> |
|
1397 |
> |
if (states[i] == zcsMoving){ |
1398 |
> |
|
1399 |
> |
zconsMols[i]->getCOMvel(COMvel); |
1400 |
> |
|
1401 |
> |
MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; |
1402 |
> |
|
1403 |
> |
totalMassOfMovingZMols_local += massOfZConsMols[i]; |
1404 |
> |
|
1405 |
> |
} |
1406 |
> |
|
1407 |
> |
} |
1408 |
> |
|
1409 |
> |
|
1410 |
> |
|
1411 |
> |
#ifndef IS_MPI |
1412 |
> |
|
1413 |
> |
MVzOfMovingMols = MVzOfMovingMols_local; |
1414 |
> |
|
1415 |
> |
totalMassOfMovingZMols = totalMassOfMovingZMols_local; |
1416 |
> |
|
1417 |
> |
#else |
1418 |
> |
|
1419 |
> |
MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE, |
1420 |
> |
|
1421 |
> |
MPI_SUM, MPI_COMM_WORLD); |
1422 |
> |
|
1423 |
> |
MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, |
1424 |
> |
|
1425 |
> |
MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
1426 |
> |
|
1427 |
> |
#endif |
1428 |
> |
|
1429 |
> |
|
1430 |
> |
|
1431 |
> |
double vzOfMovingMols; |
1432 |
> |
|
1433 |
> |
vzOfMovingMols = MVzOfMovingMols / |
1434 |
> |
|
1435 |
> |
(totalMassOfUncons + totalMassOfMovingZMols); |
1436 |
> |
|
1437 |
> |
|
1438 |
> |
|
1439 |
> |
//modify the velocites of unconstrained molecules |
1440 |
> |
|
1441 |
> |
Atom** unconsAtoms; |
1442 |
> |
|
1443 |
> |
for (int i = 0; i < (int) (unconsMols.size()); i++){ |
1444 |
> |
|
1445 |
> |
unconsAtoms = unconsMols[i]->getMyAtoms(); |
1446 |
> |
|
1447 |
> |
for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
1448 |
> |
|
1449 |
> |
vel = unconsAtoms[j]->getVel(); |
1450 |
> |
|
1451 |
> |
vel[whichDirection] -= vzOfMovingMols; |
1452 |
> |
|
1453 |
> |
unconsAtoms[j]->setVel(vel); |
1454 |
> |
|
1455 |
> |
} |
1456 |
> |
|
1457 |
> |
} |
1458 |
> |
|
1459 |
> |
|
1460 |
> |
|
1461 |
> |
//modify the velocities of moving z-constrained molecuels |
1462 |
> |
|
1463 |
> |
Atom** movingZAtoms; |
1464 |
> |
|
1465 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
1466 |
> |
|
1467 |
> |
if (states[i] == zcsMoving){ |
1468 |
> |
|
1469 |
> |
movingZAtoms = zconsMols[i]->getMyAtoms(); |
1470 |
> |
|
1471 |
> |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
1472 |
> |
|
1473 |
> |
vel = movingZAtoms[j]->getVel(); |
1474 |
> |
|
1475 |
> |
vel[whichDirection] -= vzOfMovingMols; |
1476 |
> |
|
1477 |
> |
movingZAtoms[j]->setVel(vel); |
1478 |
> |
|
1479 |
> |
} |
1480 |
> |
|
1481 |
> |
} |
1482 |
> |
|
1483 |
> |
} |
1484 |
> |
|
1485 |
> |
|
1486 |
> |
|
1487 |
> |
|
1488 |
> |
|
1489 |
> |
zSysCOMVel = calcSysCOMVel(); |
1490 |
> |
|
1491 |
> |
#ifdef IS_MPI |
1492 |
> |
|
1493 |
> |
if (worldRank == 0){ |
1494 |
> |
|
1495 |
> |
#endif |
1496 |
> |
|
1497 |
> |
#ifdef IS_MPI |
1498 |
> |
|
1499 |
> |
} |
1500 |
> |
|
1501 |
> |
#endif |
1502 |
> |
|
1503 |
> |
} |
1504 |
> |
|
1505 |
> |
|
1506 |
> |
|
1507 |
> |
|
1508 |
> |
|
1509 |
> |
template<typename T> void ZConstraint<T>::doZconstraintForce(){ |
1510 |
> |
|
1511 |
> |
Atom** zconsAtoms; |
1512 |
> |
|
1513 |
> |
double totalFZ; |
1514 |
> |
|
1515 |
> |
double totalFZ_local; |
1516 |
> |
|
1517 |
> |
double COM[3]; |
1518 |
> |
|
1519 |
> |
double force[3]; |
1520 |
> |
|
1521 |
> |
|
1522 |
> |
|
1523 |
> |
//constrain the molecules which do not reach the specified positions |
1524 |
> |
|
1525 |
> |
|
1526 |
> |
|
1527 |
> |
//Zero Out the force of z-contrained molecules |
1528 |
> |
|
1529 |
> |
totalFZ_local = 0; |
1530 |
> |
|
1531 |
> |
|
1532 |
> |
|
1533 |
> |
//calculate the total z-contrained force of fixed z-contrained molecules |
1534 |
> |
|
1535 |
> |
|
1536 |
> |
|
1537 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
1538 |
> |
|
1539 |
> |
if (states[i] == zcsFixed){ |
1540 |
> |
|
1541 |
> |
zconsMols[i]->getCOM(COM); |
1542 |
> |
|
1543 |
> |
zconsAtoms = zconsMols[i]->getMyAtoms(); |
1544 |
> |
|
1545 |
> |
|
1546 |
> |
|
1547 |
> |
fz[i] = 0; |
1548 |
> |
|
1549 |
> |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
1550 |
> |
|
1551 |
> |
zconsAtoms[j]->getFrc(force); |
1552 |
> |
|
1553 |
> |
fz[i] += force[whichDirection]; |
1554 |
> |
|
1555 |
> |
} |
1556 |
> |
|
1557 |
> |
totalFZ_local += fz[i]; |
1558 |
> |
|
1559 |
> |
|
1560 |
> |
|
1561 |
> |
} |
1562 |
> |
|
1563 |
> |
} |
1564 |
> |
|
1565 |
> |
|
1566 |
> |
|
1567 |
> |
//calculate total z-constraint force |
1568 |
> |
|
1569 |
> |
#ifdef IS_MPI |
1570 |
> |
|
1571 |
> |
MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
1572 |
> |
|
1573 |
> |
#else |
1574 |
> |
|
1575 |
> |
totalFZ = totalFZ_local; |
1576 |
> |
|
1577 |
> |
#endif |
1578 |
> |
|
1579 |
> |
|
1580 |
> |
|
1581 |
> |
|
1582 |
> |
|
1583 |
> |
// apply negative to fixed z-constrained molecues; |
1584 |
> |
|
1585 |
> |
force[0] = 0; |
1586 |
> |
|
1587 |
> |
force[1] = 0; |
1588 |
> |
|
1589 |
> |
force[2] = 0; |
1590 |
> |
|
1591 |
> |
|
1592 |
> |
|
1593 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
1594 |
> |
|
1595 |
> |
if (states[i] == zcsFixed){ |
1596 |
> |
|
1597 |
> |
int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms(); |
1598 |
> |
|
1599 |
> |
zconsAtoms = zconsMols[i]->getMyAtoms(); |
1600 |
> |
|
1601 |
> |
|
1602 |
> |
|
1603 |
> |
for (int j = 0; j < nAtomOfCurZConsMol; j++){ |
1604 |
> |
|
1605 |
> |
//force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; |
1606 |
> |
|
1607 |
> |
force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i], |
1608 |
> |
|
1609 |
> |
zconsAtoms[j], |
1610 |
> |
|
1611 |
> |
fz[i]); |
1612 |
> |
|
1613 |
> |
zconsAtoms[j]->addFrc(force); |
1614 |
> |
|
1615 |
> |
} |
1616 |
> |
|
1617 |
> |
} |
1618 |
> |
|
1619 |
> |
} |
1620 |
> |
|
1621 |
> |
|
1622 |
> |
|
1623 |
> |
force[0] = 0; |
1624 |
> |
|
1625 |
> |
force[1] = 0; |
1626 |
> |
|
1627 |
> |
force[2] = 0; |
1628 |
> |
|
1629 |
> |
|
1630 |
> |
|
1631 |
> |
//modify the forces of unconstrained molecules |
1632 |
> |
|
1633 |
> |
for (int i = 0; i < (int) (unconsMols.size()); i++){ |
1634 |
> |
|
1635 |
> |
Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
1636 |
> |
|
1637 |
> |
|
1638 |
> |
|
1639 |
> |
for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
1640 |
> |
|
1641 |
> |
//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); |
1642 |
> |
|
1643 |
> |
force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j], |
1644 |
> |
|
1645 |
> |
totalFZ); |
1646 |
> |
|
1647 |
> |
unconsAtoms[j]->addFrc(force); |
1648 |
> |
|
1649 |
> |
} |
1650 |
> |
|
1651 |
> |
} |
1652 |
> |
|
1653 |
> |
|
1654 |
> |
|
1655 |
> |
//modify the forces of moving z-constrained molecules |
1656 |
> |
|
1657 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
1658 |
> |
|
1659 |
> |
if (states[i] == zcsMoving){ |
1660 |
> |
|
1661 |
> |
Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
1662 |
> |
|
1663 |
> |
|
1664 |
> |
|
1665 |
> |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
1666 |
> |
|
1667 |
> |
//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); |
1668 |
> |
|
1669 |
> |
force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j], |
1670 |
> |
|
1671 |
> |
totalFZ); |
1672 |
> |
|
1673 |
> |
movingZAtoms[j]->addFrc(force); |
1674 |
> |
|
1675 |
> |
} |
1676 |
> |
|
1677 |
> |
} |
1678 |
> |
|
1679 |
> |
} |
1680 |
> |
|
1681 |
> |
} |
1682 |
> |
|
1683 |
> |
|
1684 |
> |
|
1685 |
> |
|
1686 |
> |
|
1687 |
> |
template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){ |
1688 |
> |
|
1689 |
> |
double force[3]; |
1690 |
> |
|
1691 |
> |
double harmonicU; |
1692 |
> |
|
1693 |
> |
double harmonicF; |
1694 |
> |
|
1695 |
> |
double COM[3]; |
1696 |
> |
|
1697 |
> |
double diff; |
1698 |
> |
|
1699 |
> |
double totalFZ_local; |
1700 |
> |
|
1701 |
> |
double totalFZ; |
1702 |
> |
|
1703 |
> |
|
1704 |
> |
|
1705 |
> |
force[0] = 0; |
1706 |
> |
|
1707 |
> |
force[1] = 0; |
1708 |
> |
|
1709 |
> |
force[2] = 0; |
1710 |
> |
|
1711 |
> |
|
1712 |
> |
|
1713 |
> |
totalFZ_local = 0; |
1714 |
> |
|
1715 |
> |
|
1716 |
> |
|
1717 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
1718 |
> |
|
1719 |
> |
if (states[i] == zcsMoving){ |
1720 |
> |
|
1721 |
> |
zconsMols[i]->getCOM(COM); |
1722 |
> |
|
1723 |
> |
|
1724 |
> |
|
1725 |
> |
diff = COM[whichDirection] - resPos[i]; |
1726 |
> |
|
1727 |
> |
|
1728 |
> |
|
1729 |
> |
harmonicU = 0.5 * kz[i] * diff * diff; |
1730 |
> |
|
1731 |
> |
info->lrPot += harmonicU; |
1732 |
> |
|
1733 |
> |
|
1734 |
> |
|
1735 |
> |
harmonicF = -kz[i] * diff; |
1736 |
> |
|
1737 |
> |
totalFZ_local += harmonicF; |
1738 |
> |
|
1739 |
> |
|
1740 |
> |
|
1741 |
> |
//adjust force |
1742 |
> |
|
1743 |
> |
|
1744 |
> |
|
1745 |
> |
Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
1746 |
> |
|
1747 |
> |
|
1748 |
> |
|
1749 |
> |
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
1750 |
> |
|
1751 |
> |
force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], |
1752 |
> |
|
1753 |
> |
movingZAtoms[j], |
1754 |
> |
|
1755 |
> |
harmonicF); |
1756 |
> |
|
1757 |
> |
movingZAtoms[j]->addFrc(force); |
1758 |
> |
|
1759 |
> |
} |
1760 |
> |
|
1761 |
> |
} |
1762 |
> |
|
1763 |
> |
} |
1764 |
> |
|
1765 |
> |
|
1766 |
> |
|
1767 |
> |
#ifndef IS_MPI |
1768 |
> |
|
1769 |
> |
totalFZ = totalFZ_local; |
1770 |
> |
|
1771 |
> |
#else |
1772 |
> |
|
1773 |
> |
MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
1774 |
> |
|
1775 |
> |
#endif |
1776 |
> |
|
1777 |
> |
|
1778 |
> |
|
1779 |
> |
force[0] = 0; |
1780 |
> |
|
1781 |
> |
force[1] = 0; |
1782 |
> |
|
1783 |
> |
force[2] = 0; |
1784 |
> |
|
1785 |
> |
|
1786 |
> |
|
1787 |
> |
//modify the forces of unconstrained molecules |
1788 |
> |
|
1789 |
> |
for (int i = 0; i < (int) (unconsMols.size()); i++){ |
1790 |
> |
|
1791 |
> |
Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
1792 |
> |
|
1793 |
> |
|
1794 |
> |
|
1795 |
> |
for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ |
1796 |
> |
|
1797 |
> |
//force[whichDirection] = - totalFZ /totNumOfUnconsAtoms; |
1798 |
> |
|
1799 |
> |
force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j], |
1800 |
> |
|
1801 |
> |
totalFZ); |
1802 |
> |
|
1803 |
> |
unconsAtoms[j]->addFrc(force); |
1804 |
> |
|
1805 |
> |
} |
1806 |
> |
|
1807 |
> |
} |
1808 |
> |
|
1809 |
> |
|
1810 |
> |
|
1811 |
> |
} |
1812 |
> |
|
1813 |
> |
|
1814 |
> |
|
1815 |
> |
template<typename T> bool ZConstraint<T>::checkZConsState(){ |
1816 |
> |
|
1817 |
> |
double COM[3]; |
1818 |
> |
|
1819 |
> |
double diff; |
1820 |
> |
|
1821 |
> |
|
1822 |
> |
|
1823 |
> |
int changed_local; |
1824 |
> |
|
1825 |
> |
int changed; |
1826 |
> |
|
1827 |
> |
|
1828 |
> |
|
1829 |
> |
changed_local = 0; |
1830 |
> |
|
1831 |
> |
|
1832 |
> |
|
1833 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
1834 |
> |
|
1835 |
> |
zconsMols[i]->getCOM(COM); |
1836 |
> |
|
1837 |
> |
diff = fabs(COM[whichDirection] - zPos[i]); |
1838 |
> |
|
1839 |
> |
if (diff <= zconsTol && states[i] == zcsMoving){ |
1840 |
> |
|
1841 |
> |
states[i] = zcsFixed; |
1842 |
> |
|
1843 |
> |
changed_local = 1; |
1844 |
> |
|
1845 |
> |
|
1846 |
> |
|
1847 |
> |
if(usingSMD) |
1848 |
> |
|
1849 |
> |
prevCantPos = cantPos; |
1850 |
> |
|
1851 |
> |
|
1852 |
> |
|
1853 |
> |
if (hasZConsGap) |
1854 |
> |
|
1855 |
> |
endFixTime[i] = info->getTime() + zconsFixTime; |
1856 |
> |
|
1857 |
> |
} |
1858 |
> |
|
1859 |
> |
else if (diff > zconsTol && states[i] == zcsFixed){ |
1860 |
> |
|
1861 |
> |
states[i] = zcsMoving; |
1862 |
> |
|
1863 |
> |
changed_local = 1; |
1864 |
> |
|
1865 |
> |
|
1866 |
> |
|
1867 |
> |
if(usingSMD) |
1868 |
> |
|
1869 |
> |
cantPos = prevCantPos; |
1870 |
> |
|
1871 |
> |
|
1872 |
> |
|
1873 |
> |
if (hasZConsGap) |
1874 |
> |
|
1875 |
> |
endFixTime[i] = INFINITE_TIME; |
1876 |
> |
|
1877 |
> |
} |
1878 |
> |
|
1879 |
> |
} |
1880 |
> |
|
1881 |
> |
|
1882 |
> |
|
1883 |
> |
#ifndef IS_MPI |
1884 |
> |
|
1885 |
> |
changed = changed_local; |
1886 |
> |
|
1887 |
> |
#else |
1888 |
> |
|
1889 |
> |
MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
1890 |
> |
|
1891 |
> |
#endif |
1892 |
> |
|
1893 |
> |
|
1894 |
> |
|
1895 |
> |
return (changed > 0); |
1896 |
> |
|
1897 |
> |
} |
1898 |
> |
|
1899 |
> |
|
1900 |
> |
|
1901 |
> |
template<typename T> bool ZConstraint<T>::haveFixedZMols(){ |
1902 |
> |
|
1903 |
> |
int havingFixed_local; |
1904 |
> |
|
1905 |
> |
int havingFixed; |
1906 |
> |
|
1907 |
> |
|
1908 |
> |
|
1909 |
> |
havingFixed_local = 0; |
1910 |
> |
|
1911 |
> |
|
1912 |
> |
|
1913 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++) |
1914 |
> |
|
1915 |
> |
if (states[i] == zcsFixed){ |
1916 |
> |
|
1917 |
> |
havingFixed_local = 1; |
1918 |
> |
|
1919 |
> |
break; |
1920 |
> |
|
1921 |
> |
} |
1922 |
> |
|
1923 |
> |
|
1924 |
> |
|
1925 |
> |
#ifndef IS_MPI |
1926 |
> |
|
1927 |
> |
havingFixed = havingFixed_local; |
1928 |
> |
|
1929 |
> |
#else |
1930 |
> |
|
1931 |
> |
MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM, |
1932 |
> |
|
1933 |
> |
MPI_COMM_WORLD); |
1934 |
> |
|
1935 |
> |
#endif |
1936 |
> |
|
1937 |
> |
|
1938 |
> |
|
1939 |
> |
return (havingFixed > 0); |
1940 |
> |
|
1941 |
> |
} |
1942 |
> |
|
1943 |
> |
|
1944 |
> |
|
1945 |
> |
|
1946 |
> |
|
1947 |
> |
template<typename T> bool ZConstraint<T>::haveMovingZMols(){ |
1948 |
> |
|
1949 |
> |
int havingMoving_local; |
1950 |
> |
|
1951 |
> |
int havingMoving; |
1952 |
> |
|
1953 |
> |
|
1954 |
> |
|
1955 |
> |
havingMoving_local = 0; |
1956 |
> |
|
1957 |
> |
|
1958 |
> |
|
1959 |
> |
for (int i = 0; i < (int) (zconsMols.size()); i++) |
1960 |
> |
|
1961 |
> |
if (states[i] == zcsMoving){ |
1962 |
> |
|
1963 |
> |
havingMoving_local = 1; |
1964 |
> |
|
1965 |
> |
break; |
1966 |
> |
|
1967 |
> |
} |
1968 |
> |
|
1969 |
> |
|
1970 |
> |
|
1971 |
> |
#ifndef IS_MPI |
1972 |
> |
|
1973 |
> |
havingMoving = havingMoving_local; |
1974 |
> |
|
1975 |
> |
#else |
1976 |
> |
|
1977 |
> |
MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM, |
1978 |
> |
|
1979 |
> |
MPI_COMM_WORLD); |
1980 |
> |
|
1981 |
> |
#endif |
1982 |
> |
|
1983 |
> |
|
1984 |
> |
|
1985 |
> |
return (havingMoving > 0); |
1986 |
> |
|
1987 |
> |
} |
1988 |
> |
|
1989 |
> |
|
1990 |
> |
|
1991 |
> |
|
1992 |
> |
|
1993 |
> |
template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){ |
1994 |
> |
|
1995 |
> |
double MVzOfMovingMols_local; |
1996 |
> |
|
1997 |
> |
double MVzOfMovingMols; |
1998 |
> |
|
1999 |
> |
double totalMassOfMovingZMols_local; |
2000 |
> |
|
2001 |
> |
double totalMassOfMovingZMols; |
2002 |
> |
|
2003 |
> |
double COMvel[3]; |
2004 |
> |
|
2005 |
> |
|
2006 |
> |
|
2007 |
> |
MVzOfMovingMols_local = 0; |
2008 |
> |
|
2009 |
> |
totalMassOfMovingZMols_local = 0; |
2010 |
> |
|
2011 |
> |
|
2012 |
> |
|
2013 |
> |
for (int i = 0; i < unconsMols.size(); i++){ |
2014 |
> |
|
2015 |
> |
unconsMols[i]->getCOMvel(COMvel); |
2016 |
> |
|
2017 |
> |
MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
2018 |
> |
|
2019 |
> |
} |
2020 |
> |
|
2021 |
> |
|
2022 |
> |
|
2023 |
> |
for (int i = 0; i < zconsMols.size(); i++){ |
2024 |
> |
|
2025 |
> |
if (states[i] == zcsMoving){ |
2026 |
> |
|
2027 |
> |
zconsMols[i]->getCOMvel(COMvel); |
2028 |
> |
|
2029 |
> |
MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; |
2030 |
> |
|
2031 |
> |
totalMassOfMovingZMols_local += massOfZConsMols[i]; |
2032 |
> |
|
2033 |
> |
} |
2034 |
> |
|
2035 |
> |
} |
2036 |
> |
|
2037 |
> |
|
2038 |
> |
|
2039 |
> |
#ifndef IS_MPI |
2040 |
> |
|
2041 |
> |
MVzOfMovingMols = MVzOfMovingMols_local; |
2042 |
> |
|
2043 |
> |
totalMassOfMovingZMols = totalMassOfMovingZMols_local; |
2044 |
> |
|
2045 |
> |
#else |
2046 |
> |
|
2047 |
> |
MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE, |
2048 |
> |
|
2049 |
> |
MPI_SUM, MPI_COMM_WORLD); |
2050 |
> |
|
2051 |
> |
MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, |
2052 |
> |
|
2053 |
> |
MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
2054 |
> |
|
2055 |
> |
#endif |
2056 |
> |
|
2057 |
> |
|
2058 |
> |
|
2059 |
> |
double vzOfMovingMols; |
2060 |
> |
|
2061 |
> |
vzOfMovingMols = MVzOfMovingMols / |
2062 |
> |
|
2063 |
> |
(totalMassOfUncons + totalMassOfMovingZMols); |
2064 |
> |
|
2065 |
> |
|
2066 |
> |
|
2067 |
> |
return vzOfMovingMols; |
2068 |
> |
|
2069 |
> |
} |
2070 |
> |
|
2071 |
> |
|
2072 |
> |
|
2073 |
> |
template<typename T> double ZConstraint<T>::calcSysCOMVel(){ |
2074 |
> |
|
2075 |
> |
double COMvel[3]; |
2076 |
> |
|
2077 |
> |
double tempMVz_local; |
2078 |
> |
|
2079 |
> |
double tempMVz; |
2080 |
> |
|
2081 |
> |
double massOfZCons_local; |
2082 |
> |
|
2083 |
> |
double massOfZCons; |
2084 |
> |
|
2085 |
> |
|
2086 |
> |
|
2087 |
> |
|
2088 |
> |
|
2089 |
> |
tempMVz_local = 0; |
2090 |
> |
|
2091 |
> |
|
2092 |
> |
|
2093 |
> |
for (int i = 0 ; i < nMols; i++){ |
2094 |
> |
|
2095 |
> |
molecules[i].getCOMvel(COMvel); |
2096 |
> |
|
2097 |
> |
tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection]; |
2098 |
> |
|
2099 |
> |
} |
2100 |
> |
|
2101 |
> |
|
2102 |
> |
|
2103 |
> |
massOfZCons_local = 0; |
2104 |
> |
|
2105 |
> |
|
2106 |
> |
|
2107 |
> |
for (int i = 0; i < (int) (massOfZConsMols.size()); i++){ |
2108 |
> |
|
2109 |
> |
massOfZCons_local += massOfZConsMols[i]; |
2110 |
> |
|
2111 |
> |
} |
2112 |
> |
|
2113 |
> |
#ifndef IS_MPI |
2114 |
> |
|
2115 |
> |
massOfZCons = massOfZCons_local; |
2116 |
> |
|
2117 |
> |
tempMVz = tempMVz_local; |
2118 |
> |
|
2119 |
> |
#else |
2120 |
> |
|
2121 |
> |
MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM, |
2122 |
> |
|
2123 |
> |
MPI_COMM_WORLD); |
2124 |
> |
|
2125 |
> |
MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
2126 |
> |
|
2127 |
> |
#endif |
2128 |
> |
|
2129 |
> |
|
2130 |
> |
|
2131 |
> |
return tempMVz / (totalMassOfUncons + massOfZCons); |
2132 |
> |
|
2133 |
> |
} |
2134 |
> |
|
2135 |
> |
|
2136 |
> |
|
2137 |
> |
template<typename T> double ZConstraint<T>::calcTotalForce(){ |
2138 |
> |
|
2139 |
> |
double force[3]; |
2140 |
> |
|
2141 |
> |
double totalForce_local; |
2142 |
> |
|
2143 |
> |
double totalForce; |
2144 |
> |
|
2145 |
> |
|
2146 |
> |
|
2147 |
> |
totalForce_local = 0; |
2148 |
> |
|
2149 |
> |
|
2150 |
> |
|
2151 |
> |
for (int i = 0; i < nAtoms; i++){ |
2152 |
> |
|
2153 |
> |
atoms[i]->getFrc(force); |
2154 |
> |
|
2155 |
> |
totalForce_local += force[whichDirection]; |
2156 |
> |
|
2157 |
> |
} |
2158 |
> |
|
2159 |
> |
|
2160 |
> |
|
2161 |
> |
#ifndef IS_MPI |
2162 |
> |
|
2163 |
> |
totalForce = totalForce_local; |
2164 |
> |
|
2165 |
> |
#else |
2166 |
> |
|
2167 |
> |
MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM, |
2168 |
> |
|
2169 |
> |
MPI_COMM_WORLD); |
2170 |
> |
|
2171 |
> |
#endif |
2172 |
> |
|
2173 |
> |
|
2174 |
> |
|
2175 |
> |
return totalForce; |
2176 |
> |
|
2177 |
> |
} |
2178 |
> |
|
2179 |
> |
|
2180 |
> |
|
2181 |
> |
template<typename T> void ZConstraint<T>::PolicyByNumber::update(){ |
2182 |
> |
|
2183 |
> |
//calculate the number of atoms of moving z-constrained molecules |
2184 |
> |
|
2185 |
> |
int nMovingZAtoms_local; |
2186 |
> |
|
2187 |
> |
int nMovingZAtoms; |
2188 |
> |
|
2189 |
> |
|
2190 |
> |
|
2191 |
> |
nMovingZAtoms_local = 0; |
2192 |
> |
|
2193 |
> |
for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++) |
2194 |
> |
|
2195 |
> |
if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){ |
2196 |
> |
|
2197 |
> |
nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms(); |
2198 |
> |
|
2199 |
> |
} |
2200 |
> |
|
2201 |
> |
|
2202 |
> |
|
2203 |
> |
#ifdef IS_MPI |
2204 |
> |
|
2205 |
> |
MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, |
2206 |
> |
|
2207 |
> |
MPI_COMM_WORLD); |
2208 |
> |
|
2209 |
> |
#else |
2210 |
> |
|
2211 |
> |
nMovingZAtoms = nMovingZAtoms_local; |
2212 |
> |
|
2213 |
> |
#endif |
2214 |
> |
|
2215 |
> |
totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms; |
2216 |
> |
|
2217 |
> |
} |
2218 |
> |
|
2219 |
> |
|
2220 |
> |
|
2221 |
> |
template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, |
2222 |
> |
|
2223 |
> |
Atom* atom, |
2224 |
> |
|
2225 |
> |
double totalForce){ |
2226 |
> |
|
2227 |
> |
return totalForce / mol->getNAtoms(); |
2228 |
> |
|
2229 |
> |
} |
2230 |
> |
|
2231 |
> |
|
2232 |
> |
|
2233 |
> |
template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, |
2234 |
> |
|
2235 |
> |
double totalForce){ |
2236 |
> |
|
2237 |
> |
return totalForce / totNumOfMovingAtoms; |
2238 |
> |
|
2239 |
> |
} |
2240 |
> |
|
2241 |
> |
|
2242 |
> |
|
2243 |
> |
template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, |
2244 |
> |
|
2245 |
> |
Atom* atom, |
2246 |
> |
|
2247 |
> |
double totalForce){ |
2248 |
> |
|
2249 |
> |
return totalForce / mol->getNAtoms(); |
2250 |
> |
|
2251 |
> |
} |
2252 |
> |
|
2253 |
> |
|
2254 |
> |
|
2255 |
> |
template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, |
2256 |
> |
|
2257 |
> |
double totalForce){ |
2258 |
> |
|
2259 |
> |
return totalForce / zconsIntegrator->totNumOfUnconsAtoms; |
2260 |
> |
|
2261 |
> |
} |
2262 |
> |
|
2263 |
> |
|
2264 |
> |
|
2265 |
> |
|
2266 |
> |
|
2267 |
> |
template<typename T> void ZConstraint<T>::PolicyByMass::update(){ |
2268 |
> |
|
2269 |
> |
//calculate the number of atoms of moving z-constrained molecules |
2270 |
> |
|
2271 |
> |
double massOfMovingZAtoms_local; |
2272 |
> |
|
2273 |
> |
double massOfMovingZAtoms; |
2274 |
> |
|
2275 |
> |
|
2276 |
> |
|
2277 |
> |
massOfMovingZAtoms_local = 0; |
2278 |
> |
|
2279 |
> |
for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++) |
2280 |
> |
|
2281 |
> |
if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){ |
2282 |
> |
|
2283 |
> |
massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass(); |
2284 |
> |
|
2285 |
> |
} |
2286 |
> |
|
2287 |
> |
|
2288 |
> |
|
2289 |
> |
#ifdef IS_MPI |
2290 |
> |
|
2291 |
> |
MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE, |
2292 |
> |
|
2293 |
> |
MPI_SUM, MPI_COMM_WORLD); |
2294 |
> |
|
2295 |
> |
#else |
2296 |
> |
|
2297 |
> |
massOfMovingZAtoms = massOfMovingZAtoms_local; |
2298 |
> |
|
2299 |
> |
#endif |
2300 |
> |
|
2301 |
> |
totMassOfMovingAtoms = massOfMovingZAtoms + |
2302 |
> |
|
2303 |
> |
zconsIntegrator->totalMassOfUncons; |
2304 |
> |
|
2305 |
> |
} |
2306 |
> |
|
2307 |
> |
|
2308 |
> |
|
2309 |
> |
template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, |
2310 |
> |
|
2311 |
> |
Atom* atom, |
2312 |
> |
|
2313 |
> |
double totalForce){ |
2314 |
> |
|
2315 |
> |
return totalForce * atom->getMass() / mol->getTotalMass(); |
2316 |
> |
|
2317 |
> |
} |
2318 |
> |
|
2319 |
> |
|
2320 |
> |
|
2321 |
> |
template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom, |
2322 |
> |
|
2323 |
> |
double totalForce){ |
2324 |
> |
|
2325 |
> |
return totalForce * atom->getMass() / totMassOfMovingAtoms; |
2326 |
> |
|
2327 |
> |
} |
2328 |
> |
|
2329 |
> |
|
2330 |
> |
|
2331 |
> |
template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, |
2332 |
> |
|
2333 |
> |
Atom* atom, |
2334 |
> |
|
2335 |
> |
double totalForce){ |
2336 |
> |
|
2337 |
> |
return totalForce * atom->getMass() / mol->getTotalMass(); |
2338 |
> |
|
2339 |
> |
} |
2340 |
> |
|
2341 |
> |
|
2342 |
> |
|
2343 |
> |
template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, |
2344 |
> |
|
2345 |
> |
double totalForce){ |
2346 |
> |
|
2347 |
> |
return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons; |
2348 |
> |
|
2349 |
> |
} |
2350 |
> |
|
2351 |
> |
|
2352 |
> |
|
2353 |
> |
template<typename T> void ZConstraint<T>::updateZPos(){ |
2354 |
> |
|
2355 |
> |
double curTime; |
2356 |
> |
|
2357 |
> |
double COM[3]; |
2358 |
> |
|
2359 |
> |
|
2360 |
> |
|
2361 |
> |
curTime = info->getTime(); |
2362 |
> |
|
2363 |
> |
|
2364 |
> |
|
2365 |
> |
for (size_t i = 0; i < zconsMols.size(); i++){ |
2366 |
> |
|
2367 |
> |
|
2368 |
> |
|
2369 |
> |
if (states[i] == zcsFixed && curTime >= endFixTime[i]){ |
2370 |
> |
|
2371 |
> |
zPos[i] += zconsGap; |
2372 |
> |
|
2373 |
> |
|
2374 |
> |
|
2375 |
> |
if (usingSMD){ |
2376 |
> |
|
2377 |
> |
zconsMols[i]->getCOM(COM); |
2378 |
> |
|
2379 |
> |
cantPos[i] = COM[whichDirection]; |
2380 |
> |
|
2381 |
> |
} |
2382 |
> |
|
2383 |
> |
|
2384 |
> |
|
2385 |
> |
} |
2386 |
> |
|
2387 |
> |
|
2388 |
> |
|
2389 |
> |
} |
2390 |
> |
|
2391 |
> |
|
2392 |
> |
|
2393 |
> |
} |
2394 |
> |
|
2395 |
> |
|
2396 |
> |
|
2397 |
> |
template<typename T> void ZConstraint<T>::updateCantPos(){ |
2398 |
> |
|
2399 |
> |
double curTime; |
2400 |
> |
|
2401 |
> |
double dt; |
2402 |
> |
|
2403 |
> |
|
2404 |
> |
|
2405 |
> |
curTime = info->getTime(); |
2406 |
> |
|
2407 |
> |
dt = info->dt; |
2408 |
> |
|
2409 |
> |
|
2410 |
> |
|
2411 |
> |
for (size_t i = 0; i < zconsMols.size(); i++){ |
2412 |
> |
|
2413 |
> |
if (states[i] == zcsMoving){ |
2414 |
> |
|
2415 |
> |
cantPos[i] += cantVel[i] * dt; |
2416 |
> |
|
2417 |
> |
} |
2418 |
> |
|
2419 |
> |
} |
2420 |
> |
|
2421 |
> |
|
2422 |
> |
|
2423 |
> |
} |
2424 |
> |
|