456 |
|
* |
457 |
|
* |
458 |
|
* |
459 |
< |
*/ |
460 |
< |
|
461 |
< |
|
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(); |
466 |
> |
zeroOutVel(); |
467 |
|
|
468 |
|
zsys = calcZSys(); |
469 |
|
cout << "---------------------------------------------------------------------" <<endl; |
470 |
< |
cout << "current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl; |
471 |
< |
cout << "before calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl; |
472 |
< |
|
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 |
< |
|
480 |
> |
|
481 |
|
//use harmonical poteintial to move the molecules to the specified positions |
482 |
|
if (haveMovingZMols()) |
483 |
|
this->doHarmonic(); |
484 |
< |
|
484 |
> |
|
485 |
> |
//cout << "after doHarmonic, totalForce is " << calcTotalForce() << endl; |
486 |
> |
|
487 |
|
fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); |
488 |
< |
cout << "after calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl; |
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; |
497 |
– |
double massOfUncons_local; |
499 |
|
double massOfCurMol; |
500 |
|
double COM[3]; |
501 |
< |
|
501 |
> |
|
502 |
|
totalMass_local = 0; |
502 |
– |
totalMass = 0; |
503 |
|
totalMZ_local = 0; |
504 |
– |
totalMZ = 0; |
505 |
– |
massOfUncons_local = 0; |
506 |
– |
|
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 |
< |
|
515 |
< |
if(isZConstraintMol(&molecules[i]) == -1){ |
516 |
< |
|
517 |
< |
massOfUncons_local += massOfCurMol; |
518 |
< |
} |
519 |
< |
|
511 |
> |
|
512 |
|
} |
513 |
+ |
|
514 |
|
|
522 |
– |
|
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); |
527 |
< |
#else |
518 |
> |
#else |
519 |
|
totalMass = totalMass_local; |
520 |
|
totalMZ = totalMZ_local; |
521 |
< |
totalMassOfUncons = massOfUncons_local; |
531 |
< |
#endif |
521 |
> |
#endif |
522 |
|
|
523 |
|
double zsys; |
524 |
|
zsys = totalMZ / totalMass; |
547 |
|
double COMvel[3]; |
548 |
|
double vel[3]; |
549 |
|
|
560 |
– |
double tempMass = 0; |
561 |
– |
double tempMVz =0; |
562 |
– |
double tempVz = 0; |
563 |
– |
for(int i = 0; i < nMols; i++){ |
564 |
– |
molecules[i].getCOMvel(COMvel); |
565 |
– |
tempMass +=molecules[i].getTotalMass(); |
566 |
– |
tempMVz += molecules[i].getTotalMass() * COMvel[whichDirection]; |
567 |
– |
tempVz += COMvel[whichDirection]; |
568 |
– |
} |
569 |
– |
cout << "initial velocity of center of mass is " << tempMVz / tempMass << endl; |
570 |
– |
|
550 |
|
//zero out the velocities of center of mass of fixed z-constrained molecules |
551 |
|
|
552 |
|
for(int i = 0; i < zconsMols.size(); i++){ |
554 |
|
if (states[i] == zcsFixed){ |
555 |
|
|
556 |
|
zconsMols[i]->getCOMvel(COMvel); |
557 |
< |
cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
557 |
> |
//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
558 |
|
|
559 |
|
fixedZAtoms = zconsMols[i]->getMyAtoms(); |
560 |
|
|
565 |
|
} |
566 |
|
|
567 |
|
zconsMols[i]->getCOMvel(COMvel); |
568 |
< |
cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
568 |
> |
//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
569 |
|
} |
570 |
|
|
571 |
|
} |
572 |
|
|
573 |
< |
cout << "before resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl; |
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; |
637 |
|
|
638 |
|
} |
639 |
|
|
640 |
< |
cout << "after resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl; |
640 |
> |
cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() <<endl; |
641 |
|
|
642 |
|
} |
643 |
|
|
650 |
|
double COM[3]; |
651 |
|
double force[3]; |
652 |
|
|
673 |
– |
int nMovingZMols_local; |
674 |
– |
int nMovingZMols; |
653 |
|
|
654 |
+ |
|
655 |
|
//constrain the molecules which do not reach the specified positions |
656 |
|
|
657 |
|
//Zero Out the force of z-contrained molecules |
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; |
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]; |
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 forces of unconstrained molecules |
730 |
+ |
int accessCount = 0; |
731 |
|
for(int i = 0; i < unconsMols.size(); i++){ |
732 |
|
|
733 |
|
Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
739 |
|
|
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 |
|
} |
742 |
– |
|
743 |
– |
// apply negative to fixed z-constrained molecues; |
744 |
– |
force[0]= 0; |
745 |
– |
force[1]= 0; |
746 |
– |
force[2]= 0; |
747 |
– |
|
748 |
– |
for(int i = 0; i < zconsMols.size(); i++){ |
749 |
– |
|
750 |
– |
if (states[i] == zcsFixed){ |
750 |
|
|
751 |
< |
int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms(); |
752 |
< |
zconsAtoms = zconsMols[i]->getMyAtoms(); |
754 |
< |
|
755 |
< |
for(int j =0; j < nAtomOfCurZConsMol; j++) { |
756 |
< |
force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; |
757 |
< |
zconsAtoms[j]->addFrc(force); |
758 |
< |
} |
759 |
< |
|
760 |
< |
} |
761 |
< |
|
762 |
< |
} |
751 |
> |
cout << "after substracting z-constraint force from moving molecuels " |
752 |
> |
<< "total force is " << calcTotalForce() << endl; |
753 |
|
|
754 |
|
} |
755 |
|
|
831 |
|
harmonicU = 0.5 * kz[i] * diff * diff; |
832 |
|
info->lrPot += harmonicU; |
833 |
|
|
834 |
< |
harmonicF = - kz[i] * diff / zconsMols[i]->getNAtoms(); |
845 |
< |
force[whichDirection] = harmonicF; |
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 |
< |
|
847 |
> |
|
848 |
|
force[0]= 0; |
849 |
|
force[1]= 0; |
850 |
|
force[2]= 0; |
861 |
|
|
862 |
|
} |
863 |
|
|
864 |
< |
template<typename T> double ZConstraint<T>::calcCOMVel() |
864 |
> |
template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel() |
865 |
|
{ |
866 |
|
double MVzOfMovingMols_local; |
867 |
|
double MVzOfMovingMols; |
902 |
|
} |
903 |
|
|
904 |
|
|
905 |
< |
template<typename T> double ZConstraint<T>::calcCOMVel2() |
905 |
> |
template<typename T> double ZConstraint<T>::calcSysCOMVel() |
906 |
|
{ |
907 |
|
double COMvel[3]; |
908 |
|
double tempMVz = 0; |
917 |
– |
int index; |
909 |
|
|
910 |
|
for(int i =0 ; i < nMols; i++){ |
911 |
< |
index = isZConstraintMol(&molecules[i]); |
912 |
< |
if( index == -1){ |
922 |
< |
molecules[i].getCOMvel(COMvel); |
923 |
< |
tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection]; |
924 |
< |
} |
925 |
< |
else if(states[i] == zcsMoving){ |
926 |
< |
zconsMols[index]->getCOMvel(COMvel); |
927 |
< |
tempMVz += massOfZConsMols[index]*COMvel[whichDirection]; |
928 |
< |
} |
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 |
< |
return tempMVz /totalMassOfUncons; |
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 |
+ |
} |