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 676 by tim, Mon Aug 11 19:40:06 2003 UTC vs.
Revision 693 by tim, Wed Aug 13 19:21:53 2003 UTC

# Line 8 | Line 8 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
8  
9    //get properties from SimInfo
10    GenericData* data;
11 <  IndexData* index;
11 >  ZConsParaData* zConsParaData;
12    DoubleData* sampleTime;
13 +  DoubleData* tolerance;
14    StringData* filename;
15 <  
15 <  //retrieve index of z-constraint molecules
16 <  data = info->getProperty("zconsindex");
17 <  if(!data) {
15 >  double COM[3];
16  
17 <    sprintf( painCave.errMsg,
18 <               "ZConstraint error: If you use an ZConstraint\n"
19 <               " , you must set index of z-constraint molecules.\n");
20 <    painCave.isFatal = 1;
21 <    simError();  
24 <  }
25 <  else{
26 <    index = dynamic_cast<IndexData*>(data);
27 <    
28 <    if(!index){
17 >  //by default, the direction of constraint is z
18 >  // 0 --> x
19 >  // 1 --> y
20 >  // 2 --> z
21 >  whichDirection = 2;
22  
23 <      sprintf( painCave.errMsg,
24 <                 "ZConstraint error: Can not get property from SimInfo\n");
32 <      painCave.isFatal = 1;
33 <      simError();  
34 <    
35 <    }
36 <    else{
37 <          
38 <      indexOfAllZConsMols = index->getIndexData();
39 <      
40 <      //the maximum value of index is the last one(we sorted the index data in SimSetup.cpp)
41 <      int maxIndex;
42 <        int minIndex;
43 <      int totalNumMol;
44 <
45 <      minIndex = indexOfAllZConsMols[0];
46 <      if(minIndex < 0){
47 <        sprintf( painCave.errMsg,
48 <               "ZConstraint error: index is out of range\n");
49 <        painCave.isFatal = 1;
50 <        simError();
51 <        }
52 <          
53 <      maxIndex = indexOfAllZConsMols[indexOfAllZConsMols.size() - 1];
54 <
55 < #ifndef IS_MPI
56 <      totalNumMol = nMols;
57 < #else
58 <      totalNumMol = mpiSim->getTotNmol();  
59 < #endif      
60 <      
61 <      if(maxIndex > totalNumMol - 1){
62 <        sprintf( painCave.errMsg,
63 <               "ZConstraint error: index is out of range\n");
64 <        painCave.isFatal = 1;
65 <        simError();
66 <                
67 <      }
68 <      
69 <    }
70 <        
71 <  }
23 >  //estimate the force constant of harmonical potential
24 >  double Kb = 1.986E-3 ; //in kcal/K
25    
26 +  double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2;
27 +  zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox);
28 +
29    //retrieve sample time of z-contraint
30 <  data = info->getProperty("zconstime");
30 >  data = info->getProperty(ZCONSTIME_ID);
31    
32    if(!data) {
33        
# Line 99 | Line 55 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
55  
56    }
57    
102 –  
58    //retrieve output filename of z force
59 <  data = info->getProperty("zconsfilename");
59 >  data = info->getProperty(ZCONSFILENAME_ID);
60    if(!data) {
61  
62        
# Line 130 | Line 85 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
85      
86  
87    }
88 +
89 +  //retrieve tolerance for z-constraint molecuels
90 +  data = info->getProperty(ZCONSTOL_ID);
91    
92 +  if(!data) {
93 +      
94 +    sprintf( painCave.errMsg,
95 +               "ZConstraint error: can not get tolerance \n");
96 +    painCave.isFatal = 1;
97 +    simError();      
98 +  }
99 +  else{
100 +  
101 +    tolerance = dynamic_cast<DoubleData*>(data);
102 +    
103 +    if(!tolerance){
104 +
105 +      sprintf( painCave.errMsg,
106 +                 "ZConstraint error: Can not get property from SimInfo\n");
107 +      painCave.isFatal = 1;
108 +      simError();  
109 +      
110 +    }
111 +    else{
112 +      this->zconsTol = tolerance->getData();
113 +    }
114 +
115 +  }
116 +        
117 +  //retrieve index of z-constraint molecules
118 +  data = info->getProperty(ZCONSPARADATA_ID);
119 +  if(!data) {
120 +
121 +    sprintf( painCave.errMsg,
122 +               "ZConstraint error: If you use an ZConstraint\n"
123 +               " , you must set index of z-constraint molecules.\n");
124 +    painCave.isFatal = 1;
125 +    simError();  
126 +  }
127 +  else{
128 +  
129 +    zConsParaData = dynamic_cast<ZConsParaData*>(data);
130 +    
131 +    if(!zConsParaData){
132 +
133 +      sprintf( painCave.errMsg,
134 +                 "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
135 +      painCave.isFatal = 1;
136 +      simError();  
137 +    
138 +    }
139 +    else{
140 +      
141 +      parameters = zConsParaData->getData();
142 +
143 +      //check the range of zconsIndex
144 +      //and the minimum value of index is the first one (we already sorted the data)
145 +      //the maximum value of index is the last one
146 +
147 +      int maxIndex;
148 +      int minIndex;
149 +      int totalNumMol;
150 +
151 +      minIndex = (*parameters)[0].zconsIndex;
152 +      if(minIndex < 0){
153 +        sprintf( painCave.errMsg,
154 +               "ZConstraint error: index is out of range\n");
155 +        painCave.isFatal = 1;
156 +        simError();
157 +        }
158 +
159 +      maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
160 +
161 + #ifndef IS_MPI
162 +      totalNumMol = nMols;
163 + #else
164 +      totalNumMol = mpiSim->getTotNmol();  
165 + #endif      
166 +      
167 +      if(maxIndex > totalNumMol - 1){
168 +        sprintf( painCave.errMsg,
169 +               "ZConstraint error: index is out of range\n");
170 +        painCave.isFatal = 1;
171 +        simError();                  
172 +      }
173 +
174 +      //if user does not specify the zpos for the zconstraint molecule
175 +      //its initial z coordinate  will be used as default
176 +      for(int i = 0; i < parameters->size(); i++){
177 +
178 +              if(!(*parameters)[i].havingZPos){
179 +
180 + #ifndef IS_MPI
181 +            for(int j = 0; j < nMols; j++){
182 +              if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
183 +                 molecules[j].getCOM(COM);
184 +                          break;
185 +              }
186 +            }
187 + #else
188 +            //query which processor current zconstraint molecule belongs to
189 +           int *MolToProcMap;
190 +           int whichNode;
191 +                         double initZPos;
192 +           MolToProcMap = mpiSim->getMolToProcMap();
193 +           whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
194 +                          
195 +           //broadcast the zpos of current z-contraint molecule
196 +           //the node which contain this
197 +          
198 +           if (worldRank == whichNode ){
199 +                                                
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 +                                
206 +           }
207 +
208 +            MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);                          
209 + #endif
210 +            
211 +                 (*parameters)[i].zPos = COM[whichDirection];
212 +
213 +            sprintf( painCave.errMsg,
214 +                     "ZConstraint warningr: Does not specify zpos for z-constraint molecule "
215 +                     "initial z coornidate will be used \n");
216 +            painCave.isFatal = 0;
217 +            simError();  
218 +          
219 +              }
220 +            }
221 +                        
222 +    }//end if (!zConsParaData)
223 +  }//end  if (!data)
224 +            
225 + //  
226   #ifdef IS_MPI
227    update();
228   #else  
229    int searchResult;
138 –  double COM[3];
230        
231    for(int i = 0; i < nMols; i++){
232      
# Line 145 | Line 236 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
236      
237        zconsMols.push_back(&molecules[i]);      
238        massOfZConsMols.push_back(molecules[i].getTotalMass());  
239 +
240 +      zPos.push_back((*parameters)[searchResult].zPos);
241 +           kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
242        
243        molecules[i].getCOM(COM);
244      }
# Line 167 | 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 182 | 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 <
188 <  
303 >  // creat zconsWriter  
304    fzOut = new ZConsWriter(zconsOutput.c_str());  
305    
306    if(!fzOut){
# Line 217 | Line 332 | template<typename T> void ZConstraint<T>::update()
332    
333    zconsMols.clear();
334    massOfZConsMols.clear();
335 +  zPos.clear();
336 +  kz.clear();
337    
338    unconsMols.clear();
339    massOfUnconsMols.clear();
# Line 230 | Line 347 | template<typename T> void ZConstraint<T>::update()
347      if(index > -1){
348      
349        zconsMols.push_back(&molecules[i]);      
350 +      zPos.push_back((*parameters)[index].zPos);
351 +        kz.push_back((*parameters)[index].kRatio * zForceConst);
352 +                        
353        massOfZConsMols.push_back(molecules[i].getTotalMass());  
354        
355        molecules[i].getCOM(COM);
# Line 242 | 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 294 | Line 424 | template<typename T> int ZConstraint<T>::isZConstraint
424    index = mol->getGlobalIndex();
425    
426    low = 0;
427 <  high = indexOfAllZConsMols.size() - 1;
427 >  high = parameters->size() - 1;
428    
429    //Binary Search (we have sorted the array)  
430    while(low <= high){
431      mid = (low + high) /2;
432 <    if (indexOfAllZConsMols[mid] == index)
432 >    if ((*parameters)[mid].zconsIndex == index)
433        return mid;
434 <    else if (indexOfAllZConsMols[mid] > index )
434 >    else if ((*parameters)[mid].zconsIndex > index )
435         high = mid -1;
436      else    
437        low = mid + 1;
# Line 310 | Line 440 | template<typename T> int ZConstraint<T>::isZConstraint
440    return -1;
441   }
442  
313 – /**
314 – * Description:
315 – *  Reset the z coordinates
316 – */
443   template<typename T> void ZConstraint<T>::integrate(){
444    
445    //zero out the velocities of center of mass of unconstrained molecules
# Line 334 | 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 421 | 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 455 | 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++){
459 <
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 491 | 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 516 | Line 669 | template<typename T> void ZConstraint<T>::doZconstrain
669    double COMvel[3];  
670    double COM[3];
671    double force[3];
519 –  double zsys;
672  
673    int nMovingZMols_local;
674    int nMovingZMols;
675  
676    //constrain the molecules which do not reach the specified positions  
525 –
526 –   zsys = calcZSys();
527 –   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;
680  
681    //calculate the total z-contrained force of fixed z-contrained molecules
682 +  cout << "Fixed Molecules" << endl;
683    for(int i = 0; i < zconsMols.size(); i++){
684                  
685      if (states[i] == zcsFixed){
# Line 552 | 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 569 | 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 579 | 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 624 | Line 774 | template<typename T> bool ZConstraint<T>::checkZConsSt
774    for(int i =0; i < zconsMols.size(); i++){
775  
776      zconsMols[i]->getCOM(COM);
777 <    diff = fabs(COM[whichDirection] - ZPos[i]);  
778 <    if (  diff <= ztol && states[i] == zcsMoving){
777 >    diff = fabs(COM[whichDirection] - zPos[i]);  
778 >    if (  diff <= zconsTol && states[i] == zcsMoving){
779        states[i] = zcsFixed;
780          changed = true;
781      }
782 <    else if ( diff > ztol && states[i] == zcsFixed){
782 >    else if ( diff > zconsTol && states[i] == zcsFixed){
783        states[i] = zcsMoving;
784          changed = true;  
785      }
# Line 658 | Line 808 | template<typename T> bool ZConstraint<T>::haveMovingZM
808  
809    return false;
810    
811 < }
811 > }
812 >
813 > /**
814 >  *
815 >  *
816 >  */
817 >
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);
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->lrPot += harmonicU;
843 >
844 >      harmonicF =  - kz[i] * diff / zconsMols[i]->getNAtoms();
845 >                force[whichDirection] = harmonicF;
846 >      totalFZ += harmonicF;
847 >                
848 >      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
849 >
850 >       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
851 >         movingZAtoms[j]->addFrc(force);
852 >    }
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