156 |
|
simError(); |
157 |
|
} |
158 |
|
|
159 |
< |
maxIndex = (*parameters)[parameters->size()].zconsIndex; |
159 |
> |
maxIndex = (*parameters)[parameters->size() - 1].zconsIndex; |
160 |
|
|
161 |
|
#ifndef IS_MPI |
162 |
|
totalNumMol = nMols; |
179 |
|
|
180 |
|
#ifndef IS_MPI |
181 |
|
for(int j = 0; j < nMols; j++){ |
182 |
< |
if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){ |
183 |
< |
molecules[i].getCOM(COM); |
182 |
> |
if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ |
183 |
> |
molecules[j].getCOM(COM); |
184 |
|
break; |
185 |
|
} |
186 |
|
} |
197 |
|
|
198 |
|
if (worldRank == whichNode ){ |
199 |
|
|
200 |
< |
for(int i = 0; i < nMols; i++) |
201 |
< |
if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){ |
202 |
< |
molecules[i].getCOM(COM); |
200 |
> |
for(int j = 0; j < nMols; j++) |
201 |
> |
if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ |
202 |
> |
molecules[j].getCOM(COM); |
203 |
|
break; |
204 |
|
} |
205 |
|
|
261 |
|
simError(); |
262 |
|
} |
263 |
|
|
264 |
< |
for(int i = 0; i < zconsMols.size(); i++) |
264 |
> |
//determine the states of z-constraint molecules |
265 |
> |
for(int i = 0; i < zconsMols.size(); i++){ |
266 |
|
indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); |
267 |
+ |
|
268 |
+ |
zconsMols[i]->getCOM(COM); |
269 |
+ |
if (fabs(zPos[i] - COM[whichDirection]) < zconsTol) |
270 |
+ |
states.push_back(zcsFixed); |
271 |
+ |
else |
272 |
+ |
states.push_back(zcsMoving); |
273 |
+ |
} |
274 |
|
|
275 |
|
#endif |
276 |
+ |
|
277 |
+ |
//get total masss of unconstraint molecules |
278 |
+ |
double totalMassOfUncons_local; |
279 |
+ |
totalMassOfUncons_local = 0; |
280 |
+ |
|
281 |
+ |
for(int i = 0; i < unconsMols.size(); i++) |
282 |
+ |
totalMassOfUncons_local += unconsMols[i]->getTotalMass(); |
283 |
+ |
|
284 |
+ |
#ifndef IS_MPI |
285 |
+ |
totalMassOfUncons = totalMassOfUncons_local; |
286 |
+ |
#else |
287 |
+ |
MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
288 |
+ |
#endif |
289 |
|
|
290 |
+ |
|
291 |
|
//get total number of unconstrained atoms |
292 |
|
int nUnconsAtoms_local; |
293 |
|
nUnconsAtoms_local = 0; |
298 |
|
totNumOfUnconsAtoms = nUnconsAtoms_local; |
299 |
|
#else |
300 |
|
MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
301 |
< |
#endif |
301 |
> |
#endif |
302 |
|
|
303 |
< |
checkZConsState(); |
282 |
< |
|
283 |
< |
// |
303 |
> |
// creat zconsWriter |
304 |
|
fzOut = new ZConsWriter(zconsOutput.c_str()); |
305 |
|
|
306 |
|
if(!fzOut){ |
348 |
|
|
349 |
|
zconsMols.push_back(&molecules[i]); |
350 |
|
zPos.push_back((*parameters)[index].zPos); |
351 |
< |
kz.push_back((*parameters)[index].kRatio * zForceConst); |
352 |
< |
|
351 |
> |
kz.push_back((*parameters)[index].kRatio * zForceConst); |
352 |
> |
|
353 |
|
massOfZConsMols.push_back(molecules[i].getTotalMass()); |
354 |
|
|
355 |
|
molecules[i].getCOM(COM); |
361 |
|
massOfUnconsMols.push_back(molecules[i].getTotalMass()); |
362 |
|
|
363 |
|
} |
364 |
+ |
} |
365 |
+ |
|
366 |
+ |
//determine the states of z-constraint molecules |
367 |
+ |
for(int i = 0; i < zconsMols.size(); i++){ |
368 |
+ |
zconsMols[i]->getCOM(COM); |
369 |
+ |
if (fabs(zPos[i] - COM[whichDirection]) < zconsTol) |
370 |
+ |
states.push_back(zcsFixed); |
371 |
+ |
else |
372 |
+ |
states.push_back(zcsMoving); |
373 |
|
} |
374 |
+ |
|
375 |
|
|
376 |
|
//The reason to declare fz and indexOfZconsMols as pointer to array is |
377 |
|
// that we want to make the MPI communication simple |
440 |
|
return -1; |
441 |
|
} |
442 |
|
|
413 |
– |
/** |
414 |
– |
* Description: |
415 |
– |
* Reset the z coordinates |
416 |
– |
*/ |
443 |
|
template<typename T> void ZConstraint<T>::integrate(){ |
444 |
|
|
445 |
|
//zero out the velocities of center of mass of unconstrained molecules |
456 |
|
* |
457 |
|
* |
458 |
|
* |
459 |
< |
*/ |
434 |
< |
|
435 |
< |
|
459 |
> |
*/ |
460 |
|
template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){ |
461 |
+ |
double zsys; |
462 |
|
|
463 |
|
T::calcForce(calcPot, calcStress); |
464 |
|
|
465 |
|
if (checkZConsState()) |
466 |
|
zeroOutVel(); |
467 |
+ |
|
468 |
+ |
zsys = calcZSys(); |
469 |
+ |
cout << "---------------------------------------------------------------------" <<endl; |
470 |
+ |
cout << "current time: " << info->getTime() << endl; |
471 |
+ |
cout << "center of mass at z: " << zsys << endl; |
472 |
+ |
//cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; |
473 |
+ |
cout << "before calcForce, the COMVel of system is " << calcSysCOMVel() <<endl; |
474 |
|
|
475 |
+ |
//cout << "before doZConstraintForce, totalForce is " << calcTotalForce() << endl; |
476 |
+ |
|
477 |
|
//do zconstraint force; |
478 |
|
if (haveFixedZMols()) |
479 |
|
this->doZconstraintForce(); |
480 |
< |
|
480 |
> |
|
481 |
|
//use harmonical poteintial to move the molecules to the specified positions |
482 |
|
if (haveMovingZMols()) |
483 |
< |
//this->doHarmonic(); |
484 |
< |
|
483 |
> |
this->doHarmonic(); |
484 |
> |
|
485 |
> |
//cout << "after doHarmonic, totalForce is " << calcTotalForce() << endl; |
486 |
> |
|
487 |
|
fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); |
488 |
< |
|
488 |
> |
//cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; |
489 |
> |
cout << "after calcForce, the COMVel of system is " << calcSysCOMVel() <<endl; |
490 |
|
} |
491 |
|
|
492 |
|
template<typename T> double ZConstraint<T>::calcZSys() |
496 |
|
double totalMass; |
497 |
|
double totalMZ_local; |
498 |
|
double totalMZ; |
462 |
– |
double massOfUncons_local; |
499 |
|
double massOfCurMol; |
500 |
|
double COM[3]; |
501 |
< |
|
501 |
> |
|
502 |
|
totalMass_local = 0; |
467 |
– |
totalMass = 0; |
503 |
|
totalMZ_local = 0; |
469 |
– |
totalMZ = 0; |
470 |
– |
massOfUncons_local = 0; |
471 |
– |
|
504 |
|
|
505 |
|
for(int i = 0; i < nMols; i++){ |
506 |
|
massOfCurMol = molecules[i].getTotalMass(); |
508 |
|
|
509 |
|
totalMass_local += massOfCurMol; |
510 |
|
totalMZ_local += massOfCurMol * COM[whichDirection]; |
511 |
< |
|
480 |
< |
if(isZConstraintMol(&molecules[i]) == -1){ |
481 |
< |
|
482 |
< |
massOfUncons_local += massOfCurMol; |
483 |
< |
} |
484 |
< |
|
511 |
> |
|
512 |
|
} |
513 |
+ |
|
514 |
|
|
487 |
– |
|
515 |
|
#ifdef IS_MPI |
516 |
|
MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
517 |
|
MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
518 |
< |
MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
492 |
< |
#else |
518 |
> |
#else |
519 |
|
totalMass = totalMass_local; |
520 |
|
totalMZ = totalMZ_local; |
521 |
< |
totalMassOfUncons = massOfUncons_local; |
496 |
< |
#endif |
521 |
> |
#endif |
522 |
|
|
523 |
|
double zsys; |
524 |
|
zsys = totalMZ / totalMass; |
546 |
|
Atom** fixedZAtoms; |
547 |
|
double COMvel[3]; |
548 |
|
double vel[3]; |
549 |
< |
|
549 |
> |
|
550 |
|
//zero out the velocities of center of mass of fixed z-constrained molecules |
551 |
|
|
552 |
|
for(int i = 0; i < zconsMols.size(); i++){ |
553 |
|
|
554 |
|
if (states[i] == zcsFixed){ |
555 |
|
|
556 |
< |
zconsMols[i]->getCOMvel(COMvel); |
556 |
> |
zconsMols[i]->getCOMvel(COMvel); |
557 |
> |
//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
558 |
> |
|
559 |
|
fixedZAtoms = zconsMols[i]->getMyAtoms(); |
560 |
|
|
561 |
|
for(int j =0; j < zconsMols[i]->getNAtoms(); j++){ |
562 |
|
fixedZAtoms[j]->getVel(vel); |
563 |
< |
vel[whichDirection] -= COMvel[whichDirection]; |
564 |
< |
fixedZAtoms[j]->setVel(vel); |
563 |
> |
vel[whichDirection] -= COMvel[whichDirection]; |
564 |
> |
fixedZAtoms[j]->setVel(vel); |
565 |
|
} |
566 |
< |
|
566 |
> |
|
567 |
> |
zconsMols[i]->getCOMvel(COMvel); |
568 |
> |
//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
569 |
|
} |
570 |
|
|
571 |
|
} |
572 |
< |
|
572 |
> |
|
573 |
> |
//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; |
574 |
> |
cout << "before resetting the COMVel of sytem is " << calcSysCOMVel() << endl; |
575 |
> |
|
576 |
|
// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules |
577 |
|
double MVzOfMovingMols_local; |
578 |
|
double MVzOfMovingMols; |
587 |
|
MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
588 |
|
} |
589 |
|
|
590 |
< |
for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){ |
559 |
< |
|
590 |
> |
for(int i = 0; i < zconsMols.size(); i++){ |
591 |
|
if (states[i] == zcsMoving){ |
592 |
|
zconsMols[i]->getCOMvel(COMvel); |
593 |
|
MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; |
622 |
|
|
623 |
|
//modify the velocities of moving z-constrained molecuels |
624 |
|
Atom** movingZAtoms; |
625 |
< |
for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){ |
625 |
> |
for(int i = 0; i < zconsMols.size(); i++){ |
626 |
|
|
627 |
|
if (states[i] ==zcsMoving){ |
628 |
|
|
629 |
|
movingZAtoms = zconsMols[i]->getMyAtoms(); |
630 |
< |
for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
630 |
> |
for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
631 |
|
movingZAtoms[j]->getVel(vel); |
632 |
|
vel[whichDirection] -= vzOfMovingMols; |
633 |
< |
movingZAtoms[j]->setVel(vel); |
634 |
< |
} |
633 |
> |
movingZAtoms[j]->setVel(vel); |
634 |
> |
} |
635 |
|
|
636 |
< |
} |
636 |
> |
} |
637 |
|
|
638 |
< |
} |
638 |
> |
} |
639 |
|
|
640 |
+ |
cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() <<endl; |
641 |
+ |
|
642 |
|
} |
643 |
|
|
644 |
|
template<typename T> void ZConstraint<T>::doZconstraintForce(){ |
649 |
|
double COMvel[3]; |
650 |
|
double COM[3]; |
651 |
|
double force[3]; |
619 |
– |
double zsys; |
652 |
|
|
621 |
– |
int nMovingZMols_local; |
622 |
– |
int nMovingZMols; |
653 |
|
|
624 |
– |
//constrain the molecules which do not reach the specified positions |
654 |
|
|
655 |
< |
zsys = calcZSys(); |
627 |
< |
cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl; |
655 |
> |
//constrain the molecules which do not reach the specified positions |
656 |
|
|
657 |
|
//Zero Out the force of z-contrained molecules |
658 |
|
totalFZ_local = 0; |
673 |
|
} |
674 |
|
totalFZ_local += fz[i]; |
675 |
|
|
676 |
< |
cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; |
676 |
> |
cout << "index: " << indexOfZConsMols[i] |
677 |
> |
<<"\tcurrent zpos: " << COM[whichDirection] |
678 |
> |
<< "\tcurrent fz: " <<fz[i] << endl; |
679 |
|
|
680 |
|
} |
681 |
|
|
682 |
|
} |
683 |
+ |
|
684 |
+ |
// apply negative to fixed z-constrained molecues; |
685 |
+ |
force[0]= 0; |
686 |
+ |
force[1]= 0; |
687 |
+ |
force[2]= 0; |
688 |
|
|
689 |
+ |
for(int i = 0; i < zconsMols.size(); i++){ |
690 |
+ |
|
691 |
+ |
if (states[i] == zcsFixed){ |
692 |
+ |
|
693 |
+ |
int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms(); |
694 |
+ |
zconsAtoms = zconsMols[i]->getMyAtoms(); |
695 |
+ |
|
696 |
+ |
for(int j =0; j < nAtomOfCurZConsMol; j++) { |
697 |
+ |
force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; |
698 |
+ |
zconsAtoms[j]->addFrc(force); |
699 |
+ |
} |
700 |
+ |
|
701 |
+ |
} |
702 |
+ |
|
703 |
+ |
} |
704 |
+ |
|
705 |
+ |
cout << "after zero out z-constraint force on fixed z-constraint molecuels " |
706 |
+ |
<< "total force is " << calcTotalForce() << endl; |
707 |
|
//calculate the number of atoms of moving z-constrained molecules |
708 |
< |
nMovingZMols_local = 0; |
709 |
< |
for(int i = 0; zconsMols.size(); i++){ |
708 |
> |
int nMovingZAtoms_local; |
709 |
> |
int nMovingZAtoms; |
710 |
> |
|
711 |
> |
nMovingZAtoms_local = 0; |
712 |
> |
for(int i = 0; i < zconsMols.size(); i++) |
713 |
|
if(states[i] == zcsMoving) |
714 |
< |
nMovingZMols_local += massOfZConsMols[i]; |
715 |
< |
} |
714 |
> |
nMovingZAtoms_local += zconsMols[i]->getNAtoms(); |
715 |
> |
|
716 |
|
#ifdef IS_MPI |
717 |
|
MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
718 |
< |
MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
718 |
> |
MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
719 |
|
#else |
720 |
|
totalFZ = totalFZ_local; |
721 |
< |
nMovingZMols = nMovingZMols_local; |
721 |
> |
nMovingZAtoms = nMovingZAtoms_local; |
722 |
|
#endif |
723 |
|
|
724 |
|
force[0]= 0; |
725 |
|
force[1]= 0; |
726 |
|
force[2]= 0; |
727 |
< |
force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols); |
727 |
> |
force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); |
728 |
|
|
729 |
< |
//modify the velocites of unconstrained molecules |
729 |
> |
//modify the forces of unconstrained molecules |
730 |
> |
int accessCount = 0; |
731 |
|
for(int i = 0; i < unconsMols.size(); i++){ |
732 |
|
|
733 |
|
Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
737 |
|
|
738 |
|
} |
739 |
|
|
740 |
< |
//modify the velocities of moving z-constrained molecules |
740 |
> |
//modify the forces of moving z-constrained molecules |
741 |
|
for(int i = 0; i < zconsMols.size(); i++) { |
742 |
< |
if (states[i] == zcsMoving){ |
742 |
> |
if (states[i] == zcsMoving){ |
743 |
|
|
744 |
< |
Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
744 |
> |
Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
745 |
|
|
746 |
< |
for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) |
747 |
< |
movingZAtoms[j]->addFrc(force); |
748 |
< |
} |
746 |
> |
for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) |
747 |
> |
movingZAtoms[j]->addFrc(force); |
748 |
> |
} |
749 |
|
} |
693 |
– |
|
694 |
– |
// apply negative to fixed z-constrained molecues; |
695 |
– |
force[0]= 0; |
696 |
– |
force[1]= 0; |
697 |
– |
force[2]= 0; |
698 |
– |
|
699 |
– |
for(int i = 0; i < zconsMols.size(); i++){ |
700 |
– |
|
701 |
– |
if (states[i] == zcsFixed){ |
750 |
|
|
751 |
< |
int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms(); |
752 |
< |
zconsAtoms = zconsMols[i]->getMyAtoms(); |
705 |
< |
|
706 |
< |
for(int j =0; j < nAtomOfCurZConsMol; j++) { |
707 |
< |
force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; |
708 |
< |
zconsAtoms[j]->addFrc(force); |
709 |
< |
} |
710 |
< |
|
711 |
< |
} |
712 |
< |
|
713 |
< |
} |
751 |
> |
cout << "after substracting z-constraint force from moving molecuels " |
752 |
> |
<< "total force is " << calcTotalForce() << endl; |
753 |
|
|
754 |
|
} |
755 |
|
|
808 |
|
template<typename T> void ZConstraint<T>::doHarmonic(){ |
809 |
|
double force[3]; |
810 |
|
double harmonicU; |
811 |
+ |
double harmonicF; |
812 |
|
double COM[3]; |
813 |
|
double diff; |
814 |
+ |
double totalFZ; |
815 |
|
|
816 |
|
force[0] = 0; |
817 |
|
force[1] = 0; |
818 |
|
force[2] = 0; |
819 |
|
|
820 |
+ |
totalFZ = 0; |
821 |
+ |
|
822 |
|
cout << "Moving Molecules" << endl; |
823 |
|
for(int i = 0; i < zconsMols.size(); i++) { |
824 |
|
|
825 |
|
if (states[i] == zcsMoving){ |
826 |
< |
zconsMols[i]->getCOM(COM): |
826 |
> |
zconsMols[i]->getCOM(COM); |
827 |
|
cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; |
828 |
|
|
829 |
|
diff = COM[whichDirection] -zPos[i]; |
830 |
|
|
831 |
|
harmonicU = 0.5 * kz[i] * diff * diff; |
832 |
< |
info->ltPot += harmonicU; |
832 |
> |
info->lrPot += harmonicU; |
833 |
|
|
834 |
< |
force[whichDirection] = - kz[i] * diff / zconsMols[i]->getNAtoms(); |
834 |
> |
harmonicF = - kz[i] * diff; |
835 |
> |
totalFZ += harmonicF; |
836 |
> |
|
837 |
> |
// |
838 |
> |
force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms(); |
839 |
|
|
840 |
|
Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
841 |
|
|
844 |
|
} |
845 |
|
|
846 |
|
} |
847 |
+ |
|
848 |
+ |
force[0]= 0; |
849 |
+ |
force[1]= 0; |
850 |
+ |
force[2]= 0; |
851 |
+ |
force[whichDirection] = -totalFZ /totNumOfUnconsAtoms; |
852 |
|
|
853 |
+ |
//modify the forces of unconstrained molecules |
854 |
+ |
for(int i = 0; i < unconsMols.size(); i++){ |
855 |
+ |
|
856 |
+ |
Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
857 |
+ |
|
858 |
+ |
for(int j = 0; j < unconsMols[i]->getNAtoms(); j++) |
859 |
+ |
unconsAtoms[j]->addFrc(force); |
860 |
+ |
} |
861 |
+ |
|
862 |
|
} |
863 |
+ |
|
864 |
+ |
template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel() |
865 |
+ |
{ |
866 |
+ |
double MVzOfMovingMols_local; |
867 |
+ |
double MVzOfMovingMols; |
868 |
+ |
double totalMassOfMovingZMols_local; |
869 |
+ |
double totalMassOfMovingZMols; |
870 |
+ |
double COMvel[3]; |
871 |
+ |
|
872 |
+ |
MVzOfMovingMols_local = 0; |
873 |
+ |
totalMassOfMovingZMols_local = 0; |
874 |
+ |
|
875 |
+ |
for(int i =0; i < unconsMols.size(); i++){ |
876 |
+ |
unconsMols[i]->getCOMvel(COMvel); |
877 |
+ |
MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
878 |
+ |
} |
879 |
+ |
|
880 |
+ |
for(int i = 0; i < zconsMols.size(); i++){ |
881 |
+ |
|
882 |
+ |
if (states[i] == zcsMoving){ |
883 |
+ |
zconsMols[i]->getCOMvel(COMvel); |
884 |
+ |
MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; |
885 |
+ |
totalMassOfMovingZMols_local += massOfZConsMols[i]; |
886 |
+ |
} |
887 |
+ |
|
888 |
+ |
} |
889 |
+ |
|
890 |
+ |
#ifndef IS_MPI |
891 |
+ |
MVzOfMovingMols = MVzOfMovingMols_local; |
892 |
+ |
totalMassOfMovingZMols = totalMassOfMovingZMols_local; |
893 |
+ |
#else |
894 |
+ |
MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
895 |
+ |
MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
896 |
+ |
#endif |
897 |
+ |
|
898 |
+ |
double vzOfMovingMols; |
899 |
+ |
vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols); |
900 |
+ |
|
901 |
+ |
return vzOfMovingMols; |
902 |
+ |
} |
903 |
+ |
|
904 |
+ |
|
905 |
+ |
template<typename T> double ZConstraint<T>::calcSysCOMVel() |
906 |
+ |
{ |
907 |
+ |
double COMvel[3]; |
908 |
+ |
double tempMVz = 0; |
909 |
+ |
|
910 |
+ |
for(int i =0 ; i < nMols; i++){ |
911 |
+ |
molecules[i].getCOMvel(COMvel); |
912 |
+ |
tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection]; |
913 |
+ |
} |
914 |
+ |
|
915 |
+ |
double massOfZCons_local; |
916 |
+ |
double massOfZCons; |
917 |
+ |
|
918 |
+ |
massOfZCons_local = 0; |
919 |
+ |
|
920 |
+ |
for(int i = 0; i < massOfZConsMols.size(); i++){ |
921 |
+ |
massOfZCons_local += massOfZConsMols[i]; |
922 |
+ |
} |
923 |
+ |
#ifndef IS_MPI |
924 |
+ |
massOfZCons = massOfZCons_local; |
925 |
+ |
#else |
926 |
+ |
MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
927 |
+ |
#endif |
928 |
+ |
|
929 |
+ |
return tempMVz /(totalMassOfUncons + massOfZCons); |
930 |
+ |
} |
931 |
+ |
|
932 |
+ |
template<typename T> double ZConstraint<T>::calcTotalForce(){ |
933 |
+ |
|
934 |
+ |
double force[3]; |
935 |
+ |
double totalForce_local; |
936 |
+ |
double totalForce; |
937 |
+ |
|
938 |
+ |
totalForce_local = 0; |
939 |
+ |
|
940 |
+ |
for(int i = 0; i < nAtoms; i++){ |
941 |
+ |
atoms[i]->getFrc(force); |
942 |
+ |
totalForce_local += force[whichDirection]; |
943 |
+ |
} |
944 |
+ |
|
945 |
+ |
#ifndef IS_MPI |
946 |
+ |
totalForce = totalForce_local; |
947 |
+ |
#else |
948 |
+ |
MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
949 |
+ |
#endif |
950 |
+ |
|
951 |
+ |
return totalForce; |
952 |
+ |
|
953 |
+ |
} |