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 682 by tim, Tue Aug 12 17:51:33 2003 UTC vs.
Revision 693 by tim, Wed Aug 13 19:21:53 2003 UTC

# Line 156 | Line 156 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
156          simError();
157          }
158  
159 <      maxIndex = (*parameters)[parameters->size()].zconsIndex;
159 >      maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
160  
161   #ifndef IS_MPI
162        totalNumMol = nMols;
# Line 179 | Line 179 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
179  
180   #ifndef IS_MPI
181              for(int j = 0; j < nMols; j++){
182 <              if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){
183 <                 molecules[i].getCOM(COM);
182 >              if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
183 >                 molecules[j].getCOM(COM);
184                            break;
185                }
186              }
# Line 197 | Line 197 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
197            
198             if (worldRank == whichNode ){
199                                                  
200 <             for(int i = 0; i < nMols; i++)
201 <               if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){
202 <                 molecules[i].getCOM(COM);
200 >             for(int j = 0; j < nMols; j++)
201 >               if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
202 >                 molecules[j].getCOM(COM);
203                                           break;
204                 }
205                                  
# Line 261 | Line 261 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
261      simError();
262    }
263  
264 <  for(int i = 0; i < zconsMols.size(); i++)
264 >  //determine the states of z-constraint molecules
265 >  for(int i = 0; i < zconsMols.size(); i++){
266      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
267 +
268 +         zconsMols[i]->getCOM(COM);
269 +    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
270 +                states.push_back(zcsFixed);
271 +         else
272 +                states.push_back(zcsMoving);
273 +  }
274    
275   #endif
276 +
277 +  //get total masss of unconstraint molecules
278 +  double totalMassOfUncons_local;
279 +  totalMassOfUncons_local = 0;
280 +  
281 +  for(int i = 0; i < unconsMols.size(); i++)
282 +    totalMassOfUncons_local += unconsMols[i]->getTotalMass();
283 +    
284 + #ifndef IS_MPI
285 +  totalMassOfUncons = totalMassOfUncons_local;
286 + #else
287 +  MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
288 + #endif
289  
290 +
291    //get total number of unconstrained atoms
292    int nUnconsAtoms_local;
293    nUnconsAtoms_local = 0;
# Line 276 | Line 298 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
298    totNumOfUnconsAtoms = nUnconsAtoms_local;
299   #else
300    MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
301 < #endif
301 > #endif  
302  
303 <  checkZConsState();
282 <
283 <  //  
303 >  // creat zconsWriter  
304    fzOut = new ZConsWriter(zconsOutput.c_str());  
305    
306    if(!fzOut){
# Line 328 | Line 348 | template<typename T> void ZConstraint<T>::update()
348      
349        zconsMols.push_back(&molecules[i]);      
350        zPos.push_back((*parameters)[index].zPos);
351 <        kz.push_back((*parameters)[index].kRatio * zForceConst);
352 <
351 >        kz.push_back((*parameters)[index].kRatio * zForceConst);
352 >                        
353        massOfZConsMols.push_back(molecules[i].getTotalMass());  
354        
355        molecules[i].getCOM(COM);
# Line 342 | Line 362 | template<typename T> void ZConstraint<T>::update()
362  
363      }
364    }
365 +
366 +  //determine the states of z-constraint molecules
367 +  for(int i = 0; i < zconsMols.size(); i++){
368 +           zconsMols[i]->getCOM(COM);
369 +      if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
370 +                  states.push_back(zcsFixed);
371 +           else
372 +                  states.push_back(zcsMoving);
373 +  }
374 +
375      
376    //The reason to declare fz and indexOfZconsMols as pointer to array is
377    // that we want to make the MPI communication simple
# Line 410 | Line 440 | template<typename T> int ZConstraint<T>::isZConstraint
440    return -1;
441   }
442  
413 – /**
414 – * Description:
415 – *  Reset the z coordinates
416 – */
443   template<typename T> void ZConstraint<T>::integrate(){
444    
445    //zero out the velocities of center of mass of unconstrained molecules
# Line 434 | Line 460 | template<typename T> void ZConstraint<T>::calcForce(in
460  
461  
462   template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
463 +  double zsys;
464  
465    T::calcForce(calcPot, calcStress);
466  
467    if (checkZConsState())
468 <    zeroOutVel();
468 >  zeroOutVel();
469 >  
470 >  zsys = calcZSys();
471 >  cout << "---------------------------------------------------------------------" <<endl;
472 >  cout << "current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;      
473 >  cout << "before calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
474 >        
475  
476    //do zconstraint force;
477    if (haveFixedZMols())
478      this->doZconstraintForce();
479 <  
479 >
480 >
481 >      
482    //use harmonical poteintial to move the molecules to the specified positions
483    if (haveMovingZMols())
484 <    //this->doHarmonic();
484 >    this->doHarmonic();
485    
486    fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
487 <      
487 >  cout << "after calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
488   }
489  
490   template<typename T> double ZConstraint<T>::calcZSys()
# Line 521 | Line 556 | template<typename T> void ZConstraint<T>::zeroOutVel()
556    Atom** fixedZAtoms;  
557    double COMvel[3];
558    double vel[3];
559 <  
559 >
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 >
571    //zero out the velocities of center of mass of fixed z-constrained molecules
572    
573    for(int i = 0; i < zconsMols.size(); i++){
574  
575      if (states[i] == zcsFixed){
576  
577 <        zconsMols[i]->getCOMvel(COMvel);      
577 >           zconsMols[i]->getCOMvel(COMvel);      
578 >                cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
579 >
580        fixedZAtoms = zconsMols[i]->getMyAtoms();
581            
582        for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
583          fixedZAtoms[j]->getVel(vel);
584 <          vel[whichDirection] -= COMvel[whichDirection];
585 <          fixedZAtoms[j]->setVel(vel);
584 >             vel[whichDirection] -= COMvel[whichDirection];
585 >             fixedZAtoms[j]->setVel(vel);
586        }
587 <          
587 >
588 >                zconsMols[i]->getCOMvel(COMvel);
589 >                cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
590      }
591          
592    }
593 <  
593 >
594 >        cout << "before resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;    
595 >                  
596    // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
597    double MVzOfMovingMols_local;
598    double MVzOfMovingMols;
# Line 555 | Line 607 | template<typename T> void ZConstraint<T>::zeroOutVel()
607      MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
608    }
609  
610 <  for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
559 <
610 >  for(int i = 0; i < zconsMols.size(); i++){
611      if (states[i] == zcsMoving){
612        zconsMols[i]->getCOMvel(COMvel);
613        MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
# Line 591 | Line 642 | template<typename T> void ZConstraint<T>::zeroOutVel()
642  
643    //modify the velocities of moving z-constrained molecuels
644    Atom** movingZAtoms;
645 <  for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
645 >  for(int i = 0; i < zconsMols.size(); i++){
646  
647      if (states[i] ==zcsMoving){
648    
649        movingZAtoms = zconsMols[i]->getMyAtoms();
650 <        for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
650 >           for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
651          movingZAtoms[j]->getVel(vel);
652          vel[whichDirection] -= vzOfMovingMols;
653 <          movingZAtoms[j]->setVel(vel);
654 <        }
653 >             movingZAtoms[j]->setVel(vel);
654 >          }
655            
656 <    }
656 >   }
657  
658 <  }
658 > }
659  
660 +        cout << "after resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
661 +
662   }
663  
664   template<typename T> void ZConstraint<T>::doZconstraintForce(){
# Line 616 | Line 669 | template<typename T> void ZConstraint<T>::doZconstrain
669    double COMvel[3];  
670    double COM[3];
671    double force[3];
619 –  double zsys;
672  
673    int nMovingZMols_local;
674    int nMovingZMols;
675  
676    //constrain the molecules which do not reach the specified positions  
625 –
626 –   zsys = calcZSys();
627 –   cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;  
677      
678    //Zero Out the force of z-contrained molecules    
679    totalFZ_local = 0;
# Line 653 | Line 702 | template<typename T> void ZConstraint<T>::doZconstrain
702  
703    //calculate the number of atoms of moving z-constrained molecules
704    nMovingZMols_local = 0;
705 <  for(int i = 0; zconsMols.size(); i++){
705 >  for(int i = 0; i < zconsMols.size(); i++)
706      if(states[i] == zcsMoving)
707 <        nMovingZMols_local += massOfZConsMols[i];
708 <  }
707 >           nMovingZMols_local += massOfZConsMols[i];
708 >  
709   #ifdef IS_MPI
710    MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
711    MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
# Line 670 | Line 719 | template<typename T> void ZConstraint<T>::doZconstrain
719    force[2]= 0;
720    force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);
721  
722 <  //modify the velocites of unconstrained molecules
722 >  //modify the forces of unconstrained molecules
723    for(int i = 0; i < unconsMols.size(); i++){
724      
725       Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
# Line 680 | Line 729 | template<typename T> void ZConstraint<T>::doZconstrain
729      
730    }      
731  
732 < //modify the velocities of moving z-constrained molecules
732 > //modify the forces of moving z-constrained molecules
733    for(int i = 0; i < zconsMols.size(); i++) {
734     if (states[i] == zcsMoving){
735                  
# Line 769 | Line 818 | template<typename T> void ZConstraint<T>::doHarmonic()
818   template<typename T> void ZConstraint<T>::doHarmonic(){
819    double force[3];
820    double harmonicU;
821 +  double harmonicF;
822    double COM[3];
823    double diff;
824 +  double totalFZ;
825          
826    force[0] = 0;
827    force[1] = 0;
828    force[2] = 0;
829  
830 +  totalFZ = 0;
831 +
832    cout << "Moving Molecules" << endl;  
833    for(int i = 0; i < zconsMols.size(); i++) {
834  
835      if (states[i] == zcsMoving){
836 <      zconsMols[i]->getCOM(COM):
836 >      zconsMols[i]->getCOM(COM);
837        cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
838                  
839                  diff = COM[whichDirection] -zPos[i];
840                  
841        harmonicU = 0.5 * kz[i] * diff * diff;  
842 <                info->ltPot += harmonicU;
842 >                info->lrPot += harmonicU;
843  
844 <                force[whichDirection] = - kz[i] * diff / zconsMols[i]->getNAtoms();
844 >      harmonicF =  - kz[i] * diff / zconsMols[i]->getNAtoms();
845 >                force[whichDirection] = harmonicF;
846 >      totalFZ += harmonicF;
847                  
848        Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
849  
# Line 798 | Line 853 | template<typename T> void ZConstraint<T>::doHarmonic()
853  
854    }
855  
856 +  force[0]= 0;
857 +  force[1]= 0;
858 +  force[2]= 0;
859 +  force[whichDirection] = -totalFZ /totNumOfUnconsAtoms;
860 +
861 +  //modify the forces of unconstrained molecules
862 +  for(int i = 0; i < unconsMols.size(); i++){
863 +    
864 +     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
865 +    
866 +     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)          
867 +       unconsAtoms[j]->addFrc(force);    
868 +  }  
869 +
870   }
871 +
872 + template<typename T> double ZConstraint<T>::calcCOMVel()
873 + {
874 +  double MVzOfMovingMols_local;
875 +  double MVzOfMovingMols;
876 +  double totalMassOfMovingZMols_local;
877 +  double totalMassOfMovingZMols;
878 +  double COMvel[3];
879 +      
880 +  MVzOfMovingMols_local = 0;
881 +  totalMassOfMovingZMols_local = 0;
882 +
883 +  for(int i =0; i < unconsMols.size(); i++){
884 +    unconsMols[i]->getCOMvel(COMvel);
885 +    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
886 +  }
887 +
888 +  for(int i = 0; i < zconsMols.size(); i++){
889 +
890 +    if (states[i] == zcsMoving){
891 +      zconsMols[i]->getCOMvel(COMvel);
892 +      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
893 +      totalMassOfMovingZMols_local += massOfZConsMols[i];              
894 +    }
895 +                
896 +  }
897 +
898 + #ifndef IS_MPI
899 +  MVzOfMovingMols = MVzOfMovingMols_local;
900 +  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
901 + #else
902 +  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
903 +  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
904 + #endif
905 +
906 +  double vzOfMovingMols;
907 +  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
908 +
909 +  return vzOfMovingMols;
910 + }
911 +
912 +
913 + template<typename T> double ZConstraint<T>::calcCOMVel2()
914 + {
915 +  double COMvel[3];
916 +  double tempMVz = 0;
917 +  int index;
918 +                
919 +  for(int i =0 ; i < nMols; i++){
920 +         index = isZConstraintMol(&molecules[i]);
921 +    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 +         }
929 +  }
930 +        
931 +  return tempMVz /totalMassOfUncons;
932 + }

Diff Legend

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