ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/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 696 by tim, Thu Aug 14 16:16:39 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 341 | Line 361 | template<typename T> void ZConstraint<T>::update()
361        massOfUnconsMols.push_back(molecules[i].getTotalMass());
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 430 | Line 456 | template<typename T> void ZConstraint<T>::integrate(){
456   *
457   *
458   *
459 < */
434 <
435 <
459 > */
460   template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
461 +  double zsys;
462  
463    T::calcForce(calcPot, calcStress);
464  
465    if (checkZConsState())
466      zeroOutVel();
467 +  
468 +  zsys = calcZSys();
469 +  cout << "---------------------------------------------------------------------" <<endl;
470 +  cout << "current time: " << info->getTime() << endl;
471 +  cout << "center of mass at z: " << zsys << endl;      
472 +  //cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
473 +  cout << "before calcForce, the COMVel of system is " << calcSysCOMVel() <<endl;
474  
475 +  //cout <<      "before doZConstraintForce, totalForce is " << calcTotalForce() << endl;
476 +
477    //do zconstraint force;
478    if (haveFixedZMols())
479      this->doZconstraintForce();
480 <  
480 >    
481    //use harmonical poteintial to move the molecules to the specified positions
482    if (haveMovingZMols())
483 <    //this->doHarmonic();
484 <  
483 >    this->doHarmonic();
484 >
485 >  //cout <<      "after doHarmonic, totalForce is " << calcTotalForce() << endl;
486 >
487    fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
488 <      
488 >  //cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
489 >  cout << "after calcForce, the COMVel of system is " << calcSysCOMVel() <<endl;
490   }
491  
492   template<typename T> double ZConstraint<T>::calcZSys()
# Line 459 | Line 496 | template<typename T> double ZConstraint<T>::calcZSys()
496    double totalMass;
497    double totalMZ_local;
498    double totalMZ;
462 –  double massOfUncons_local;
499    double massOfCurMol;
500    double COM[3];
501 <  
501 >        
502    totalMass_local = 0;
467 –  totalMass = 0;
503    totalMZ_local = 0;
469 –  totalMZ = 0;
470 –  massOfUncons_local = 0;
471 –    
504    
505    for(int i = 0; i < nMols; i++){
506      massOfCurMol = molecules[i].getTotalMass();
# Line 476 | Line 508 | template<typename T> double ZConstraint<T>::calcZSys()
508      
509      totalMass_local += massOfCurMol;
510      totalMZ_local += massOfCurMol * COM[whichDirection];
511 <    
480 <    if(isZConstraintMol(&molecules[i]) == -1){
481 <    
482 <      massOfUncons_local += massOfCurMol;
483 <    }  
484 <    
511 >
512    }
513 +
514    
487 –  
515   #ifdef IS_MPI  
516    MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
517    MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
518 <  MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
492 < #else
518 > #else
519    totalMass = totalMass_local;
520    totalMZ = totalMZ_local;
521 <  totalMassOfUncons = massOfUncons_local;
496 < #endif  
521 > #endif  
522  
523    double zsys;
524    zsys = totalMZ / totalMass;
# Line 521 | Line 546 | template<typename T> void ZConstraint<T>::zeroOutVel()
546    Atom** fixedZAtoms;  
547    double COMvel[3];
548    double vel[3];
549 <  
549 >
550    //zero out the velocities of center of mass of fixed z-constrained molecules
551    
552    for(int i = 0; i < zconsMols.size(); i++){
553  
554      if (states[i] == zcsFixed){
555  
556 <        zconsMols[i]->getCOMvel(COMvel);      
556 >           zconsMols[i]->getCOMvel(COMvel);      
557 >                //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
558 >
559        fixedZAtoms = zconsMols[i]->getMyAtoms();
560            
561        for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
562          fixedZAtoms[j]->getVel(vel);
563 <          vel[whichDirection] -= COMvel[whichDirection];
564 <          fixedZAtoms[j]->setVel(vel);
563 >             vel[whichDirection] -= COMvel[whichDirection];
564 >             fixedZAtoms[j]->setVel(vel);
565        }
566 <          
566 >
567 >                zconsMols[i]->getCOMvel(COMvel);
568 >                //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
569      }
570          
571    }
572 <  
572 >
573 >        //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;      
574 >        cout << "before resetting the COMVel of sytem is " << calcSysCOMVel() << endl;  
575 >                  
576    // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
577    double MVzOfMovingMols_local;
578    double MVzOfMovingMols;
# Line 555 | Line 587 | template<typename T> void ZConstraint<T>::zeroOutVel()
587      MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
588    }
589  
590 <  for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
559 <
590 >  for(int i = 0; i < zconsMols.size(); i++){
591      if (states[i] == zcsMoving){
592        zconsMols[i]->getCOMvel(COMvel);
593        MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
# Line 591 | Line 622 | template<typename T> void ZConstraint<T>::zeroOutVel()
622  
623    //modify the velocities of moving z-constrained molecuels
624    Atom** movingZAtoms;
625 <  for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
625 >  for(int i = 0; i < zconsMols.size(); i++){
626  
627      if (states[i] ==zcsMoving){
628    
629        movingZAtoms = zconsMols[i]->getMyAtoms();
630 <        for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
630 >           for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
631          movingZAtoms[j]->getVel(vel);
632          vel[whichDirection] -= vzOfMovingMols;
633 <          movingZAtoms[j]->setVel(vel);
634 <        }
633 >             movingZAtoms[j]->setVel(vel);
634 >          }
635            
636 <    }
636 >   }
637  
638 <  }
638 > }
639  
640 +        cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() <<endl;
641 +
642   }
643  
644   template<typename T> void ZConstraint<T>::doZconstraintForce(){
# Line 616 | Line 649 | template<typename T> void ZConstraint<T>::doZconstrain
649    double COMvel[3];  
650    double COM[3];
651    double force[3];
619 –  double zsys;
652  
621 –  int nMovingZMols_local;
622 –  int nMovingZMols;
653  
624 –  //constrain the molecules which do not reach the specified positions  
654  
655 <   zsys = calcZSys();
627 <   cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;  
655 >  //constrain the molecules which do not reach the specified positions  
656      
657    //Zero Out the force of z-contrained molecules    
658    totalFZ_local = 0;
# Line 645 | Line 673 | template<typename T> void ZConstraint<T>::doZconstrain
673        }
674        totalFZ_local += fz[i];
675  
676 <      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
676 >      cout << "index: " << indexOfZConsMols[i]
677 >                                <<"\tcurrent zpos: " << COM[whichDirection]
678 >                                << "\tcurrent fz: " <<fz[i] << endl;
679  
680      }
681            
682    }
683 +        
684 +  // apply negative to fixed z-constrained molecues;
685 +  force[0]= 0;
686 +  force[1]= 0;
687 +  force[2]= 0;
688  
689 +  for(int i = 0; i < zconsMols.size(); i++){
690 +
691 +    if (states[i] == zcsFixed){  
692 +        
693 +      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
694 +      zconsAtoms = zconsMols[i]->getMyAtoms();  
695 +    
696 +      for(int j =0; j < nAtomOfCurZConsMol; j++) {
697 +        force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
698 +        zconsAtoms[j]->addFrc(force);
699 +      }
700 +                
701 +    }
702 +        
703 +  }
704 +
705 +  cout << "after zero out z-constraint force on fixed z-constraint molecuels "
706 +                  << "total force is " << calcTotalForce() << endl;
707    //calculate the number of atoms of moving z-constrained molecules
708 <  nMovingZMols_local = 0;
709 <  for(int i = 0; zconsMols.size(); i++){
708 >  int nMovingZAtoms_local;
709 >  int nMovingZAtoms;
710 >        
711 >  nMovingZAtoms_local = 0;
712 >  for(int i = 0; i < zconsMols.size(); i++)
713      if(states[i] == zcsMoving)
714 <        nMovingZMols_local += massOfZConsMols[i];
715 <  }
714 >           nMovingZAtoms_local += zconsMols[i]->getNAtoms();
715 >  
716   #ifdef IS_MPI
717    MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
718 <  MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
718 >  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
719   #else
720    totalFZ = totalFZ_local;
721 <  nMovingZMols = nMovingZMols_local;
721 >  nMovingZAtoms = nMovingZAtoms_local;
722   #endif
723  
724    force[0]= 0;
725    force[1]= 0;
726    force[2]= 0;
727 <  force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);
727 >  force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
728  
729 <  //modify the velocites of unconstrained molecules
729 >  //modify the forces of unconstrained molecules
730 >  int accessCount = 0;
731    for(int i = 0; i < unconsMols.size(); i++){
732      
733       Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
# Line 680 | Line 737 | template<typename T> void ZConstraint<T>::doZconstrain
737      
738    }      
739  
740 < //modify the velocities of moving z-constrained molecules
740 > //modify the forces of moving z-constrained molecules
741    for(int i = 0; i < zconsMols.size(); i++) {
742 <   if (states[i] == zcsMoving){
742 >    if (states[i] == zcsMoving){
743                  
744 <     Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
744 >      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
745  
746 <     for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
747 <       movingZAtoms[j]->addFrc(force);
748 <     }
746 >      for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)
747 >        movingZAtoms[j]->addFrc(force);
748 >    }
749    }
693 –
694 –  // apply negative to fixed z-constrained molecues;
695 –  force[0]= 0;
696 –  force[1]= 0;
697 –  force[2]= 0;
698 –
699 –  for(int i = 0; i < zconsMols.size(); i++){
700 –
701 –    if (states[i] == zcsFixed){  
750          
751 <      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
752 <      zconsAtoms = zconsMols[i]->getMyAtoms();  
705 <    
706 <      for(int j =0; j < nAtomOfCurZConsMol; j++) {
707 <        force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
708 <        zconsAtoms[j]->addFrc(force);
709 <      }
710 <                
711 <    }
712 <        
713 <  }
751 >  cout << "after substracting z-constraint force from moving molecuels "
752 >                  << "total force is " << calcTotalForce()  << endl;
753  
754   }
755  
# Line 769 | Line 808 | template<typename T> void ZConstraint<T>::doHarmonic()
808   template<typename T> void ZConstraint<T>::doHarmonic(){
809    double force[3];
810    double harmonicU;
811 +  double harmonicF;
812    double COM[3];
813    double diff;
814 +  double totalFZ;
815          
816    force[0] = 0;
817    force[1] = 0;
818    force[2] = 0;
819  
820 +  totalFZ = 0;
821 +
822    cout << "Moving Molecules" << endl;  
823    for(int i = 0; i < zconsMols.size(); i++) {
824  
825      if (states[i] == zcsMoving){
826 <      zconsMols[i]->getCOM(COM):
826 >      zconsMols[i]->getCOM(COM);
827        cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
828                  
829                  diff = COM[whichDirection] -zPos[i];
830                  
831        harmonicU = 0.5 * kz[i] * diff * diff;  
832 <                info->ltPot += harmonicU;
832 >                info->lrPot += harmonicU;
833  
834 <                force[whichDirection] = - kz[i] * diff / zconsMols[i]->getNAtoms();
834 >      harmonicF =  - kz[i] * diff;
835 >      totalFZ += harmonicF;
836 >
837 >       //
838 >                force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
839                  
840        Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
841  
# Line 797 | Line 844 | template<typename T> void ZConstraint<T>::doHarmonic()
844      }
845  
846    }
847 +        
848 +  force[0]= 0;
849 +  force[1]= 0;
850 +  force[2]= 0;
851 +  force[whichDirection] = -totalFZ /totNumOfUnconsAtoms;
852  
853 +  //modify the forces of unconstrained molecules
854 +  for(int i = 0; i < unconsMols.size(); i++){
855 +    
856 +     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
857 +    
858 +     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)          
859 +       unconsAtoms[j]->addFrc(force);    
860 +  }  
861 +
862   }
863 +
864 + template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
865 + {
866 +  double MVzOfMovingMols_local;
867 +  double MVzOfMovingMols;
868 +  double totalMassOfMovingZMols_local;
869 +  double totalMassOfMovingZMols;
870 +  double COMvel[3];
871 +      
872 +  MVzOfMovingMols_local = 0;
873 +  totalMassOfMovingZMols_local = 0;
874 +
875 +  for(int i =0; i < unconsMols.size(); i++){
876 +    unconsMols[i]->getCOMvel(COMvel);
877 +    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
878 +  }
879 +
880 +  for(int i = 0; i < zconsMols.size(); i++){
881 +
882 +    if (states[i] == zcsMoving){
883 +      zconsMols[i]->getCOMvel(COMvel);
884 +      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
885 +      totalMassOfMovingZMols_local += massOfZConsMols[i];              
886 +    }
887 +                
888 +  }
889 +
890 + #ifndef IS_MPI
891 +  MVzOfMovingMols = MVzOfMovingMols_local;
892 +  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
893 + #else
894 +  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
895 +  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
896 + #endif
897 +
898 +  double vzOfMovingMols;
899 +  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
900 +
901 +  return vzOfMovingMols;
902 + }
903 +
904 +
905 + template<typename T> double ZConstraint<T>::calcSysCOMVel()
906 + {
907 +  double COMvel[3];
908 +  double tempMVz = 0;
909 +                
910 +  for(int i =0 ; i < nMols; i++){
911 +    molecules[i].getCOMvel(COMvel);
912 +         tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection];
913 +  }
914 +
915 +  double massOfZCons_local;
916 +  double massOfZCons;
917 +
918 +  massOfZCons_local = 0;
919 +        
920 +  for(int i = 0; i < massOfZConsMols.size(); i++){
921 +    massOfZCons_local += massOfZConsMols[i];
922 +  }
923 + #ifndef IS_MPI
924 +  massOfZCons = massOfZCons_local;
925 + #else
926 +  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
927 + #endif
928 +
929 +  return tempMVz /(totalMassOfUncons + massOfZCons);
930 + }
931 +
932 + template<typename T> double ZConstraint<T>::calcTotalForce(){
933 +
934 +  double force[3];  
935 +  double totalForce_local;
936 +  double totalForce;
937 +
938 +  totalForce_local = 0;
939 +
940 +  for(int i = 0; i < nAtoms; i++){
941 +    atoms[i]->getFrc(force);
942 +    totalForce_local += force[whichDirection];
943 +  }
944 +
945 + #ifndef IS_MPI
946 +  totalForce = totalForce_local;
947 + #else
948 +  MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
949 + #endif
950 +
951 +  return totalForce;
952 +
953 + }

Diff Legend

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