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 693 by tim, Wed Aug 13 19:21:53 2003 UTC vs.
Revision 696 by tim, Thu Aug 14 16:16:39 2003 UTC

# Line 456 | Line 456 | template<typename T> void ZConstraint<T>::integrate(){
456   *
457   *
458   *
459 < */
460 <
461 <
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();
466 >    zeroOutVel();
467    
468    zsys = calcZSys();
469    cout << "---------------------------------------------------------------------" <<endl;
470 <  cout << "current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;      
471 <  cout << "before calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
472 <        
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 <      
480 >    
481    //use harmonical poteintial to move the molecules to the specified positions
482    if (haveMovingZMols())
483      this->doHarmonic();
484 <  
484 >
485 >  //cout <<      "after doHarmonic, totalForce is " << calcTotalForce() << endl;
486 >
487    fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
488 <  cout << "after calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
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 494 | Line 496 | template<typename T> double ZConstraint<T>::calcZSys()
496    double totalMass;
497    double totalMZ_local;
498    double totalMZ;
497 –  double massOfUncons_local;
499    double massOfCurMol;
500    double COM[3];
501 <  
501 >        
502    totalMass_local = 0;
502 –  totalMass = 0;
503    totalMZ_local = 0;
504 –  totalMZ = 0;
505 –  massOfUncons_local = 0;
506 –    
504    
505    for(int i = 0; i < nMols; i++){
506      massOfCurMol = molecules[i].getTotalMass();
# Line 511 | Line 508 | template<typename T> double ZConstraint<T>::calcZSys()
508      
509      totalMass_local += massOfCurMol;
510      totalMZ_local += massOfCurMol * COM[whichDirection];
511 <    
515 <    if(isZConstraintMol(&molecules[i]) == -1){
516 <    
517 <      massOfUncons_local += massOfCurMol;
518 <    }  
519 <    
511 >
512    }
513 +
514    
522 –  
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);
527 < #else
518 > #else
519    totalMass = totalMass_local;
520    totalMZ = totalMZ_local;
521 <  totalMassOfUncons = massOfUncons_local;
531 < #endif  
521 > #endif  
522  
523    double zsys;
524    zsys = totalMZ / totalMass;
# Line 557 | Line 547 | template<typename T> void ZConstraint<T>::zeroOutVel()
547    double COMvel[3];
548    double vel[3];
549  
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 –
550    //zero out the velocities of center of mass of fixed z-constrained molecules
551    
552    for(int i = 0; i < zconsMols.size(); i++){
# Line 575 | Line 554 | template<typename T> void ZConstraint<T>::zeroOutVel()
554      if (states[i] == zcsFixed){
555  
556             zconsMols[i]->getCOMvel(COMvel);      
557 <                cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
557 >                //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
558  
559        fixedZAtoms = zconsMols[i]->getMyAtoms();
560            
# Line 586 | Line 565 | template<typename T> void ZConstraint<T>::zeroOutVel()
565        }
566  
567                  zconsMols[i]->getCOMvel(COMvel);
568 <                cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
568 >                //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
569      }
570          
571    }
572  
573 <        cout << "before resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;    
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;
# Line 657 | Line 637 | template<typename T> void ZConstraint<T>::zeroOutVel()
637  
638   }
639  
640 <        cout << "after resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
640 >        cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() <<endl;
641  
642   }
643  
# Line 670 | Line 650 | template<typename T> void ZConstraint<T>::doZconstrain
650    double COM[3];
651    double force[3];
652  
673 –  int nMovingZMols_local;
674 –  int nMovingZMols;
653  
654 +
655    //constrain the molecules which do not reach the specified positions  
656      
657    //Zero Out the force of z-contrained molecules    
# Line 694 | 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;
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];
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 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 731 | Line 739 | template<typename T> void ZConstraint<T>::doZconstrain
739  
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    }
742 –
743 –  // apply negative to fixed z-constrained molecues;
744 –  force[0]= 0;
745 –  force[1]= 0;
746 –  force[2]= 0;
747 –
748 –  for(int i = 0; i < zconsMols.size(); i++){
749 –
750 –    if (states[i] == zcsFixed){  
750          
751 <      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
752 <      zconsAtoms = zconsMols[i]->getMyAtoms();  
754 <    
755 <      for(int j =0; j < nAtomOfCurZConsMol; j++) {
756 <        force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
757 <        zconsAtoms[j]->addFrc(force);
758 <      }
759 <                
760 <    }
761 <        
762 <  }
751 >  cout << "after substracting z-constraint force from moving molecuels "
752 >                  << "total force is " << calcTotalForce()  << endl;
753  
754   }
755  
# Line 841 | Line 831 | template<typename T> void ZConstraint<T>::doHarmonic()
831        harmonicU = 0.5 * kz[i] * diff * diff;  
832                  info->lrPot += harmonicU;
833  
834 <      harmonicF =  - kz[i] * diff / zconsMols[i]->getNAtoms();
845 <                force[whichDirection] = harmonicF;
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 852 | Line 844 | template<typename T> void ZConstraint<T>::doHarmonic()
844      }
845  
846    }
847 <
847 >        
848    force[0]= 0;
849    force[1]= 0;
850    force[2]= 0;
# Line 869 | Line 861 | template<typename T> double ZConstraint<T>::calcCOMVel
861  
862   }
863  
864 < template<typename T> double ZConstraint<T>::calcCOMVel()
864 > template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
865   {
866    double MVzOfMovingMols_local;
867    double MVzOfMovingMols;
# Line 910 | Line 902 | template<typename T> double ZConstraint<T>::calcCOMVel
902   }
903  
904  
905 < template<typename T> double ZConstraint<T>::calcCOMVel2()
905 > template<typename T> double ZConstraint<T>::calcSysCOMVel()
906   {
907    double COMvel[3];
908    double tempMVz = 0;
917 –  int index;
909                  
910    for(int i =0 ; i < nMols; i++){
911 <         index = isZConstraintMol(&molecules[i]);
912 <    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 <         }
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 <  return tempMVz /totalMassOfUncons;
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