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 736 by tim, Thu Aug 28 21:09:47 2003 UTC vs.
Revision 1093 by tim, Wed Mar 17 14:22:59 2004 UTC

# Line 1 | Line 1
1 < #include "Integrator.hpp"
2 < #include "simError.h"
3 < #include <cmath>
4 < template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff)
5 <                                    : T(theInfo, the_ff), fz(NULL), curZPos(NULL),
6 <                                 indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0)
7 < {
8 <
9 <  //get properties from SimInfo
10 <  GenericData* data;
11 <  ZConsParaData* zConsParaData;
12 <  DoubleData* sampleTime;
13 <  DoubleData* tolerance;
14 <  StringData* policy;
15 <  StringData* filename;
16 <  double COM[3];
17 <
18 <  //by default, the direction of constraint is z
19 <  // 0 --> x
20 <  // 1 --> y
21 <  // 2 --> z
22 <  whichDirection = 2;
23 <
24 <  //estimate the force constant of harmonical potential
25 <  double Kb = 1.986E-3 ; //in kcal/K
26 <  
27 <  double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2;
28 <  zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox);
29 <
30 <  //creat force substraction policy
31 <  data = info->getProperty(ZCONSFORCEPOLICY_ID);
32 <  if(!data){
33 <    sprintf( painCave.errMsg,
34 <               "ZConstraint Warning: User does not set force substraction policy, "
35 <               "PolicyByMass is used\n");
36 <    painCave.isFatal = 0;
37 <    simError();      
38 <
39 <    forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
40 <  }
41 <  else{
42 <    policy = dynamic_cast<StringData*>(data);
43 <    
44 <    if(!policy){
45 <      sprintf( painCave.errMsg,
46 <                 "ZConstraint Error: Convertion from GenericData to StringData failure, "
47 <                 "PolicyByMass is used\n");
48 <      painCave.isFatal = 0;
49 <      simError();      
50 <
51 <      forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
52 <    }
53 <    else{
54 <      if(policy->getData() == "BYNUMBER")
55 <        forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
56 <      else if(policy->getData() == "BYMASS")
57 <        forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
58 <      else{
59 <        sprintf( painCave.errMsg,
60 <                  "ZConstraint Warning: unknown force substraction policy, "
61 <                  "PolicyByMass is used\n");
62 <        painCave.isFatal = 0;
63 <        simError();      
64 <        forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
65 <      }  
66 <    }
67 <  }
68 <  
69 <  
70 <  //retrieve sample time of z-contraint
71 <  data = info->getProperty(ZCONSTIME_ID);
72 <  
73 <  if(!data) {
74 <      
75 <    sprintf( painCave.errMsg,
76 <               "ZConstraint error: If you use an ZConstraint\n"
77 <               " , you must set sample time.\n");
78 <    painCave.isFatal = 1;
79 <    simError();      
80 <  }
81 <  else{
82 <  
83 <    sampleTime = dynamic_cast<DoubleData*>(data);
84 <    
85 <    if(!sampleTime){
86 <
87 <      sprintf( painCave.errMsg,
88 <                 "ZConstraint error: Can not get property from SimInfo\n");
89 <      painCave.isFatal = 1;
90 <      simError();  
91 <      
92 <    }
93 <    else{
94 <      this->zconsTime = sampleTime->getData();
95 <    }
96 <
97 <  }
98 <  
99 <  //retrieve output filename of z force
100 <  data = info->getProperty(ZCONSFILENAME_ID);
101 <  if(!data) {
102 <
103 <      
104 <    sprintf( painCave.errMsg,
105 <               "ZConstraint error: If you use an ZConstraint\n"
106 <               " , you must set output filename of z-force.\n");
107 <    painCave.isFatal = 1;
108 <    simError();  
109 <
110 <  }
111 <  else{
112 <
113 <    filename = dynamic_cast<StringData*>(data);
114 <    
115 <    if(!filename){
116 <
117 <      sprintf( painCave.errMsg,
118 <                 "ZConstraint error: Can not get property from SimInfo\n");
119 <      painCave.isFatal = 1;
120 <      simError();  
121 <        
122 <    }
123 <    else{
124 <      this->zconsOutput = filename->getData();
125 <    }
126 <    
127 <  }
128 <
129 <  //retrieve tolerance for z-constraint molecuels
130 <  data = info->getProperty(ZCONSTOL_ID);
131 <  
132 <  if(!data) {
133 <      
134 <    sprintf( painCave.errMsg,
135 <               "ZConstraint error: can not get tolerance \n");
136 <    painCave.isFatal = 1;
137 <    simError();      
138 <  }
139 <  else{
140 <  
141 <    tolerance = dynamic_cast<DoubleData*>(data);
142 <    
143 <    if(!tolerance){
144 <
145 <      sprintf( painCave.errMsg,
146 <                 "ZConstraint error: Can not get property from SimInfo\n");
147 <      painCave.isFatal = 1;
148 <      simError();  
149 <      
150 <    }
151 <    else{
152 <      this->zconsTol = tolerance->getData();
153 <    }
154 <
155 <  }
156 <  
157 <  //retrieve index of z-constraint molecules
158 <  data = info->getProperty(ZCONSPARADATA_ID);
159 <  if(!data) {
160 <
161 <    sprintf( painCave.errMsg,
162 <               "ZConstraint error: If you use an ZConstraint\n"
163 <               " , you must set index of z-constraint molecules.\n");
164 <    painCave.isFatal = 1;
165 <    simError();  
166 <  }
167 <  else{
168 <  
169 <    zConsParaData = dynamic_cast<ZConsParaData*>(data);
170 <    
171 <    if(!zConsParaData){
172 <
173 <      sprintf( painCave.errMsg,
174 <                 "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
175 <      painCave.isFatal = 1;
176 <      simError();  
177 <    
178 <    }
179 <    else{
180 <      
181 <      parameters = zConsParaData->getData();
182 <
183 <      //check the range of zconsIndex
184 <      //and the minimum value of index is the first one (we already sorted the data)
185 <      //the maximum value of index is the last one
186 <
187 <      int maxIndex;
188 <      int minIndex;
189 <      int totalNumMol;
190 <
191 <      minIndex = (*parameters)[0].zconsIndex;
192 <      if(minIndex < 0){
193 <        sprintf( painCave.errMsg,
194 <               "ZConstraint error: index is out of range\n");
195 <        painCave.isFatal = 1;
196 <        simError();
197 <        }
198 <
199 <      maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
200 <
201 < #ifndef IS_MPI
202 <      totalNumMol = nMols;
203 < #else
204 <      totalNumMol = mpiSim->getTotNmol();  
205 < #endif      
206 <      
207 <      if(maxIndex > totalNumMol - 1){
208 <        sprintf( painCave.errMsg,
209 <               "ZConstraint error: index is out of range\n");
210 <        painCave.isFatal = 1;
211 <        simError();                  
212 <      }
213 <
214 <      //if user does not specify the zpos for the zconstraint molecule
215 <      //its initial z coordinate  will be used as default
216 <      for(int i = 0; i < parameters->size(); i++){
217 <
218 <        if(!(*parameters)[i].havingZPos){
219 < #ifndef IS_MPI
220 <          for(int j = 0; j < nMols; j++){
221 <            if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
222 <              molecules[j].getCOM(COM);
223 <        break;
224 <            }
225 <          }
226 < #else
227 <            //query which processor current zconstraint molecule belongs to
228 <          int *MolToProcMap;
229 <          int whichNode;
230 <    double initZPos;
231 <          MolToProcMap = mpiSim->getMolToProcMap();
232 <          whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
233 <        
234 <          //broadcast the zpos of current z-contraint molecule
235 <          //the node which contain this
236 <          
237 <          if (worldRank == whichNode ){
238 <            
239 <            for(int j = 0; j < nMols; j++)
240 <              if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
241 <                molecules[j].getCOM(COM);
242 <          break;
243 <              }
244 <        
245 <          }
246 <
247 <          MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);        
248 < #endif
249 <            
250 <          (*parameters)[i].zPos = COM[whichDirection];
251 <
252 <          sprintf( painCave.errMsg,
253 <                     "ZConstraint warning: Does not specify zpos for z-constraint molecule "
254 <                     "initial z coornidate will be used \n");
255 <           painCave.isFatal = 0;
256 <           simError();  
257 <    
258 <  }
259 <      }
260 <      
261 <    }//end if (!zConsParaData)
262 <  }//end  if (!data)
263 <            
264 < //  
265 < #ifdef IS_MPI
266 <  update();
267 < #else  
268 <  int searchResult;
269 <      
270 <  for(int i = 0; i < nMols; i++){
271 <    
272 <    searchResult = isZConstraintMol(&molecules[i]);
273 <    
274 <    if(searchResult > -1){
275 <    
276 <      zconsMols.push_back(&molecules[i]);      
277 <      massOfZConsMols.push_back(molecules[i].getTotalMass());  
278 <
279 <      zPos.push_back((*parameters)[searchResult].zPos);
280 < //       cout << "index: "<< (*parameters)[searchResult].zconsIndex
281 < //              <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
282 <      kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
283 <      
284 <      molecules[i].getCOM(COM);
285 <    }
286 <    else
287 <    {
288 <    
289 <      unconsMols.push_back(&molecules[i]);
290 <      massOfUnconsMols.push_back(molecules[i].getTotalMass());
291 <
292 <    }
293 <  }
294 <
295 <  fz = new double[zconsMols.size()];
296 <  curZPos = new double[zconsMols.size()];
297 <  indexOfZConsMols = new int [zconsMols.size()];
298 <
299 <  if(!fz || !curZPos || !indexOfZConsMols){
300 <    sprintf( painCave.errMsg,
301 <             "Memory allocation failure in class Zconstraint\n");
302 <    painCave.isFatal = 1;
303 <    simError();
304 <  }
305 <
306 <  //determine the states of z-constraint molecules
307 <  for(int i = 0; i < zconsMols.size(); i++){
308 <    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
309 <
310 <    zconsMols[i]->getCOM(COM);
311 <    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
312 <      states.push_back(zcsFixed);
313 <    else
314 <      states.push_back(zcsMoving);
315 <  }
316 <  
317 < #endif
318 <
319 <  //get total masss of unconstraint molecules
320 <  double totalMassOfUncons_local;
321 <  totalMassOfUncons_local = 0;
322 <  
323 <  for(int i = 0; i < unconsMols.size(); i++)
324 <    totalMassOfUncons_local += unconsMols[i]->getTotalMass();
325 <    
326 < #ifndef IS_MPI
327 <  totalMassOfUncons = totalMassOfUncons_local;
328 < #else
329 <  MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1,
330 <                      MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
331 < #endif
332 <
333 <
334 <  //get total number of unconstrained atoms
335 <  int nUnconsAtoms_local;
336 <  nUnconsAtoms_local = 0;
337 <  for(int i = 0; i < unconsMols.size(); i++)
338 <    nUnconsAtoms_local += unconsMols[i]->getNAtoms();
339 <    
340 < #ifndef IS_MPI
341 <  totNumOfUnconsAtoms = nUnconsAtoms_local;
342 < #else
343 <  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1,
344 <                      MPI_INT,MPI_SUM, MPI_COMM_WORLD);  
345 < #endif  
346 <
347 <  // creat zconsWriter  
348 <  fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);  
349 <  
350 <  if(!fzOut){
351 <    sprintf( painCave.errMsg,
352 <             "Memory allocation failure in class Zconstraint\n");
353 <    painCave.isFatal = 1;
354 <    simError();
355 <  }
356 <
357 <  forcePolicy->update();
358 < }
359 <
360 < template<typename T> ZConstraint<T>::~ZConstraint()
361 < {
362 <  if(fz)
363 <    delete[] fz;
364 <
365 <  if(curZPos)
366 <    delete[] curZPos;
367 <  
368 <  if(indexOfZConsMols)
369 <    delete[] indexOfZConsMols;
370 <  
371 <  if(fzOut)
372 <    delete fzOut;
373 <  
374 <  if(forcePolicy)
375 <    delete forcePolicy;
376 < }
377 <
378 <
379 < /**
380 < *
381 < */
382 <
383 < #ifdef IS_MPI
384 < template<typename T> void ZConstraint<T>::update()
385 < {
386 <  double COM[3];
387 <  int index;
388 <  
389 <  zconsMols.clear();
390 <  massOfZConsMols.clear();
391 <  zPos.clear();
392 <  kz.clear();
393 <  
394 <  unconsMols.clear();
395 <  massOfUnconsMols.clear();
396 <
397 <
398 <  //creat zconsMol and unconsMol lists
399 <  for(int i = 0; i < nMols; i++){
400 <    
401 <    index = isZConstraintMol(&molecules[i]);
402 <    
403 <    if(index > -1){
404 <    
405 <      zconsMols.push_back(&molecules[i]);      
406 <      zPos.push_back((*parameters)[index].zPos);
407 <      kz.push_back((*parameters)[index].kRatio * zForceConst);
408 <      
409 <      massOfZConsMols.push_back(molecules[i].getTotalMass());  
410 <      
411 <      molecules[i].getCOM(COM);
412 <    }
413 <    else
414 <    {
415 <    
416 <      unconsMols.push_back(&molecules[i]);
417 <      massOfUnconsMols.push_back(molecules[i].getTotalMass());
418 <
419 <    }
420 <  }
421 <
422 <  //determine the states of z-constraint molecules
423 <  for(int i = 0; i < zconsMols.size(); i++){
424 <     zconsMols[i]->getCOM(COM);
425 <      if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
426 <      states.push_back(zcsFixed);
427 <     else
428 <      states.push_back(zcsMoving);
429 <  }
430 <
431 <    
432 <  //The reason to declare fz and indexOfZconsMols as pointer to array is
433 <  // that we want to make the MPI communication simple
434 <  if(fz)
435 <    delete[] fz;
436 <  
437 <  if(curZPos)
438 <    delete[] curZPos;
439 <    
440 <  if(indexOfZConsMols)
441 <    delete[] indexOfZConsMols;
442 <    
443 <  if (zconsMols.size() > 0){
444 <    fz = new double[zconsMols.size()];
445 <   curZPos = new double[zconsMols.size()];
446 <    indexOfZConsMols =  new int[zconsMols.size()];
447 <    
448 <    if(!fz || !curZPos || !indexOfZConsMols){
449 <      sprintf( painCave.errMsg,
450 <               "Memory allocation failure in class Zconstraint\n");
451 <      painCave.isFatal = 1;
452 <      simError();
453 <    }
454 <        
455 <    for(int i = 0; i < zconsMols.size(); i++){
456 <      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
457 <    }
458 <
459 <  }
460 <  else{
461 <    fz = NULL;
462 <   curZPos = NULL;
463 <    indexOfZConsMols = NULL;
464 <  }
465 <  
466 <  //
467 <  forcePolicy->update();
468 <  
469 < }
470 <
471 < #endif
472 <
473 < /**
474 < *  Function Name: isZConstraintMol
475 < *  Parameter
476 < *    Molecule* mol
477 < *  Return value:
478 < *    -1, if the molecule is not z-constraint molecule,
479 < *    other non-negative values, its index in indexOfAllZConsMols vector
480 < */
481 <
482 < template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
483 < {
484 <  int index;
485 <  int low;
486 <  int high;
487 <  int mid;
488 <
489 <  index = mol->getGlobalIndex();
490 <  
491 <  low = 0;
492 <  high = parameters->size() - 1;
493 <  
494 <  //Binary Search (we have sorted the array)  
495 <  while(low <= high){
496 <    mid = (low + high) /2;
497 <    if ((*parameters)[mid].zconsIndex == index)
498 <      return mid;
499 <    else if ((*parameters)[mid].zconsIndex > index )
500 <       high = mid -1;
501 <    else    
502 <      low = mid + 1;
503 <  }
504 <  
505 <  return -1;
506 < }
507 <
508 < template<typename T> void ZConstraint<T>::integrate(){
509 <  
510 <  //zero out the velocities of center of mass of unconstrained molecules
511 <  //and the velocities of center of mass of every single z-constrained molecueles
512 <  zeroOutVel();
513 <
514 <  curZconsTime = zconsTime + info->getTime();
515 <  
516 <  T::integrate();
517 <
518 < }
519 <
520 <
521 < /**
522 < *
523 < *
524 < *
525 < *
526 < */
527 < template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
528 <  double zsys;
529 <  double COM[3];
530 <  double force[3];
531 <  double zSysCOMVel;
532 <
533 <  T::calcForce(calcPot, calcStress);
534 <
535 <  if (checkZConsState()){
536 <    
537 < #ifdef IS_MPI
538 <    if(worldRank == 0){
539 < #endif
540 < //       std::cerr << "\n"
541 < //              << "*******************************************\n"
542 < //              << " about to call zeroOutVel()\n"
543 < //              << "*******************************************\n"
544 < //              << "\n";
545 < #ifdef IS_MPI
546 <    }
547 < #endif
548 <    zeroOutVel();
549 <
550 < #ifdef IS_MPI
551 <    if(worldRank == 0){
552 < #endif
553 < //       std::cerr << "\n"
554 < //              << "*******************************************\n"
555 < //              << " finished zeroOutVel()\n"
556 < //              << "*******************************************\n"
557 < //              << "\n";
558 < #ifdef IS_MPI
559 <    }
560 < #endif
561 <    
562 <    forcePolicy->update();
563 <  }  
564 <  
565 <  zsys = calcZSys();
566 <  zSysCOMVel = calcSysCOMVel();
567 < #ifdef IS_MPI
568 <  if(worldRank == 0){
569 < #endif
570 < //     cout << "---------------------------------------------------------------------" <<endl;
571 < //     cout << "current time: " << info->getTime() << endl;
572 < //     cout << "center of mass at z: " << zsys << endl;    
573 < //     cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
574 <
575 < #ifdef IS_MPI
576 <  }
577 < #endif
578 <
579 <  //do zconstraint force;
580 <  if (haveFixedZMols())
581 <    this->doZconstraintForce();
582 <
583 <  //use harmonical poteintial to move the molecules to the specified positions
584 <  if (haveMovingZMols())
585 <    this->doHarmonic();
586 <
587 <  //write out forces and current positions of z-constraint molecules
588 <  if(info->getTime() >= curZconsTime){    
589 <   for(int i = 0; i < zconsMols.size(); i++){
590 <      zconsMols[i]->getCOM(COM);
591 <    curZPos[i] = COM[whichDirection];
592 <
593 <    //if the z-constraint molecule is still moving, just record its force
594 <    if(states[i] == zcsMoving){
595 <         fz[i] = 0;
596 <      Atom** movingZAtoms;
597 <      movingZAtoms = zconsMols[i]->getMyAtoms();
598 <      for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
599 <           movingZAtoms[j]->getFrc(force);
600 <           fz[i] += force[whichDirection];
601 <      }
602 <     }
603 <   }
604 <    fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos);
605 <   curZconsTime += zconsTime;
606 <  }
607 <
608 <  zSysCOMVel = calcSysCOMVel();  
609 < #ifdef IS_MPI
610 <  if(worldRank == 0){
611 < #endif
612 < //    cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
613 < #ifdef IS_MPI
614 <  }
615 < #endif
616 < }
617 <
618 <
619 < /**
620 < *
621 < */
622 <
623 < template<typename T> double ZConstraint<T>::calcZSys()
624 < {
625 <  //calculate reference z coordinate for z-constraint molecules
626 <  double totalMass_local;
627 <  double totalMass;
628 <  double totalMZ_local;
629 <  double totalMZ;
630 <  double massOfCurMol;
631 <  double COM[3];
632 <  
633 <  totalMass_local = 0;
634 <  totalMZ_local = 0;
635 <  
636 <  for(int i = 0; i < nMols; i++){
637 <    massOfCurMol = molecules[i].getTotalMass();
638 <    molecules[i].getCOM(COM);
639 <    
640 <    totalMass_local += massOfCurMol;
641 <    totalMZ_local += massOfCurMol * COM[whichDirection];
642 <
643 <  }
644 <
645 <  
646 < #ifdef IS_MPI  
647 <  MPI_Allreduce(&totalMass_local, &totalMass, 1,
648 <                      MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
649 <  MPI_Allreduce(&totalMZ_local, &totalMZ, 1,
650 <                      MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
651 < #else
652 <  totalMass = totalMass_local;
653 <  totalMZ = totalMZ_local;
654 < #endif  
655 <
656 <  double zsys;
657 <  zsys = totalMZ / totalMass;
658 <
659 <  return zsys;
660 < }
661 <
662 < /**
663 < *
664 < */
665 < template<typename T> void ZConstraint<T>::thermalize( void ){
666 <
667 <  T::thermalize();
668 <  zeroOutVel();
669 < }
670 <
671 < /**
672 < *
673 < */
674 <
675 < template<typename T> void ZConstraint<T>::zeroOutVel(){
676 <
677 <  Atom** fixedZAtoms;  
678 <  double COMvel[3];
679 <  double vel[3];
680 <  double zSysCOMVel;
681 <
682 <  //zero out the velocities of center of mass of fixed z-constrained molecules
683 <  
684 <  for(int i = 0; i < zconsMols.size(); i++){
685 <
686 <    if (states[i] == zcsFixed){
687 <
688 <     zconsMols[i]->getCOMvel(COMvel);      
689 <    //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
690 <
691 <      fixedZAtoms = zconsMols[i]->getMyAtoms();
692 <    
693 <      for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
694 <        fixedZAtoms[j]->getVel(vel);
695 <       vel[whichDirection] -= COMvel[whichDirection];
696 <       fixedZAtoms[j]->setVel(vel);
697 <      }
698 <
699 <    zconsMols[i]->getCOMvel(COMvel);
700 <    //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
701 <    }
702 <  
703 <  }
704 <
705 <    //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
706 <
707 <  zSysCOMVel = calcSysCOMVel();
708 < #ifdef IS_MPI
709 <  if(worldRank == 0){
710 < #endif
711 < //     cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
712 < #ifdef IS_MPI
713 <  }
714 < #endif
715 <      
716 <  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
717 <  double MVzOfMovingMols_local;
718 <  double MVzOfMovingMols;
719 <  double totalMassOfMovingZMols_local;
720 <  double totalMassOfMovingZMols;
721 <      
722 <  MVzOfMovingMols_local = 0;
723 <  totalMassOfMovingZMols_local = 0;
724 <
725 <  for(int i =0; i < unconsMols.size(); i++){
726 <    unconsMols[i]->getCOMvel(COMvel);
727 <    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
728 <  }
729 <
730 <  for(int i = 0; i < zconsMols.size(); i++){
731 <    if (states[i] == zcsMoving){
732 <      zconsMols[i]->getCOMvel(COMvel);
733 <      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
734 <      totalMassOfMovingZMols_local += massOfZConsMols[i];      
735 <    }
736 <    
737 <  }
738 <
739 < #ifndef IS_MPI
740 <  MVzOfMovingMols = MVzOfMovingMols_local;
741 <  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
742 < #else
743 <  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
744 <  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
745 < #endif
746 <
747 <  double vzOfMovingMols;
748 <  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
749 <
750 <  //modify the velocites of unconstrained molecules  
751 <  Atom** unconsAtoms;
752 <  for(int i = 0; i < unconsMols.size(); i++){
753 <  
754 <    unconsAtoms = unconsMols[i]->getMyAtoms();
755 <    for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
756 <      unconsAtoms[j]->getVel(vel);
757 <      vel[whichDirection] -= vzOfMovingMols;
758 <      unconsAtoms[j]->setVel(vel);
759 <    }
760 <  
761 <  }  
762 <
763 <  //modify the velocities of moving z-constrained molecuels
764 <  Atom** movingZAtoms;
765 <  for(int i = 0; i < zconsMols.size(); i++){
766 <
767 <    if (states[i] ==zcsMoving){
768 <  
769 <      movingZAtoms = zconsMols[i]->getMyAtoms();
770 <     for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
771 <        movingZAtoms[j]->getVel(vel);
772 <        vel[whichDirection] -= vzOfMovingMols;
773 <       movingZAtoms[j]->setVel(vel);
774 <    }
775 <    
776 <   }
777 <
778 < }
779 <
780 <
781 <  zSysCOMVel = calcSysCOMVel();
782 < #ifdef IS_MPI
783 <  if(worldRank == 0){
784 < #endif
785 < //     cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
786 < #ifdef IS_MPI
787 <  }
788 < #endif
789 <
790 < }
791 <
792 < /**
793 < *
794 < */
795 <
796 < template<typename T> void ZConstraint<T>::doZconstraintForce(){
797 <
798 <  Atom** zconsAtoms;
799 <  double totalFZ;
800 <  double totalFZ_local;
801 <  double COMvel[3];  
802 <  double COM[3];
803 <  double force[3];
804 <
805 <  //constrain the molecules which do not reach the specified positions  
806 <    
807 <  //Zero Out the force of z-contrained molecules    
808 <  totalFZ_local = 0;
809 <
810 <  //calculate the total z-contrained force of fixed z-contrained molecules
811 <
812 <  for(int i = 0; i < zconsMols.size(); i++){
813 <    
814 <    if (states[i] == zcsFixed){
815 <      
816 <      zconsMols[i]->getCOM(COM);
817 <      zconsAtoms = zconsMols[i]->getMyAtoms();  
818 <
819 <      fz[i] = 0;      
820 <      for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
821 <        zconsAtoms[j]->getFrc(force);
822 <        fz[i] += force[whichDirection];      
823 <      }
824 <      totalFZ_local += fz[i];
825 <
826 <      //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
827 <      //      <<"\tcurrent zpos: " << COM[whichDirection]
828 <      //       << "\tcurrent fz: " <<fz[i] << endl;
829 <
830 <
831 <    }
832 <    
833 <  }
834 <
835 <  //calculate total z-constraint force
836 < #ifdef IS_MPI
837 <  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
838 < #else
839 <  totalFZ = totalFZ_local;
840 < #endif
841 <
842 <  
843 <  // apply negative to fixed z-constrained molecues;
844 <  force[0]= 0;
845 <  force[1]= 0;
846 <  force[2]= 0;
847 <
848 <  for(int i = 0; i < zconsMols.size(); i++){
849 <
850 <    if (states[i] == zcsFixed){  
851 <  
852 <      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
853 <      zconsAtoms = zconsMols[i]->getMyAtoms();  
854 <    
855 <      for(int j =0; j < nAtomOfCurZConsMol; j++) {
856 <        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
857 <        force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
858 <        zconsAtoms[j]->addFrc(force);
859 <      }
860 <    
861 <    }
862 <  
863 <  }
864 <
865 < //   cout << "after zero out z-constraint force on fixed z-constraint molecuels "
866 < //        << "total force is " << calcTotalForce() << endl;
867 <
868 <  force[0]= 0;
869 <  force[1]= 0;
870 <  force[2]= 0;
871 <
872 <  //modify the forces of unconstrained molecules
873 <  for(int i = 0; i < unconsMols.size(); i++){
874 <    
875 <     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
876 <    
877 <     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){          
878 <       //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
879 <       force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
880 <       unconsAtoms[j]->addFrc(force);
881 <     }
882 <    
883 <  }      
884 <
885 < //modify the forces of moving z-constrained molecules
886 <  for(int i = 0; i < zconsMols.size(); i++) {
887 <    if (states[i] == zcsMoving){
888 <    
889 <      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
890 <
891 <      for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
892 <        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
893 <        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
894 <        movingZAtoms[j]->addFrc(force);
895 <      }
896 <    }
897 <  }
898 <
899 <  //cout << "after substracting z-constraint force from moving molecuels "
900 <  //      << "total force is " << calcTotalForce()  << endl;
901 <
902 < }
903 <
904 < /**
905 <  *
906 <  *
907 <  */
908 <
909 < template<typename T> void ZConstraint<T>::doHarmonic(){
910 <  double force[3];
911 <  double harmonicU;
912 <  double harmonicF;
913 <  double COM[3];
914 <  double diff;
915 <  double totalFZ_local;
916 <  double totalFZ;
917 <  
918 <  force[0] = 0;
919 <  force[1] = 0;
920 <  force[2] = 0;
921 <
922 <  totalFZ_local = 0;
923 <
924 <  for(int i = 0; i < zconsMols.size(); i++) {
925 <
926 <    if (states[i] == zcsMoving){
927 <      zconsMols[i]->getCOM(COM);
928 < //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
929 < //         << "\tcurrent zpos: " << COM[whichDirection] << endl;
930 <
931 <      diff = COM[whichDirection] -zPos[i];
932 <    
933 <      harmonicU = 0.5 * kz[i] * diff * diff;  
934 <      info->lrPot += harmonicU;
935 <
936 <      harmonicF =  - kz[i] * diff;
937 <      totalFZ_local += harmonicF;
938 <
939 <       //adjust force
940 <    
941 <      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
942 <
943 <       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){          
944 <        //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
945 <        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
946 <         movingZAtoms[j]->addFrc(force);
947 <       }
948 <    }
949 <
950 <  }
951 <
952 < #ifndef IS_MPI
953 <  totalFZ = totalFZ_local;
954 < #else
955 <  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
956 < #endif
957 <
958 <  force[0]= 0;
959 <  force[1]= 0;
960 <  force[2]= 0;
961 <
962 <  //modify the forces of unconstrained molecules
963 <  for(int i = 0; i < unconsMols.size(); i++){
964 <    
965 <     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
966 <    
967 <     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){          
968 <       //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
969 <       force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
970 <       unconsAtoms[j]->addFrc(force);    
971 <     }
972 <  }  
973 <
974 < }
975 <
976 < /**
977 < *
978 < */
979 <
980 < template<typename T> bool ZConstraint<T>::checkZConsState(){
981 <  double COM[3];
982 <  double diff;
983 <  
984 <  int changed_local;
985 <  int changed;
986 <  
987 <  changed_local = 0;
988 <  
989 <  for(int i =0; i < zconsMols.size(); i++){
990 <
991 <    zconsMols[i]->getCOM(COM);
992 <    diff = fabs(COM[whichDirection] - zPos[i]);  
993 <    if (  diff <= zconsTol && states[i] == zcsMoving){
994 <      states[i] = zcsFixed;
995 <     changed_local = 1;
996 <    }
997 <    else if ( diff > zconsTol && states[i] == zcsFixed){
998 <      states[i] = zcsMoving;
999 <     changed_local = 1;  
1000 <    }
1001 <  
1002 <  }
1003 <
1004 < #ifndef IS_MPI
1005 <  changed =changed_local;
1006 < #else
1007 <  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1008 < #endif
1009 <
1010 <  return (changed > 0);
1011 <
1012 < }
1013 <
1014 < template<typename T> bool ZConstraint<T>::haveFixedZMols(){
1015 <
1016 <  int havingFixed_local;
1017 <  int havingFixed;
1018 <
1019 <  havingFixed_local = 0;
1020 <
1021 <  for(int i = 0; i < zconsMols.size(); i++)
1022 <    if (states[i] == zcsFixed){
1023 <      havingFixed_local = 1;
1024 <    break;
1025 <    }
1026 <
1027 < #ifndef IS_MPI
1028 <  havingFixed = havingFixed_local;
1029 < #else
1030 <  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1031 < #endif
1032 <
1033 <  return (havingFixed > 0);
1034 < }
1035 <
1036 <
1037 < /**
1038 < *
1039 < */
1040 < template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1041 <
1042 <  int havingMoving_local;
1043 <  int havingMoving;
1044 <
1045 <  havingMoving_local = 0;
1046 <
1047 <  for(int i = 0; i < zconsMols.size(); i++)
1048 <    if (states[i] == zcsMoving){
1049 <      havingMoving_local = 1;
1050 <    break;
1051 <    }
1052 <
1053 < #ifndef IS_MPI
1054 <  havingMoving = havingMoving_local;
1055 < #else
1056 <  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1057 < #endif
1058 <
1059 <  return (havingMoving > 0);
1060 <  
1061 < }
1062 <
1063 < /**
1064 < *
1065 < */
1066 <
1067 < template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
1068 < {
1069 <  double MVzOfMovingMols_local;
1070 <  double MVzOfMovingMols;
1071 <  double totalMassOfMovingZMols_local;
1072 <  double totalMassOfMovingZMols;
1073 <  double COMvel[3];
1074 <      
1075 <  MVzOfMovingMols_local = 0;
1076 <  totalMassOfMovingZMols_local = 0;
1077 <
1078 <  for(int i =0; i < unconsMols.size(); i++){
1079 <    unconsMols[i]->getCOMvel(COMvel);
1080 <    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
1081 <  }
1082 <
1083 <  for(int i = 0; i < zconsMols.size(); i++){
1084 <
1085 <    if (states[i] == zcsMoving){
1086 <      zconsMols[i]->getCOMvel(COMvel);
1087 <      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
1088 <      totalMassOfMovingZMols_local += massOfZConsMols[i];      
1089 <    }
1090 <    
1091 <  }
1092 <
1093 < #ifndef IS_MPI
1094 <  MVzOfMovingMols = MVzOfMovingMols_local;
1095 <  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1096 < #else
1097 <  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1098 <  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
1099 < #endif
1100 <
1101 <  double vzOfMovingMols;
1102 <  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
1103 <
1104 <  return vzOfMovingMols;
1105 < }
1106 <
1107 < /**
1108 < *
1109 < */
1110 <
1111 < template<typename T> double ZConstraint<T>::calcSysCOMVel()
1112 < {
1113 <  double COMvel[3];
1114 <  double tempMVz_local;
1115 <  double tempMVz;
1116 <  double massOfZCons_local;
1117 <  double massOfZCons;
1118 <
1119 <
1120 < tempMVz_local = 0;
1121 <
1122 <  for(int i =0 ; i < nMols; i++){
1123 <    molecules[i].getCOMvel(COMvel);
1124 <   tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
1125 <  }
1126 <
1127 <  massOfZCons_local = 0;
1128 <  
1129 <  for(int i = 0; i < massOfZConsMols.size(); i++){
1130 <    massOfZCons_local += massOfZConsMols[i];
1131 <  }
1132 < #ifndef IS_MPI
1133 <  massOfZCons = massOfZCons_local;
1134 <  tempMVz = tempMVz_local;
1135 < #else
1136 <  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1137 <  MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1138 < #endif
1139 <
1140 <  return tempMVz /(totalMassOfUncons + massOfZCons);
1141 < }
1142 <
1143 < /**
1144 < *
1145 < */
1146 <
1147 < template<typename T> double ZConstraint<T>::calcTotalForce(){
1148 <
1149 <  double force[3];  
1150 <  double totalForce_local;
1151 <  double totalForce;
1152 <
1153 <  totalForce_local = 0;
1154 <
1155 <  for(int i = 0; i < nAtoms; i++){
1156 <    atoms[i]->getFrc(force);
1157 <    totalForce_local += force[whichDirection];
1158 <  }
1159 <
1160 < #ifndef IS_MPI
1161 <  totalForce = totalForce_local;
1162 < #else
1163 <  MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1164 < #endif
1165 <
1166 <  return totalForce;
1167 <
1168 < }
1169 <
1170 < /**
1171 < *
1172 < */
1173 <
1174 < template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1175 <  //calculate the number of atoms of moving z-constrained molecules
1176 <  int nMovingZAtoms_local;
1177 <  int nMovingZAtoms;
1178 <  
1179 <  nMovingZAtoms_local = 0;
1180 <  for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1181 <    if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1182 <     nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1183 <  
1184 < #ifdef IS_MPI
1185 <  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1186 < #else
1187 <  nMovingZAtoms = nMovingZAtoms_local;
1188 < #endif
1189 <  totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1190 <
1191 < #ifdef IS_MPI
1192 <  if(worldRank == 0){
1193 < #endif
1194 < //    std::cerr << "\n"
1195 < //            << "*******************************************\n"
1196 < //            << " fiished Policy by numbr()\n"
1197 < //            << "*******************************************\n"
1198 < //            << "\n";
1199 < #ifdef IS_MPI
1200 <  }
1201 < #endif
1202 < }
1203 <
1204 < template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1205 <  return totalForce / mol->getNAtoms();
1206 < }
1207 <
1208 < template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, double totalForce){
1209 <  return totalForce / totNumOfMovingAtoms;
1210 < }
1211 <
1212 < template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1213 <    return totalForce / mol->getNAtoms();
1214 < }
1215 <
1216 < template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, double totalForce){
1217 <  return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1218 < }
1219 <
1220 < /**
1221 < *
1222 < */
1223 <
1224 < template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1225 <  //calculate the number of atoms of moving z-constrained molecules
1226 <  double massOfMovingZAtoms_local;
1227 <  double massOfMovingZAtoms;
1228 <  
1229 <  massOfMovingZAtoms_local = 0;
1230 <  for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1231 <    if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1232 <     massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1233 <  
1234 < #ifdef IS_MPI
1235 <  MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1236 < #else
1237 <  massOfMovingZAtoms = massOfMovingZAtoms_local;
1238 < #endif
1239 <  totMassOfMovingAtoms = massOfMovingZAtoms_local + zconsIntegrator->totalMassOfUncons;
1240 < }
1241 <
1242 < template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1243 <  return totalForce * atom->getMass() / mol->getTotalMass();
1244 < }
1245 <
1246 < template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols( Atom* atom, double totalForce){
1247 <    return totalForce * atom->getMass() / totMassOfMovingAtoms;
1248 < }
1249 <
1250 < template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1251 <  return totalForce * atom->getMass() / mol->getTotalMass();
1252 < }
1253 <
1254 < template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, double totalForce){
1255 <    return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1256 < }
1257 <
1 > #include "Integrator.hpp"
2 > #include "simError.h"
3 > #include <math.h>
4 >
5 > const double INFINITE_TIME = 10e30;
6 > template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo,
7 >                                                 ForceFields* the_ff): T(theInfo, the_ff),
8 >                                                                       fzOut(NULL),
9 >                                                                       curZconsTime(0),
10 >                                                                       forcePolicy(NULL),
11 >                                                                       usingSMD(false),
12 >                                                                       hasZConsGap(false){
13 >  //get properties from SimInfo
14 >  GenericData* data;
15 >  ZConsParaData* zConsParaData;
16 >  DoubleData* sampleTime;
17 >  DoubleData* tolerance;
18 >  DoubleData* gap;
19 >  DoubleData* fixtime;
20 >  StringData* policy;
21 >  StringData* filename;
22 >  IntData* smdFlag;
23 >  double COM[3];
24 >
25 >  //by default, the direction of constraint is z
26 >  // 0 --> x
27 >  // 1 --> y
28 >  // 2 --> z
29 >  whichDirection = 2;
30 >
31 >  //estimate the force constant of harmonical potential
32 >  double Kb = 1.986E-3 ; //in kcal/K
33 >
34 >  double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /
35 >                            2;
36 >  zForceConst = Kb * info->target_temp / (halfOfLargestBox * halfOfLargestBox);
37 >
38 >  //creat force Subtraction policy
39 >  data = info->getProperty(ZCONSFORCEPOLICY_ID);
40 >  if (!data){
41 >    sprintf(painCave.errMsg,
42 >            "ZConstraint Warning: User does not set force Subtraction policy, "
43 >            "PolicyByMass is used\n");
44 >    painCave.isFatal = 0;
45 >    simError();      
46 >
47 >    forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
48 >  }
49 >  else{
50 >    policy = dynamic_cast<StringData*>(data);
51 >
52 >    if (!policy){
53 >      sprintf(painCave.errMsg,
54 >              "ZConstraint Error: Convertion from GenericData to StringData failure, "
55 >              "PolicyByMass is used\n");
56 >      painCave.isFatal = 0;
57 >      simError();      
58 >
59 >      forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
60 >    }
61 >    else{
62 >      if (policy->getData() == "BYNUMBER")
63 >        forcePolicy = (ForceSubtractionPolicy *) new PolicyByNumber(this);
64 >      else if (policy->getData() == "BYMASS")
65 >        forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
66 >      else{
67 >        sprintf(painCave.errMsg,
68 >                "ZConstraint Warning: unknown force Subtraction policy, "
69 >                "PolicyByMass is used\n");
70 >        painCave.isFatal = 0;
71 >        simError();      
72 >        forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
73 >      }
74 >    }
75 >  }
76 >
77 >
78 >  //retrieve sample time of z-contraint
79 >  data = info->getProperty(ZCONSTIME_ID);
80 >
81 >  if (!data){
82 >    sprintf(painCave.errMsg,
83 >            "ZConstraint error: If you use an ZConstraint\n"
84 >            " , you must set sample time.\n");
85 >    painCave.isFatal = 1;
86 >    simError();
87 >  }
88 >  else{
89 >    sampleTime = dynamic_cast<DoubleData*>(data);
90 >
91 >    if (!sampleTime){
92 >      sprintf(painCave.errMsg,
93 >              "ZConstraint error: Can not get property from SimInfo\n");
94 >      painCave.isFatal = 1;
95 >      simError();
96 >    }
97 >    else{
98 >      this->zconsTime = sampleTime->getData();
99 >    }
100 >  }
101 >
102 >  //retrieve output filename of z force
103 >  data = info->getProperty(ZCONSFILENAME_ID);
104 >  if (!data){
105 >    sprintf(painCave.errMsg,
106 >            "ZConstraint error: If you use an ZConstraint\n"
107 >            " , you must set output filename of z-force.\n");
108 >    painCave.isFatal = 1;
109 >    simError();
110 >  }
111 >  else{
112 >    filename = dynamic_cast<StringData*>(data);
113 >
114 >    if (!filename){
115 >      sprintf(painCave.errMsg,
116 >              "ZConstraint error: Can not get property from SimInfo\n");
117 >      painCave.isFatal = 1;
118 >      simError();
119 >    }
120 >    else{
121 >      this->zconsOutput = filename->getData();
122 >    }
123 >  }
124 >
125 >  //retrieve tolerance for z-constraint molecuels
126 >  data = info->getProperty(ZCONSTOL_ID);
127 >
128 >  if (!data){
129 >    sprintf(painCave.errMsg, "ZConstraint error: can not get tolerance \n");
130 >    painCave.isFatal = 1;
131 >    simError();
132 >  }
133 >  else{
134 >    tolerance = dynamic_cast<DoubleData*>(data);
135 >
136 >    if (!tolerance){
137 >      sprintf(painCave.errMsg,
138 >              "ZConstraint error: Can not get property from SimInfo\n");
139 >      painCave.isFatal = 1;
140 >      simError();
141 >    }
142 >    else{
143 >      this->zconsTol = tolerance->getData();
144 >    }
145 >  }
146 >
147 >  //quick hack here
148 >  data = info->getProperty(ZCONSGAP_ID);
149 >
150 >  if (data){
151 >    gap = dynamic_cast<DoubleData*>(data);
152 >
153 >    if (!gap){
154 >      sprintf(painCave.errMsg,
155 >              "ZConstraint error: Can not get property from SimInfo\n");
156 >      painCave.isFatal = 1;
157 >      simError();
158 >    }
159 >    else{
160 >      this->hasZConsGap = true;
161 >      this->zconsGap = gap->getData();
162 >    }
163 >  }
164 >
165 >
166 >
167 >  data = info->getProperty(ZCONSFIXTIME_ID);
168 >
169 >  if (data){
170 >    fixtime = dynamic_cast<DoubleData*>(data);
171 >    if (!fixtime){
172 >      sprintf(painCave.errMsg,
173 >              "ZConstraint error: Can not get zconsFixTime from SimInfo\n");
174 >      painCave.isFatal = 1;
175 >      simError();
176 >    }
177 >    else{
178 >      this->zconsFixTime = fixtime->getData();
179 >    }
180 >  }
181 >  else if(hasZConsGap){
182 >      sprintf(painCave.errMsg,
183 >              "ZConstraint error: must set fixtime if already set zconsGap\n");
184 >      painCave.isFatal = 1;
185 >      simError();
186 >  }
187 >
188 >
189 >
190 >  data = info->getProperty(ZCONSUSINGSMD_ID);
191 >
192 >  if (data){
193 >    smdFlag = dynamic_cast<IntData*>(data);
194 >
195 >    if (!smdFlag){
196 >      sprintf(painCave.errMsg,
197 >              "ZConstraint error: Can not get property from SimInfo\n");
198 >      painCave.isFatal = 1;
199 >      simError();
200 >    }
201 >    else{
202 >      this->usingSMD= smdFlag->getData() ? true : false;
203 >    }
204 >
205 >  }
206 >
207 >
208 >
209 >  //retrieve index of z-constraint molecules
210 >  data = info->getProperty(ZCONSPARADATA_ID);
211 >  if (!data){
212 >    sprintf(painCave.errMsg,
213 >            "ZConstraint error: If you use an ZConstraint\n"
214 >            " , you must set index of z-constraint molecules.\n");
215 >    painCave.isFatal = 1;
216 >    simError();
217 >  }
218 >  else{
219 >    zConsParaData = dynamic_cast<ZConsParaData*>(data);
220 >
221 >    if (!zConsParaData){
222 >      sprintf(painCave.errMsg,
223 >              "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
224 >      painCave.isFatal = 1;
225 >      simError();
226 >    }
227 >    else{
228 >      parameters = zConsParaData->getData();
229 >
230 >      //check the range of zconsIndex
231 >      //and the minimum value of index is the first one (we already sorted the data)
232 >      //the maximum value of index is the last one
233 >
234 >      int maxIndex;
235 >      int minIndex;
236 >      int totalNumMol;
237 >
238 >      minIndex = (*parameters)[0].zconsIndex;
239 >      if (minIndex < 0){
240 >        sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n");
241 >        painCave.isFatal = 1;
242 >        simError();
243 >      }
244 >
245 >      maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
246 >
247 > #ifndef IS_MPI
248 >      totalNumMol = nMols;
249 > #else
250 >      totalNumMol = mpiSim->getTotNmol();  
251 > #endif      
252 >
253 >      if (maxIndex > totalNumMol - 1){
254 >        sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n");
255 >        painCave.isFatal = 1;
256 >        simError();
257 >      }
258 >
259 >      //if user does not specify the zpos for the zconstraint molecule
260 >      //its initial z coordinate  will be used as default
261 >      for (int i = 0; i < (int) (parameters->size()); i++){
262 >        if (!(*parameters)[i].havingZPos){
263 > #ifndef IS_MPI
264 >          for (int j = 0; j < nMols; j++){
265 >            if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
266 >              molecules[j].getCOM(COM);
267 >              break;
268 >            }
269 >          }
270 > #else
271 >          //query which processor current zconstraint molecule belongs to
272 >          int* MolToProcMap;
273 >          int whichNode;
274 >
275 >          MolToProcMap = mpiSim->getMolToProcMap();
276 >          whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
277 >
278 >          //broadcast the zpos of current z-contraint molecule
279 >          //the node which contain this
280 >
281 >          if (worldRank == whichNode){
282 >            for (int j = 0; j < nMols; j++)
283 >              if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
284 >                molecules[j].getCOM(COM);
285 >                break;
286 >              }
287 >          }
288 >
289 >          MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE, whichNode,
290 >                    MPI_COMM_WORLD);        
291 > #endif
292 >
293 >          (*parameters)[i].zPos = COM[whichDirection];
294 >
295 >          sprintf(painCave.errMsg,
296 >                  "ZConstraint warning: Does not specify zpos for z-constraint molecule "
297 >                  "initial z coornidate will be used \n");
298 >          painCave.isFatal = 0;
299 >          simError();
300 >        }
301 >      }
302 >    }//end if (!zConsParaData)
303 >
304 >  }//end  if (!data)
305 >
306 >  //  
307 > #ifdef IS_MPI
308 >  update();
309 > #else  
310 >  int searchResult;
311 >
312 >  for (int i = 0; i < nMols; i++){
313 >    searchResult = isZConstraintMol(&molecules[i]);
314 >
315 >    if (searchResult > -1){
316 >      zconsMols.push_back(&molecules[i]);      
317 >      massOfZConsMols.push_back(molecules[i].getTotalMass());  
318 >
319 >      zPos.push_back((*parameters)[searchResult].zPos);
320 >      kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
321 >      
322 >      if(usingSMD)
323 >        cantVel.push_back((*parameters)[searchResult].cantVel);
324 >
325 >    }
326 >    else{
327 >      unconsMols.push_back(&molecules[i]);
328 >      massOfUnconsMols.push_back(molecules[i].getTotalMass());
329 >    }
330 >  }
331 >
332 >  fz.resize(zconsMols.size());
333 >  curZPos.resize(zconsMols.size());
334 >  indexOfZConsMols.resize(zconsMols.size());  
335 >
336 >  //determine the states of z-constraint molecules
337 >  for (size_t i = 0; i < zconsMols.size(); i++){
338 >    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
339 >
340 >    zconsMols[i]->getCOM(COM);
341 >    
342 >    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
343 >      states.push_back(zcsFixed);
344 >
345 >      if (hasZConsGap)
346 >        endFixTime.push_back(info->getTime() + zconsFixTime);
347 >    }
348 >    else{
349 >      states.push_back(zcsMoving);
350 >
351 >      if (hasZConsGap)
352 >        endFixTime.push_back(INFINITE_TIME);
353 >    }
354 >
355 >    if(usingSMD)
356 >      cantPos.push_back(COM[whichDirection]);    
357 >  }
358 >
359 > #endif
360 >
361 >  
362 >  //get total masss of unconstraint molecules
363 >  double totalMassOfUncons_local;
364 >  totalMassOfUncons_local = 0;
365 >
366 >  for (size_t i = 0; i < unconsMols.size(); i++)
367 >    totalMassOfUncons_local += unconsMols[i]->getTotalMass();
368 >
369 > #ifndef IS_MPI
370 >  totalMassOfUncons = totalMassOfUncons_local;
371 > #else
372 >  MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,
373 >                MPI_SUM, MPI_COMM_WORLD);  
374 > #endif
375 >
376 >  //get total number of unconstrained atoms
377 >  int nUnconsAtoms_local;
378 >  nUnconsAtoms_local = 0;
379 >  for (int i = 0; i < (int) (unconsMols.size()); i++)
380 >    nUnconsAtoms_local += unconsMols[i]->getNAtoms();
381 >
382 > #ifndef IS_MPI
383 >  totNumOfUnconsAtoms = nUnconsAtoms_local;
384 > #else
385 >  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM,
386 >                MPI_COMM_WORLD);  
387 > #endif  
388 >
389 >  forcePolicy->update();
390 > }
391 >
392 > template<typename T> ZConstraint<T>::~ZConstraint(){
393 >
394 >  if (fzOut){
395 >    delete fzOut;
396 >  }
397 >
398 >  if (forcePolicy){
399 >    delete forcePolicy;
400 >  }
401 > }
402 >
403 >
404 > /**
405 > *
406 > */
407 >
408 > #ifdef IS_MPI
409 > template<typename T> void ZConstraint<T>::update(){
410 >  double COM[3];
411 >  int index;
412 >
413 >  zconsMols.clear();
414 >  massOfZConsMols.clear();
415 >  zPos.clear();
416 >  kz.clear();
417 >  cantPos.clear();
418 >  cantVel.clear();
419 >
420 >  unconsMols.clear();
421 >  massOfUnconsMols.clear();
422 >
423 >
424 >  //creat zconsMol and unconsMol lists
425 >  for (int i = 0; i < nMols; i++){
426 >    index = isZConstraintMol(&molecules[i]);
427 >
428 >    if (index > -1){
429 >      zconsMols.push_back(&molecules[i]);      
430 >      zPos.push_back((*parameters)[index].zPos);
431 >      kz.push_back((*parameters)[index].kRatio * zForceConst);
432 >      massOfZConsMols.push_back(molecules[i].getTotalMass());
433 >      
434 >      if(usingSMD)
435 >        cantVel.push_back((*parameters)[index].cantVel);
436 >
437 >    }
438 >    else{
439 >      unconsMols.push_back(&molecules[i]);
440 >      massOfUnconsMols.push_back(molecules[i].getTotalMass());
441 >    }
442 >  }
443 >
444 >  fz.resize(zconsMols.size());
445 >  curZPos.resize(zconsMols.size());
446 >  indexOfZConsMols.resize(zconsMols.size());  
447 >
448 >  for (size_t i = 0; i < zconsMols.size(); i++){
449 >    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
450 >  }
451 >    
452 >  //determine the states of z-constraint molecules
453 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
454 >
455 >    zconsMols[i]->getCOM(COM);
456 >    
457 >    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
458 >      states.push_back(zcsFixed);
459 >
460 >      if (hasZConsGap)
461 >        endFixTime.push_back(info->getTime() + zconsFixTime);
462 >    }
463 >    else{
464 >      states.push_back(zcsMoving);
465 >
466 >      if (hasZConsGap)
467 >        endFixTime.push_back(INFINITE_TIME);
468 >    }
469 >
470 >    if(usingSMD)
471 >      cantPos.push_back(COM[whichDirection]);        
472 >  }
473 >  //
474 >  forcePolicy->update();
475 > }
476 >
477 > #endif
478 >
479 > /**
480 > *  Function Name: isZConstraintMol
481 > *  Parameter
482 > *    Molecule* mol
483 > *  Return value:
484 > *    -1, if the molecule is not z-constraint molecule,
485 > *    other non-negative values, its index in indexOfAllZConsMols vector
486 > */
487 >
488 > template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){
489 >  int index;
490 >  int low;
491 >  int high;
492 >  int mid;
493 >
494 >  index = mol->getGlobalIndex();
495 >
496 >  low = 0;
497 >  high = parameters->size() - 1;
498 >
499 >  //Binary Search (we have sorted the array)  
500 >  while (low <= high){
501 >    mid = (low + high) / 2;
502 >    if ((*parameters)[mid].zconsIndex == index)
503 >      return mid;
504 >    else if ((*parameters)[mid].zconsIndex > index)
505 >      high = mid - 1;
506 >    else
507 >      low = mid + 1;
508 >  }
509 >
510 >  return -1;
511 > }
512 >
513 > template<typename T> void ZConstraint<T>::integrate(){
514 >  // creat zconsWriter  
515 >  fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);  
516 >
517 >  if (!fzOut){
518 >    sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
519 >    painCave.isFatal = 1;
520 >    simError();
521 >  }
522 >
523 >  //zero out the velocities of center of mass of unconstrained molecules
524 >  //and the velocities of center of mass of every single z-constrained molecueles
525 >  zeroOutVel();
526 >
527 >  curZconsTime = zconsTime + info->getTime();
528 >
529 >  T::integrate();
530 > }
531 >
532 >
533 > /**
534 > *
535 > *
536 > *
537 > *
538 > */
539 > template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
540 >  double zsys;
541 >  double COM[3];
542 >  double force[3];
543 >  double zSysCOMVel;
544 >
545 >  T::calcForce(calcPot, calcStress);
546 >
547 >
548 >  if (hasZConsGap){
549 >    updateZPos();
550 >  }
551 >
552 >  if (checkZConsState()){
553 >    zeroOutVel();    
554 >    forcePolicy->update();
555 >  }  
556 >
557 >  zsys = calcZSys();
558 >  zSysCOMVel = calcSysCOMVel();
559 > #ifdef IS_MPI
560 >  if (worldRank == 0){
561 > #endif
562 >    //cout << "---------------------------------------------------------------------" <<endl;
563 >    //cout << "current time: " << info->getTime() << endl;
564 >    //cout << "center of mass at z: " << zsys << endl;    
565 >    //cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
566 >
567 > #ifdef IS_MPI
568 >  }
569 > #endif
570 >
571 >  //do zconstraint force;
572 >  if (haveFixedZMols()){
573 >    this->doZconstraintForce();
574 >  }
575 >
576 >  //use external force to move the molecules to the specified positions
577 >  if (haveMovingZMols()){
578 >    if (usingSMD)
579 >      this->doHarmonic(cantPos);
580 >    else
581 >      this->doHarmonic(zPos);      
582 >  }
583 >
584 >  //write out forces and current positions of z-constraint molecules
585 >  if (info->getTime() >= curZconsTime){
586 >    for (int i = 0; i < (int) (zconsMols.size()); i++){
587 >      zconsMols[i]->getCOM(COM);
588 >      curZPos[i] = COM[whichDirection];
589 >
590 >      //if the z-constraint molecule is still moving, just record its force
591 >      if (states[i] == zcsMoving){
592 >        fz[i] = 0;
593 >        Atom** movingZAtoms;
594 >        movingZAtoms = zconsMols[i]->getMyAtoms();
595 >        for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
596 >          movingZAtoms[j]->getFrc(force);
597 >          fz[i] += force[whichDirection];
598 >        }
599 >      }
600 >    }
601 >    fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
602 >                   &curZPos[0], &zPos[0]);
603 >    curZconsTime += zconsTime;
604 >  }
605 >
606 >  zSysCOMVel = calcSysCOMVel();  
607 > #ifdef IS_MPI
608 >  if (worldRank == 0){
609 > #endif
610 >    //cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
611 > #ifdef IS_MPI
612 >  }
613 > #endif
614 > }
615 >
616 >
617 > /**
618 > *
619 > */
620 >
621 > template<typename T> double ZConstraint<T>::calcZSys(){
622 >  //calculate reference z coordinate for z-constraint molecules
623 >  double totalMass_local;
624 >  double totalMass;
625 >  double totalMZ_local;
626 >  double totalMZ;
627 >  double massOfCurMol;
628 >  double COM[3];
629 >
630 >  totalMass_local = 0;
631 >  totalMZ_local = 0;
632 >
633 >  for (int i = 0; i < nMols; i++){
634 >    massOfCurMol = molecules[i].getTotalMass();
635 >    molecules[i].getCOM(COM);
636 >
637 >    totalMass_local += massOfCurMol;
638 >    totalMZ_local += massOfCurMol * COM[whichDirection];
639 >  }
640 >
641 >
642 > #ifdef IS_MPI  
643 >  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM,
644 >                MPI_COMM_WORLD);
645 >  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
646 > #else
647 >  totalMass = totalMass_local;
648 >  totalMZ = totalMZ_local;
649 > #endif  
650 >
651 >  double zsys;
652 >  zsys = totalMZ / totalMass;
653 >
654 >  return zsys;
655 > }
656 >
657 > /**
658 > *
659 > */
660 > template<typename T> void ZConstraint<T>::thermalize(void){
661 >  T::thermalize();
662 >  zeroOutVel();
663 > }
664 >
665 > /**
666 > *
667 > */
668 >
669 > template<typename T> void ZConstraint<T>::zeroOutVel(){
670 >  Atom** fixedZAtoms;  
671 >  double COMvel[3];
672 >  double vel[3];
673 >  double zSysCOMVel;
674 >
675 >  //zero out the velocities of center of mass of fixed z-constrained molecules
676 >
677 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
678 >    if (states[i] == zcsFixed){
679 >      zconsMols[i]->getCOMvel(COMvel);      
680 >      //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
681 >
682 >      fixedZAtoms = zconsMols[i]->getMyAtoms();
683 >
684 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
685 >        fixedZAtoms[j]->getVel(vel);
686 >        vel[whichDirection] -= COMvel[whichDirection];
687 >        fixedZAtoms[j]->setVel(vel);
688 >      }
689 >
690 >      zconsMols[i]->getCOMvel(COMvel);
691 >      //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
692 >    }
693 >  }
694 >
695 >  //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
696 >
697 >  zSysCOMVel = calcSysCOMVel();
698 > #ifdef IS_MPI
699 >  if (worldRank == 0){
700 > #endif
701 >    //cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
702 > #ifdef IS_MPI
703 >  }
704 > #endif
705 >
706 >  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
707 >  double MVzOfMovingMols_local;
708 >  double MVzOfMovingMols;
709 >  double totalMassOfMovingZMols_local;
710 >  double totalMassOfMovingZMols;
711 >
712 >  MVzOfMovingMols_local = 0;
713 >  totalMassOfMovingZMols_local = 0;
714 >
715 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
716 >    unconsMols[i]->getCOMvel(COMvel);
717 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
718 >  }
719 >
720 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
721 >    if (states[i] == zcsMoving){
722 >      zconsMols[i]->getCOMvel(COMvel);
723 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
724 >      totalMassOfMovingZMols_local += massOfZConsMols[i];
725 >    }
726 >  }
727 >
728 > #ifndef IS_MPI
729 >  MVzOfMovingMols = MVzOfMovingMols_local;
730 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
731 > #else
732 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
733 >                MPI_SUM, MPI_COMM_WORLD);
734 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
735 >                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
736 > #endif
737 >
738 >  double vzOfMovingMols;
739 >  vzOfMovingMols = MVzOfMovingMols /
740 >                   (totalMassOfUncons + totalMassOfMovingZMols);
741 >
742 >  //modify the velocites of unconstrained molecules  
743 >  Atom** unconsAtoms;
744 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
745 >    unconsAtoms = unconsMols[i]->getMyAtoms();
746 >    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
747 >      unconsAtoms[j]->getVel(vel);
748 >      vel[whichDirection] -= vzOfMovingMols;
749 >      unconsAtoms[j]->setVel(vel);
750 >    }
751 >  }  
752 >
753 >  //modify the velocities of moving z-constrained molecuels
754 >  Atom** movingZAtoms;
755 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
756 >    if (states[i] == zcsMoving){
757 >      movingZAtoms = zconsMols[i]->getMyAtoms();
758 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
759 >        movingZAtoms[j]->getVel(vel);
760 >        vel[whichDirection] -= vzOfMovingMols;
761 >        movingZAtoms[j]->setVel(vel);
762 >      }
763 >    }
764 >  }
765 >
766 >
767 >  zSysCOMVel = calcSysCOMVel();
768 > #ifdef IS_MPI
769 >  if (worldRank == 0){
770 > #endif
771 >    //cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
772 > #ifdef IS_MPI
773 >  }
774 > #endif
775 > }
776 >
777 > /**
778 > *
779 > */
780 >
781 > template<typename T> void ZConstraint<T>::doZconstraintForce(){
782 >  Atom** zconsAtoms;
783 >  double totalFZ;
784 >  double totalFZ_local;
785 >  double COM[3];
786 >  double force[3];
787 >
788 >  //constrain the molecules which do not reach the specified positions  
789 >
790 >  //Zero Out the force of z-contrained molecules    
791 >  totalFZ_local = 0;
792 >
793 >  //calculate the total z-contrained force of fixed z-contrained molecules
794 >
795 >  //cout << "before zero out z-constraint force on fixed z-constraint molecuels "
796 >  //       << "total force is " << calcTotalForce() << endl;
797 >
798 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
799 >    if (states[i] == zcsFixed){
800 >      zconsMols[i]->getCOM(COM);
801 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
802 >
803 >      fz[i] = 0;      
804 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
805 >        zconsAtoms[j]->getFrc(force);
806 >        fz[i] += force[whichDirection];
807 >      }
808 >      totalFZ_local += fz[i];
809 >
810 >      //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
811 >      //      <<"\tcurrent zpos: " << COM[whichDirection]
812 >      //      << "\tcurrent fz: " <<fz[i] << endl;
813 >    }
814 >  }
815 >
816 >  //calculate total z-constraint force
817 > #ifdef IS_MPI
818 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
819 > #else
820 >  totalFZ = totalFZ_local;
821 > #endif
822 >
823 >
824 >  // apply negative to fixed z-constrained molecues;
825 >  force[0] = 0;
826 >  force[1] = 0;
827 >  force[2] = 0;
828 >
829 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
830 >    if (states[i] == zcsFixed){
831 >      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
832 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
833 >
834 >      for (int j = 0; j < nAtomOfCurZConsMol; j++){
835 >        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
836 >        force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i],
837 >                                                                zconsAtoms[j],
838 >                                                                fz[i]);
839 >        zconsAtoms[j]->addFrc(force);
840 >      }
841 >    }
842 >  }
843 >
844 >  //cout << "after zero out z-constraint force on fixed z-constraint molecuels "
845 >  //      << "total force is " << calcTotalForce() << endl;
846 >
847 >
848 >  force[0] = 0;
849 >  force[1] = 0;
850 >  force[2] = 0;
851 >
852 >  //modify the forces of unconstrained molecules
853 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
854 >    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
855 >
856 >    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
857 >      //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
858 >      force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],
859 >                                                             totalFZ);
860 >      unconsAtoms[j]->addFrc(force);
861 >    }
862 >  }      
863 >
864 >  //modify the forces of moving z-constrained molecules
865 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
866 >    if (states[i] == zcsMoving){
867 >      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
868 >
869 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
870 >        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
871 >        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],
872 >                                                               totalFZ);
873 >        movingZAtoms[j]->addFrc(force);
874 >      }
875 >    }
876 >  }
877 >  //  cout << "after substracting z-constraint force from moving molecuels "
878 >  //        << "total force is " << calcTotalForce()  << endl;
879 > }
880 >
881 > /**
882 >  *
883 >  *
884 >  */
885 >
886 > template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
887 >  double force[3];
888 >  double harmonicU;
889 >  double harmonicF;
890 >  double COM[3];
891 >  double diff;
892 >  double totalFZ_local;
893 >  double totalFZ;
894 >
895 >  force[0] = 0;
896 >  force[1] = 0;
897 >  force[2] = 0;
898 >
899 >  totalFZ_local = 0;
900 >
901 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
902 >    if (states[i] == zcsMoving){
903 >      zconsMols[i]->getCOM(COM);
904 >      //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
905 >      //     << "\tcurrent zpos: " << COM[whichDirection] << endl;
906 >
907 >      diff = COM[whichDirection] - resPos[i];
908 >
909 >      harmonicU = 0.5 * kz[i] * diff * diff;  
910 >      info->lrPot += harmonicU;
911 >
912 >      harmonicF = -kz[i] * diff;
913 >      totalFZ_local += harmonicF;
914 >
915 >      //adjust force
916 >
917 >      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
918 >
919 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
920 >        //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
921 >        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i],
922 >                                                               movingZAtoms[j],
923 >                                                               harmonicF);
924 >        movingZAtoms[j]->addFrc(force);
925 >      }
926 >    }
927 >  }
928 >
929 > #ifndef IS_MPI
930 >  totalFZ = totalFZ_local;
931 > #else
932 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
933 > #endif
934 >
935 >  //cout << "before substracting harmonic force from moving molecuels "
936 >  //      << "total force is " << calcTotalForce()  << endl;
937 >
938 >  force[0] = 0;
939 >  force[1] = 0;
940 >  force[2] = 0;
941 >
942 >  //modify the forces of unconstrained molecules
943 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
944 >    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
945 >
946 >    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
947 >      //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
948 >      force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j],
949 >                                                              totalFZ);
950 >      unconsAtoms[j]->addFrc(force);
951 >    }
952 >  }  
953 >
954 >  //cout << "after substracting harmonic force from moving molecuels "
955 >  //      << "total force is " << calcTotalForce()  << endl;
956 > }
957 >
958 > /**
959 > *
960 > */
961 >
962 > template<typename T> bool ZConstraint<T>::checkZConsState(){
963 >  double COM[3];
964 >  double diff;
965 >
966 >  int changed_local;
967 >  int changed;
968 >
969 >  changed_local = 0;
970 >
971 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
972 >    zconsMols[i]->getCOM(COM);
973 >    diff = fabs(COM[whichDirection] - zPos[i]);  
974 >    if (diff <= zconsTol && states[i] == zcsMoving){
975 >      states[i] = zcsFixed;
976 >      changed_local = 1;
977 >
978 >      if (hasZConsGap)
979 >        endFixTime[i] = info->getTime() + zconsFixTime;
980 >    }
981 >    else if (diff > zconsTol && states[i] == zcsFixed){
982 >      states[i] = zcsMoving;
983 >      changed_local = 1;  
984 >
985 >      if (hasZConsGap)
986 >        endFixTime[i] = INFINITE_TIME;
987 >    }
988 >  }
989 >
990 > #ifndef IS_MPI
991 >  changed = changed_local;
992 > #else
993 >  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
994 > #endif
995 >
996 >  return (changed > 0);
997 > }
998 >
999 > template<typename T> bool ZConstraint<T>::haveFixedZMols(){
1000 >  int havingFixed_local;
1001 >  int havingFixed;
1002 >
1003 >  havingFixed_local = 0;
1004 >
1005 >  for (int i = 0; i < (int) (zconsMols.size()); i++)
1006 >    if (states[i] == zcsFixed){
1007 >      havingFixed_local = 1;
1008 >      break;
1009 >    }
1010 >
1011 > #ifndef IS_MPI
1012 >  havingFixed = havingFixed_local;
1013 > #else
1014 >  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
1015 >                MPI_COMM_WORLD);
1016 > #endif
1017 >
1018 >  return (havingFixed > 0);
1019 > }
1020 >
1021 >
1022 > /**
1023 > *
1024 > */
1025 > template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1026 >  int havingMoving_local;
1027 >  int havingMoving;
1028 >
1029 >  havingMoving_local = 0;
1030 >
1031 >  for (int i = 0; i < (int) (zconsMols.size()); i++)
1032 >    if (states[i] == zcsMoving){
1033 >      havingMoving_local = 1;
1034 >      break;
1035 >    }
1036 >
1037 > #ifndef IS_MPI
1038 >  havingMoving = havingMoving_local;
1039 > #else
1040 >  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
1041 >                MPI_COMM_WORLD);
1042 > #endif
1043 >
1044 >  return (havingMoving > 0);
1045 > }
1046 >
1047 > /**
1048 > *
1049 > */
1050 >
1051 > template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){
1052 >  double MVzOfMovingMols_local;
1053 >  double MVzOfMovingMols;
1054 >  double totalMassOfMovingZMols_local;
1055 >  double totalMassOfMovingZMols;
1056 >  double COMvel[3];
1057 >
1058 >  MVzOfMovingMols_local = 0;
1059 >  totalMassOfMovingZMols_local = 0;
1060 >
1061 >  for (int i = 0; i < unconsMols.size(); i++){
1062 >    unconsMols[i]->getCOMvel(COMvel);
1063 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1064 >  }
1065 >
1066 >  for (int i = 0; i < zconsMols.size(); i++){
1067 >    if (states[i] == zcsMoving){
1068 >      zconsMols[i]->getCOMvel(COMvel);
1069 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
1070 >      totalMassOfMovingZMols_local += massOfZConsMols[i];
1071 >    }
1072 >  }
1073 >
1074 > #ifndef IS_MPI
1075 >  MVzOfMovingMols = MVzOfMovingMols_local;
1076 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1077 > #else
1078 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
1079 >                MPI_SUM, MPI_COMM_WORLD);
1080 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
1081 >                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
1082 > #endif
1083 >
1084 >  double vzOfMovingMols;
1085 >  vzOfMovingMols = MVzOfMovingMols /
1086 >                   (totalMassOfUncons + totalMassOfMovingZMols);
1087 >
1088 >  return vzOfMovingMols;
1089 > }
1090 >
1091 > /**
1092 > *
1093 > */
1094 >
1095 > template<typename T> double ZConstraint<T>::calcSysCOMVel(){
1096 >  double COMvel[3];
1097 >  double tempMVz_local;
1098 >  double tempMVz;
1099 >  double massOfZCons_local;
1100 >  double massOfZCons;
1101 >
1102 >
1103 >  tempMVz_local = 0;
1104 >
1105 >  for (int i = 0 ; i < nMols; i++){
1106 >    molecules[i].getCOMvel(COMvel);
1107 >    tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection];
1108 >  }
1109 >
1110 >  massOfZCons_local = 0;
1111 >
1112 >  for (int i = 0; i < (int) (massOfZConsMols.size()); i++){
1113 >    massOfZCons_local += massOfZConsMols[i];
1114 >  }
1115 > #ifndef IS_MPI
1116 >  massOfZCons = massOfZCons_local;
1117 >  tempMVz = tempMVz_local;
1118 > #else
1119 >  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM,
1120 >                MPI_COMM_WORLD);
1121 >  MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
1122 > #endif
1123 >
1124 >  return tempMVz / (totalMassOfUncons + massOfZCons);
1125 > }
1126 >
1127 > /**
1128 > *
1129 > */
1130 >
1131 > template<typename T> double ZConstraint<T>::calcTotalForce(){
1132 >  double force[3];  
1133 >  double totalForce_local;
1134 >  double totalForce;
1135 >
1136 >  totalForce_local = 0;
1137 >
1138 >  for (int i = 0; i < nAtoms; i++){
1139 >    atoms[i]->getFrc(force);
1140 >    totalForce_local += force[whichDirection];
1141 >  }
1142 >
1143 > #ifndef IS_MPI
1144 >  totalForce = totalForce_local;
1145 > #else
1146 >  MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM,
1147 >                MPI_COMM_WORLD);
1148 > #endif
1149 >
1150 >  return totalForce;
1151 > }
1152 >
1153 > /**
1154 > *
1155 > */
1156 >
1157 > template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1158 >  //calculate the number of atoms of moving z-constrained molecules
1159 >  int nMovingZAtoms_local;
1160 >  int nMovingZAtoms;
1161 >
1162 >  nMovingZAtoms_local = 0;
1163 >  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1164 >    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1165 >      nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1166 >    }
1167 >
1168 > #ifdef IS_MPI
1169 >  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM,
1170 >                MPI_COMM_WORLD);
1171 > #else
1172 >  nMovingZAtoms = nMovingZAtoms_local;
1173 > #endif
1174 >  totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1175 > }
1176 >
1177 > template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol,
1178 >                                                                              Atom* atom,
1179 >                                                                              double totalForce){
1180 >  return totalForce / mol->getNAtoms();
1181 > }
1182 >
1183 > template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom,
1184 >                                                                              double totalForce){
1185 >  return totalForce / totNumOfMovingAtoms;
1186 > }
1187 >
1188 > template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol,
1189 >                                                                              Atom* atom,
1190 >                                                                              double totalForce){
1191 >  return totalForce / mol->getNAtoms();
1192 > }
1193 >
1194 > template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom,
1195 >                                                                              double totalForce){
1196 >  return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1197 > }
1198 >
1199 > /**
1200 > *
1201 > */
1202 >
1203 > template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1204 >  //calculate the number of atoms of moving z-constrained molecules
1205 >  double massOfMovingZAtoms_local;
1206 >  double massOfMovingZAtoms;
1207 >
1208 >  massOfMovingZAtoms_local = 0;
1209 >  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1210 >    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1211 >      massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1212 >    }
1213 >
1214 > #ifdef IS_MPI
1215 >  MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,
1216 >                MPI_SUM, MPI_COMM_WORLD);
1217 > #else
1218 >  massOfMovingZAtoms = massOfMovingZAtoms_local;
1219 > #endif
1220 >  totMassOfMovingAtoms = massOfMovingZAtoms +
1221 >                         zconsIntegrator->totalMassOfUncons;
1222 > }
1223 >
1224 > template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol,
1225 >                                                                            Atom* atom,
1226 >                                                                            double totalForce){
1227 >  return totalForce * atom->getMass() / mol->getTotalMass();
1228 > }
1229 >
1230 > template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom,
1231 >                                                                            double totalForce){
1232 >  return totalForce * atom->getMass() / totMassOfMovingAtoms;
1233 > }
1234 >
1235 > template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol,
1236 >                                                                            Atom* atom,
1237 >                                                                            double totalForce){
1238 >  return totalForce * atom->getMass() / mol->getTotalMass();
1239 > }
1240 >
1241 > template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom,
1242 >                                                                            double totalForce){
1243 >  return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1244 > }
1245 >
1246 > template<typename T> void ZConstraint<T>::updateZPos(){
1247 >  double curTime;
1248 >  double COM[3];
1249 >  
1250 >  curTime = info->getTime();
1251 >
1252 >  for (size_t i = 0; i < zconsMols.size(); i++){
1253 >
1254 >    if (states[i] == zcsFixed && curTime >= endFixTime[i]){
1255 >      zPos[i] += zconsGap;
1256 >
1257 >      if (usingSMD){
1258 >        zconsMols[i]->getCOM(COM);
1259 >        cantPos[i] = COM[whichDirection];
1260 >      }
1261 >      
1262 >    }
1263 >    
1264 >  }
1265 >  
1266 > }
1267 >
1268 > template<typename T> void ZConstraint<T>::updateCantPos(){
1269 >  double curTime;
1270 >  double dt;
1271 >
1272 >  curTime = info->getTime();
1273 >  dt = info->dt;
1274 >
1275 >  for (size_t i = 0; i < zconsMols.size(); i++){
1276 >    if (states[i] == zcsMoving){
1277 >      cantPos[i] += cantVel[i] * dt;
1278 >    }
1279 >  }
1280 >
1281 > }

Diff Legend

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