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 696 by tim, Thu Aug 14 16:16:39 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 330 | Line 456 | template<typename T> void ZConstraint<T>::integrate(){
456   *
457   *
458   *
459 < */
334 <
335 <
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 359 | Line 496 | template<typename T> double ZConstraint<T>::calcZSys()
496    double totalMass;
497    double totalMZ_local;
498    double totalMZ;
362 –  double massOfUncons_local;
499    double massOfCurMol;
500    double COM[3];
501 <  
501 >        
502    totalMass_local = 0;
367 –  totalMass = 0;
503    totalMZ_local = 0;
369 –  totalMZ = 0;
370 –  massOfUncons_local = 0;
371 –    
504    
505    for(int i = 0; i < nMols; i++){
506      massOfCurMol = molecules[i].getTotalMass();
# Line 376 | Line 508 | template<typename T> double ZConstraint<T>::calcZSys()
508      
509      totalMass_local += massOfCurMol;
510      totalMZ_local += massOfCurMol * COM[whichDirection];
511 <    
380 <    if(isZConstraintMol(&molecules[i]) == -1){
381 <    
382 <      massOfUncons_local += massOfCurMol;
383 <    }  
384 <    
511 >
512    }
513 +
514    
387 –  
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);
392 < #else
518 > #else
519    totalMass = totalMass_local;
520    totalMZ = totalMZ_local;
521 <  totalMassOfUncons = massOfUncons_local;
396 < #endif  
521 > #endif  
522  
523    double zsys;
524    zsys = totalMZ / totalMass;
# Line 421 | 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 455 | 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++){
459 <
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 491 | 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 516 | Line 649 | template<typename T> void ZConstraint<T>::doZconstrain
649    double COMvel[3];  
650    double COM[3];
651    double force[3];
519 –  double zsys;
652  
521 –  int nMovingZMols_local;
522 –  int nMovingZMols;
653  
524 –  //constrain the molecules which do not reach the specified positions  
654  
655 <   zsys = calcZSys();
527 <   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;
659  
660    //calculate the total z-contrained force of fixed z-contrained molecules
661 +  cout << "Fixed Molecules" << endl;
662    for(int i = 0; i < zconsMols.size(); i++){
663                  
664      if (states[i] == zcsFixed){
# Line 544 | 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 579 | 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    }
592 –
593 –  // apply negative to fixed z-constrained molecues;
594 –  force[0]= 0;
595 –  force[1]= 0;
596 –  force[2]= 0;
597 –
598 –  for(int i = 0; i < zconsMols.size(); i++){
599 –
600 –    if (states[i] == zcsFixed){  
750          
751 <      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
752 <      zconsAtoms = zconsMols[i]->getMyAtoms();  
604 <    
605 <      for(int j =0; j < nAtomOfCurZConsMol; j++) {
606 <        force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
607 <        zconsAtoms[j]->addFrc(force);
608 <      }
609 <                
610 <    }
611 <        
612 <  }
751 >  cout << "after substracting z-constraint force from moving molecuels "
752 >                  << "total force is " << calcTotalForce()  << endl;
753  
754   }
755  
# Line 624 | Line 764 | template<typename T> bool ZConstraint<T>::checkZConsSt
764    for(int i =0; i < zconsMols.size(); i++){
765  
766      zconsMols[i]->getCOM(COM);
767 <    diff = fabs(COM[whichDirection] - ZPos[i]);  
768 <    if (  diff <= ztol && states[i] == zcsMoving){
767 >    diff = fabs(COM[whichDirection] - zPos[i]);  
768 >    if (  diff <= zconsTol && states[i] == zcsMoving){
769        states[i] = zcsFixed;
770          changed = true;
771      }
772 <    else if ( diff > ztol && states[i] == zcsFixed){
772 >    else if ( diff > zconsTol && states[i] == zcsFixed){
773        states[i] = zcsMoving;
774          changed = true;  
775      }
# Line 658 | Line 798 | template<typename T> bool ZConstraint<T>::haveMovingZM
798  
799    return false;
800    
801 < }
801 > }
802 >
803 > /**
804 >  *
805 >  *
806 >  */
807 >
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);
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->lrPot += harmonicU;
833 >
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 >
842 >       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
843 >         movingZAtoms[j]->addFrc(force);
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