476 |
|
if(usingSMD) |
477 |
|
prevCantPos = cantPos; |
478 |
|
|
479 |
– |
// |
479 |
|
forcePolicy->update(); |
480 |
|
} |
481 |
|
|
534 |
|
T::integrate(); |
535 |
|
} |
536 |
|
|
538 |
– |
|
539 |
– |
/** |
540 |
– |
* |
541 |
– |
* |
542 |
– |
* |
543 |
– |
* |
544 |
– |
*/ |
537 |
|
template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){ |
538 |
|
double zsys; |
539 |
|
double COM[3]; |
557 |
|
#ifdef IS_MPI |
558 |
|
if (worldRank == 0){ |
559 |
|
#endif |
568 |
– |
//cout << "---------------------------------------------------------------------" <<endl; |
569 |
– |
//cout << "current time: " << info->getTime() << endl; |
570 |
– |
//cout << "center of mass at z: " << zsys << endl; |
571 |
– |
//cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl; |
560 |
|
|
561 |
|
#ifdef IS_MPI |
562 |
|
} |
601 |
|
#ifdef IS_MPI |
602 |
|
if (worldRank == 0){ |
603 |
|
#endif |
616 |
– |
//cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl; |
604 |
|
#ifdef IS_MPI |
605 |
|
} |
606 |
|
#endif |
607 |
|
} |
608 |
|
|
609 |
|
|
623 |
– |
/** |
624 |
– |
* |
625 |
– |
*/ |
626 |
– |
|
610 |
|
template<typename T> double ZConstraint<T>::calcZSys(){ |
611 |
|
//calculate reference z coordinate for z-constraint molecules |
612 |
|
double totalMass_local; |
643 |
|
return zsys; |
644 |
|
} |
645 |
|
|
663 |
– |
/** |
664 |
– |
* |
665 |
– |
*/ |
646 |
|
template<typename T> void ZConstraint<T>::thermalize(void){ |
647 |
|
T::thermalize(); |
648 |
|
zeroOutVel(); |
649 |
|
} |
650 |
|
|
671 |
– |
/** |
672 |
– |
* |
673 |
– |
*/ |
674 |
– |
|
651 |
|
template<typename T> void ZConstraint<T>::zeroOutVel(){ |
652 |
|
Atom** fixedZAtoms; |
653 |
|
double COMvel[3]; |
670 |
|
} |
671 |
|
|
672 |
|
zconsMols[i]->getCOMvel(COMvel); |
697 |
– |
//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; |
673 |
|
} |
674 |
|
} |
675 |
|
|
701 |
– |
//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; |
702 |
– |
|
676 |
|
zSysCOMVel = calcSysCOMVel(); |
677 |
|
#ifdef IS_MPI |
678 |
|
if (worldRank == 0){ |
679 |
|
#endif |
707 |
– |
//cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl; |
680 |
|
#ifdef IS_MPI |
681 |
|
} |
682 |
|
#endif |
746 |
|
#ifdef IS_MPI |
747 |
|
if (worldRank == 0){ |
748 |
|
#endif |
777 |
– |
//cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl; |
749 |
|
#ifdef IS_MPI |
750 |
|
} |
751 |
|
#endif |
752 |
|
} |
753 |
|
|
783 |
– |
/** |
784 |
– |
* |
785 |
– |
*/ |
754 |
|
|
755 |
|
template<typename T> void ZConstraint<T>::doZconstraintForce(){ |
756 |
|
Atom** zconsAtoms; |
766 |
|
|
767 |
|
//calculate the total z-contrained force of fixed z-contrained molecules |
768 |
|
|
801 |
– |
//cout << "before zero out z-constraint force on fixed z-constraint molecuels " |
802 |
– |
// << "total force is " << calcTotalForce() << endl; |
803 |
– |
|
769 |
|
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
770 |
|
if (states[i] == zcsFixed){ |
771 |
|
zconsMols[i]->getCOM(COM); |
778 |
|
} |
779 |
|
totalFZ_local += fz[i]; |
780 |
|
|
816 |
– |
//cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i] |
817 |
– |
// <<"\tcurrent zpos: " << COM[whichDirection] |
818 |
– |
// << "\tcurrent fz: " <<fz[i] << endl; |
781 |
|
} |
782 |
|
} |
783 |
|
|
809 |
|
} |
810 |
|
} |
811 |
|
|
850 |
– |
//cout << "after zero out z-constraint force on fixed z-constraint molecuels " |
851 |
– |
// << "total force is " << calcTotalForce() << endl; |
852 |
– |
|
853 |
– |
|
812 |
|
force[0] = 0; |
813 |
|
force[1] = 0; |
814 |
|
force[2] = 0; |
838 |
|
} |
839 |
|
} |
840 |
|
} |
883 |
– |
// cout << "after substracting z-constraint force from moving molecuels " |
884 |
– |
// << "total force is " << calcTotalForce() << endl; |
841 |
|
} |
842 |
|
|
887 |
– |
/** |
888 |
– |
* |
889 |
– |
* |
890 |
– |
*/ |
843 |
|
|
844 |
|
template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){ |
845 |
|
double force[3]; |
859 |
|
for (int i = 0; i < (int) (zconsMols.size()); i++){ |
860 |
|
if (states[i] == zcsMoving){ |
861 |
|
zconsMols[i]->getCOM(COM); |
910 |
– |
// cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] |
911 |
– |
// << "\tcurrent zpos: " << COM[whichDirection] << endl; |
862 |
|
|
863 |
|
diff = COM[whichDirection] - resPos[i]; |
864 |
|
|
873 |
|
Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
874 |
|
|
875 |
|
for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
926 |
– |
//force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms(); |
876 |
|
force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], |
877 |
|
movingZAtoms[j], |
878 |
|
harmonicF); |
887 |
|
MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); |
888 |
|
#endif |
889 |
|
|
941 |
– |
//cout << "before substracting harmonic force from moving molecuels " |
942 |
– |
// << "total force is " << calcTotalForce() << endl; |
943 |
– |
|
890 |
|
force[0] = 0; |
891 |
|
force[1] = 0; |
892 |
|
force[2] = 0; |
903 |
|
} |
904 |
|
} |
905 |
|
|
960 |
– |
//cout << "after substracting harmonic force from moving molecuels " |
961 |
– |
// << "total force is " << calcTotalForce() << endl; |
906 |
|
} |
907 |
|
|
964 |
– |
/** |
965 |
– |
* |
966 |
– |
*/ |
967 |
– |
|
908 |
|
template<typename T> bool ZConstraint<T>::checkZConsState(){ |
909 |
|
double COM[3]; |
910 |
|
double diff; |
971 |
|
} |
972 |
|
|
973 |
|
|
1034 |
– |
/** |
1035 |
– |
* |
1036 |
– |
*/ |
974 |
|
template<typename T> bool ZConstraint<T>::haveMovingZMols(){ |
975 |
|
int havingMoving_local; |
976 |
|
int havingMoving; |
993 |
|
return (havingMoving > 0); |
994 |
|
} |
995 |
|
|
1059 |
– |
/** |
1060 |
– |
* |
1061 |
– |
*/ |
996 |
|
|
997 |
|
template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){ |
998 |
|
double MVzOfMovingMols_local; |
1034 |
|
return vzOfMovingMols; |
1035 |
|
} |
1036 |
|
|
1103 |
– |
/** |
1104 |
– |
* |
1105 |
– |
*/ |
1106 |
– |
|
1037 |
|
template<typename T> double ZConstraint<T>::calcSysCOMVel(){ |
1038 |
|
double COMvel[3]; |
1039 |
|
double tempMVz_local; |
1066 |
|
return tempMVz / (totalMassOfUncons + massOfZCons); |
1067 |
|
} |
1068 |
|
|
1139 |
– |
/** |
1140 |
– |
* |
1141 |
– |
*/ |
1142 |
– |
|
1069 |
|
template<typename T> double ZConstraint<T>::calcTotalForce(){ |
1070 |
|
double force[3]; |
1071 |
|
double totalForce_local; |
1088 |
|
return totalForce; |
1089 |
|
} |
1090 |
|
|
1165 |
– |
/** |
1166 |
– |
* |
1167 |
– |
*/ |
1168 |
– |
|
1091 |
|
template<typename T> void ZConstraint<T>::PolicyByNumber::update(){ |
1092 |
|
//calculate the number of atoms of moving z-constrained molecules |
1093 |
|
int nMovingZAtoms_local; |
1130 |
|
return totalForce / zconsIntegrator->totNumOfUnconsAtoms; |
1131 |
|
} |
1132 |
|
|
1211 |
– |
/** |
1212 |
– |
* |
1213 |
– |
*/ |
1133 |
|
|
1134 |
|
template<typename T> void ZConstraint<T>::PolicyByMass::update(){ |
1135 |
|
//calculate the number of atoms of moving z-constrained molecules |