ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/constraints/ZConstraint.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-4/src/constraints/ZConstraint.cpp (file contents):
Revision 1829 by tim, Wed Nov 3 16:08:43 2004 UTC vs.
Revision 1830 by tim, Thu Dec 2 05:17:10 2004 UTC

# Line 1 | Line 1
1 < #include "integrators/Integrator.hpp"
2 < #include "utils/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 <  DoubleGenericData* sampleTime;
17 <  DoubleGenericData* tolerance;
18 <  DoubleGenericData* gap;
19 <  DoubleGenericData* fixtime;
20 <  StringGenericData* policy;
21 <  StringGenericData* filename;
22 <  IntGenericData* 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->getPropertyByName(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<StringGenericData*>(data);
51 <
52 <    if (!policy){
53 <      sprintf(painCave.errMsg,
54 <              "ZConstraint Error: Convertion from GenericData to StringGenericData 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->getPropertyByName(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<DoubleGenericData*>(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->getPropertyByName(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<StringGenericData*>(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->getPropertyByName(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<DoubleGenericData*>(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->getPropertyByName(ZCONSGAP_ID);
149 <
150 <  if (data){
151 <    gap = dynamic_cast<DoubleGenericData*>(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->getPropertyByName(ZCONSFIXTIME_ID);
168 <
169 <  if (data){
170 <    fixtime = dynamic_cast<DoubleGenericData*>(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->getPropertyByName(ZCONSUSINGSMD_ID);
191 <
192 <  if (data){
193 <    smdFlag = dynamic_cast<IntGenericData*>(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->getPropertyByName(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->getNMolGlobal();  
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 <  if(usingSMD)
360 <    prevCantPos = cantPos;
361 < #endif
362 <
363 <  
364 <  //get total masss of unconstraint molecules
365 <  double totalMassOfUncons_local;
366 <  totalMassOfUncons_local = 0;
367 <
368 <  for (size_t i = 0; i < unconsMols.size(); i++)
369 <    totalMassOfUncons_local += unconsMols[i]->getTotalMass();
370 <
371 < #ifndef IS_MPI
372 <  totalMassOfUncons = totalMassOfUncons_local;
373 < #else
374 <  MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,
375 <                MPI_SUM, MPI_COMM_WORLD);  
376 < #endif
377 <
378 <  //get total number of unconstrained atoms
379 <  int nUnconsAtoms_local;
380 <  nUnconsAtoms_local = 0;
381 <  for (int i = 0; i < (int) (unconsMols.size()); i++)
382 <    nUnconsAtoms_local += unconsMols[i]->getNAtoms();
383 <
384 < #ifndef IS_MPI
385 <  totNumOfUnconsAtoms = nUnconsAtoms_local;
386 < #else
387 <  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM,
388 <                MPI_COMM_WORLD);  
389 < #endif  
390 <
391 <  forcePolicy->update();
392 < }
393 <
394 < template<typename T> ZConstraint<T>::~ZConstraint(){
395 <
396 <  if (fzOut){
397 <    delete fzOut;
398 <  }
399 <
400 <  if (forcePolicy){
401 <    delete forcePolicy;
402 <  }
403 < }
404 <
405 <
406 < /**
407 < *
408 < */
409 <
410 < #ifdef IS_MPI
411 < template<typename T> void ZConstraint<T>::update(){
412 <  double COM[3];
413 <  int index;
414 <
415 <  zconsMols.clear();
416 <  massOfZConsMols.clear();
417 <  zPos.clear();
418 <  kz.clear();
419 <  cantPos.clear();
420 <  cantVel.clear();
421 <
422 <  unconsMols.clear();
423 <  massOfUnconsMols.clear();
424 <
425 <
426 <  //creat zconsMol and unconsMol lists
427 <  for (int i = 0; i < nMols; i++){
428 <    index = isZConstraintMol(&molecules[i]);
429 <
430 <    if (index > -1){
431 <      zconsMols.push_back(&molecules[i]);      
432 <      zPos.push_back((*parameters)[index].zPos);
433 <      kz.push_back((*parameters)[index].kRatio * zForceConst);
434 <      massOfZConsMols.push_back(molecules[i].getTotalMass());
435 <      
436 <      if(usingSMD)
437 <        cantVel.push_back((*parameters)[index].cantVel);
438 <
439 <    }
440 <    else{
441 <      unconsMols.push_back(&molecules[i]);
442 <      massOfUnconsMols.push_back(molecules[i].getTotalMass());
443 <    }
444 <  }
445 <
446 <  fz.resize(zconsMols.size());
447 <  curZPos.resize(zconsMols.size());
448 <  indexOfZConsMols.resize(zconsMols.size());  
449 <
450 <  for (size_t i = 0; i < zconsMols.size(); i++){
451 <    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
452 <  }
453 <    
454 <  //determine the states of z-constraint molecules
455 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
456 <
457 <    zconsMols[i]->getCOM(COM);
458 <    
459 <    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
460 <      states.push_back(zcsFixed);
461 <
462 <      if (hasZConsGap)
463 <        endFixTime.push_back(info->getTime() + zconsFixTime);
464 <    }
465 <    else{
466 <      states.push_back(zcsMoving);
467 <
468 <      if (hasZConsGap)
469 <        endFixTime.push_back(INFINITE_TIME);
470 <    }
471 <
472 <    if(usingSMD)
473 <      cantPos.push_back(COM[whichDirection]);        
474 <  }
475 <
476 <  if(usingSMD)
477 <  prevCantPos = cantPos;
478 <
479 <  forcePolicy->update();
480 < }
481 <
482 < #endif
483 <
484 < /**
485 < *  Function Name: isZConstraintMol
486 < *  Parameter
487 < *    Molecule* mol
488 < *  Return value:
489 < *    -1, if the molecule is not z-constraint molecule,
490 < *    other non-negative values, its index in indexOfAllZConsMols vector
491 < */
492 <
493 < template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){
494 <  int index;
495 <  int low;
496 <  int high;
497 <  int mid;
498 <
499 <  index = mol->getGlobalIndex();
500 <
501 <  low = 0;
502 <  high = parameters->size() - 1;
503 <
504 <  //Binary Search (we have sorted the array)  
505 <  while (low <= high){
506 <    mid = (low + high) / 2;
507 <    if ((*parameters)[mid].zconsIndex == index)
508 <      return mid;
509 <    else if ((*parameters)[mid].zconsIndex > index)
510 <      high = mid - 1;
511 <    else
512 <      low = mid + 1;
513 <  }
514 <
515 <  return -1;
516 < }
517 <
518 < template<typename T> void ZConstraint<T>::integrate(){
519 <  // creat zconsWriter  
520 <  fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);  
521 <
522 <  if (!fzOut){
523 <    sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
524 <    painCave.isFatal = 1;
525 <    simError();
526 <  }
527 <
528 <  //zero out the velocities of center of mass of unconstrained molecules
529 <  //and the velocities of center of mass of every single z-constrained molecueles
530 <  zeroOutVel();
531 <
532 <  curZconsTime = zconsTime + info->getTime();
533 <
534 <  T::integrate();
535 < }
536 <
537 < template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
538 <  double zsys;
539 <  double COM[3];
540 <  double force[3];
541 <  double zSysCOMVel;
542 <
543 <  T::calcForce(calcPot, calcStress);
544 <
545 <
546 <  if (hasZConsGap){
547 <    updateZPos();
548 <  }
549 <
550 <  if (checkZConsState()){
551 <    zeroOutVel();    
552 <    forcePolicy->update();
553 <  }  
554 <
555 <  zsys = calcZSys();
556 <  zSysCOMVel = calcSysCOMVel();
557 < #ifdef IS_MPI
558 <  if (worldRank == 0){
559 < #endif
560 <
561 < #ifdef IS_MPI
562 <  }
563 < #endif
564 <
565 <  //do zconstraint force;
566 <  if (haveFixedZMols()){
567 <    this->doZconstraintForce();
568 <  }
569 <
570 <  //use external force to move the molecules to the specified positions
571 <  if (haveMovingZMols()){
572 <    if (usingSMD)
573 <      this->doHarmonic(cantPos);
574 <    else
575 <      this->doHarmonic(zPos);      
576 <  }
577 <
578 <  //write out forces and current positions of z-constraint molecules
579 <  if (info->getTime() >= curZconsTime){
580 <    for (int i = 0; i < (int) (zconsMols.size()); i++){
581 <      zconsMols[i]->getCOM(COM);
582 <      curZPos[i] = COM[whichDirection];
583 <
584 <      //if the z-constraint molecule is still moving, just record its force
585 <      if (states[i] == zcsMoving){
586 <        fz[i] = 0;
587 <        Atom** movingZAtoms;
588 <        movingZAtoms = zconsMols[i]->getMyAtoms();
589 <        for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
590 <          movingZAtoms[j]->getFrc(force);
591 <          fz[i] += force[whichDirection];
592 <        }
593 <      }
594 <    }
595 <    fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
596 <                   &curZPos[0], &zPos[0]);
597 <    curZconsTime += zconsTime;
598 <  }
599 <
600 <  zSysCOMVel = calcSysCOMVel();  
601 < #ifdef IS_MPI
602 <  if (worldRank == 0){
603 < #endif
604 < #ifdef IS_MPI
605 <  }
606 < #endif
607 < }
608 <
609 <
610 < template<typename T> double ZConstraint<T>::calcZSys(){
611 <  //calculate reference z coordinate for z-constraint molecules
612 <  double totalMass_local;
613 <  double totalMass;
614 <  double totalMZ_local;
615 <  double totalMZ;
616 <  double massOfCurMol;
617 <  double COM[3];
618 <
619 <  totalMass_local = 0;
620 <  totalMZ_local = 0;
621 <
622 <  for (int i = 0; i < nMols; i++){
623 <    massOfCurMol = molecules[i].getTotalMass();
624 <    molecules[i].getCOM(COM);
625 <
626 <    totalMass_local += massOfCurMol;
627 <    totalMZ_local += massOfCurMol * COM[whichDirection];
628 <  }
629 <
630 <
631 < #ifdef IS_MPI  
632 <  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM,
633 <                MPI_COMM_WORLD);
634 <  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
635 < #else
636 <  totalMass = totalMass_local;
637 <  totalMZ = totalMZ_local;
638 < #endif  
639 <
640 <  double zsys;
641 <  zsys = totalMZ / totalMass;
642 <
643 <  return zsys;
644 < }
645 <
646 < template<typename T> void ZConstraint<T>::thermalize(void){
647 <  T::thermalize();
648 <  zeroOutVel();
649 < }
650 <
651 < template<typename T> void ZConstraint<T>::zeroOutVel(){
652 <  Atom** fixedZAtoms;  
653 <  double COMvel[3];
654 <  double vel[3];
655 <  double zSysCOMVel;
656 <
657 <  //zero out the velocities of center of mass of fixed z-constrained molecules
658 <
659 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
660 <    if (states[i] == zcsFixed){
661 <      zconsMols[i]->getCOMvel(COMvel);      
662 <      //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
663 <
664 <      fixedZAtoms = zconsMols[i]->getMyAtoms();
665 <
666 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
667 <        vel = fixedZAtoms[j]->getVel();
668 <        vel[whichDirection] -= COMvel[whichDirection];
669 <        fixedZAtoms[j]->setVel(vel);
670 <      }
671 <
672 <      zconsMols[i]->getCOMvel(COMvel);
673 <    }
674 <  }
675 <
676 <  zSysCOMVel = calcSysCOMVel();
677 < #ifdef IS_MPI
678 <  if (worldRank == 0){
679 < #endif
680 < #ifdef IS_MPI
681 <  }
682 < #endif
683 <
684 <  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
685 <  double MVzOfMovingMols_local;
686 <  double MVzOfMovingMols;
687 <  double totalMassOfMovingZMols_local;
688 <  double totalMassOfMovingZMols;
689 <
690 <  MVzOfMovingMols_local = 0;
691 <  totalMassOfMovingZMols_local = 0;
692 <
693 <  for (int i = 0; i < (int) (unconsMols.size()); i++){
694 <    unconsMols[i]->getCOMvel(COMvel);
695 <    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
696 <  }
697 <
698 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
699 <    if (states[i] == zcsMoving){
700 <      zconsMols[i]->getCOMvel(COMvel);
701 <      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
702 <      totalMassOfMovingZMols_local += massOfZConsMols[i];
703 <    }
704 <  }
705 <
706 < #ifndef IS_MPI
707 <  MVzOfMovingMols = MVzOfMovingMols_local;
708 <  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
709 < #else
710 <  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
711 <                MPI_SUM, MPI_COMM_WORLD);
712 <  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
713 <                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
714 < #endif
715 <
716 <  double vzOfMovingMols;
717 <  vzOfMovingMols = MVzOfMovingMols /
718 <                   (totalMassOfUncons + totalMassOfMovingZMols);
719 <
720 <  //modify the velocites of unconstrained molecules  
721 <  Atom** unconsAtoms;
722 <  for (int i = 0; i < (int) (unconsMols.size()); i++){
723 <    unconsAtoms = unconsMols[i]->getMyAtoms();
724 <    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
725 <      vel = unconsAtoms[j]->getVel();
726 <      vel[whichDirection] -= vzOfMovingMols;
727 <      unconsAtoms[j]->setVel(vel);
728 <    }
729 <  }  
730 <
731 <  //modify the velocities of moving z-constrained molecuels
732 <  Atom** movingZAtoms;
733 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
734 <    if (states[i] == zcsMoving){
735 <      movingZAtoms = zconsMols[i]->getMyAtoms();
736 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
737 <        vel = movingZAtoms[j]->getVel();
738 <        vel[whichDirection] -= vzOfMovingMols;
739 <        movingZAtoms[j]->setVel(vel);
740 <      }
741 <    }
742 <  }
743 <
744 <
745 <  zSysCOMVel = calcSysCOMVel();
746 < #ifdef IS_MPI
747 <  if (worldRank == 0){
748 < #endif
749 < #ifdef IS_MPI
750 <  }
751 < #endif
752 < }
753 <
754 <
755 < template<typename T> void ZConstraint<T>::doZconstraintForce(){
756 <  Atom** zconsAtoms;
757 <  double totalFZ;
758 <  double totalFZ_local;
759 <  double COM[3];
760 <  double force[3];
761 <
762 <  //constrain the molecules which do not reach the specified positions  
763 <
764 <  //Zero Out the force of z-contrained molecules    
765 <  totalFZ_local = 0;
766 <
767 <  //calculate the total z-contrained force of fixed z-contrained molecules
768 <
769 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
770 <    if (states[i] == zcsFixed){
771 <      zconsMols[i]->getCOM(COM);
772 <      zconsAtoms = zconsMols[i]->getMyAtoms();  
773 <
774 <      fz[i] = 0;      
775 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
776 <        zconsAtoms[j]->getFrc(force);
777 <        fz[i] += force[whichDirection];
778 <      }
779 <      totalFZ_local += fz[i];
780 <
781 <    }
782 <  }
783 <
784 <  //calculate total z-constraint force
785 < #ifdef IS_MPI
786 <  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
787 < #else
788 <  totalFZ = totalFZ_local;
789 < #endif
790 <
791 <
792 <  // apply negative to fixed z-constrained molecues;
793 <  force[0] = 0;
794 <  force[1] = 0;
795 <  force[2] = 0;
796 <
797 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
798 <    if (states[i] == zcsFixed){
799 <      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
800 <      zconsAtoms = zconsMols[i]->getMyAtoms();  
801 <
802 <      for (int j = 0; j < nAtomOfCurZConsMol; j++){
803 <        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
804 <        force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i],
805 <                                                                zconsAtoms[j],
806 <                                                                fz[i]);
807 <        zconsAtoms[j]->addFrc(force);
808 <      }
809 <    }
810 <  }
811 <
812 <  force[0] = 0;
813 <  force[1] = 0;
814 <  force[2] = 0;
815 <
816 <  //modify the forces of unconstrained molecules
817 <  for (int i = 0; i < (int) (unconsMols.size()); i++){
818 <    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
819 <
820 <    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
821 <      //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
822 <      force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],
823 <                                                             totalFZ);
824 <      unconsAtoms[j]->addFrc(force);
825 <    }
826 <  }      
827 <
828 <  //modify the forces of moving z-constrained molecules
829 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
830 <    if (states[i] == zcsMoving){
831 <      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
832 <
833 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
834 <        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
835 <        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],
836 <                                                               totalFZ);
837 <        movingZAtoms[j]->addFrc(force);
838 <      }
839 <    }
840 <  }
841 < }
842 <
843 <
844 < template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
845 <  double force[3];
846 <  double harmonicU;
847 <  double harmonicF;
848 <  double COM[3];
849 <  double diff;
850 <  double totalFZ_local;
851 <  double totalFZ;
852 <
853 <  force[0] = 0;
854 <  force[1] = 0;
855 <  force[2] = 0;
856 <
857 <  totalFZ_local = 0;
858 <
859 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
860 <    if (states[i] == zcsMoving){
861 <      zconsMols[i]->getCOM(COM);
862 <
863 <      diff = COM[whichDirection] - resPos[i];
864 <
865 <      harmonicU = 0.5 * kz[i] * diff * diff;  
866 <      info->lrPot += harmonicU;
867 <
868 <      harmonicF = -kz[i] * diff;
869 <      totalFZ_local += harmonicF;
870 <
871 <      //adjust force
872 <
873 <      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
874 <
875 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
876 <        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i],
877 <                                                               movingZAtoms[j],
878 <                                                               harmonicF);
879 <        movingZAtoms[j]->addFrc(force);
880 <      }
881 <    }
882 <  }
883 <
884 < #ifndef IS_MPI
885 <  totalFZ = totalFZ_local;
886 < #else
887 <  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
888 < #endif
889 <
890 <  force[0] = 0;
891 <  force[1] = 0;
892 <  force[2] = 0;
893 <
894 <  //modify the forces of unconstrained molecules
895 <  for (int i = 0; i < (int) (unconsMols.size()); i++){
896 <    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
897 <
898 <    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
899 <      //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
900 <      force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j],
901 <                                                              totalFZ);
902 <      unconsAtoms[j]->addFrc(force);
903 <    }
904 <  }  
905 <
906 < }
907 <
908 < template<typename T> bool ZConstraint<T>::checkZConsState(){
909 <  double COM[3];
910 <  double diff;
911 <
912 <  int changed_local;
913 <  int changed;
914 <
915 <  changed_local = 0;
916 <
917 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
918 <    zconsMols[i]->getCOM(COM);
919 <    diff = fabs(COM[whichDirection] - zPos[i]);  
920 <    if (diff <= zconsTol && states[i] == zcsMoving){
921 <      states[i] = zcsFixed;
922 <      changed_local = 1;
923 <
924 <      if(usingSMD)
925 <        prevCantPos = cantPos;
926 <
927 <      if (hasZConsGap)
928 <        endFixTime[i] = info->getTime() + zconsFixTime;
929 <    }
930 <    else if (diff > zconsTol && states[i] == zcsFixed){
931 <      states[i] = zcsMoving;
932 <      changed_local = 1;  
933 <
934 <      if(usingSMD)
935 <         cantPos = prevCantPos;
936 <      
937 <      if (hasZConsGap)
938 <        endFixTime[i] = INFINITE_TIME;
939 <    }
940 <  }
941 <
942 < #ifndef IS_MPI
943 <  changed = changed_local;
944 < #else
945 <  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
946 < #endif
947 <
948 <  return (changed > 0);
949 < }
950 <
951 < template<typename T> bool ZConstraint<T>::haveFixedZMols(){
952 <  int havingFixed_local;
953 <  int havingFixed;
954 <
955 <  havingFixed_local = 0;
956 <
957 <  for (int i = 0; i < (int) (zconsMols.size()); i++)
958 <    if (states[i] == zcsFixed){
959 <      havingFixed_local = 1;
960 <      break;
961 <    }
962 <
963 < #ifndef IS_MPI
964 <  havingFixed = havingFixed_local;
965 < #else
966 <  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
967 <                MPI_COMM_WORLD);
968 < #endif
969 <
970 <  return (havingFixed > 0);
971 < }
972 <
973 <
974 < template<typename T> bool ZConstraint<T>::haveMovingZMols(){
975 <  int havingMoving_local;
976 <  int havingMoving;
977 <
978 <  havingMoving_local = 0;
979 <
980 <  for (int i = 0; i < (int) (zconsMols.size()); i++)
981 <    if (states[i] == zcsMoving){
982 <      havingMoving_local = 1;
983 <      break;
984 <    }
985 <
986 < #ifndef IS_MPI
987 <  havingMoving = havingMoving_local;
988 < #else
989 <  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
990 <                MPI_COMM_WORLD);
991 < #endif
992 <
993 <  return (havingMoving > 0);
994 < }
995 <
996 <
997 < template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){
998 <  double MVzOfMovingMols_local;
999 <  double MVzOfMovingMols;
1000 <  double totalMassOfMovingZMols_local;
1001 <  double totalMassOfMovingZMols;
1002 <  double COMvel[3];
1003 <
1004 <  MVzOfMovingMols_local = 0;
1005 <  totalMassOfMovingZMols_local = 0;
1006 <
1007 <  for (int i = 0; i < unconsMols.size(); i++){
1008 <    unconsMols[i]->getCOMvel(COMvel);
1009 <    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1010 <  }
1011 <
1012 <  for (int i = 0; i < zconsMols.size(); i++){
1013 <    if (states[i] == zcsMoving){
1014 <      zconsMols[i]->getCOMvel(COMvel);
1015 <      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
1016 <      totalMassOfMovingZMols_local += massOfZConsMols[i];
1017 <    }
1018 <  }
1019 <
1020 < #ifndef IS_MPI
1021 <  MVzOfMovingMols = MVzOfMovingMols_local;
1022 <  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1023 < #else
1024 <  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
1025 <                MPI_SUM, MPI_COMM_WORLD);
1026 <  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
1027 <                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
1028 < #endif
1029 <
1030 <  double vzOfMovingMols;
1031 <  vzOfMovingMols = MVzOfMovingMols /
1032 <                   (totalMassOfUncons + totalMassOfMovingZMols);
1033 <
1034 <  return vzOfMovingMols;
1035 < }
1036 <
1037 < template<typename T> double ZConstraint<T>::calcSysCOMVel(){
1038 <  double COMvel[3];
1039 <  double tempMVz_local;
1040 <  double tempMVz;
1041 <  double massOfZCons_local;
1042 <  double massOfZCons;
1043 <
1044 <
1045 <  tempMVz_local = 0;
1046 <
1047 <  for (int i = 0 ; i < nMols; i++){
1048 <    molecules[i].getCOMvel(COMvel);
1049 <    tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection];
1050 <  }
1051 <
1052 <  massOfZCons_local = 0;
1053 <
1054 <  for (int i = 0; i < (int) (massOfZConsMols.size()); i++){
1055 <    massOfZCons_local += massOfZConsMols[i];
1056 <  }
1057 < #ifndef IS_MPI
1058 <  massOfZCons = massOfZCons_local;
1059 <  tempMVz = tempMVz_local;
1060 < #else
1061 <  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM,
1062 <                MPI_COMM_WORLD);
1063 <  MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
1064 < #endif
1065 <
1066 <  return tempMVz / (totalMassOfUncons + massOfZCons);
1067 < }
1068 <
1069 < template<typename T> double ZConstraint<T>::calcTotalForce(){
1070 <  double force[3];  
1071 <  double totalForce_local;
1072 <  double totalForce;
1073 <
1074 <  totalForce_local = 0;
1075 <
1076 <  for (int i = 0; i < nAtoms; i++){
1077 <    atoms[i]->getFrc(force);
1078 <    totalForce_local += force[whichDirection];
1079 <  }
1080 <
1081 < #ifndef IS_MPI
1082 <  totalForce = totalForce_local;
1083 < #else
1084 <  MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM,
1085 <                MPI_COMM_WORLD);
1086 < #endif
1087 <
1088 <  return totalForce;
1089 < }
1090 <
1091 < template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1092 <  //calculate the number of atoms of moving z-constrained molecules
1093 <  int nMovingZAtoms_local;
1094 <  int nMovingZAtoms;
1095 <
1096 <  nMovingZAtoms_local = 0;
1097 <  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1098 <    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1099 <      nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1100 <    }
1101 <
1102 < #ifdef IS_MPI
1103 <  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM,
1104 <                MPI_COMM_WORLD);
1105 < #else
1106 <  nMovingZAtoms = nMovingZAtoms_local;
1107 < #endif
1108 <  totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1109 < }
1110 <
1111 < template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol,
1112 <                                                                              Atom* atom,
1113 <                                                                              double totalForce){
1114 <  return totalForce / mol->getNAtoms();
1115 < }
1116 <
1117 < template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom,
1118 <                                                                              double totalForce){
1119 <  return totalForce / totNumOfMovingAtoms;
1120 < }
1121 <
1122 < template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol,
1123 <                                                                              Atom* atom,
1124 <                                                                              double totalForce){
1125 <  return totalForce / mol->getNAtoms();
1126 < }
1127 <
1128 < template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom,
1129 <                                                                              double totalForce){
1130 <  return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1131 < }
1132 <
1133 <
1134 < template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1135 <  //calculate the number of atoms of moving z-constrained molecules
1136 <  double massOfMovingZAtoms_local;
1137 <  double massOfMovingZAtoms;
1138 <
1139 <  massOfMovingZAtoms_local = 0;
1140 <  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1141 <    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1142 <      massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1143 <    }
1144 <
1145 < #ifdef IS_MPI
1146 <  MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,
1147 <                MPI_SUM, MPI_COMM_WORLD);
1148 < #else
1149 <  massOfMovingZAtoms = massOfMovingZAtoms_local;
1150 < #endif
1151 <  totMassOfMovingAtoms = massOfMovingZAtoms +
1152 <                         zconsIntegrator->totalMassOfUncons;
1153 < }
1154 <
1155 < template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol,
1156 <                                                                            Atom* atom,
1157 <                                                                            double totalForce){
1158 <  return totalForce * atom->getMass() / mol->getTotalMass();
1159 < }
1160 <
1161 < template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom,
1162 <                                                                            double totalForce){
1163 <  return totalForce * atom->getMass() / totMassOfMovingAtoms;
1164 < }
1165 <
1166 < template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol,
1167 <                                                                            Atom* atom,
1168 <                                                                            double totalForce){
1169 <  return totalForce * atom->getMass() / mol->getTotalMass();
1170 < }
1171 <
1172 < template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom,
1173 <                                                                            double totalForce){
1174 <  return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1175 < }
1176 <
1177 < template<typename T> void ZConstraint<T>::updateZPos(){
1178 <  double curTime;
1179 <  double COM[3];
1180 <  
1181 <  curTime = info->getTime();
1182 <
1183 <  for (size_t i = 0; i < zconsMols.size(); i++){
1184 <
1185 <    if (states[i] == zcsFixed && curTime >= endFixTime[i]){
1186 <      zPos[i] += zconsGap;
1187 <
1188 <      if (usingSMD){
1189 <        zconsMols[i]->getCOM(COM);
1190 <        cantPos[i] = COM[whichDirection];
1191 <      }
1192 <      
1193 <    }
1194 <    
1195 <  }
1196 <  
1197 < }
1198 <
1199 < template<typename T> void ZConstraint<T>::updateCantPos(){
1200 <  double curTime;
1201 <  double dt;
1202 <
1203 <  curTime = info->getTime();
1204 <  dt = info->dt;
1205 <
1206 <  for (size_t i = 0; i < zconsMols.size(); i++){
1207 <    if (states[i] == zcsMoving){
1208 <      cantPos[i] += cantVel[i] * dt;
1209 <    }
1210 <  }
1211 <
1212 < }
1 > #include "integrators/Integrator.hpp"
2 >
3 > #include "utils/simError.h"
4 >
5 > #include <math.h>
6 >
7 >
8 >
9 > const double INFINITE_TIME = 10e30;
10 >
11 > template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo,
12 >
13 >                                                 ForceFields* the_ff): T(theInfo, the_ff),
14 >
15 >                                                                       fzOut(NULL),
16 >
17 >                                                                       curZconsTime(0),
18 >
19 >                                                                       forcePolicy(NULL),
20 >
21 >                                                                       usingSMD(false),
22 >
23 >                                                                       hasZConsGap(false){
24 >
25 >  //get properties from SimInfo
26 >
27 >  GenericData* data;
28 >
29 >  ZConsParaData* zConsParaData;
30 >
31 >  DoubleGenericData* sampleTime;
32 >
33 >  DoubleGenericData* tolerance;
34 >
35 >  DoubleGenericData* gap;
36 >
37 >  DoubleGenericData* fixtime;
38 >
39 >  StringGenericData* policy;
40 >
41 >  StringGenericData* filename;
42 >
43 >  IntGenericData* smdFlag;
44 >
45 >  double COM[3];
46 >
47 >
48 >
49 >  //by default, the direction of constraint is z
50 >
51 >  // 0 --> x
52 >
53 >  // 1 --> y
54 >
55 >  // 2 --> z
56 >
57 >  whichDirection = 2;
58 >
59 >
60 >
61 >  //estimate the force constant of harmonical potential
62 >
63 >  double Kb = 1.986E-3 ; //in kcal/K
64 >
65 >
66 >
67 >  double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /
68 >
69 >                            2;
70 >
71 >  zForceConst = Kb * info->target_temp / (halfOfLargestBox * halfOfLargestBox);
72 >
73 >
74 >
75 >  //creat force Subtraction policy
76 >
77 >  data = info->getPropertyByName(ZCONSFORCEPOLICY_ID);
78 >
79 >  if (!data){
80 >
81 >    sprintf(painCave.errMsg,
82 >
83 >            "ZConstraint Warning: User does not set force Subtraction policy, "
84 >
85 >            "PolicyByMass is used\n");
86 >
87 >    painCave.isFatal = 0;
88 >
89 >    simError();      
90 >
91 >
92 >
93 >    forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
94 >
95 >  }
96 >
97 >  else{
98 >
99 >    policy = dynamic_cast<StringGenericData*>(data);
100 >
101 >
102 >
103 >    if (!policy){
104 >
105 >      sprintf(painCave.errMsg,
106 >
107 >              "ZConstraint Error: Convertion from GenericData to StringGenericData failure, "
108 >
109 >              "PolicyByMass is used\n");
110 >
111 >      painCave.isFatal = 0;
112 >
113 >      simError();      
114 >
115 >
116 >
117 >      forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
118 >
119 >    }
120 >
121 >    else{
122 >
123 >      if (policy->getData() == "BYNUMBER")
124 >
125 >        forcePolicy = (ForceSubtractionPolicy *) new PolicyByNumber(this);
126 >
127 >      else if (policy->getData() == "BYMASS")
128 >
129 >        forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
130 >
131 >      else{
132 >
133 >        sprintf(painCave.errMsg,
134 >
135 >                "ZConstraint Warning: unknown force Subtraction policy, "
136 >
137 >                "PolicyByMass is used\n");
138 >
139 >        painCave.isFatal = 0;
140 >
141 >        simError();      
142 >
143 >        forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
144 >
145 >      }
146 >
147 >    }
148 >
149 >  }
150 >
151 >
152 >
153 >
154 >
155 >  //retrieve sample time of z-contraint
156 >
157 >  data = info->getPropertyByName(ZCONSTIME_ID);
158 >
159 >
160 >
161 >  if (!data){
162 >
163 >    sprintf(painCave.errMsg,
164 >
165 >            "ZConstraint error: If you use an ZConstraint\n"
166 >
167 >            " , you must set sample time.\n");
168 >
169 >    painCave.isFatal = 1;
170 >
171 >    simError();
172 >
173 >  }
174 >
175 >  else{
176 >
177 >    sampleTime = dynamic_cast<DoubleGenericData*>(data);
178 >
179 >
180 >
181 >    if (!sampleTime){
182 >
183 >      sprintf(painCave.errMsg,
184 >
185 >              "ZConstraint error: Can not get property from SimInfo\n");
186 >
187 >      painCave.isFatal = 1;
188 >
189 >      simError();
190 >
191 >    }
192 >
193 >    else{
194 >
195 >      this->zconsTime = sampleTime->getData();
196 >
197 >    }
198 >
199 >  }
200 >
201 >
202 >
203 >  //retrieve output filename of z force
204 >
205 >  data = info->getPropertyByName(ZCONSFILENAME_ID);
206 >
207 >  if (!data){
208 >
209 >    sprintf(painCave.errMsg,
210 >
211 >            "ZConstraint error: If you use an ZConstraint\n"
212 >
213 >            " , you must set output filename of z-force.\n");
214 >
215 >    painCave.isFatal = 1;
216 >
217 >    simError();
218 >
219 >  }
220 >
221 >  else{
222 >
223 >    filename = dynamic_cast<StringGenericData*>(data);
224 >
225 >
226 >
227 >    if (!filename){
228 >
229 >      sprintf(painCave.errMsg,
230 >
231 >              "ZConstraint error: Can not get property from SimInfo\n");
232 >
233 >      painCave.isFatal = 1;
234 >
235 >      simError();
236 >
237 >    }
238 >
239 >    else{
240 >
241 >      this->zconsOutput = filename->getData();
242 >
243 >    }
244 >
245 >  }
246 >
247 >
248 >
249 >  //retrieve tolerance for z-constraint molecuels
250 >
251 >  data = info->getPropertyByName(ZCONSTOL_ID);
252 >
253 >
254 >
255 >  if (!data){
256 >
257 >    sprintf(painCave.errMsg, "ZConstraint error: can not get tolerance \n");
258 >
259 >    painCave.isFatal = 1;
260 >
261 >    simError();
262 >
263 >  }
264 >
265 >  else{
266 >
267 >    tolerance = dynamic_cast<DoubleGenericData*>(data);
268 >
269 >
270 >
271 >    if (!tolerance){
272 >
273 >      sprintf(painCave.errMsg,
274 >
275 >              "ZConstraint error: Can not get property from SimInfo\n");
276 >
277 >      painCave.isFatal = 1;
278 >
279 >      simError();
280 >
281 >    }
282 >
283 >    else{
284 >
285 >      this->zconsTol = tolerance->getData();
286 >
287 >    }
288 >
289 >  }
290 >
291 >
292 >
293 >  //quick hack here
294 >
295 >  data = info->getPropertyByName(ZCONSGAP_ID);
296 >
297 >
298 >
299 >  if (data){
300 >
301 >    gap = dynamic_cast<DoubleGenericData*>(data);
302 >
303 >
304 >
305 >    if (!gap){
306 >
307 >      sprintf(painCave.errMsg,
308 >
309 >              "ZConstraint error: Can not get property from SimInfo\n");
310 >
311 >      painCave.isFatal = 1;
312 >
313 >      simError();
314 >
315 >    }
316 >
317 >    else{
318 >
319 >      this->hasZConsGap = true;
320 >
321 >      this->zconsGap = gap->getData();
322 >
323 >    }
324 >
325 >  }
326 >
327 >
328 >
329 >
330 >
331 >
332 >
333 >  data = info->getPropertyByName(ZCONSFIXTIME_ID);
334 >
335 >
336 >
337 >  if (data){
338 >
339 >    fixtime = dynamic_cast<DoubleGenericData*>(data);
340 >
341 >    if (!fixtime){
342 >
343 >      sprintf(painCave.errMsg,
344 >
345 >              "ZConstraint error: Can not get zconsFixTime from SimInfo\n");
346 >
347 >      painCave.isFatal = 1;
348 >
349 >      simError();
350 >
351 >    }
352 >
353 >    else{
354 >
355 >      this->zconsFixTime = fixtime->getData();
356 >
357 >    }
358 >
359 >  }
360 >
361 >  else if(hasZConsGap){
362 >
363 >      sprintf(painCave.errMsg,
364 >
365 >              "ZConstraint error: must set fixtime if already set zconsGap\n");
366 >
367 >      painCave.isFatal = 1;
368 >
369 >      simError();
370 >
371 >  }
372 >
373 >
374 >
375 >
376 >
377 >
378 >
379 >  data = info->getPropertyByName(ZCONSUSINGSMD_ID);
380 >
381 >
382 >
383 >  if (data){
384 >
385 >    smdFlag = dynamic_cast<IntGenericData*>(data);
386 >
387 >
388 >
389 >    if (!smdFlag){
390 >
391 >      sprintf(painCave.errMsg,
392 >
393 >              "ZConstraint error: Can not get property from SimInfo\n");
394 >
395 >      painCave.isFatal = 1;
396 >
397 >      simError();
398 >
399 >    }
400 >
401 >    else{
402 >
403 >      this->usingSMD= smdFlag->getData() ? true : false;
404 >
405 >    }
406 >
407 >
408 >
409 >  }
410 >
411 >
412 >
413 >
414 >
415 >
416 >
417 >  //retrieve index of z-constraint molecules
418 >
419 >  data = info->getPropertyByName(ZCONSPARADATA_ID);
420 >
421 >  if (!data){
422 >
423 >    sprintf(painCave.errMsg,
424 >
425 >            "ZConstraint error: If you use an ZConstraint\n"
426 >
427 >            " , you must set index of z-constraint molecules.\n");
428 >
429 >    painCave.isFatal = 1;
430 >
431 >    simError();
432 >
433 >  }
434 >
435 >  else{
436 >
437 >    zConsParaData = dynamic_cast<ZConsParaData*>(data);
438 >
439 >
440 >
441 >    if (!zConsParaData){
442 >
443 >      sprintf(painCave.errMsg,
444 >
445 >              "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
446 >
447 >      painCave.isFatal = 1;
448 >
449 >      simError();
450 >
451 >    }
452 >
453 >    else{
454 >
455 >      parameters = zConsParaData->getData();
456 >
457 >
458 >
459 >      //check the range of zconsIndex
460 >
461 >      //and the minimum value of index is the first one (we already sorted the data)
462 >
463 >      //the maximum value of index is the last one
464 >
465 >
466 >
467 >      int maxIndex;
468 >
469 >      int minIndex;
470 >
471 >      int totalNumMol;
472 >
473 >
474 >
475 >      minIndex = (*parameters)[0].zconsIndex;
476 >
477 >      if (minIndex < 0){
478 >
479 >        sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n");
480 >
481 >        painCave.isFatal = 1;
482 >
483 >        simError();
484 >
485 >      }
486 >
487 >
488 >
489 >      maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
490 >
491 >
492 >
493 > #ifndef IS_MPI
494 >
495 >      totalNumMol = nMols;
496 >
497 > #else
498 >
499 >      totalNumMol = mpiSim->getNMolGlobal();  
500 >
501 > #endif      
502 >
503 >
504 >
505 >      if (maxIndex > totalNumMol - 1){
506 >
507 >        sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n");
508 >
509 >        painCave.isFatal = 1;
510 >
511 >        simError();
512 >
513 >      }
514 >
515 >
516 >
517 >      //if user does not specify the zpos for the zconstraint molecule
518 >
519 >      //its initial z coordinate  will be used as default
520 >
521 >      for (int i = 0; i < (int) (parameters->size()); i++){
522 >
523 >        if (!(*parameters)[i].havingZPos){
524 >
525 > #ifndef IS_MPI
526 >
527 >          for (int j = 0; j < nMols; j++){
528 >
529 >            if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
530 >
531 >              molecules[j].getCOM(COM);
532 >
533 >              break;
534 >
535 >            }
536 >
537 >          }
538 >
539 > #else
540 >
541 >          //query which processor current zconstraint molecule belongs to
542 >
543 >          int* MolToProcMap;
544 >
545 >          int whichNode;
546 >
547 >
548 >
549 >          MolToProcMap = mpiSim->getMolToProcMap();
550 >
551 >          whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
552 >
553 >
554 >
555 >          //broadcast the zpos of current z-contraint molecule
556 >
557 >          //the node which contain this
558 >
559 >
560 >
561 >          if (worldRank == whichNode){
562 >
563 >            for (int j = 0; j < nMols; j++)
564 >
565 >              if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
566 >
567 >                molecules[j].getCOM(COM);
568 >
569 >                break;
570 >
571 >              }
572 >
573 >          }
574 >
575 >
576 >
577 >          MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE, whichNode,
578 >
579 >                    MPI_COMM_WORLD);        
580 >
581 > #endif
582 >
583 >
584 >
585 >          (*parameters)[i].zPos = COM[whichDirection];
586 >
587 >
588 >
589 >          sprintf(painCave.errMsg,
590 >
591 >                  "ZConstraint warning: Does not specify zpos for z-constraint molecule "
592 >
593 >                  "initial z coornidate will be used \n");
594 >
595 >          painCave.isFatal = 0;
596 >
597 >          simError();
598 >
599 >        }
600 >
601 >      }
602 >
603 >    }//end if (!zConsParaData)
604 >
605 >
606 >
607 >  }//end  if (!data)
608 >
609 >
610 >
611 >  //  
612 >
613 > #ifdef IS_MPI
614 >
615 >  update();
616 >
617 > #else  
618 >
619 >  int searchResult;
620 >
621 >
622 >
623 >  for (int i = 0; i < nMols; i++){
624 >
625 >    searchResult = isZConstraintMol(&molecules[i]);
626 >
627 >
628 >
629 >    if (searchResult > -1){
630 >
631 >      zconsMols.push_back(&molecules[i]);      
632 >
633 >      massOfZConsMols.push_back(molecules[i].getTotalMass());  
634 >
635 >
636 >
637 >      zPos.push_back((*parameters)[searchResult].zPos);
638 >
639 >      kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
640 >
641 >      
642 >
643 >      if(usingSMD)
644 >
645 >        cantVel.push_back((*parameters)[searchResult].cantVel);
646 >
647 >
648 >
649 >    }
650 >
651 >    else{
652 >
653 >      unconsMols.push_back(&molecules[i]);
654 >
655 >      massOfUnconsMols.push_back(molecules[i].getTotalMass());
656 >
657 >    }
658 >
659 >  }
660 >
661 >
662 >
663 >  fz.resize(zconsMols.size());
664 >
665 >  curZPos.resize(zconsMols.size());
666 >
667 >  indexOfZConsMols.resize(zconsMols.size());  
668 >
669 >
670 >
671 >  //determine the states of z-constraint molecules
672 >
673 >  for (size_t i = 0; i < zconsMols.size(); i++){
674 >
675 >    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
676 >
677 >
678 >
679 >    zconsMols[i]->getCOM(COM);
680 >
681 >    
682 >
683 >    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
684 >
685 >      states.push_back(zcsFixed);
686 >
687 >
688 >
689 >      if (hasZConsGap)
690 >
691 >        endFixTime.push_back(info->getTime() + zconsFixTime);
692 >
693 >    }
694 >
695 >    else{
696 >
697 >      states.push_back(zcsMoving);
698 >
699 >
700 >
701 >      if (hasZConsGap)
702 >
703 >        endFixTime.push_back(INFINITE_TIME);
704 >
705 >    }
706 >
707 >
708 >
709 >    if(usingSMD)
710 >
711 >      cantPos.push_back(COM[whichDirection]);    
712 >
713 >  }
714 >
715 >
716 >
717 >  if(usingSMD)
718 >
719 >    prevCantPos = cantPos;
720 >
721 > #endif
722 >
723 >
724 >
725 >  
726 >
727 >  //get total masss of unconstraint molecules
728 >
729 >  double totalMassOfUncons_local;
730 >
731 >  totalMassOfUncons_local = 0;
732 >
733 >
734 >
735 >  for (size_t i = 0; i < unconsMols.size(); i++)
736 >
737 >    totalMassOfUncons_local += unconsMols[i]->getTotalMass();
738 >
739 >
740 >
741 > #ifndef IS_MPI
742 >
743 >  totalMassOfUncons = totalMassOfUncons_local;
744 >
745 > #else
746 >
747 >  MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,
748 >
749 >                MPI_SUM, MPI_COMM_WORLD);  
750 >
751 > #endif
752 >
753 >
754 >
755 >  //get total number of unconstrained atoms
756 >
757 >  int nUnconsAtoms_local;
758 >
759 >  nUnconsAtoms_local = 0;
760 >
761 >  for (int i = 0; i < (int) (unconsMols.size()); i++)
762 >
763 >    nUnconsAtoms_local += unconsMols[i]->getNAtoms();
764 >
765 >
766 >
767 > #ifndef IS_MPI
768 >
769 >  totNumOfUnconsAtoms = nUnconsAtoms_local;
770 >
771 > #else
772 >
773 >  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM,
774 >
775 >                MPI_COMM_WORLD);  
776 >
777 > #endif  
778 >
779 >
780 >
781 >  forcePolicy->update();
782 >
783 > }
784 >
785 >
786 >
787 > template<typename T> ZConstraint<T>::~ZConstraint(){
788 >
789 >
790 >
791 >  if (fzOut){
792 >
793 >    delete fzOut;
794 >
795 >  }
796 >
797 >
798 >
799 >  if (forcePolicy){
800 >
801 >    delete forcePolicy;
802 >
803 >  }
804 >
805 > }
806 >
807 >
808 >
809 >
810 >
811 > /**
812 >
813 > *
814 >
815 > */
816 >
817 >
818 >
819 > #ifdef IS_MPI
820 >
821 > template<typename T> void ZConstraint<T>::update(){
822 >
823 >  double COM[3];
824 >
825 >  int index;
826 >
827 >
828 >
829 >  zconsMols.clear();
830 >
831 >  massOfZConsMols.clear();
832 >
833 >  zPos.clear();
834 >
835 >  kz.clear();
836 >
837 >  cantPos.clear();
838 >
839 >  cantVel.clear();
840 >
841 >
842 >
843 >  unconsMols.clear();
844 >
845 >  massOfUnconsMols.clear();
846 >
847 >
848 >
849 >
850 >
851 >  //creat zconsMol and unconsMol lists
852 >
853 >  for (int i = 0; i < nMols; i++){
854 >
855 >    index = isZConstraintMol(&molecules[i]);
856 >
857 >
858 >
859 >    if (index > -1){
860 >
861 >      zconsMols.push_back(&molecules[i]);      
862 >
863 >      zPos.push_back((*parameters)[index].zPos);
864 >
865 >      kz.push_back((*parameters)[index].kRatio * zForceConst);
866 >
867 >      massOfZConsMols.push_back(molecules[i].getTotalMass());
868 >
869 >      
870 >
871 >      if(usingSMD)
872 >
873 >        cantVel.push_back((*parameters)[index].cantVel);
874 >
875 >
876 >
877 >    }
878 >
879 >    else{
880 >
881 >      unconsMols.push_back(&molecules[i]);
882 >
883 >      massOfUnconsMols.push_back(molecules[i].getTotalMass());
884 >
885 >    }
886 >
887 >  }
888 >
889 >
890 >
891 >  fz.resize(zconsMols.size());
892 >
893 >  curZPos.resize(zconsMols.size());
894 >
895 >  indexOfZConsMols.resize(zconsMols.size());  
896 >
897 >
898 >
899 >  for (size_t i = 0; i < zconsMols.size(); i++){
900 >
901 >    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
902 >
903 >  }
904 >
905 >    
906 >
907 >  //determine the states of z-constraint molecules
908 >
909 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
910 >
911 >
912 >
913 >    zconsMols[i]->getCOM(COM);
914 >
915 >    
916 >
917 >    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
918 >
919 >      states.push_back(zcsFixed);
920 >
921 >
922 >
923 >      if (hasZConsGap)
924 >
925 >        endFixTime.push_back(info->getTime() + zconsFixTime);
926 >
927 >    }
928 >
929 >    else{
930 >
931 >      states.push_back(zcsMoving);
932 >
933 >
934 >
935 >      if (hasZConsGap)
936 >
937 >        endFixTime.push_back(INFINITE_TIME);
938 >
939 >    }
940 >
941 >
942 >
943 >    if(usingSMD)
944 >
945 >      cantPos.push_back(COM[whichDirection]);        
946 >
947 >  }
948 >
949 >
950 >
951 >  if(usingSMD)
952 >
953 >  prevCantPos = cantPos;
954 >
955 >
956 >
957 >  forcePolicy->update();
958 >
959 > }
960 >
961 >
962 >
963 > #endif
964 >
965 >
966 >
967 > /**
968 >
969 > *  Function Name: isZConstraintMol
970 >
971 > *  Parameter
972 >
973 > *    Molecule* mol
974 >
975 > *  Return value:
976 >
977 > *    -1, if the molecule is not z-constraint molecule,
978 >
979 > *    other non-negative values, its index in indexOfAllZConsMols vector
980 >
981 > */
982 >
983 >
984 >
985 > template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){
986 >
987 >  int index;
988 >
989 >  int low;
990 >
991 >  int high;
992 >
993 >  int mid;
994 >
995 >
996 >
997 >  index = mol->getGlobalIndex();
998 >
999 >
1000 >
1001 >  low = 0;
1002 >
1003 >  high = parameters->size() - 1;
1004 >
1005 >
1006 >
1007 >  //Binary Search (we have sorted the array)  
1008 >
1009 >  while (low <= high){
1010 >
1011 >    mid = (low + high) / 2;
1012 >
1013 >    if ((*parameters)[mid].zconsIndex == index)
1014 >
1015 >      return mid;
1016 >
1017 >    else if ((*parameters)[mid].zconsIndex > index)
1018 >
1019 >      high = mid - 1;
1020 >
1021 >    else
1022 >
1023 >      low = mid + 1;
1024 >
1025 >  }
1026 >
1027 >
1028 >
1029 >  return -1;
1030 >
1031 > }
1032 >
1033 >
1034 >
1035 > template<typename T> void ZConstraint<T>::integrate(){
1036 >
1037 >  // creat zconsWriter  
1038 >
1039 >  fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);  
1040 >
1041 >
1042 >
1043 >  if (!fzOut){
1044 >
1045 >    sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
1046 >
1047 >    painCave.isFatal = 1;
1048 >
1049 >    simError();
1050 >
1051 >  }
1052 >
1053 >
1054 >
1055 >  //zero out the velocities of center of mass of unconstrained molecules
1056 >
1057 >  //and the velocities of center of mass of every single z-constrained molecueles
1058 >
1059 >  zeroOutVel();
1060 >
1061 >
1062 >
1063 >  curZconsTime = zconsTime + info->getTime();
1064 >
1065 >
1066 >
1067 >  T::integrate();
1068 >
1069 > }
1070 >
1071 >
1072 >
1073 > template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
1074 >
1075 >  double zsys;
1076 >
1077 >  double COM[3];
1078 >
1079 >  double force[3];
1080 >
1081 >  double zSysCOMVel;
1082 >
1083 >
1084 >
1085 >  T::calcForce(calcPot, calcStress);
1086 >
1087 >
1088 >
1089 >
1090 >
1091 >  if (hasZConsGap){
1092 >
1093 >    updateZPos();
1094 >
1095 >  }
1096 >
1097 >
1098 >
1099 >  if (checkZConsState()){
1100 >
1101 >    zeroOutVel();    
1102 >
1103 >    forcePolicy->update();
1104 >
1105 >  }  
1106 >
1107 >
1108 >
1109 >  zsys = calcZSys();
1110 >
1111 >  zSysCOMVel = calcSysCOMVel();
1112 >
1113 > #ifdef IS_MPI
1114 >
1115 >  if (worldRank == 0){
1116 >
1117 > #endif
1118 >
1119 >
1120 >
1121 > #ifdef IS_MPI
1122 >
1123 >  }
1124 >
1125 > #endif
1126 >
1127 >
1128 >
1129 >  //do zconstraint force;
1130 >
1131 >  if (haveFixedZMols()){
1132 >
1133 >    this->doZconstraintForce();
1134 >
1135 >  }
1136 >
1137 >
1138 >
1139 >  //use external force to move the molecules to the specified positions
1140 >
1141 >  if (haveMovingZMols()){
1142 >
1143 >    if (usingSMD)
1144 >
1145 >      this->doHarmonic(cantPos);
1146 >
1147 >    else
1148 >
1149 >      this->doHarmonic(zPos);      
1150 >
1151 >  }
1152 >
1153 >
1154 >
1155 >  //write out forces and current positions of z-constraint molecules
1156 >
1157 >  if (info->getTime() >= curZconsTime){
1158 >
1159 >    for (int i = 0; i < (int) (zconsMols.size()); i++){
1160 >
1161 >      zconsMols[i]->getCOM(COM);
1162 >
1163 >      curZPos[i] = COM[whichDirection];
1164 >
1165 >
1166 >
1167 >      //if the z-constraint molecule is still moving, just record its force
1168 >
1169 >      if (states[i] == zcsMoving){
1170 >
1171 >        fz[i] = 0;
1172 >
1173 >        Atom** movingZAtoms;
1174 >
1175 >        movingZAtoms = zconsMols[i]->getMyAtoms();
1176 >
1177 >        for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
1178 >
1179 >          movingZAtoms[j]->getFrc(force);
1180 >
1181 >          fz[i] += force[whichDirection];
1182 >
1183 >        }
1184 >
1185 >      }
1186 >
1187 >    }
1188 >
1189 >    fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
1190 >
1191 >                   &curZPos[0], &zPos[0]);
1192 >
1193 >    curZconsTime += zconsTime;
1194 >
1195 >  }
1196 >
1197 >
1198 >
1199 >  zSysCOMVel = calcSysCOMVel();  
1200 >
1201 > #ifdef IS_MPI
1202 >
1203 >  if (worldRank == 0){
1204 >
1205 > #endif
1206 >
1207 > #ifdef IS_MPI
1208 >
1209 >  }
1210 >
1211 > #endif
1212 >
1213 > }
1214 >
1215 >
1216 >
1217 >
1218 >
1219 > template<typename T> double ZConstraint<T>::calcZSys(){
1220 >
1221 >  //calculate reference z coordinate for z-constraint molecules
1222 >
1223 >  double totalMass_local;
1224 >
1225 >  double totalMass;
1226 >
1227 >  double totalMZ_local;
1228 >
1229 >  double totalMZ;
1230 >
1231 >  double massOfCurMol;
1232 >
1233 >  double COM[3];
1234 >
1235 >
1236 >
1237 >  totalMass_local = 0;
1238 >
1239 >  totalMZ_local = 0;
1240 >
1241 >
1242 >
1243 >  for (int i = 0; i < nMols; i++){
1244 >
1245 >    massOfCurMol = molecules[i].getTotalMass();
1246 >
1247 >    molecules[i].getCOM(COM);
1248 >
1249 >
1250 >
1251 >    totalMass_local += massOfCurMol;
1252 >
1253 >    totalMZ_local += massOfCurMol * COM[whichDirection];
1254 >
1255 >  }
1256 >
1257 >
1258 >
1259 >
1260 >
1261 > #ifdef IS_MPI  
1262 >
1263 >  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM,
1264 >
1265 >                MPI_COMM_WORLD);
1266 >
1267 >  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
1268 >
1269 > #else
1270 >
1271 >  totalMass = totalMass_local;
1272 >
1273 >  totalMZ = totalMZ_local;
1274 >
1275 > #endif  
1276 >
1277 >
1278 >
1279 >  double zsys;
1280 >
1281 >  zsys = totalMZ / totalMass;
1282 >
1283 >
1284 >
1285 >  return zsys;
1286 >
1287 > }
1288 >
1289 >
1290 >
1291 > template<typename T> void ZConstraint<T>::thermalize(void){
1292 >
1293 >  T::thermalize();
1294 >
1295 >  zeroOutVel();
1296 >
1297 > }
1298 >
1299 >
1300 >
1301 > template<typename T> void ZConstraint<T>::zeroOutVel(){
1302 >
1303 >  Atom** fixedZAtoms;  
1304 >
1305 >  double COMvel[3];
1306 >
1307 >  double vel[3];
1308 >
1309 >  double zSysCOMVel;
1310 >
1311 >
1312 >
1313 >  //zero out the velocities of center of mass of fixed z-constrained molecules
1314 >
1315 >
1316 >
1317 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
1318 >
1319 >    if (states[i] == zcsFixed){
1320 >
1321 >      zconsMols[i]->getCOMvel(COMvel);      
1322 >
1323 >      //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << std::endl;
1324 >
1325 >
1326 >
1327 >      fixedZAtoms = zconsMols[i]->getMyAtoms();
1328 >
1329 >
1330 >
1331 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
1332 >
1333 >        vel = fixedZAtoms[j]->getVel();
1334 >
1335 >        vel[whichDirection] -= COMvel[whichDirection];
1336 >
1337 >        fixedZAtoms[j]->setVel(vel);
1338 >
1339 >      }
1340 >
1341 >
1342 >
1343 >      zconsMols[i]->getCOMvel(COMvel);
1344 >
1345 >    }
1346 >
1347 >  }
1348 >
1349 >
1350 >
1351 >  zSysCOMVel = calcSysCOMVel();
1352 >
1353 > #ifdef IS_MPI
1354 >
1355 >  if (worldRank == 0){
1356 >
1357 > #endif
1358 >
1359 > #ifdef IS_MPI
1360 >
1361 >  }
1362 >
1363 > #endif
1364 >
1365 >
1366 >
1367 >  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
1368 >
1369 >  double MVzOfMovingMols_local;
1370 >
1371 >  double MVzOfMovingMols;
1372 >
1373 >  double totalMassOfMovingZMols_local;
1374 >
1375 >  double totalMassOfMovingZMols;
1376 >
1377 >
1378 >
1379 >  MVzOfMovingMols_local = 0;
1380 >
1381 >  totalMassOfMovingZMols_local = 0;
1382 >
1383 >
1384 >
1385 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
1386 >
1387 >    unconsMols[i]->getCOMvel(COMvel);
1388 >
1389 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1390 >
1391 >  }
1392 >
1393 >
1394 >
1395 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
1396 >
1397 >    if (states[i] == zcsMoving){
1398 >
1399 >      zconsMols[i]->getCOMvel(COMvel);
1400 >
1401 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
1402 >
1403 >      totalMassOfMovingZMols_local += massOfZConsMols[i];
1404 >
1405 >    }
1406 >
1407 >  }
1408 >
1409 >
1410 >
1411 > #ifndef IS_MPI
1412 >
1413 >  MVzOfMovingMols = MVzOfMovingMols_local;
1414 >
1415 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1416 >
1417 > #else
1418 >
1419 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
1420 >
1421 >                MPI_SUM, MPI_COMM_WORLD);
1422 >
1423 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
1424 >
1425 >                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
1426 >
1427 > #endif
1428 >
1429 >
1430 >
1431 >  double vzOfMovingMols;
1432 >
1433 >  vzOfMovingMols = MVzOfMovingMols /
1434 >
1435 >                   (totalMassOfUncons + totalMassOfMovingZMols);
1436 >
1437 >
1438 >
1439 >  //modify the velocites of unconstrained molecules  
1440 >
1441 >  Atom** unconsAtoms;
1442 >
1443 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
1444 >
1445 >    unconsAtoms = unconsMols[i]->getMyAtoms();
1446 >
1447 >    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
1448 >
1449 >      vel = unconsAtoms[j]->getVel();
1450 >
1451 >      vel[whichDirection] -= vzOfMovingMols;
1452 >
1453 >      unconsAtoms[j]->setVel(vel);
1454 >
1455 >    }
1456 >
1457 >  }  
1458 >
1459 >
1460 >
1461 >  //modify the velocities of moving z-constrained molecuels
1462 >
1463 >  Atom** movingZAtoms;
1464 >
1465 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
1466 >
1467 >    if (states[i] == zcsMoving){
1468 >
1469 >      movingZAtoms = zconsMols[i]->getMyAtoms();
1470 >
1471 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
1472 >
1473 >        vel = movingZAtoms[j]->getVel();
1474 >
1475 >        vel[whichDirection] -= vzOfMovingMols;
1476 >
1477 >        movingZAtoms[j]->setVel(vel);
1478 >
1479 >      }
1480 >
1481 >    }
1482 >
1483 >  }
1484 >
1485 >
1486 >
1487 >
1488 >
1489 >  zSysCOMVel = calcSysCOMVel();
1490 >
1491 > #ifdef IS_MPI
1492 >
1493 >  if (worldRank == 0){
1494 >
1495 > #endif
1496 >
1497 > #ifdef IS_MPI
1498 >
1499 >  }
1500 >
1501 > #endif
1502 >
1503 > }
1504 >
1505 >
1506 >
1507 >
1508 >
1509 > template<typename T> void ZConstraint<T>::doZconstraintForce(){
1510 >
1511 >  Atom** zconsAtoms;
1512 >
1513 >  double totalFZ;
1514 >
1515 >  double totalFZ_local;
1516 >
1517 >  double COM[3];
1518 >
1519 >  double force[3];
1520 >
1521 >
1522 >
1523 >  //constrain the molecules which do not reach the specified positions  
1524 >
1525 >
1526 >
1527 >  //Zero Out the force of z-contrained molecules    
1528 >
1529 >  totalFZ_local = 0;
1530 >
1531 >
1532 >
1533 >  //calculate the total z-contrained force of fixed z-contrained molecules
1534 >
1535 >
1536 >
1537 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
1538 >
1539 >    if (states[i] == zcsFixed){
1540 >
1541 >      zconsMols[i]->getCOM(COM);
1542 >
1543 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
1544 >
1545 >
1546 >
1547 >      fz[i] = 0;      
1548 >
1549 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
1550 >
1551 >        zconsAtoms[j]->getFrc(force);
1552 >
1553 >        fz[i] += force[whichDirection];
1554 >
1555 >      }
1556 >
1557 >      totalFZ_local += fz[i];
1558 >
1559 >
1560 >
1561 >    }
1562 >
1563 >  }
1564 >
1565 >
1566 >
1567 >  //calculate total z-constraint force
1568 >
1569 > #ifdef IS_MPI
1570 >
1571 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
1572 >
1573 > #else
1574 >
1575 >  totalFZ = totalFZ_local;
1576 >
1577 > #endif
1578 >
1579 >
1580 >
1581 >
1582 >
1583 >  // apply negative to fixed z-constrained molecues;
1584 >
1585 >  force[0] = 0;
1586 >
1587 >  force[1] = 0;
1588 >
1589 >  force[2] = 0;
1590 >
1591 >
1592 >
1593 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
1594 >
1595 >    if (states[i] == zcsFixed){
1596 >
1597 >      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
1598 >
1599 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
1600 >
1601 >
1602 >
1603 >      for (int j = 0; j < nAtomOfCurZConsMol; j++){
1604 >
1605 >        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
1606 >
1607 >        force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i],
1608 >
1609 >                                                                zconsAtoms[j],
1610 >
1611 >                                                                fz[i]);
1612 >
1613 >        zconsAtoms[j]->addFrc(force);
1614 >
1615 >      }
1616 >
1617 >    }
1618 >
1619 >  }
1620 >
1621 >
1622 >
1623 >  force[0] = 0;
1624 >
1625 >  force[1] = 0;
1626 >
1627 >  force[2] = 0;
1628 >
1629 >
1630 >
1631 >  //modify the forces of unconstrained molecules
1632 >
1633 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
1634 >
1635 >    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
1636 >
1637 >
1638 >
1639 >    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
1640 >
1641 >      //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
1642 >
1643 >      force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],
1644 >
1645 >                                                             totalFZ);
1646 >
1647 >      unconsAtoms[j]->addFrc(force);
1648 >
1649 >    }
1650 >
1651 >  }      
1652 >
1653 >
1654 >
1655 >  //modify the forces of moving z-constrained molecules
1656 >
1657 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
1658 >
1659 >    if (states[i] == zcsMoving){
1660 >
1661 >      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
1662 >
1663 >
1664 >
1665 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
1666 >
1667 >        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
1668 >
1669 >        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],
1670 >
1671 >                                                               totalFZ);
1672 >
1673 >        movingZAtoms[j]->addFrc(force);
1674 >
1675 >      }
1676 >
1677 >    }
1678 >
1679 >  }
1680 >
1681 > }
1682 >
1683 >
1684 >
1685 >
1686 >
1687 > template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
1688 >
1689 >  double force[3];
1690 >
1691 >  double harmonicU;
1692 >
1693 >  double harmonicF;
1694 >
1695 >  double COM[3];
1696 >
1697 >  double diff;
1698 >
1699 >  double totalFZ_local;
1700 >
1701 >  double totalFZ;
1702 >
1703 >
1704 >
1705 >  force[0] = 0;
1706 >
1707 >  force[1] = 0;
1708 >
1709 >  force[2] = 0;
1710 >
1711 >
1712 >
1713 >  totalFZ_local = 0;
1714 >
1715 >
1716 >
1717 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
1718 >
1719 >    if (states[i] == zcsMoving){
1720 >
1721 >      zconsMols[i]->getCOM(COM);
1722 >
1723 >
1724 >
1725 >      diff = COM[whichDirection] - resPos[i];
1726 >
1727 >
1728 >
1729 >      harmonicU = 0.5 * kz[i] * diff * diff;  
1730 >
1731 >      info->lrPot += harmonicU;
1732 >
1733 >
1734 >
1735 >      harmonicF = -kz[i] * diff;
1736 >
1737 >      totalFZ_local += harmonicF;
1738 >
1739 >
1740 >
1741 >      //adjust force
1742 >
1743 >
1744 >
1745 >      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
1746 >
1747 >
1748 >
1749 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
1750 >
1751 >        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i],
1752 >
1753 >                                                               movingZAtoms[j],
1754 >
1755 >                                                               harmonicF);
1756 >
1757 >        movingZAtoms[j]->addFrc(force);
1758 >
1759 >      }
1760 >
1761 >    }
1762 >
1763 >  }
1764 >
1765 >
1766 >
1767 > #ifndef IS_MPI
1768 >
1769 >  totalFZ = totalFZ_local;
1770 >
1771 > #else
1772 >
1773 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
1774 >
1775 > #endif
1776 >
1777 >
1778 >
1779 >  force[0] = 0;
1780 >
1781 >  force[1] = 0;
1782 >
1783 >  force[2] = 0;
1784 >
1785 >
1786 >
1787 >  //modify the forces of unconstrained molecules
1788 >
1789 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
1790 >
1791 >    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
1792 >
1793 >
1794 >
1795 >    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
1796 >
1797 >      //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
1798 >
1799 >      force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j],
1800 >
1801 >                                                              totalFZ);
1802 >
1803 >      unconsAtoms[j]->addFrc(force);
1804 >
1805 >    }
1806 >
1807 >  }  
1808 >
1809 >
1810 >
1811 > }
1812 >
1813 >
1814 >
1815 > template<typename T> bool ZConstraint<T>::checkZConsState(){
1816 >
1817 >  double COM[3];
1818 >
1819 >  double diff;
1820 >
1821 >
1822 >
1823 >  int changed_local;
1824 >
1825 >  int changed;
1826 >
1827 >
1828 >
1829 >  changed_local = 0;
1830 >
1831 >
1832 >
1833 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
1834 >
1835 >    zconsMols[i]->getCOM(COM);
1836 >
1837 >    diff = fabs(COM[whichDirection] - zPos[i]);  
1838 >
1839 >    if (diff <= zconsTol && states[i] == zcsMoving){
1840 >
1841 >      states[i] = zcsFixed;
1842 >
1843 >      changed_local = 1;
1844 >
1845 >
1846 >
1847 >      if(usingSMD)
1848 >
1849 >        prevCantPos = cantPos;
1850 >
1851 >
1852 >
1853 >      if (hasZConsGap)
1854 >
1855 >        endFixTime[i] = info->getTime() + zconsFixTime;
1856 >
1857 >    }
1858 >
1859 >    else if (diff > zconsTol && states[i] == zcsFixed){
1860 >
1861 >      states[i] = zcsMoving;
1862 >
1863 >      changed_local = 1;  
1864 >
1865 >
1866 >
1867 >      if(usingSMD)
1868 >
1869 >         cantPos = prevCantPos;
1870 >
1871 >      
1872 >
1873 >      if (hasZConsGap)
1874 >
1875 >        endFixTime[i] = INFINITE_TIME;
1876 >
1877 >    }
1878 >
1879 >  }
1880 >
1881 >
1882 >
1883 > #ifndef IS_MPI
1884 >
1885 >  changed = changed_local;
1886 >
1887 > #else
1888 >
1889 >  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1890 >
1891 > #endif
1892 >
1893 >
1894 >
1895 >  return (changed > 0);
1896 >
1897 > }
1898 >
1899 >
1900 >
1901 > template<typename T> bool ZConstraint<T>::haveFixedZMols(){
1902 >
1903 >  int havingFixed_local;
1904 >
1905 >  int havingFixed;
1906 >
1907 >
1908 >
1909 >  havingFixed_local = 0;
1910 >
1911 >
1912 >
1913 >  for (int i = 0; i < (int) (zconsMols.size()); i++)
1914 >
1915 >    if (states[i] == zcsFixed){
1916 >
1917 >      havingFixed_local = 1;
1918 >
1919 >      break;
1920 >
1921 >    }
1922 >
1923 >
1924 >
1925 > #ifndef IS_MPI
1926 >
1927 >  havingFixed = havingFixed_local;
1928 >
1929 > #else
1930 >
1931 >  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
1932 >
1933 >                MPI_COMM_WORLD);
1934 >
1935 > #endif
1936 >
1937 >
1938 >
1939 >  return (havingFixed > 0);
1940 >
1941 > }
1942 >
1943 >
1944 >
1945 >
1946 >
1947 > template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1948 >
1949 >  int havingMoving_local;
1950 >
1951 >  int havingMoving;
1952 >
1953 >
1954 >
1955 >  havingMoving_local = 0;
1956 >
1957 >
1958 >
1959 >  for (int i = 0; i < (int) (zconsMols.size()); i++)
1960 >
1961 >    if (states[i] == zcsMoving){
1962 >
1963 >      havingMoving_local = 1;
1964 >
1965 >      break;
1966 >
1967 >    }
1968 >
1969 >
1970 >
1971 > #ifndef IS_MPI
1972 >
1973 >  havingMoving = havingMoving_local;
1974 >
1975 > #else
1976 >
1977 >  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
1978 >
1979 >                MPI_COMM_WORLD);
1980 >
1981 > #endif
1982 >
1983 >
1984 >
1985 >  return (havingMoving > 0);
1986 >
1987 > }
1988 >
1989 >
1990 >
1991 >
1992 >
1993 > template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){
1994 >
1995 >  double MVzOfMovingMols_local;
1996 >
1997 >  double MVzOfMovingMols;
1998 >
1999 >  double totalMassOfMovingZMols_local;
2000 >
2001 >  double totalMassOfMovingZMols;
2002 >
2003 >  double COMvel[3];
2004 >
2005 >
2006 >
2007 >  MVzOfMovingMols_local = 0;
2008 >
2009 >  totalMassOfMovingZMols_local = 0;
2010 >
2011 >
2012 >
2013 >  for (int i = 0; i < unconsMols.size(); i++){
2014 >
2015 >    unconsMols[i]->getCOMvel(COMvel);
2016 >
2017 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
2018 >
2019 >  }
2020 >
2021 >
2022 >
2023 >  for (int i = 0; i < zconsMols.size(); i++){
2024 >
2025 >    if (states[i] == zcsMoving){
2026 >
2027 >      zconsMols[i]->getCOMvel(COMvel);
2028 >
2029 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
2030 >
2031 >      totalMassOfMovingZMols_local += massOfZConsMols[i];
2032 >
2033 >    }
2034 >
2035 >  }
2036 >
2037 >
2038 >
2039 > #ifndef IS_MPI
2040 >
2041 >  MVzOfMovingMols = MVzOfMovingMols_local;
2042 >
2043 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
2044 >
2045 > #else
2046 >
2047 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
2048 >
2049 >                MPI_SUM, MPI_COMM_WORLD);
2050 >
2051 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
2052 >
2053 >                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
2054 >
2055 > #endif
2056 >
2057 >
2058 >
2059 >  double vzOfMovingMols;
2060 >
2061 >  vzOfMovingMols = MVzOfMovingMols /
2062 >
2063 >                   (totalMassOfUncons + totalMassOfMovingZMols);
2064 >
2065 >
2066 >
2067 >  return vzOfMovingMols;
2068 >
2069 > }
2070 >
2071 >
2072 >
2073 > template<typename T> double ZConstraint<T>::calcSysCOMVel(){
2074 >
2075 >  double COMvel[3];
2076 >
2077 >  double tempMVz_local;
2078 >
2079 >  double tempMVz;
2080 >
2081 >  double massOfZCons_local;
2082 >
2083 >  double massOfZCons;
2084 >
2085 >
2086 >
2087 >
2088 >
2089 >  tempMVz_local = 0;
2090 >
2091 >
2092 >
2093 >  for (int i = 0 ; i < nMols; i++){
2094 >
2095 >    molecules[i].getCOMvel(COMvel);
2096 >
2097 >    tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection];
2098 >
2099 >  }
2100 >
2101 >
2102 >
2103 >  massOfZCons_local = 0;
2104 >
2105 >
2106 >
2107 >  for (int i = 0; i < (int) (massOfZConsMols.size()); i++){
2108 >
2109 >    massOfZCons_local += massOfZConsMols[i];
2110 >
2111 >  }
2112 >
2113 > #ifndef IS_MPI
2114 >
2115 >  massOfZCons = massOfZCons_local;
2116 >
2117 >  tempMVz = tempMVz_local;
2118 >
2119 > #else
2120 >
2121 >  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM,
2122 >
2123 >                MPI_COMM_WORLD);
2124 >
2125 >  MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
2126 >
2127 > #endif
2128 >
2129 >
2130 >
2131 >  return tempMVz / (totalMassOfUncons + massOfZCons);
2132 >
2133 > }
2134 >
2135 >
2136 >
2137 > template<typename T> double ZConstraint<T>::calcTotalForce(){
2138 >
2139 >  double force[3];  
2140 >
2141 >  double totalForce_local;
2142 >
2143 >  double totalForce;
2144 >
2145 >
2146 >
2147 >  totalForce_local = 0;
2148 >
2149 >
2150 >
2151 >  for (int i = 0; i < nAtoms; i++){
2152 >
2153 >    atoms[i]->getFrc(force);
2154 >
2155 >    totalForce_local += force[whichDirection];
2156 >
2157 >  }
2158 >
2159 >
2160 >
2161 > #ifndef IS_MPI
2162 >
2163 >  totalForce = totalForce_local;
2164 >
2165 > #else
2166 >
2167 >  MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM,
2168 >
2169 >                MPI_COMM_WORLD);
2170 >
2171 > #endif
2172 >
2173 >
2174 >
2175 >  return totalForce;
2176 >
2177 > }
2178 >
2179 >
2180 >
2181 > template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
2182 >
2183 >  //calculate the number of atoms of moving z-constrained molecules
2184 >
2185 >  int nMovingZAtoms_local;
2186 >
2187 >  int nMovingZAtoms;
2188 >
2189 >
2190 >
2191 >  nMovingZAtoms_local = 0;
2192 >
2193 >  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
2194 >
2195 >    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
2196 >
2197 >      nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
2198 >
2199 >    }
2200 >
2201 >
2202 >
2203 > #ifdef IS_MPI
2204 >
2205 >  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM,
2206 >
2207 >                MPI_COMM_WORLD);
2208 >
2209 > #else
2210 >
2211 >  nMovingZAtoms = nMovingZAtoms_local;
2212 >
2213 > #endif
2214 >
2215 >  totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
2216 >
2217 > }
2218 >
2219 >
2220 >
2221 > template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol,
2222 >
2223 >                                                                              Atom* atom,
2224 >
2225 >                                                                              double totalForce){
2226 >
2227 >  return totalForce / mol->getNAtoms();
2228 >
2229 > }
2230 >
2231 >
2232 >
2233 > template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom,
2234 >
2235 >                                                                              double totalForce){
2236 >
2237 >  return totalForce / totNumOfMovingAtoms;
2238 >
2239 > }
2240 >
2241 >
2242 >
2243 > template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol,
2244 >
2245 >                                                                              Atom* atom,
2246 >
2247 >                                                                              double totalForce){
2248 >
2249 >  return totalForce / mol->getNAtoms();
2250 >
2251 > }
2252 >
2253 >
2254 >
2255 > template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom,
2256 >
2257 >                                                                              double totalForce){
2258 >
2259 >  return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
2260 >
2261 > }
2262 >
2263 >
2264 >
2265 >
2266 >
2267 > template<typename T> void ZConstraint<T>::PolicyByMass::update(){
2268 >
2269 >  //calculate the number of atoms of moving z-constrained molecules
2270 >
2271 >  double massOfMovingZAtoms_local;
2272 >
2273 >  double massOfMovingZAtoms;
2274 >
2275 >
2276 >
2277 >  massOfMovingZAtoms_local = 0;
2278 >
2279 >  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
2280 >
2281 >    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
2282 >
2283 >      massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
2284 >
2285 >    }
2286 >
2287 >
2288 >
2289 > #ifdef IS_MPI
2290 >
2291 >  MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,
2292 >
2293 >                MPI_SUM, MPI_COMM_WORLD);
2294 >
2295 > #else
2296 >
2297 >  massOfMovingZAtoms = massOfMovingZAtoms_local;
2298 >
2299 > #endif
2300 >
2301 >  totMassOfMovingAtoms = massOfMovingZAtoms +
2302 >
2303 >                         zconsIntegrator->totalMassOfUncons;
2304 >
2305 > }
2306 >
2307 >
2308 >
2309 > template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol,
2310 >
2311 >                                                                            Atom* atom,
2312 >
2313 >                                                                            double totalForce){
2314 >
2315 >  return totalForce * atom->getMass() / mol->getTotalMass();
2316 >
2317 > }
2318 >
2319 >
2320 >
2321 > template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom,
2322 >
2323 >                                                                            double totalForce){
2324 >
2325 >  return totalForce * atom->getMass() / totMassOfMovingAtoms;
2326 >
2327 > }
2328 >
2329 >
2330 >
2331 > template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol,
2332 >
2333 >                                                                            Atom* atom,
2334 >
2335 >                                                                            double totalForce){
2336 >
2337 >  return totalForce * atom->getMass() / mol->getTotalMass();
2338 >
2339 > }
2340 >
2341 >
2342 >
2343 > template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom,
2344 >
2345 >                                                                            double totalForce){
2346 >
2347 >  return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
2348 >
2349 > }
2350 >
2351 >
2352 >
2353 > template<typename T> void ZConstraint<T>::updateZPos(){
2354 >
2355 >  double curTime;
2356 >
2357 >  double COM[3];
2358 >
2359 >  
2360 >
2361 >  curTime = info->getTime();
2362 >
2363 >
2364 >
2365 >  for (size_t i = 0; i < zconsMols.size(); i++){
2366 >
2367 >
2368 >
2369 >    if (states[i] == zcsFixed && curTime >= endFixTime[i]){
2370 >
2371 >      zPos[i] += zconsGap;
2372 >
2373 >
2374 >
2375 >      if (usingSMD){
2376 >
2377 >        zconsMols[i]->getCOM(COM);
2378 >
2379 >        cantPos[i] = COM[whichDirection];
2380 >
2381 >      }
2382 >
2383 >      
2384 >
2385 >    }
2386 >
2387 >    
2388 >
2389 >  }
2390 >
2391 >  
2392 >
2393 > }
2394 >
2395 >
2396 >
2397 > template<typename T> void ZConstraint<T>::updateCantPos(){
2398 >
2399 >  double curTime;
2400 >
2401 >  double dt;
2402 >
2403 >
2404 >
2405 >  curTime = info->getTime();
2406 >
2407 >  dt = info->dt;
2408 >
2409 >
2410 >
2411 >  for (size_t i = 0; i < zconsMols.size(); i++){
2412 >
2413 >    if (states[i] == zcsMoving){
2414 >
2415 >      cantPos[i] += cantVel[i] * dt;
2416 >
2417 >    }
2418 >
2419 >  }
2420 >
2421 >
2422 >
2423 > }
2424 >

Diff Legend

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