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 702 by tim, Wed Aug 20 14:50:32 2003 UTC vs.
Revision 736 by tim, Thu Aug 28 21:09:47 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 58 | Line 58 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
58        else{
59          sprintf( painCave.errMsg,
60                    "ZConstraint Warning: unknown force substraction policy, "
61 <                  "average force substraction policy is used\n");
61 >                  "PolicyByMass is used\n");
62          painCave.isFatal = 0;
63          simError();      
64 +        forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
65        }  
66      }
67    }
# Line 249 | Line 250 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
250            (*parameters)[i].zPos = COM[whichDirection];
251  
252            sprintf( painCave.errMsg,
253 <                     "ZConstraint warningr: Does not specify zpos for z-constraint molecule "
253 >                     "ZConstraint warning: Does not specify zpos for z-constraint molecule "
254                       "initial z coornidate will be used \n");
255             painCave.isFatal = 0;
256             simError();  
# Line 276 | Line 277 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
277        massOfZConsMols.push_back(molecules[i].getTotalMass());  
278  
279        zPos.push_back((*parameters)[searchResult].zPos);
280 <      cout << "index: "<< (*parameters)[searchResult].zconsIndex
281 <             <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
280 > //       cout << "index: "<< (*parameters)[searchResult].zconsIndex
281 > //              <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
282        kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
283        
284        molecules[i].getCOM(COM);
# Line 509 | Line 510 | template<typename T> void ZConstraint<T>::integrate(){
510    //zero out the velocities of center of mass of unconstrained molecules
511    //and the velocities of center of mass of every single z-constrained molecueles
512    zeroOutVel();
513 +
514 +  curZconsTime = zconsTime + info->getTime();
515    
516    T::integrate();
517  
# Line 530 | Line 533 | template<typename T> void ZConstraint<T>::calcForce(in
533    T::calcForce(calcPot, calcStress);
534  
535    if (checkZConsState()){
536 +    
537 + #ifdef IS_MPI
538 +    if(worldRank == 0){
539 + #endif
540 + //       std::cerr << "\n"
541 + //              << "*******************************************\n"
542 + //              << " about to call zeroOutVel()\n"
543 + //              << "*******************************************\n"
544 + //              << "\n";
545 + #ifdef IS_MPI
546 +    }
547 + #endif
548      zeroOutVel();
534 –   forcePolicy->update();
535 –  }  
549  
550 + #ifdef IS_MPI
551 +    if(worldRank == 0){
552 + #endif
553 + //       std::cerr << "\n"
554 + //              << "*******************************************\n"
555 + //              << " finished zeroOutVel()\n"
556 + //              << "*******************************************\n"
557 + //              << "\n";
558 + #ifdef IS_MPI
559 +    }
560 + #endif
561 +    
562 +    forcePolicy->update();
563 +  }  
564 +  
565    zsys = calcZSys();
566    zSysCOMVel = calcSysCOMVel();
567   #ifdef IS_MPI
568    if(worldRank == 0){
569   #endif
570 <    cout << "---------------------------------------------------------------------" <<endl;
571 <    cout << "current time: " << info->getTime() << endl;
572 <    cout << "center of mass at z: " << zsys << endl;    
573 <    cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
570 > //     cout << "---------------------------------------------------------------------" <<endl;
571 > //     cout << "current time: " << info->getTime() << endl;
572 > //     cout << "center of mass at z: " << zsys << endl;    
573 > //     cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
574  
575   #ifdef IS_MPI
576    }
# Line 551 | Line 579 | template<typename T> void ZConstraint<T>::calcForce(in
579    //do zconstraint force;
580    if (haveFixedZMols())
581      this->doZconstraintForce();
582 <    
582 >
583    //use harmonical poteintial to move the molecules to the specified positions
584    if (haveMovingZMols())
585      this->doHarmonic();
# Line 581 | Line 609 | template<typename T> void ZConstraint<T>::calcForce(in
609   #ifdef IS_MPI
610    if(worldRank == 0){
611   #endif
612 <    cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
612 > //    cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
613   #ifdef IS_MPI
614    }
615   #endif
588 –
616   }
617  
618  
# Line 681 | Line 708 | template<typename T> void ZConstraint<T>::zeroOutVel()
708   #ifdef IS_MPI
709    if(worldRank == 0){
710   #endif
711 <    cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
711 > //     cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
712   #ifdef IS_MPI
713    }
714   #endif
# Line 755 | Line 782 | template<typename T> void ZConstraint<T>::zeroOutVel()
782   #ifdef IS_MPI
783    if(worldRank == 0){
784   #endif
785 <    cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
785 > //     cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
786   #ifdef IS_MPI
787    }
788   #endif
# Line 775 | Line 802 | template<typename T> void ZConstraint<T>::doZconstrain
802    double COM[3];
803    double force[3];
804  
778 –
779 –
805    //constrain the molecules which do not reach the specified positions  
806      
807    //Zero Out the force of z-contrained molecules    
# Line 798 | Line 823 | template<typename T> void ZConstraint<T>::doZconstrain
823        }
824        totalFZ_local += fz[i];
825  
826 <      cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
827 <             <<"\tcurrent zpos: " << COM[whichDirection]
828 <             << "\tcurrent fz: " <<fz[i] << endl;
826 >      //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
827 >      //      <<"\tcurrent zpos: " << COM[whichDirection]
828 >      //       << "\tcurrent fz: " <<fz[i] << endl;
829  
830 +
831      }
832      
833    }
# Line 827 | Line 853 | template<typename T> void ZConstraint<T>::doZconstrain
853        zconsAtoms = zconsMols[i]->getMyAtoms();  
854      
855        for(int j =0; j < nAtomOfCurZConsMol; j++) {
856 <      force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
857 <        //force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
856 >        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
857 >        force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
858          zconsAtoms[j]->addFrc(force);
859        }
860      
# Line 836 | Line 862 | template<typename T> void ZConstraint<T>::doZconstrain
862    
863    }
864  
865 <  //cout << "after zero out z-constraint force on fixed z-constraint molecuels "
866 <  //       << "total force is " << calcTotalForce() << endl;
865 > //   cout << "after zero out z-constraint force on fixed z-constraint molecuels "
866 > //        << "total force is " << calcTotalForce() << endl;
867  
842 –  //calculate the number of atoms of moving z-constrained molecules
843 –  int nMovingZAtoms_local;
844 –  int nMovingZAtoms;
845 –  
846 –  nMovingZAtoms_local = 0;
847 –  for(int i = 0; i < zconsMols.size(); i++)
848 –    if(states[i] == zcsMoving)
849 –     nMovingZAtoms_local += zconsMols[i]->getNAtoms();
850 –  
851 – #ifdef IS_MPI
852 –  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1,
853 –                      MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
854 – #else
855 –  nMovingZAtoms = nMovingZAtoms_local;
856 – #endif
857 –
868    force[0]= 0;
869    force[1]= 0;
870    force[2]= 0;
# Line 865 | Line 875 | template<typename T> void ZConstraint<T>::doZconstrain
875       Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
876      
877       for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){          
878 <       force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
879 <       //force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
878 >       //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
879 >       force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
880         unconsAtoms[j]->addFrc(force);
881       }
882      
# Line 879 | Line 889 | template<typename T> void ZConstraint<T>::doZconstrain
889        Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
890  
891        for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
892 <        force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
893 <        //force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
892 >        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
893 >        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
894          movingZAtoms[j]->addFrc(force);
895        }
896      }
# Line 915 | Line 925 | template<typename T> void ZConstraint<T>::doHarmonic()
925  
926      if (states[i] == zcsMoving){
927        zconsMols[i]->getCOM(COM);
928 <      cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
928 > //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
929 > //         << "\tcurrent zpos: " << COM[whichDirection] << endl;
930 >
931 >      diff = COM[whichDirection] -zPos[i];
932      
920 –    diff = COM[whichDirection] -zPos[i];
921 –    
933        harmonicU = 0.5 * kz[i] * diff * diff;  
934 <    info->lrPot += harmonicU;
934 >      info->lrPot += harmonicU;
935  
936        harmonicF =  - kz[i] * diff;
937        totalFZ_local += harmonicF;
# Line 930 | Line 941 | template<typename T> void ZConstraint<T>::doHarmonic()
941        Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
942  
943         for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){          
944 <        force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
945 <         //force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
944 >        //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
945 >        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
946           movingZAtoms[j]->addFrc(force);
947         }
948      }
# Line 954 | Line 965 | template<typename T> void ZConstraint<T>::doHarmonic()
965       Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
966      
967       for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){          
968 <       force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
969 <       //force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
968 >       //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
969 >       force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
970         unconsAtoms[j]->addFrc(force);    
971       }
972    }  
# Line 996 | Line 1007 | template<typename T> bool ZConstraint<T>::checkZConsSt
1007    MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1008   #endif
1009  
1010 <  return changed > 0 ? true : false;
1010 >  return (changed > 0);
1011 >
1012   }
1013  
1014   template<typename T> bool ZConstraint<T>::haveFixedZMols(){
# Line 1018 | Line 1030 | template<typename T> bool ZConstraint<T>::haveFixedZMo
1030    MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1031   #endif
1032  
1033 <  return havingFixed > 0 ? true : false;
1033 >  return (havingFixed > 0);
1034   }
1035  
1036  
# Line 1044 | Line 1056 | template<typename T> bool ZConstraint<T>::haveMovingZM
1056    MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1057   #endif
1058  
1059 <  return havingMoving > 0 ? true : false;
1059 >  return (havingMoving > 0);
1060    
1061   }
1062  
# Line 1175 | Line 1187 | template<typename T> void ZConstraint<T>::PolicyByNumb
1187    nMovingZAtoms = nMovingZAtoms_local;
1188   #endif
1189    totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1190 +
1191 + #ifdef IS_MPI
1192 +  if(worldRank == 0){
1193 + #endif
1194 + //    std::cerr << "\n"
1195 + //            << "*******************************************\n"
1196 + //            << " fiished Policy by numbr()\n"
1197 + //            << "*******************************************\n"
1198 + //            << "\n";
1199 + #ifdef IS_MPI
1200 +  }
1201 + #endif
1202   }
1203  
1204   template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){

Diff Legend

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