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 |
|