ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/ZConstraint.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 701 by tim, Wed Aug 20 14:34:04 2003 UTC vs.
Revision 733 by tim, Wed Aug 27 19:23:29 2003 UTC

# Line 32 | Line 32 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
32    if(!data){
33      sprintf( painCave.errMsg,
34                 "ZConstraint Warning: User does not set force substraction policy, "
35 <               "average force substraction policy is used\n");
35 >               "PolicyByMass is used\n");
36      painCave.isFatal = 0;
37      simError();      
38  
39 <    forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
39 >    forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
40    }
41    else{
42      policy = dynamic_cast<StringData*>(data);
# Line 44 | Line 44 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
44      if(!policy){
45        sprintf( painCave.errMsg,
46                   "ZConstraint Error: Convertion from GenericData to StringData failure, "
47 <                 "average force substraction policy is used\n");
47 >                 "PolicyByMass is used\n");
48        painCave.isFatal = 0;
49        simError();      
50  
51 <      forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
51 >      forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
52      }
53      else{
54        if(policy->getData() == "BYNUMBER")
# Line 276 | Line 276 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
276        massOfZConsMols.push_back(molecules[i].getTotalMass());  
277  
278        zPos.push_back((*parameters)[searchResult].zPos);
279 <      cout << "index: "<< (*parameters)[searchResult].zconsIndex
280 <             <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
279 > //       cout << "index: "<< (*parameters)[searchResult].zconsIndex
280 > //              <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
281        kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
282        
283        molecules[i].getCOM(COM);
# Line 509 | Line 509 | template<typename T> void ZConstraint<T>::integrate(){
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  
# Line 525 | Line 527 | template<typename T> void ZConstraint<T>::calcForce(in
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())
580      this->doZconstraintForce();
581 <    
581 >
582    //use harmonical poteintial to move the molecules to the specified positions
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++){
# Line 571 | Line 603 | template<typename T> void ZConstraint<T>::calcForce(in
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  
# Line 669 | Line 707 | template<typename T> void ZConstraint<T>::zeroOutVel()
707   #ifdef IS_MPI
708    if(worldRank == 0){
709   #endif
710 <    cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
710 > //     cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
711   #ifdef IS_MPI
712    }
713   #endif
# Line 743 | Line 781 | template<typename T> void ZConstraint<T>::zeroOutVel()
781   #ifdef IS_MPI
782    if(worldRank == 0){
783   #endif
784 <    cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
784 > //     cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
785   #ifdef IS_MPI
786    }
787   #endif
# Line 763 | Line 801 | template<typename T> void ZConstraint<T>::doZconstrain
801    double COM[3];
802    double force[3];
803  
766 –
767 –
804    //constrain the molecules which do not reach the specified positions  
805      
806    //Zero Out the force of z-contrained molecules    
# Line 786 | Line 822 | template<typename T> void ZConstraint<T>::doZconstrain
822        }
823        totalFZ_local += fz[i];
824  
825 <      cout << "Fixed Molecule --\tindex: " << indexOfZConsMols[i]
826 <             <<"\tcurrent zpos: " << COM[whichDirection]
827 <             << "\tcurrent fz: " <<fz[i] << endl;
825 >      //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
826 >      //      <<"\tcurrent zpos: " << COM[whichDirection]
827 >      //       << "\tcurrent fz: " <<fz[i] << endl;
828  
829 +
830      }
831      
832    }
# Line 815 | Line 852 | template<typename T> void ZConstraint<T>::doZconstrain
852        zconsAtoms = zconsMols[i]->getMyAtoms();  
853      
854        for(int j =0; j < nAtomOfCurZConsMol; j++) {
855 <      force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
856 <        //force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
855 >        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
856 >        force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
857          zconsAtoms[j]->addFrc(force);
858        }
859      
# Line 824 | Line 861 | template<typename T> void ZConstraint<T>::doZconstrain
861    
862    }
863  
864 <  //cout << "after zero out z-constraint force on fixed z-constraint molecuels "
865 <  //       << "total force is " << calcTotalForce() << endl;
864 > //   cout << "after zero out z-constraint force on fixed z-constraint molecuels "
865 > //        << "total force is " << calcTotalForce() << endl;
866  
830 –  //calculate the number of atoms of moving z-constrained molecules
831 –  int nMovingZAtoms_local;
832 –  int nMovingZAtoms;
833 –  
834 –  nMovingZAtoms_local = 0;
835 –  for(int i = 0; i < zconsMols.size(); i++)
836 –    if(states[i] == zcsMoving)
837 –     nMovingZAtoms_local += zconsMols[i]->getNAtoms();
838 –  
839 – #ifdef IS_MPI
840 –  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1,
841 –                      MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
842 – #else
843 –  nMovingZAtoms = nMovingZAtoms_local;
844 – #endif
845 –
867    force[0]= 0;
868    force[1]= 0;
869    force[2]= 0;
# Line 853 | Line 874 | template<typename T> void ZConstraint<T>::doZconstrain
874       Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
875      
876       for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){          
877 <       force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
878 <       //force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
877 >       //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
878 >       force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
879         unconsAtoms[j]->addFrc(force);
880       }
881      
# Line 867 | Line 888 | template<typename T> void ZConstraint<T>::doZconstrain
888        Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
889  
890        for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
891 <        force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
892 <        //force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
891 >        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
892 >        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
893          movingZAtoms[j]->addFrc(force);
894        }
895      }
# Line 903 | Line 924 | template<typename T> void ZConstraint<T>::doHarmonic()
924  
925      if (states[i] == zcsMoving){
926        zconsMols[i]->getCOM(COM);
927 <      cout << "Moving Molecule --\tindex: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
927 > //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
928 > //         << "\tcurrent zpos: " << COM[whichDirection] << endl;
929 >
930 >      diff = COM[whichDirection] -zPos[i];
931      
908 –    diff = COM[whichDirection] -zPos[i];
909 –    
932        harmonicU = 0.5 * kz[i] * diff * diff;  
933 <    info->lrPot += harmonicU;
933 >      info->lrPot += harmonicU;
934  
935        harmonicF =  - kz[i] * diff;
936        totalFZ_local += harmonicF;
# Line 918 | Line 940 | template<typename T> void ZConstraint<T>::doHarmonic()
940        Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
941  
942         for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){          
943 <        force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
944 <         //force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
943 >        //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
944 >        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
945           movingZAtoms[j]->addFrc(force);
946         }
947      }
# Line 942 | Line 964 | template<typename T> void ZConstraint<T>::doHarmonic()
964       Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
965      
966       for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){          
967 <       force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
968 <       //force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
967 >       //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
968 >       force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
969         unconsAtoms[j]->addFrc(force);    
970       }
971    }  
# Line 984 | Line 1006 | template<typename T> bool ZConstraint<T>::checkZConsSt
1006    MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1007   #endif
1008  
1009 <  return changed > 0 ? true : false;
1009 >  return (changed > 0);
1010 >
1011   }
1012  
1013   template<typename T> bool ZConstraint<T>::haveFixedZMols(){
# Line 1006 | Line 1029 | template<typename T> bool ZConstraint<T>::haveFixedZMo
1029    MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1030   #endif
1031  
1032 <  return havingFixed > 0 ? true : false;
1032 >  return (havingFixed > 0);
1033   }
1034  
1035  
# Line 1032 | Line 1055 | template<typename T> bool ZConstraint<T>::haveMovingZM
1055    MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1056   #endif
1057  
1058 <  return havingMoving > 0 ? true : false;
1058 >  return (havingMoving > 0);
1059    
1060   }
1061  
# Line 1163 | Line 1186 | template<typename T> void ZConstraint<T>::PolicyByNumb
1186    nMovingZAtoms = nMovingZAtoms_local;
1187   #endif
1188    totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1189 +
1190 + #ifdef IS_MPI
1191 +  if(worldRank == 0){
1192 + #endif
1193 + //    std::cerr << "\n"
1194 + //            << "*******************************************\n"
1195 + //            << " fiished Policy by numbr()\n"
1196 + //            << "*******************************************\n"
1197 + //            << "\n";
1198 + #ifdef IS_MPI
1199 +  }
1200 + #endif
1201   }
1202  
1203   template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines