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 1141 by tim, Wed Apr 28 23:09:32 2004 UTC vs.
Revision 1203 by gezelter, Thu May 27 18:59:17 2004 UTC

# Line 1 | Line 1
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 <  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 <  //
480 <  forcePolicy->update();
481 < }
482 <
483 < #endif
484 <
485 < /**
486 < *  Function Name: isZConstraintMol
487 < *  Parameter
488 < *    Molecule* mol
489 < *  Return value:
490 < *    -1, if the molecule is not z-constraint molecule,
491 < *    other non-negative values, its index in indexOfAllZConsMols vector
492 < */
493 <
494 < template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){
495 <  int index;
496 <  int low;
497 <  int high;
498 <  int mid;
499 <
500 <  index = mol->getGlobalIndex();
501 <
502 <  low = 0;
503 <  high = parameters->size() - 1;
504 <
505 <  //Binary Search (we have sorted the array)  
506 <  while (low <= high){
507 <    mid = (low + high) / 2;
508 <    if ((*parameters)[mid].zconsIndex == index)
509 <      return mid;
510 <    else if ((*parameters)[mid].zconsIndex > index)
511 <      high = mid - 1;
512 <    else
513 <      low = mid + 1;
514 <  }
515 <
516 <  return -1;
517 < }
518 <
519 < template<typename T> void ZConstraint<T>::integrate(){
520 <  // creat zconsWriter  
521 <  fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);  
522 <
523 <  if (!fzOut){
524 <    sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
525 <    painCave.isFatal = 1;
526 <    simError();
527 <  }
528 <
529 <  //zero out the velocities of center of mass of unconstrained molecules
530 <  //and the velocities of center of mass of every single z-constrained molecueles
531 <  zeroOutVel();
532 <
533 <  curZconsTime = zconsTime + info->getTime();
534 <
535 <  T::integrate();
536 < }
537 <
538 <
539 < /**
540 < *
541 < *
542 < *
543 < *
544 < */
545 < template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
546 <  double zsys;
547 <  double COM[3];
548 <  double force[3];
549 <  double zSysCOMVel;
550 <
551 <  T::calcForce(calcPot, calcStress);
552 <
553 <
554 <  if (hasZConsGap){
555 <    updateZPos();
556 <  }
557 <
558 <  if (checkZConsState()){
559 <    zeroOutVel();    
560 <    forcePolicy->update();
561 <  }  
562 <
563 <  zsys = calcZSys();
564 <  zSysCOMVel = calcSysCOMVel();
565 < #ifdef IS_MPI
566 <  if (worldRank == 0){
567 < #endif
568 <    //cout << "---------------------------------------------------------------------" <<endl;
569 <    //cout << "current time: " << info->getTime() << endl;
570 <    //cout << "center of mass at z: " << zsys << endl;    
571 <    //cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
572 <
573 < #ifdef IS_MPI
574 <  }
575 < #endif
576 <
577 <  //do zconstraint force;
578 <  if (haveFixedZMols()){
579 <    this->doZconstraintForce();
580 <  }
581 <
582 <  //use external force to move the molecules to the specified positions
583 <  if (haveMovingZMols()){
584 <    if (usingSMD)
585 <      this->doHarmonic(cantPos);
586 <    else
587 <      this->doHarmonic(zPos);      
588 <  }
589 <
590 <  //write out forces and current positions of z-constraint molecules
591 <  if (info->getTime() >= curZconsTime){
592 <    for (int i = 0; i < (int) (zconsMols.size()); i++){
593 <      zconsMols[i]->getCOM(COM);
594 <      curZPos[i] = COM[whichDirection];
595 <
596 <      //if the z-constraint molecule is still moving, just record its force
597 <      if (states[i] == zcsMoving){
598 <        fz[i] = 0;
599 <        Atom** movingZAtoms;
600 <        movingZAtoms = zconsMols[i]->getMyAtoms();
601 <        for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
602 <          movingZAtoms[j]->getFrc(force);
603 <          fz[i] += force[whichDirection];
604 <        }
605 <      }
606 <    }
607 <    fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
608 <                   &curZPos[0], &zPos[0]);
609 <    curZconsTime += zconsTime;
610 <  }
611 <
612 <  zSysCOMVel = calcSysCOMVel();  
613 < #ifdef IS_MPI
614 <  if (worldRank == 0){
615 < #endif
616 <    //cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
617 < #ifdef IS_MPI
618 <  }
619 < #endif
620 < }
621 <
622 <
623 < /**
624 < *
625 < */
626 <
627 < template<typename T> double ZConstraint<T>::calcZSys(){
628 <  //calculate reference z coordinate for z-constraint molecules
629 <  double totalMass_local;
630 <  double totalMass;
631 <  double totalMZ_local;
632 <  double totalMZ;
633 <  double massOfCurMol;
634 <  double COM[3];
635 <
636 <  totalMass_local = 0;
637 <  totalMZ_local = 0;
638 <
639 <  for (int i = 0; i < nMols; i++){
640 <    massOfCurMol = molecules[i].getTotalMass();
641 <    molecules[i].getCOM(COM);
642 <
643 <    totalMass_local += massOfCurMol;
644 <    totalMZ_local += massOfCurMol * COM[whichDirection];
645 <  }
646 <
647 <
648 < #ifdef IS_MPI  
649 <  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM,
650 <                MPI_COMM_WORLD);
651 <  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
652 < #else
653 <  totalMass = totalMass_local;
654 <  totalMZ = totalMZ_local;
655 < #endif  
656 <
657 <  double zsys;
658 <  zsys = totalMZ / totalMass;
659 <
660 <  return zsys;
661 < }
662 <
663 < /**
664 < *
665 < */
666 < template<typename T> void ZConstraint<T>::thermalize(void){
667 <  T::thermalize();
668 <  zeroOutVel();
669 < }
670 <
671 < /**
672 < *
673 < */
674 <
675 < template<typename T> void ZConstraint<T>::zeroOutVel(){
676 <  Atom** fixedZAtoms;  
677 <  double COMvel[3];
678 <  double vel[3];
679 <  double zSysCOMVel;
680 <
681 <  //zero out the velocities of center of mass of fixed z-constrained molecules
682 <
683 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
684 <    if (states[i] == zcsFixed){
685 <      zconsMols[i]->getCOMvel(COMvel);      
686 <      //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
687 <
688 <      fixedZAtoms = zconsMols[i]->getMyAtoms();
689 <
690 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
691 <        fixedZAtoms[j]->getVel(vel);
692 <        vel[whichDirection] -= COMvel[whichDirection];
693 <        fixedZAtoms[j]->setVel(vel);
694 <      }
695 <
696 <      zconsMols[i]->getCOMvel(COMvel);
697 <      //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
698 <    }
699 <  }
700 <
701 <  //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
702 <
703 <  zSysCOMVel = calcSysCOMVel();
704 < #ifdef IS_MPI
705 <  if (worldRank == 0){
706 < #endif
707 <    //cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
708 < #ifdef IS_MPI
709 <  }
710 < #endif
711 <
712 <  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
713 <  double MVzOfMovingMols_local;
714 <  double MVzOfMovingMols;
715 <  double totalMassOfMovingZMols_local;
716 <  double totalMassOfMovingZMols;
717 <
718 <  MVzOfMovingMols_local = 0;
719 <  totalMassOfMovingZMols_local = 0;
720 <
721 <  for (int i = 0; i < (int) (unconsMols.size()); i++){
722 <    unconsMols[i]->getCOMvel(COMvel);
723 <    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
724 <  }
725 <
726 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
727 <    if (states[i] == zcsMoving){
728 <      zconsMols[i]->getCOMvel(COMvel);
729 <      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
730 <      totalMassOfMovingZMols_local += massOfZConsMols[i];
731 <    }
732 <  }
733 <
734 < #ifndef IS_MPI
735 <  MVzOfMovingMols = MVzOfMovingMols_local;
736 <  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
737 < #else
738 <  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
739 <                MPI_SUM, MPI_COMM_WORLD);
740 <  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
741 <                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
742 < #endif
743 <
744 <  double vzOfMovingMols;
745 <  vzOfMovingMols = MVzOfMovingMols /
746 <                   (totalMassOfUncons + totalMassOfMovingZMols);
747 <
748 <  //modify the velocites of unconstrained molecules  
749 <  Atom** unconsAtoms;
750 <  for (int i = 0; i < (int) (unconsMols.size()); i++){
751 <    unconsAtoms = unconsMols[i]->getMyAtoms();
752 <    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
753 <      unconsAtoms[j]->getVel(vel);
754 <      vel[whichDirection] -= vzOfMovingMols;
755 <      unconsAtoms[j]->setVel(vel);
756 <    }
757 <  }  
758 <
759 <  //modify the velocities of moving z-constrained molecuels
760 <  Atom** movingZAtoms;
761 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
762 <    if (states[i] == zcsMoving){
763 <      movingZAtoms = zconsMols[i]->getMyAtoms();
764 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
765 <        movingZAtoms[j]->getVel(vel);
766 <        vel[whichDirection] -= vzOfMovingMols;
767 <        movingZAtoms[j]->setVel(vel);
768 <      }
769 <    }
770 <  }
771 <
772 <
773 <  zSysCOMVel = calcSysCOMVel();
774 < #ifdef IS_MPI
775 <  if (worldRank == 0){
776 < #endif
777 <    //cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
778 < #ifdef IS_MPI
779 <  }
780 < #endif
781 < }
782 <
783 < /**
784 < *
785 < */
786 <
787 < template<typename T> void ZConstraint<T>::doZconstraintForce(){
788 <  Atom** zconsAtoms;
789 <  double totalFZ;
790 <  double totalFZ_local;
791 <  double COM[3];
792 <  double force[3];
793 <
794 <  //constrain the molecules which do not reach the specified positions  
795 <
796 <  //Zero Out the force of z-contrained molecules    
797 <  totalFZ_local = 0;
798 <
799 <  //calculate the total z-contrained force of fixed z-contrained molecules
800 <
801 <  //cout << "before zero out z-constraint force on fixed z-constraint molecuels "
802 <  //       << "total force is " << calcTotalForce() << endl;
803 <
804 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
805 <    if (states[i] == zcsFixed){
806 <      zconsMols[i]->getCOM(COM);
807 <      zconsAtoms = zconsMols[i]->getMyAtoms();  
808 <
809 <      fz[i] = 0;      
810 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
811 <        zconsAtoms[j]->getFrc(force);
812 <        fz[i] += force[whichDirection];
813 <      }
814 <      totalFZ_local += fz[i];
815 <
816 <      //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
817 <      //      <<"\tcurrent zpos: " << COM[whichDirection]
818 <      //      << "\tcurrent fz: " <<fz[i] << endl;
819 <    }
820 <  }
821 <
822 <  //calculate total z-constraint force
823 < #ifdef IS_MPI
824 <  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
825 < #else
826 <  totalFZ = totalFZ_local;
827 < #endif
828 <
829 <
830 <  // apply negative to fixed z-constrained molecues;
831 <  force[0] = 0;
832 <  force[1] = 0;
833 <  force[2] = 0;
834 <
835 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
836 <    if (states[i] == zcsFixed){
837 <      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
838 <      zconsAtoms = zconsMols[i]->getMyAtoms();  
839 <
840 <      for (int j = 0; j < nAtomOfCurZConsMol; j++){
841 <        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
842 <        force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i],
843 <                                                                zconsAtoms[j],
844 <                                                                fz[i]);
845 <        zconsAtoms[j]->addFrc(force);
846 <      }
847 <    }
848 <  }
849 <
850 <  //cout << "after zero out z-constraint force on fixed z-constraint molecuels "
851 <  //      << "total force is " << calcTotalForce() << endl;
852 <
853 <
854 <  force[0] = 0;
855 <  force[1] = 0;
856 <  force[2] = 0;
857 <
858 <  //modify the forces of unconstrained molecules
859 <  for (int i = 0; i < (int) (unconsMols.size()); i++){
860 <    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
861 <
862 <    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
863 <      //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
864 <      force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],
865 <                                                             totalFZ);
866 <      unconsAtoms[j]->addFrc(force);
867 <    }
868 <  }      
869 <
870 <  //modify the forces of moving z-constrained molecules
871 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
872 <    if (states[i] == zcsMoving){
873 <      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
874 <
875 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
876 <        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
877 <        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],
878 <                                                               totalFZ);
879 <        movingZAtoms[j]->addFrc(force);
880 <      }
881 <    }
882 <  }
883 <  //  cout << "after substracting z-constraint force from moving molecuels "
884 <  //        << "total force is " << calcTotalForce()  << endl;
885 < }
886 <
887 < /**
888 <  *
889 <  *
890 <  */
891 <
892 < template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
893 <  double force[3];
894 <  double harmonicU;
895 <  double harmonicF;
896 <  double COM[3];
897 <  double diff;
898 <  double totalFZ_local;
899 <  double totalFZ;
900 <
901 <  force[0] = 0;
902 <  force[1] = 0;
903 <  force[2] = 0;
904 <
905 <  totalFZ_local = 0;
906 <
907 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
908 <    if (states[i] == zcsMoving){
909 <      zconsMols[i]->getCOM(COM);
910 <      //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
911 <      //     << "\tcurrent zpos: " << COM[whichDirection] << endl;
912 <
913 <      diff = COM[whichDirection] - resPos[i];
914 <
915 <      harmonicU = 0.5 * kz[i] * diff * diff;  
916 <      info->lrPot += harmonicU;
917 <
918 <      harmonicF = -kz[i] * diff;
919 <      totalFZ_local += harmonicF;
920 <
921 <      //adjust force
922 <
923 <      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
924 <
925 <      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
926 <        //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
927 <        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i],
928 <                                                               movingZAtoms[j],
929 <                                                               harmonicF);
930 <        movingZAtoms[j]->addFrc(force);
931 <      }
932 <    }
933 <  }
934 <
935 < #ifndef IS_MPI
936 <  totalFZ = totalFZ_local;
937 < #else
938 <  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
939 < #endif
940 <
941 <  //cout << "before substracting harmonic force from moving molecuels "
942 <  //      << "total force is " << calcTotalForce()  << endl;
943 <
944 <  force[0] = 0;
945 <  force[1] = 0;
946 <  force[2] = 0;
947 <
948 <  //modify the forces of unconstrained molecules
949 <  for (int i = 0; i < (int) (unconsMols.size()); i++){
950 <    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
951 <
952 <    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
953 <      //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
954 <      force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j],
955 <                                                              totalFZ);
956 <      unconsAtoms[j]->addFrc(force);
957 <    }
958 <  }  
959 <
960 <  //cout << "after substracting harmonic force from moving molecuels "
961 <  //      << "total force is " << calcTotalForce()  << endl;
962 < }
963 <
964 < /**
965 < *
966 < */
967 <
968 < template<typename T> bool ZConstraint<T>::checkZConsState(){
969 <  double COM[3];
970 <  double diff;
971 <
972 <  int changed_local;
973 <  int changed;
974 <
975 <  changed_local = 0;
976 <
977 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
978 <    zconsMols[i]->getCOM(COM);
979 <    diff = fabs(COM[whichDirection] - zPos[i]);  
980 <    if (diff <= zconsTol && states[i] == zcsMoving){
981 <      states[i] = zcsFixed;
982 <      changed_local = 1;
983 <
984 <      if(usingSMD)
985 <        prevCantPos = cantPos;
986 <
987 <      if (hasZConsGap)
988 <        endFixTime[i] = info->getTime() + zconsFixTime;
989 <    }
990 <    else if (diff > zconsTol && states[i] == zcsFixed){
991 <      states[i] = zcsMoving;
992 <      changed_local = 1;  
993 <
994 <      if(usingSMD)
995 <         cantPos = prevCantPos;
996 <      
997 <      if (hasZConsGap)
998 <        endFixTime[i] = INFINITE_TIME;
999 <    }
1000 <  }
1001 <
1002 < #ifndef IS_MPI
1003 <  changed = changed_local;
1004 < #else
1005 <  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1006 < #endif
1007 <
1008 <  return (changed > 0);
1009 < }
1010 <
1011 < template<typename T> bool ZConstraint<T>::haveFixedZMols(){
1012 <  int havingFixed_local;
1013 <  int havingFixed;
1014 <
1015 <  havingFixed_local = 0;
1016 <
1017 <  for (int i = 0; i < (int) (zconsMols.size()); i++)
1018 <    if (states[i] == zcsFixed){
1019 <      havingFixed_local = 1;
1020 <      break;
1021 <    }
1022 <
1023 < #ifndef IS_MPI
1024 <  havingFixed = havingFixed_local;
1025 < #else
1026 <  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
1027 <                MPI_COMM_WORLD);
1028 < #endif
1029 <
1030 <  return (havingFixed > 0);
1031 < }
1032 <
1033 <
1034 < /**
1035 < *
1036 < */
1037 < template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1038 <  int havingMoving_local;
1039 <  int havingMoving;
1040 <
1041 <  havingMoving_local = 0;
1042 <
1043 <  for (int i = 0; i < (int) (zconsMols.size()); i++)
1044 <    if (states[i] == zcsMoving){
1045 <      havingMoving_local = 1;
1046 <      break;
1047 <    }
1048 <
1049 < #ifndef IS_MPI
1050 <  havingMoving = havingMoving_local;
1051 < #else
1052 <  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
1053 <                MPI_COMM_WORLD);
1054 < #endif
1055 <
1056 <  return (havingMoving > 0);
1057 < }
1058 <
1059 < /**
1060 < *
1061 < */
1062 <
1063 < template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){
1064 <  double MVzOfMovingMols_local;
1065 <  double MVzOfMovingMols;
1066 <  double totalMassOfMovingZMols_local;
1067 <  double totalMassOfMovingZMols;
1068 <  double COMvel[3];
1069 <
1070 <  MVzOfMovingMols_local = 0;
1071 <  totalMassOfMovingZMols_local = 0;
1072 <
1073 <  for (int i = 0; i < unconsMols.size(); i++){
1074 <    unconsMols[i]->getCOMvel(COMvel);
1075 <    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1076 <  }
1077 <
1078 <  for (int i = 0; i < zconsMols.size(); i++){
1079 <    if (states[i] == zcsMoving){
1080 <      zconsMols[i]->getCOMvel(COMvel);
1081 <      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
1082 <      totalMassOfMovingZMols_local += massOfZConsMols[i];
1083 <    }
1084 <  }
1085 <
1086 < #ifndef IS_MPI
1087 <  MVzOfMovingMols = MVzOfMovingMols_local;
1088 <  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1089 < #else
1090 <  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
1091 <                MPI_SUM, MPI_COMM_WORLD);
1092 <  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
1093 <                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
1094 < #endif
1095 <
1096 <  double vzOfMovingMols;
1097 <  vzOfMovingMols = MVzOfMovingMols /
1098 <                   (totalMassOfUncons + totalMassOfMovingZMols);
1099 <
1100 <  return vzOfMovingMols;
1101 < }
1102 <
1103 < /**
1104 < *
1105 < */
1106 <
1107 < template<typename T> double ZConstraint<T>::calcSysCOMVel(){
1108 <  double COMvel[3];
1109 <  double tempMVz_local;
1110 <  double tempMVz;
1111 <  double massOfZCons_local;
1112 <  double massOfZCons;
1113 <
1114 <
1115 <  tempMVz_local = 0;
1116 <
1117 <  for (int i = 0 ; i < nMols; i++){
1118 <    molecules[i].getCOMvel(COMvel);
1119 <    tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection];
1120 <  }
1121 <
1122 <  massOfZCons_local = 0;
1123 <
1124 <  for (int i = 0; i < (int) (massOfZConsMols.size()); i++){
1125 <    massOfZCons_local += massOfZConsMols[i];
1126 <  }
1127 < #ifndef IS_MPI
1128 <  massOfZCons = massOfZCons_local;
1129 <  tempMVz = tempMVz_local;
1130 < #else
1131 <  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM,
1132 <                MPI_COMM_WORLD);
1133 <  MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
1134 < #endif
1135 <
1136 <  return tempMVz / (totalMassOfUncons + massOfZCons);
1137 < }
1138 <
1139 < /**
1140 < *
1141 < */
1142 <
1143 < template<typename T> double ZConstraint<T>::calcTotalForce(){
1144 <  double force[3];  
1145 <  double totalForce_local;
1146 <  double totalForce;
1147 <
1148 <  totalForce_local = 0;
1149 <
1150 <  for (int i = 0; i < nAtoms; i++){
1151 <    atoms[i]->getFrc(force);
1152 <    totalForce_local += force[whichDirection];
1153 <  }
1154 <
1155 < #ifndef IS_MPI
1156 <  totalForce = totalForce_local;
1157 < #else
1158 <  MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM,
1159 <                MPI_COMM_WORLD);
1160 < #endif
1161 <
1162 <  return totalForce;
1163 < }
1164 <
1165 < /**
1166 < *
1167 < */
1168 <
1169 < template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1170 <  //calculate the number of atoms of moving z-constrained molecules
1171 <  int nMovingZAtoms_local;
1172 <  int nMovingZAtoms;
1173 <
1174 <  nMovingZAtoms_local = 0;
1175 <  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1176 <    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1177 <      nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1178 <    }
1179 <
1180 < #ifdef IS_MPI
1181 <  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM,
1182 <                MPI_COMM_WORLD);
1183 < #else
1184 <  nMovingZAtoms = nMovingZAtoms_local;
1185 < #endif
1186 <  totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1187 < }
1188 <
1189 < template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol,
1190 <                                                                              Atom* atom,
1191 <                                                                              double totalForce){
1192 <  return totalForce / mol->getNAtoms();
1193 < }
1194 <
1195 < template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom,
1196 <                                                                              double totalForce){
1197 <  return totalForce / totNumOfMovingAtoms;
1198 < }
1199 <
1200 < template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol,
1201 <                                                                              Atom* atom,
1202 <                                                                              double totalForce){
1203 <  return totalForce / mol->getNAtoms();
1204 < }
1205 <
1206 < template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom,
1207 <                                                                              double totalForce){
1208 <  return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1209 < }
1210 <
1211 < /**
1212 < *
1213 < */
1214 <
1215 < template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1216 <  //calculate the number of atoms of moving z-constrained molecules
1217 <  double massOfMovingZAtoms_local;
1218 <  double massOfMovingZAtoms;
1219 <
1220 <  massOfMovingZAtoms_local = 0;
1221 <  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1222 <    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1223 <      massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1224 <    }
1225 <
1226 < #ifdef IS_MPI
1227 <  MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,
1228 <                MPI_SUM, MPI_COMM_WORLD);
1229 < #else
1230 <  massOfMovingZAtoms = massOfMovingZAtoms_local;
1231 < #endif
1232 <  totMassOfMovingAtoms = massOfMovingZAtoms +
1233 <                         zconsIntegrator->totalMassOfUncons;
1234 < }
1235 <
1236 < template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol,
1237 <                                                                            Atom* atom,
1238 <                                                                            double totalForce){
1239 <  return totalForce * atom->getMass() / mol->getTotalMass();
1240 < }
1241 <
1242 < template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom,
1243 <                                                                            double totalForce){
1244 <  return totalForce * atom->getMass() / totMassOfMovingAtoms;
1245 < }
1246 <
1247 < template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol,
1248 <                                                                            Atom* atom,
1249 <                                                                            double totalForce){
1250 <  return totalForce * atom->getMass() / mol->getTotalMass();
1251 < }
1252 <
1253 < template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom,
1254 <                                                                            double totalForce){
1255 <  return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1256 < }
1257 <
1258 < template<typename T> void ZConstraint<T>::updateZPos(){
1259 <  double curTime;
1260 <  double COM[3];
1261 <  
1262 <  curTime = info->getTime();
1263 <
1264 <  for (size_t i = 0; i < zconsMols.size(); i++){
1265 <
1266 <    if (states[i] == zcsFixed && curTime >= endFixTime[i]){
1267 <      zPos[i] += zconsGap;
1268 <
1269 <      if (usingSMD){
1270 <        zconsMols[i]->getCOM(COM);
1271 <        cantPos[i] = COM[whichDirection];
1272 <      }
1273 <      
1274 <    }
1275 <    
1276 <  }
1277 <  
1278 < }
1279 <
1280 < template<typename T> void ZConstraint<T>::updateCantPos(){
1281 <  double curTime;
1282 <  double dt;
1283 <
1284 <  curTime = info->getTime();
1285 <  dt = info->dt;
1286 <
1287 <  for (size_t i = 0; i < zconsMols.size(); i++){
1288 <    if (states[i] == zcsMoving){
1289 <      cantPos[i] += cantVel[i] * dt;
1290 <    }
1291 <  }
1292 <
1293 < }
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->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 >  //
480 >  forcePolicy->update();
481 > }
482 >
483 > #endif
484 >
485 > /**
486 > *  Function Name: isZConstraintMol
487 > *  Parameter
488 > *    Molecule* mol
489 > *  Return value:
490 > *    -1, if the molecule is not z-constraint molecule,
491 > *    other non-negative values, its index in indexOfAllZConsMols vector
492 > */
493 >
494 > template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){
495 >  int index;
496 >  int low;
497 >  int high;
498 >  int mid;
499 >
500 >  index = mol->getGlobalIndex();
501 >
502 >  low = 0;
503 >  high = parameters->size() - 1;
504 >
505 >  //Binary Search (we have sorted the array)  
506 >  while (low <= high){
507 >    mid = (low + high) / 2;
508 >    if ((*parameters)[mid].zconsIndex == index)
509 >      return mid;
510 >    else if ((*parameters)[mid].zconsIndex > index)
511 >      high = mid - 1;
512 >    else
513 >      low = mid + 1;
514 >  }
515 >
516 >  return -1;
517 > }
518 >
519 > template<typename T> void ZConstraint<T>::integrate(){
520 >  // creat zconsWriter  
521 >  fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);  
522 >
523 >  if (!fzOut){
524 >    sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
525 >    painCave.isFatal = 1;
526 >    simError();
527 >  }
528 >
529 >  //zero out the velocities of center of mass of unconstrained molecules
530 >  //and the velocities of center of mass of every single z-constrained molecueles
531 >  zeroOutVel();
532 >
533 >  curZconsTime = zconsTime + info->getTime();
534 >
535 >  T::integrate();
536 > }
537 >
538 >
539 > /**
540 > *
541 > *
542 > *
543 > *
544 > */
545 > template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
546 >  double zsys;
547 >  double COM[3];
548 >  double force[3];
549 >  double zSysCOMVel;
550 >
551 >  T::calcForce(calcPot, calcStress);
552 >
553 >
554 >  if (hasZConsGap){
555 >    updateZPos();
556 >  }
557 >
558 >  if (checkZConsState()){
559 >    zeroOutVel();    
560 >    forcePolicy->update();
561 >  }  
562 >
563 >  zsys = calcZSys();
564 >  zSysCOMVel = calcSysCOMVel();
565 > #ifdef IS_MPI
566 >  if (worldRank == 0){
567 > #endif
568 >    //cout << "---------------------------------------------------------------------" <<endl;
569 >    //cout << "current time: " << info->getTime() << endl;
570 >    //cout << "center of mass at z: " << zsys << endl;    
571 >    //cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
572 >
573 > #ifdef IS_MPI
574 >  }
575 > #endif
576 >
577 >  //do zconstraint force;
578 >  if (haveFixedZMols()){
579 >    this->doZconstraintForce();
580 >  }
581 >
582 >  //use external force to move the molecules to the specified positions
583 >  if (haveMovingZMols()){
584 >    if (usingSMD)
585 >      this->doHarmonic(cantPos);
586 >    else
587 >      this->doHarmonic(zPos);      
588 >  }
589 >
590 >  //write out forces and current positions of z-constraint molecules
591 >  if (info->getTime() >= curZconsTime){
592 >    for (int i = 0; i < (int) (zconsMols.size()); i++){
593 >      zconsMols[i]->getCOM(COM);
594 >      curZPos[i] = COM[whichDirection];
595 >
596 >      //if the z-constraint molecule is still moving, just record its force
597 >      if (states[i] == zcsMoving){
598 >        fz[i] = 0;
599 >        Atom** movingZAtoms;
600 >        movingZAtoms = zconsMols[i]->getMyAtoms();
601 >        for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
602 >          movingZAtoms[j]->getFrc(force);
603 >          fz[i] += force[whichDirection];
604 >        }
605 >      }
606 >    }
607 >    fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
608 >                   &curZPos[0], &zPos[0]);
609 >    curZconsTime += zconsTime;
610 >  }
611 >
612 >  zSysCOMVel = calcSysCOMVel();  
613 > #ifdef IS_MPI
614 >  if (worldRank == 0){
615 > #endif
616 >    //cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
617 > #ifdef IS_MPI
618 >  }
619 > #endif
620 > }
621 >
622 >
623 > /**
624 > *
625 > */
626 >
627 > template<typename T> double ZConstraint<T>::calcZSys(){
628 >  //calculate reference z coordinate for z-constraint molecules
629 >  double totalMass_local;
630 >  double totalMass;
631 >  double totalMZ_local;
632 >  double totalMZ;
633 >  double massOfCurMol;
634 >  double COM[3];
635 >
636 >  totalMass_local = 0;
637 >  totalMZ_local = 0;
638 >
639 >  for (int i = 0; i < nMols; i++){
640 >    massOfCurMol = molecules[i].getTotalMass();
641 >    molecules[i].getCOM(COM);
642 >
643 >    totalMass_local += massOfCurMol;
644 >    totalMZ_local += massOfCurMol * COM[whichDirection];
645 >  }
646 >
647 >
648 > #ifdef IS_MPI  
649 >  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM,
650 >                MPI_COMM_WORLD);
651 >  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
652 > #else
653 >  totalMass = totalMass_local;
654 >  totalMZ = totalMZ_local;
655 > #endif  
656 >
657 >  double zsys;
658 >  zsys = totalMZ / totalMass;
659 >
660 >  return zsys;
661 > }
662 >
663 > /**
664 > *
665 > */
666 > template<typename T> void ZConstraint<T>::thermalize(void){
667 >  T::thermalize();
668 >  zeroOutVel();
669 > }
670 >
671 > /**
672 > *
673 > */
674 >
675 > template<typename T> void ZConstraint<T>::zeroOutVel(){
676 >  Atom** fixedZAtoms;  
677 >  double COMvel[3];
678 >  double vel[3];
679 >  double zSysCOMVel;
680 >
681 >  //zero out the velocities of center of mass of fixed z-constrained molecules
682 >
683 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
684 >    if (states[i] == zcsFixed){
685 >      zconsMols[i]->getCOMvel(COMvel);      
686 >      //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
687 >
688 >      fixedZAtoms = zconsMols[i]->getMyAtoms();
689 >
690 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
691 >        fixedZAtoms[j]->getVel(vel);
692 >        vel[whichDirection] -= COMvel[whichDirection];
693 >        fixedZAtoms[j]->setVel(vel);
694 >      }
695 >
696 >      zconsMols[i]->getCOMvel(COMvel);
697 >      //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
698 >    }
699 >  }
700 >
701 >  //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
702 >
703 >  zSysCOMVel = calcSysCOMVel();
704 > #ifdef IS_MPI
705 >  if (worldRank == 0){
706 > #endif
707 >    //cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
708 > #ifdef IS_MPI
709 >  }
710 > #endif
711 >
712 >  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
713 >  double MVzOfMovingMols_local;
714 >  double MVzOfMovingMols;
715 >  double totalMassOfMovingZMols_local;
716 >  double totalMassOfMovingZMols;
717 >
718 >  MVzOfMovingMols_local = 0;
719 >  totalMassOfMovingZMols_local = 0;
720 >
721 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
722 >    unconsMols[i]->getCOMvel(COMvel);
723 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
724 >  }
725 >
726 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
727 >    if (states[i] == zcsMoving){
728 >      zconsMols[i]->getCOMvel(COMvel);
729 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
730 >      totalMassOfMovingZMols_local += massOfZConsMols[i];
731 >    }
732 >  }
733 >
734 > #ifndef IS_MPI
735 >  MVzOfMovingMols = MVzOfMovingMols_local;
736 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
737 > #else
738 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
739 >                MPI_SUM, MPI_COMM_WORLD);
740 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
741 >                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
742 > #endif
743 >
744 >  double vzOfMovingMols;
745 >  vzOfMovingMols = MVzOfMovingMols /
746 >                   (totalMassOfUncons + totalMassOfMovingZMols);
747 >
748 >  //modify the velocites of unconstrained molecules  
749 >  Atom** unconsAtoms;
750 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
751 >    unconsAtoms = unconsMols[i]->getMyAtoms();
752 >    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
753 >      unconsAtoms[j]->getVel(vel);
754 >      vel[whichDirection] -= vzOfMovingMols;
755 >      unconsAtoms[j]->setVel(vel);
756 >    }
757 >  }  
758 >
759 >  //modify the velocities of moving z-constrained molecuels
760 >  Atom** movingZAtoms;
761 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
762 >    if (states[i] == zcsMoving){
763 >      movingZAtoms = zconsMols[i]->getMyAtoms();
764 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
765 >        movingZAtoms[j]->getVel(vel);
766 >        vel[whichDirection] -= vzOfMovingMols;
767 >        movingZAtoms[j]->setVel(vel);
768 >      }
769 >    }
770 >  }
771 >
772 >
773 >  zSysCOMVel = calcSysCOMVel();
774 > #ifdef IS_MPI
775 >  if (worldRank == 0){
776 > #endif
777 >    //cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
778 > #ifdef IS_MPI
779 >  }
780 > #endif
781 > }
782 >
783 > /**
784 > *
785 > */
786 >
787 > template<typename T> void ZConstraint<T>::doZconstraintForce(){
788 >  Atom** zconsAtoms;
789 >  double totalFZ;
790 >  double totalFZ_local;
791 >  double COM[3];
792 >  double force[3];
793 >
794 >  //constrain the molecules which do not reach the specified positions  
795 >
796 >  //Zero Out the force of z-contrained molecules    
797 >  totalFZ_local = 0;
798 >
799 >  //calculate the total z-contrained force of fixed z-contrained molecules
800 >
801 >  //cout << "before zero out z-constraint force on fixed z-constraint molecuels "
802 >  //       << "total force is " << calcTotalForce() << endl;
803 >
804 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
805 >    if (states[i] == zcsFixed){
806 >      zconsMols[i]->getCOM(COM);
807 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
808 >
809 >      fz[i] = 0;      
810 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
811 >        zconsAtoms[j]->getFrc(force);
812 >        fz[i] += force[whichDirection];
813 >      }
814 >      totalFZ_local += fz[i];
815 >
816 >      //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
817 >      //      <<"\tcurrent zpos: " << COM[whichDirection]
818 >      //      << "\tcurrent fz: " <<fz[i] << endl;
819 >    }
820 >  }
821 >
822 >  //calculate total z-constraint force
823 > #ifdef IS_MPI
824 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
825 > #else
826 >  totalFZ = totalFZ_local;
827 > #endif
828 >
829 >
830 >  // apply negative to fixed z-constrained molecues;
831 >  force[0] = 0;
832 >  force[1] = 0;
833 >  force[2] = 0;
834 >
835 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
836 >    if (states[i] == zcsFixed){
837 >      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
838 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
839 >
840 >      for (int j = 0; j < nAtomOfCurZConsMol; j++){
841 >        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
842 >        force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i],
843 >                                                                zconsAtoms[j],
844 >                                                                fz[i]);
845 >        zconsAtoms[j]->addFrc(force);
846 >      }
847 >    }
848 >  }
849 >
850 >  //cout << "after zero out z-constraint force on fixed z-constraint molecuels "
851 >  //      << "total force is " << calcTotalForce() << endl;
852 >
853 >
854 >  force[0] = 0;
855 >  force[1] = 0;
856 >  force[2] = 0;
857 >
858 >  //modify the forces of unconstrained molecules
859 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
860 >    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
861 >
862 >    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
863 >      //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
864 >      force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],
865 >                                                             totalFZ);
866 >      unconsAtoms[j]->addFrc(force);
867 >    }
868 >  }      
869 >
870 >  //modify the forces of moving z-constrained molecules
871 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
872 >    if (states[i] == zcsMoving){
873 >      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
874 >
875 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
876 >        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
877 >        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],
878 >                                                               totalFZ);
879 >        movingZAtoms[j]->addFrc(force);
880 >      }
881 >    }
882 >  }
883 >  //  cout << "after substracting z-constraint force from moving molecuels "
884 >  //        << "total force is " << calcTotalForce()  << endl;
885 > }
886 >
887 > /**
888 >  *
889 >  *
890 >  */
891 >
892 > template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
893 >  double force[3];
894 >  double harmonicU;
895 >  double harmonicF;
896 >  double COM[3];
897 >  double diff;
898 >  double totalFZ_local;
899 >  double totalFZ;
900 >
901 >  force[0] = 0;
902 >  force[1] = 0;
903 >  force[2] = 0;
904 >
905 >  totalFZ_local = 0;
906 >
907 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
908 >    if (states[i] == zcsMoving){
909 >      zconsMols[i]->getCOM(COM);
910 >      //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
911 >      //     << "\tcurrent zpos: " << COM[whichDirection] << endl;
912 >
913 >      diff = COM[whichDirection] - resPos[i];
914 >
915 >      harmonicU = 0.5 * kz[i] * diff * diff;  
916 >      info->lrPot += harmonicU;
917 >
918 >      harmonicF = -kz[i] * diff;
919 >      totalFZ_local += harmonicF;
920 >
921 >      //adjust force
922 >
923 >      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
924 >
925 >      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
926 >        //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
927 >        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i],
928 >                                                               movingZAtoms[j],
929 >                                                               harmonicF);
930 >        movingZAtoms[j]->addFrc(force);
931 >      }
932 >    }
933 >  }
934 >
935 > #ifndef IS_MPI
936 >  totalFZ = totalFZ_local;
937 > #else
938 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
939 > #endif
940 >
941 >  //cout << "before substracting harmonic force from moving molecuels "
942 >  //      << "total force is " << calcTotalForce()  << endl;
943 >
944 >  force[0] = 0;
945 >  force[1] = 0;
946 >  force[2] = 0;
947 >
948 >  //modify the forces of unconstrained molecules
949 >  for (int i = 0; i < (int) (unconsMols.size()); i++){
950 >    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
951 >
952 >    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
953 >      //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
954 >      force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j],
955 >                                                              totalFZ);
956 >      unconsAtoms[j]->addFrc(force);
957 >    }
958 >  }  
959 >
960 >  //cout << "after substracting harmonic force from moving molecuels "
961 >  //      << "total force is " << calcTotalForce()  << endl;
962 > }
963 >
964 > /**
965 > *
966 > */
967 >
968 > template<typename T> bool ZConstraint<T>::checkZConsState(){
969 >  double COM[3];
970 >  double diff;
971 >
972 >  int changed_local;
973 >  int changed;
974 >
975 >  changed_local = 0;
976 >
977 >  for (int i = 0; i < (int) (zconsMols.size()); i++){
978 >    zconsMols[i]->getCOM(COM);
979 >    diff = fabs(COM[whichDirection] - zPos[i]);  
980 >    if (diff <= zconsTol && states[i] == zcsMoving){
981 >      states[i] = zcsFixed;
982 >      changed_local = 1;
983 >
984 >      if(usingSMD)
985 >        prevCantPos = cantPos;
986 >
987 >      if (hasZConsGap)
988 >        endFixTime[i] = info->getTime() + zconsFixTime;
989 >    }
990 >    else if (diff > zconsTol && states[i] == zcsFixed){
991 >      states[i] = zcsMoving;
992 >      changed_local = 1;  
993 >
994 >      if(usingSMD)
995 >         cantPos = prevCantPos;
996 >      
997 >      if (hasZConsGap)
998 >        endFixTime[i] = INFINITE_TIME;
999 >    }
1000 >  }
1001 >
1002 > #ifndef IS_MPI
1003 >  changed = changed_local;
1004 > #else
1005 >  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1006 > #endif
1007 >
1008 >  return (changed > 0);
1009 > }
1010 >
1011 > template<typename T> bool ZConstraint<T>::haveFixedZMols(){
1012 >  int havingFixed_local;
1013 >  int havingFixed;
1014 >
1015 >  havingFixed_local = 0;
1016 >
1017 >  for (int i = 0; i < (int) (zconsMols.size()); i++)
1018 >    if (states[i] == zcsFixed){
1019 >      havingFixed_local = 1;
1020 >      break;
1021 >    }
1022 >
1023 > #ifndef IS_MPI
1024 >  havingFixed = havingFixed_local;
1025 > #else
1026 >  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
1027 >                MPI_COMM_WORLD);
1028 > #endif
1029 >
1030 >  return (havingFixed > 0);
1031 > }
1032 >
1033 >
1034 > /**
1035 > *
1036 > */
1037 > template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1038 >  int havingMoving_local;
1039 >  int havingMoving;
1040 >
1041 >  havingMoving_local = 0;
1042 >
1043 >  for (int i = 0; i < (int) (zconsMols.size()); i++)
1044 >    if (states[i] == zcsMoving){
1045 >      havingMoving_local = 1;
1046 >      break;
1047 >    }
1048 >
1049 > #ifndef IS_MPI
1050 >  havingMoving = havingMoving_local;
1051 > #else
1052 >  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
1053 >                MPI_COMM_WORLD);
1054 > #endif
1055 >
1056 >  return (havingMoving > 0);
1057 > }
1058 >
1059 > /**
1060 > *
1061 > */
1062 >
1063 > template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){
1064 >  double MVzOfMovingMols_local;
1065 >  double MVzOfMovingMols;
1066 >  double totalMassOfMovingZMols_local;
1067 >  double totalMassOfMovingZMols;
1068 >  double COMvel[3];
1069 >
1070 >  MVzOfMovingMols_local = 0;
1071 >  totalMassOfMovingZMols_local = 0;
1072 >
1073 >  for (int i = 0; i < unconsMols.size(); i++){
1074 >    unconsMols[i]->getCOMvel(COMvel);
1075 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1076 >  }
1077 >
1078 >  for (int i = 0; i < zconsMols.size(); i++){
1079 >    if (states[i] == zcsMoving){
1080 >      zconsMols[i]->getCOMvel(COMvel);
1081 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
1082 >      totalMassOfMovingZMols_local += massOfZConsMols[i];
1083 >    }
1084 >  }
1085 >
1086 > #ifndef IS_MPI
1087 >  MVzOfMovingMols = MVzOfMovingMols_local;
1088 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1089 > #else
1090 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
1091 >                MPI_SUM, MPI_COMM_WORLD);
1092 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
1093 >                MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
1094 > #endif
1095 >
1096 >  double vzOfMovingMols;
1097 >  vzOfMovingMols = MVzOfMovingMols /
1098 >                   (totalMassOfUncons + totalMassOfMovingZMols);
1099 >
1100 >  return vzOfMovingMols;
1101 > }
1102 >
1103 > /**
1104 > *
1105 > */
1106 >
1107 > template<typename T> double ZConstraint<T>::calcSysCOMVel(){
1108 >  double COMvel[3];
1109 >  double tempMVz_local;
1110 >  double tempMVz;
1111 >  double massOfZCons_local;
1112 >  double massOfZCons;
1113 >
1114 >
1115 >  tempMVz_local = 0;
1116 >
1117 >  for (int i = 0 ; i < nMols; i++){
1118 >    molecules[i].getCOMvel(COMvel);
1119 >    tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection];
1120 >  }
1121 >
1122 >  massOfZCons_local = 0;
1123 >
1124 >  for (int i = 0; i < (int) (massOfZConsMols.size()); i++){
1125 >    massOfZCons_local += massOfZConsMols[i];
1126 >  }
1127 > #ifndef IS_MPI
1128 >  massOfZCons = massOfZCons_local;
1129 >  tempMVz = tempMVz_local;
1130 > #else
1131 >  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM,
1132 >                MPI_COMM_WORLD);
1133 >  MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
1134 > #endif
1135 >
1136 >  return tempMVz / (totalMassOfUncons + massOfZCons);
1137 > }
1138 >
1139 > /**
1140 > *
1141 > */
1142 >
1143 > template<typename T> double ZConstraint<T>::calcTotalForce(){
1144 >  double force[3];  
1145 >  double totalForce_local;
1146 >  double totalForce;
1147 >
1148 >  totalForce_local = 0;
1149 >
1150 >  for (int i = 0; i < nAtoms; i++){
1151 >    atoms[i]->getFrc(force);
1152 >    totalForce_local += force[whichDirection];
1153 >  }
1154 >
1155 > #ifndef IS_MPI
1156 >  totalForce = totalForce_local;
1157 > #else
1158 >  MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM,
1159 >                MPI_COMM_WORLD);
1160 > #endif
1161 >
1162 >  return totalForce;
1163 > }
1164 >
1165 > /**
1166 > *
1167 > */
1168 >
1169 > template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1170 >  //calculate the number of atoms of moving z-constrained molecules
1171 >  int nMovingZAtoms_local;
1172 >  int nMovingZAtoms;
1173 >
1174 >  nMovingZAtoms_local = 0;
1175 >  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1176 >    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1177 >      nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1178 >    }
1179 >
1180 > #ifdef IS_MPI
1181 >  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM,
1182 >                MPI_COMM_WORLD);
1183 > #else
1184 >  nMovingZAtoms = nMovingZAtoms_local;
1185 > #endif
1186 >  totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1187 > }
1188 >
1189 > template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol,
1190 >                                                                              Atom* atom,
1191 >                                                                              double totalForce){
1192 >  return totalForce / mol->getNAtoms();
1193 > }
1194 >
1195 > template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom,
1196 >                                                                              double totalForce){
1197 >  return totalForce / totNumOfMovingAtoms;
1198 > }
1199 >
1200 > template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol,
1201 >                                                                              Atom* atom,
1202 >                                                                              double totalForce){
1203 >  return totalForce / mol->getNAtoms();
1204 > }
1205 >
1206 > template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom,
1207 >                                                                              double totalForce){
1208 >  return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1209 > }
1210 >
1211 > /**
1212 > *
1213 > */
1214 >
1215 > template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1216 >  //calculate the number of atoms of moving z-constrained molecules
1217 >  double massOfMovingZAtoms_local;
1218 >  double massOfMovingZAtoms;
1219 >
1220 >  massOfMovingZAtoms_local = 0;
1221 >  for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1222 >    if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1223 >      massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1224 >    }
1225 >
1226 > #ifdef IS_MPI
1227 >  MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,
1228 >                MPI_SUM, MPI_COMM_WORLD);
1229 > #else
1230 >  massOfMovingZAtoms = massOfMovingZAtoms_local;
1231 > #endif
1232 >  totMassOfMovingAtoms = massOfMovingZAtoms +
1233 >                         zconsIntegrator->totalMassOfUncons;
1234 > }
1235 >
1236 > template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol,
1237 >                                                                            Atom* atom,
1238 >                                                                            double totalForce){
1239 >  return totalForce * atom->getMass() / mol->getTotalMass();
1240 > }
1241 >
1242 > template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom,
1243 >                                                                            double totalForce){
1244 >  return totalForce * atom->getMass() / totMassOfMovingAtoms;
1245 > }
1246 >
1247 > template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol,
1248 >                                                                            Atom* atom,
1249 >                                                                            double totalForce){
1250 >  return totalForce * atom->getMass() / mol->getTotalMass();
1251 > }
1252 >
1253 > template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom,
1254 >                                                                            double totalForce){
1255 >  return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1256 > }
1257 >
1258 > template<typename T> void ZConstraint<T>::updateZPos(){
1259 >  double curTime;
1260 >  double COM[3];
1261 >  
1262 >  curTime = info->getTime();
1263 >
1264 >  for (size_t i = 0; i < zconsMols.size(); i++){
1265 >
1266 >    if (states[i] == zcsFixed && curTime >= endFixTime[i]){
1267 >      zPos[i] += zconsGap;
1268 >
1269 >      if (usingSMD){
1270 >        zconsMols[i]->getCOM(COM);
1271 >        cantPos[i] = COM[whichDirection];
1272 >      }
1273 >      
1274 >    }
1275 >    
1276 >  }
1277 >  
1278 > }
1279 >
1280 > template<typename T> void ZConstraint<T>::updateCantPos(){
1281 >  double curTime;
1282 >  double dt;
1283 >
1284 >  curTime = info->getTime();
1285 >  dt = info->dt;
1286 >
1287 >  for (size_t i = 0; i < zconsMols.size(); i++){
1288 >    if (states[i] == zcsMoving){
1289 >      cantPos[i] += cantVel[i] * dt;
1290 >    }
1291 >  }
1292 >
1293 > }

Diff Legend

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