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