# | Line 5 | Line 5 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
---|---|---|
5 | const double INFINITE_TIME = 10e30; | |
6 | template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, | |
7 | ForceFields* the_ff): T(theInfo, the_ff), | |
8 | – | indexOfZConsMols(NULL), |
9 | – | fz(NULL), |
10 | – | curZPos(NULL), |
8 | fzOut(NULL), | |
9 | curZconsTime(0), | |
10 | forcePolicy(NULL), | |
11 | + | usingSMD(false), |
12 | hasZConsGap(false){ | |
13 | //get properties from SimInfo | |
14 | GenericData* data; | |
# | Line 21 | Line 19 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
19 | DoubleData* fixtime; | |
20 | StringData* policy; | |
21 | StringData* filename; | |
22 | + | IntData* smdFlag; |
23 | double COM[3]; | |
24 | ||
25 | //by default, the direction of constraint is z | |
# | Line 150 | Line 149 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
149 | ||
150 | if (data){ | |
151 | gap = dynamic_cast<DoubleData*>(data); | |
153 | – | } |
152 | ||
153 | < | if (!gap){ |
154 | < | sprintf(painCave.errMsg, |
155 | < | "ZConstraint error: Can not get property from SimInfo\n"); |
156 | < | painCave.isFatal = 1; |
157 | < | simError(); |
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 | } | |
161 | – | else{ |
162 | – | this->hasZConsGap = true; |
163 | – | this->zconsGap = gap->getData(); |
164 | – | } |
164 | ||
165 | + | |
166 | + | |
167 | data = info->getProperty(ZCONSFIXTIME_ID); | |
168 | ||
169 | if (data){ | |
170 | fixtime = dynamic_cast<DoubleData*>(data); | |
171 | + | if (!fixtime){ |
172 | + | sprintf(painCave.errMsg, |
173 | + | "ZConstraint error: Can not get zconsFixTime from SimInfo\n"); |
174 | + | painCave.isFatal = 1; |
175 | + | simError(); |
176 | + | } |
177 | + | else{ |
178 | + | this->zconsFixTime = fixtime->getData(); |
179 | + | } |
180 | } | |
181 | < | |
182 | < | if (!fixtime){ |
183 | < | sprintf(painCave.errMsg, |
184 | < | "ZConstraint error: Can not get property from SimInfo\n"); |
185 | < | painCave.isFatal = 1; |
176 | < | simError(); |
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 | < | else{ |
188 | < | this->zconsFixTime = fixtime->getData(); |
187 | > | |
188 | > | |
189 | > | |
190 | > | data = info->getProperty(ZCONSUSINGSMD_ID); |
191 | > | |
192 | > | if (data){ |
193 | > | smdFlag = dynamic_cast<IntData*>(data); |
194 | > | |
195 | > | if (!smdFlag){ |
196 | > | sprintf(painCave.errMsg, |
197 | > | "ZConstraint error: Can not get property from SimInfo\n"); |
198 | > | painCave.isFatal = 1; |
199 | > | simError(); |
200 | > | } |
201 | > | else{ |
202 | > | this->usingSMD= smdFlag->getData() ? true : false; |
203 | > | } |
204 | > | |
205 | } | |
206 | ||
207 | ||
# | Line 292 | Line 317 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
317 | massOfZConsMols.push_back(molecules[i].getTotalMass()); | |
318 | ||
319 | zPos.push_back((*parameters)[searchResult].zPos); | |
295 | – | // cout << "index: "<< (*parameters)[searchResult].zconsIndex |
296 | – | // <<"\tzPos = " << (*parameters)[searchResult].zPos << endl; |
297 | – | |
320 | kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); | |
321 | < | molecules[i].getCOM(COM); |
321 | > | |
322 | > | if(usingSMD) |
323 | > | cantVel.push_back((*parameters)[searchResult].cantVel); |
324 | > | |
325 | } | |
326 | else{ | |
327 | unconsMols.push_back(&molecules[i]); | |
# | Line 304 | Line 329 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
329 | } | |
330 | } | |
331 | ||
332 | < | fz = new double[zconsMols.size()]; |
333 | < | curZPos = new double[zconsMols.size()]; |
334 | < | indexOfZConsMols = new int [zconsMols.size()]; |
332 | > | fz.resize(zconsMols.size()); |
333 | > | curZPos.resize(zconsMols.size()); |
334 | > | indexOfZConsMols.resize(zconsMols.size()); |
335 | ||
311 | – | if (!fz || !curZPos || !indexOfZConsMols){ |
312 | – | sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n"); |
313 | – | painCave.isFatal = 1; |
314 | – | simError(); |
315 | – | } |
316 | – | |
336 | //determine the states of z-constraint molecules | |
337 | < | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
337 | > | for (size_t i = 0; i < zconsMols.size(); i++){ |
338 | indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); | |
339 | ||
340 | zconsMols[i]->getCOM(COM); | |
341 | + | |
342 | if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){ | |
343 | states.push_back(zcsFixed); | |
344 | ||
# | Line 331 | Line 351 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
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 (int i = 0; i < (int) (unconsMols.size()); i++) |
368 | > | for (size_t i = 0; i < unconsMols.size(); i++) |
369 | totalMassOfUncons_local += unconsMols[i]->getTotalMass(); | |
370 | ||
371 | #ifndef IS_MPI | |
# | Line 366 | Line 392 | template<typename T> ZConstraint<T>::~ZConstraint(){ | |
392 | } | |
393 | ||
394 | template<typename T> ZConstraint<T>::~ZConstraint(){ | |
369 | – | if (fz){ |
370 | – | delete[] fz; |
371 | – | } |
395 | ||
373 | – | if (curZPos){ |
374 | – | delete[] curZPos; |
375 | – | } |
376 | – | |
377 | – | if (indexOfZConsMols){ |
378 | – | delete[] indexOfZConsMols; |
379 | – | } |
380 | – | |
396 | if (fzOut){ | |
397 | delete fzOut; | |
398 | } | |
# | Line 401 | Line 416 | template<typename T> void ZConstraint<T>::update(){ | |
416 | massOfZConsMols.clear(); | |
417 | zPos.clear(); | |
418 | kz.clear(); | |
419 | + | cantPos.clear(); |
420 | + | cantVel.clear(); |
421 | ||
422 | unconsMols.clear(); | |
423 | massOfUnconsMols.clear(); | |
# | Line 414 | Line 431 | template<typename T> void ZConstraint<T>::update(){ | |
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()); |
434 | > | massOfZConsMols.push_back(molecules[i].getTotalMass()); |
435 | > | |
436 | > | if(usingSMD) |
437 | > | cantVel.push_back((*parameters)[index].cantVel); |
438 | ||
419 | – | molecules[i].getCOM(COM); |
439 | } | |
440 | else{ | |
441 | unconsMols.push_back(&molecules[i]); | |
# | Line 424 | Line 443 | template<typename T> void ZConstraint<T>::update(){ | |
443 | } | |
444 | } | |
445 | ||
446 | < | |
447 | < | //The reason to declare fz and indexOfZconsMols as pointer to array is |
448 | < | // that we want to make the MPI communication simple |
449 | < | if (fz){ |
450 | < | delete[] fz; |
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 | < | |
434 | < | if (curZPos){ |
435 | < | delete[] curZPos; |
436 | < | } |
437 | < | |
438 | < | if (indexOfZConsMols){ |
439 | < | delete[] indexOfZConsMols; |
440 | < | } |
441 | < | |
442 | < | if (zconsMols.size() > 0){ |
443 | < | fz = new double[zconsMols.size()]; |
444 | < | curZPos = new double[zconsMols.size()]; |
445 | < | indexOfZConsMols = new int[zconsMols.size()]; |
446 | < | |
447 | < | if (!fz || !curZPos || !indexOfZConsMols){ |
448 | < | sprintf(painCave.errMsg, |
449 | < | "Memory allocation failure in class Zconstraint\n"); |
450 | < | painCave.isFatal = 1; |
451 | < | simError(); |
452 | < | } |
453 | < | |
454 | < | for (int i = 0; i < (int) (zconsMols.size()); i++){ |
455 | < | indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); |
456 | < | } |
457 | < | } |
458 | < | else{ |
459 | < | fz = NULL; |
460 | < | curZPos = NULL; |
461 | < | indexOfZConsMols = NULL; |
462 | < | } |
463 | < | |
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 | ||
# | Line 476 | Line 468 | template<typename T> void ZConstraint<T>::update(){ | |
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 | // | |
480 | forcePolicy->update(); | |
481 | } | |
# | Line 580 | Line 579 | template<typename T> void ZConstraint<T>::calcForce(in | |
579 | this->doZconstraintForce(); | |
580 | } | |
581 | ||
582 | < | //use harmonical poteintial to move the molecules to the specified positions |
582 | > | //use external force to move the molecules to the specified positions |
583 | if (haveMovingZMols()){ | |
584 | < | this->doHarmonic(); |
584 | > | if (usingSMD) |
585 | > | this->doHarmonic(cantPos); |
586 | > | else |
587 | > | this->doHarmonic(zPos); |
588 | } | |
589 | ||
590 | //write out forces and current positions of z-constraint molecules | |
# | Line 602 | Line 604 | template<typename T> void ZConstraint<T>::calcForce(in | |
604 | } | |
605 | } | |
606 | } | |
607 | < | fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, |
608 | < | curZPos, &zPos[0]); |
607 | > | fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0], |
608 | > | &curZPos[0], &zPos[0]); |
609 | curZconsTime += zconsTime; | |
610 | } | |
611 | ||
# | Line 887 | Line 889 | template<typename T> void ZConstraint<T>::doZconstrain | |
889 | * | |
890 | */ | |
891 | ||
892 | < | template<typename T> void ZConstraint<T>::doHarmonic(){ |
892 | > | template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){ |
893 | double force[3]; | |
894 | double harmonicU; | |
895 | double harmonicF; | |
# | Line 908 | Line 910 | template<typename T> void ZConstraint<T>::doHarmonic() | |
910 | // cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] | |
911 | // << "\tcurrent zpos: " << COM[whichDirection] << endl; | |
912 | ||
913 | < | diff = COM[whichDirection] - zPos[i]; |
913 | > | diff = COM[whichDirection] - resPos[i]; |
914 | ||
915 | harmonicU = 0.5 * kz[i] * diff * diff; | |
916 | info->lrPot += harmonicU; | |
# | Line 978 | Line 980 | template<typename T> bool ZConstraint<T>::checkZConsSt | |
980 | if (diff <= zconsTol && states[i] == zcsMoving){ | |
981 | states[i] = zcsFixed; | |
982 | changed_local = 1; | |
983 | + | |
984 | + | if(usingSMD) |
985 | + | prevCantPos = cantPos; |
986 | ||
987 | if (hasZConsGap) | |
988 | endFixTime[i] = info->getTime() + zconsFixTime; | |
# | Line 986 | Line 991 | template<typename T> bool ZConstraint<T>::checkZConsSt | |
991 | states[i] = zcsMoving; | |
992 | changed_local = 1; | |
993 | ||
994 | + | if(usingSMD) |
995 | + | cantPos = prevCantPos; |
996 | + | |
997 | if (hasZConsGap) | |
998 | endFixTime[i] = INFINITE_TIME; | |
999 | } | |
# | Line 1249 | Line 1257 | template<typename T> void ZConstraint<T>::updateZPos() | |
1257 | ||
1258 | template<typename T> void ZConstraint<T>::updateZPos(){ | |
1259 | double curTime; | |
1260 | < | |
1260 | > | double COM[3]; |
1261 | > | |
1262 | curTime = info->getTime(); | |
1263 | ||
1264 | for (size_t i = 0; i < zconsMols.size(); i++){ | |
1265 | + | |
1266 | if (states[i] == zcsFixed && curTime >= endFixTime[i]){ | |
1267 | zPos[i] += zconsGap; | |
1268 | + | |
1269 | + | if (usingSMD){ |
1270 | + | zconsMols[i]->getCOM(COM); |
1271 | + | cantPos[i] = COM[whichDirection]; |
1272 | + | } |
1273 | + | |
1274 | } | |
1275 | + | |
1276 | } | |
1277 | + | |
1278 | } | |
1279 | + | |
1280 | + | template<typename T> void ZConstraint<T>::updateCantPos(){ |
1281 | + | double curTime; |
1282 | + | double dt; |
1283 | + | |
1284 | + | curTime = info->getTime(); |
1285 | + | dt = info->dt; |
1286 | + | |
1287 | + | for (size_t i = 0; i < zconsMols.size(); i++){ |
1288 | + | if (states[i] == zcsMoving){ |
1289 | + | cantPos[i] += cantVel[i] * dt; |
1290 | + | } |
1291 | + | } |
1292 | + | |
1293 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |