| 509 | 
  | 
  //zero out the velocities of center of mass of unconstrained molecules  | 
| 510 | 
  | 
  //and the velocities of center of mass of every single z-constrained molecueles | 
| 511 | 
  | 
  zeroOutVel(); | 
| 512 | 
+ | 
 | 
| 513 | 
+ | 
  curZconsTime = zconsTime + info->getTime(); | 
| 514 | 
  | 
   | 
| 515 | 
  | 
  T::integrate(); | 
| 516 | 
  | 
 | 
| 527 | 
  | 
  double zsys; | 
| 528 | 
  | 
  double COM[3]; | 
| 529 | 
  | 
  double force[3]; | 
| 530 | 
+ | 
  double zSysCOMVel; | 
| 531 | 
  | 
 | 
| 532 | 
  | 
  T::calcForce(calcPot, calcStress); | 
| 533 | 
  | 
 | 
| 534 | 
  | 
  if (checkZConsState()){ | 
| 535 | 
+ | 
     | 
| 536 | 
+ | 
#ifdef IS_MPI | 
| 537 | 
+ | 
    if(worldRank == 0){ | 
| 538 | 
+ | 
#endif | 
| 539 | 
+ | 
      std::cerr << "\n" | 
| 540 | 
+ | 
                << "*******************************************\n" | 
| 541 | 
+ | 
                << " about to call zeroOutVel()\n" | 
| 542 | 
+ | 
                << "*******************************************\n" | 
| 543 | 
+ | 
                << "\n"; | 
| 544 | 
+ | 
#ifdef IS_MPI | 
| 545 | 
+ | 
    } | 
| 546 | 
+ | 
#endif | 
| 547 | 
  | 
    zeroOutVel(); | 
| 548 | 
< | 
   forcePolicy->update(); | 
| 548 | 
> | 
 | 
| 549 | 
> | 
#ifdef IS_MPI | 
| 550 | 
> | 
    if(worldRank == 0){ | 
| 551 | 
> | 
#endif | 
| 552 | 
> | 
      std::cerr << "\n" | 
| 553 | 
> | 
                << "*******************************************\n" | 
| 554 | 
> | 
                << " finished zeroOutVel()\n" | 
| 555 | 
> | 
                << "*******************************************\n" | 
| 556 | 
> | 
                << "\n"; | 
| 557 | 
> | 
#ifdef IS_MPI | 
| 558 | 
> | 
    } | 
| 559 | 
> | 
#endif | 
| 560 | 
> | 
     | 
| 561 | 
> | 
    forcePolicy->update(); | 
| 562 | 
  | 
  }   | 
| 563 | 
+ | 
   | 
| 564 | 
  | 
  zsys = calcZSys(); | 
| 565 | 
< | 
  cout << "---------------------------------------------------------------------" <<endl; | 
| 566 | 
< | 
  cout << "current time: " << info->getTime() << endl; | 
| 567 | 
< | 
  cout << "center of mass at z: " << zsys << endl;      | 
| 568 | 
< | 
  //cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; | 
| 569 | 
< | 
  cout << "before calcForce, the COMVel of system is " << calcSysCOMVel() <<endl; | 
| 565 | 
> | 
  zSysCOMVel = calcSysCOMVel(); | 
| 566 | 
> | 
#ifdef IS_MPI | 
| 567 | 
> | 
  if(worldRank == 0){ | 
| 568 | 
> | 
#endif | 
| 569 | 
> | 
    cout << "---------------------------------------------------------------------" <<endl; | 
| 570 | 
> | 
    cout << "current time: " << info->getTime() << endl; | 
| 571 | 
> | 
    cout << "center of mass at z: " << zsys << endl;      | 
| 572 | 
> | 
    cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl; | 
| 573 | 
  | 
 | 
| 574 | 
< | 
  //cout <<  "before doZConstraintForce, totalForce is " << calcTotalForce() << endl; | 
| 574 | 
> | 
#ifdef IS_MPI | 
| 575 | 
> | 
  } | 
| 576 | 
> | 
#endif | 
| 577 | 
  | 
 | 
| 578 | 
  | 
  //do zconstraint force;  | 
| 579 | 
  | 
  if (haveFixedZMols()) | 
| 583 | 
  | 
  if (haveMovingZMols()) | 
| 584 | 
  | 
    this->doHarmonic(); | 
| 585 | 
  | 
 | 
| 552 | 
– | 
  //cout <<  "after doHarmonic, totalForce is " << calcTotalForce() << endl; | 
| 553 | 
– | 
 | 
| 586 | 
  | 
  //write out forces and current positions of z-constraint molecules | 
| 587 | 
  | 
  if(info->getTime() >= curZconsTime){     | 
| 588 | 
  | 
   for(int i = 0; i < zconsMols.size(); i++){ | 
| 603 | 
  | 
    fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos); | 
| 604 | 
  | 
   curZconsTime += zconsTime; | 
| 605 | 
  | 
  } | 
| 606 | 
< | 
   | 
| 607 | 
< | 
  //cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;  | 
| 608 | 
< | 
  cout << "after calcForce, the COMVel of system is " << calcSysCOMVel() <<endl;  | 
| 606 | 
> | 
 | 
| 607 | 
> | 
  zSysCOMVel = calcSysCOMVel();   | 
| 608 | 
> | 
#ifdef IS_MPI | 
| 609 | 
> | 
  if(worldRank == 0){ | 
| 610 | 
> | 
#endif | 
| 611 | 
> | 
    cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl; | 
| 612 | 
> | 
#ifdef IS_MPI | 
| 613 | 
> | 
  } | 
| 614 | 
> | 
#endif | 
| 615 | 
> | 
 | 
| 616 | 
  | 
} | 
| 617 | 
  | 
 | 
| 618 | 
  | 
 | 
| 825 | 
  | 
      }  | 
| 826 | 
  | 
      totalFZ_local += fz[i]; | 
| 827 | 
  | 
 | 
| 828 | 
< | 
      cout << "Fixed Molecule --\tindex: " << indexOfZConsMols[i]  | 
| 829 | 
< | 
             <<"\tcurrent zpos: " << COM[whichDirection]  | 
| 830 | 
< | 
             << "\tcurrent fz: " <<fz[i] << endl; | 
| 828 | 
> | 
//       cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]  | 
| 829 | 
> | 
//              <<"\tcurrent zpos: " << COM[whichDirection]  | 
| 830 | 
> | 
//              << "\tcurrent fz: " <<fz[i] << endl; | 
| 831 | 
  | 
 | 
| 832 | 
  | 
    } | 
| 833 | 
  | 
     | 
| 863 | 
  | 
   | 
| 864 | 
  | 
  }  | 
| 865 | 
  | 
 | 
| 866 | 
< | 
  //cout << "after zero out z-constraint force on fixed z-constraint molecuels " | 
| 867 | 
< | 
  //       << "total force is " << calcTotalForce() << endl; | 
| 866 | 
> | 
//   cout << "after zero out z-constraint force on fixed z-constraint molecuels " | 
| 867 | 
> | 
//        << "total force is " << calcTotalForce() << endl; | 
| 868 | 
  | 
 | 
| 869 | 
  | 
  //calculate the number of atoms of moving z-constrained molecules | 
| 870 | 
  | 
  int nMovingZAtoms_local; | 
| 942 | 
  | 
 | 
| 943 | 
  | 
    if (states[i] == zcsMoving){ | 
| 944 | 
  | 
      zconsMols[i]->getCOM(COM); | 
| 945 | 
< | 
      cout << "Moving Molecule --\tindex: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; | 
| 945 | 
> | 
      cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]  | 
| 946 | 
> | 
           << "\tcurrent zpos: " << COM[whichDirection] << endl; | 
| 947 | 
  | 
     | 
| 948 | 
  | 
    diff = COM[whichDirection] -zPos[i]; | 
| 949 | 
  | 
     | 
| 1023 | 
  | 
#else | 
| 1024 | 
  | 
  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 1025 | 
  | 
#endif | 
| 1026 | 
< | 
 | 
| 1027 | 
< | 
  return changed > 0 ? true : false; | 
| 1026 | 
> | 
   | 
| 1027 | 
> | 
  return (changed > 0); | 
| 1028 | 
  | 
} | 
| 1029 | 
  | 
 | 
| 1030 | 
  | 
template<typename T> bool ZConstraint<T>::haveFixedZMols(){ | 
| 1203 | 
  | 
  nMovingZAtoms = nMovingZAtoms_local; | 
| 1204 | 
  | 
#endif | 
| 1205 | 
  | 
  totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms; | 
| 1206 | 
+ | 
 | 
| 1207 | 
+ | 
#ifdef IS_MPI | 
| 1208 | 
+ | 
  if(worldRank == 0){ | 
| 1209 | 
+ | 
#endif | 
| 1210 | 
+ | 
    std::cerr << "\n" | 
| 1211 | 
+ | 
              << "*******************************************\n" | 
| 1212 | 
+ | 
              << " fiished Policy by numbr()\n" | 
| 1213 | 
+ | 
              << "*******************************************\n" | 
| 1214 | 
+ | 
              << "\n"; | 
| 1215 | 
+ | 
#ifdef IS_MPI | 
| 1216 | 
+ | 
  } | 
| 1217 | 
+ | 
#endif | 
| 1218 | 
  | 
} | 
| 1219 | 
  | 
 | 
| 1220 | 
  | 
template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){ |