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 660 by tim, Thu Jul 31 19:59:34 2003 UTC vs.
Revision 682 by tim, Tue Aug 12 17:51:33 2003 UTC

# Line 1 | Line 1
1   #include "Integrator.hpp"
2   #include "simError.h"
3 <
3 > #include <cmath>
4   template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff)
5 <                                    : T(theInfo, the_ff), fz(NULL), indexOfZConsMols(NULL)
5 >                                    : T(theInfo, the_ff), fz(NULL),
6 >                                      indexOfZConsMols(NULL)
7   {
8  
9    //get properties from SimInfo
10    GenericData* data;
11 <  IndexData* index;
11 >  ZConsParaData* zConsParaData;
12    DoubleData* sampleTime;
13 +  DoubleData* tolerance;
14    StringData* filename;
15 <  
14 <  
15 <  data = info->getProperty("zconsindex");
16 <  if(!data) {
15 >  double COM[3];
16  
17 <    sprintf( painCave.errMsg,
18 <               "ZConstraint error: If you use an ZConstraint\n"
19 <               " , you must set index of z-constraint molecules.\n");
20 <    painCave.isFatal = 1;
21 <    simError();  
23 <  }
24 <  else{
25 <    index = dynamic_cast<IndexData*>(data);
26 <    
27 <    if(!index){
17 >  //by default, the direction of constraint is z
18 >  // 0 --> x
19 >  // 1 --> y
20 >  // 2 --> z
21 >  whichDirection = 2;
22  
23 <      sprintf( painCave.errMsg,
24 <                 "ZConstraint error: Can not get property from SimInfo\n");
25 <      painCave.isFatal = 1;
26 <      simError();  
27 <    
34 <    }
35 <    else{
36 <          
37 <      indexOfAllZConsMols = index->getIndexData();
38 <      
39 <      //the maximum value of index is the last one(we sorted the index data in SimSetup.cpp)
40 <      int maxIndex;
41 <      int totalNumMol;
42 <      
43 <      maxIndex = indexOfAllZConsMols[indexOfAllZConsMols.size() - 1];
23 >  //estimate the force constant of harmonical potential
24 >  double Kb = 1.986E-3 ; //in kcal/K
25 >  
26 >  double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2;
27 >  zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox);
28  
29 < #ifndef IS_MPI
30 <      totalNumMol = nMols;
47 < #else
48 <      totalNumMol = mpiSim->getTotNmol();  
49 < #endif      
50 <      
51 <      if(maxIndex > totalNumMol - 1){
52 <        sprintf( painCave.errMsg,
53 <               "ZConstraint error: index is out of range\n");
54 <        painCave.isFatal = 1;
55 <        simError();
56 <                
57 <      }
58 <      
59 <    }
60 <        
61 <  }
29 >  //retrieve sample time of z-contraint
30 >  data = info->getProperty(ZCONSTIME_ID);
31    
63 –  //retrive sample time of z-contraint
64 –  data = info->getProperty("zconstime");
65 –  
32    if(!data) {
33        
34      sprintf( painCave.errMsg,
# Line 89 | Line 55 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
55  
56    }
57    
58 <  
59 <  //retrive output filename of z force
94 <  data = info->getProperty("zconsfilename");
58 >  //retrieve output filename of z force
59 >  data = info->getProperty(ZCONSFILENAME_ID);
60    if(!data) {
61  
62        
# Line 121 | Line 86 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
86  
87    }
88  
89 <
90 <  //calculate reference z coordinate for z-constraint molecules
126 <  double totalMass_local;
127 <  double totalMass;
128 <  double totalMZ_local;
129 <  double totalMZ;
130 <  double massOfUncons_local;
131 <  double massOfCurMol;
132 <  double COM[3];
89 >  //retrieve tolerance for z-constraint molecuels
90 >  data = info->getProperty(ZCONSTOL_ID);
91    
92 <  totalMass_local = 0;
93 <  totalMass = 0;
94 <  totalMZ_local = 0;
95 <  totalMZ = 0;
96 <  massOfUncons_local = 0;
97 <    
140 <  
141 <  for(int i = 0; i < nMols; i++){
142 <    massOfCurMol = molecules[i].getTotalMass();
143 <    molecules[i].getCOM(COM);
144 <    
145 <    totalMass_local += massOfCurMol;
146 <    totalMZ_local += massOfCurMol * COM[2];
147 <    
148 <    if(isZConstraintMol(&molecules[i]) == -1){
149 <    
150 <      massOfUncons_local += massOfCurMol;
151 <    }  
152 <    
92 >  if(!data) {
93 >      
94 >    sprintf( painCave.errMsg,
95 >               "ZConstraint error: can not get tolerance \n");
96 >    painCave.isFatal = 1;
97 >    simError();      
98    }
99 +  else{
100    
101 <  
102 < #ifdef IS_MPI  
103 <  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
158 <  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
159 <  MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
160 < #else
161 <  totalMass = totalMass_local;
162 <  totalMZ = totalMZ_local;
163 <  totalMassOfUncons = massOfUncons_local;
164 < #endif  
101 >    tolerance = dynamic_cast<DoubleData*>(data);
102 >    
103 >    if(!tolerance){
104  
105 <  double zsys;
106 <  zsys = totalMZ / totalMass;
107 <
108 < #ifndef IS_MPI  
109 <  for(int i = 0; i < nMols; i++){
171 <    
172 <    if(isZConstraintMol(&molecules[i]) > -1 ){
173 <      molecules[i].getCOM(COM);
174 <      allRefZ.push_back(COM[2] - zsys);  
105 >      sprintf( painCave.errMsg,
106 >                 "ZConstraint error: Can not get property from SimInfo\n");
107 >      painCave.isFatal = 1;
108 >      simError();  
109 >      
110      }
111 <    
111 >    else{
112 >      this->zconsTol = tolerance->getData();
113 >    }
114 >
115    }
116 < #else
117 <
118 <  int whichNode;
119 <  enum CommType { RequestMolZPos, EndOfRequest} status;
120 <  //int status;
121 <  double zpos;
122 <  int localIndex;
123 <  MPI_Status ierr;
124 <  int tag = 0;
125 <  
188 <  if(worldRank == 0){
189 <    
190 <    int globalIndexOfCurMol;
191 <    int *MolToProcMap;
192 <    MolToProcMap = mpiSim->getMolToProcMap();
193 <    
194 <    for(int i = 0; i < indexOfAllZConsMols.size(); i++){
195 <      
196 <      whichNode = MolToProcMap[indexOfAllZConsMols[i]];
197 <      globalIndexOfCurMol = indexOfAllZConsMols[i];
198 <      
199 <      if(whichNode == 0){
200 <        
201 <        for(int j = 0; j < nMols; j++)
202 <          if(molecules[j].getGlobalIndex() == globalIndexOfCurMol){
203 <            localIndex = j;
204 <            break;
205 <          }
206 <                  
207 <        molecules[localIndex].getCOM(COM);
208 <        allRefZ.push_back(COM[2] - zsys);  
209 <              
210 <      }
211 <      else{
212 <        status = RequestMolZPos;
213 <        MPI_Send(&status, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
214 <        MPI_Send(&globalIndexOfCurMol, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
215 <        MPI_Recv(&zpos, 1, MPI_DOUBLE_PRECISION, whichNode, tag, MPI_COMM_WORLD, &ierr);
216 <        
217 <        allRefZ.push_back(zpos - zsys);
218 <      
219 <      }
220 <              
221 <    } //End of Request Loop
222 <    
223 <    //Send ending request message to slave nodes    
224 <    status = EndOfRequest;
225 <    for(int i =1; i < mpiSim->getNumberProcessors(); i++)
226 <      MPI_Send(&status, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
227 <    
116 >        
117 >  //retrieve index of z-constraint molecules
118 >  data = info->getProperty(ZCONSPARADATA_ID);
119 >  if(!data) {
120 >
121 >    sprintf( painCave.errMsg,
122 >               "ZConstraint error: If you use an ZConstraint\n"
123 >               " , you must set index of z-constraint molecules.\n");
124 >    painCave.isFatal = 1;
125 >    simError();  
126    }
127    else{
128    
129 <    int whichMol;
130 <    bool done = false;
129 >    zConsParaData = dynamic_cast<ZConsParaData*>(data);
130 >    
131 >    if(!zConsParaData){
132  
133 <    while (!done){  
134 <      
135 <      MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
133 >      sprintf( painCave.errMsg,
134 >                 "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
135 >      painCave.isFatal = 1;
136 >      simError();  
137      
238 –      switch (status){
239 –          
240 –        case RequestMolZPos :
241 –          
242 –          MPI_Recv(&whichMol, 1, MPI_INT, 0, tag, MPI_COMM_WORLD,&ierr);
243 –          
244 –          for(int i = 0; i < nMols; i++)
245 –            if(molecules[i].getGlobalIndex() == whichMol){
246 –              localIndex = i;
247 –              break;
248 –            }
249 –          
250 –          molecules[localIndex].getCOM(COM);
251 –          zpos = COM[2];          
252 –          MPI_Send(&zpos, 1, MPI_DOUBLE_PRECISION, 0, tag, MPI_COMM_WORLD);      
253 –          
254 –          break;
255 –            
256 –        case EndOfRequest :
257 –        
258 –         done = true;
259 –         break;
260 –      }
261 –      
138      }
139 <          
140 <  }
139 >    else{
140 >      
141 >      parameters = zConsParaData->getData();
142  
143 <  //Brocast the allRefZ to slave nodes;
144 <  double* allRefZBuf;
145 <  int nZConsMols;
269 <  nZConsMols = indexOfAllZConsMols.size();
270 <  
271 <  allRefZBuf = new double[nZConsMols];
272 <  
273 <  if(worldRank == 0){
143 >      //check the range of zconsIndex
144 >      //and the minimum value of index is the first one (we already sorted the data)
145 >      //the maximum value of index is the last one
146  
147 <    for(int i = 0; i < nZConsMols; i++)
148 <      allRefZBuf[i] = allRefZ[i];
149 <  }    
150 <  
151 <    MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
152 <  
153 <  if(worldRank != 0){
154 <    
155 <    for(int i = 0; i < nZConsMols; i++)
156 <      allRefZ.push_back(allRefZBuf[i]);  
157 <  }
158 <  
159 <  delete[] allRefZBuf;
147 >      int maxIndex;
148 >      int minIndex;
149 >      int totalNumMol;
150 >
151 >      minIndex = (*parameters)[0].zconsIndex;
152 >      if(minIndex < 0){
153 >        sprintf( painCave.errMsg,
154 >               "ZConstraint error: index is out of range\n");
155 >        painCave.isFatal = 1;
156 >        simError();
157 >        }
158 >
159 >      maxIndex = (*parameters)[parameters->size()].zconsIndex;
160 >
161 > #ifndef IS_MPI
162 >      totalNumMol = nMols;
163 > #else
164 >      totalNumMol = mpiSim->getTotNmol();  
165 > #endif      
166 >      
167 >      if(maxIndex > totalNumMol - 1){
168 >        sprintf( painCave.errMsg,
169 >               "ZConstraint error: index is out of range\n");
170 >        painCave.isFatal = 1;
171 >        simError();                  
172 >      }
173 >
174 >      //if user does not specify the zpos for the zconstraint molecule
175 >      //its initial z coordinate  will be used as default
176 >      for(int i = 0; i < parameters->size(); i++){
177 >
178 >              if(!(*parameters)[i].havingZPos){
179 >
180 > #ifndef IS_MPI
181 >            for(int j = 0; j < nMols; j++){
182 >              if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){
183 >                 molecules[i].getCOM(COM);
184 >                          break;
185 >              }
186 >            }
187 > #else
188 >            //query which processor current zconstraint molecule belongs to
189 >           int *MolToProcMap;
190 >           int whichNode;
191 >                         double initZPos;
192 >           MolToProcMap = mpiSim->getMolToProcMap();
193 >           whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
194 >                          
195 >           //broadcast the zpos of current z-contraint molecule
196 >           //the node which contain this
197 >          
198 >           if (worldRank == whichNode ){
199 >                                                
200 >             for(int i = 0; i < nMols; i++)
201 >               if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){
202 >                 molecules[i].getCOM(COM);
203 >                                         break;
204 >               }
205 >                                
206 >           }
207 >
208 >            MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);                          
209   #endif
210 +            
211 +                 (*parameters)[i].zPos = COM[whichDirection];
212  
213 <  
213 >            sprintf( painCave.errMsg,
214 >                     "ZConstraint warningr: Does not specify zpos for z-constraint molecule "
215 >                     "initial z coornidate will be used \n");
216 >            painCave.isFatal = 0;
217 >            simError();  
218 >          
219 >              }
220 >            }
221 >                        
222 >    }//end if (!zConsParaData)
223 >  }//end  if (!data)
224 >            
225 > //  
226   #ifdef IS_MPI
227    update();
228   #else  
229    int searchResult;
230 <  
296 <  refZ = allRefZ;
297 <
230 >      
231    for(int i = 0; i < nMols; i++){
232      
233      searchResult = isZConstraintMol(&molecules[i]);
# Line 303 | Line 236 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
236      
237        zconsMols.push_back(&molecules[i]);      
238        massOfZConsMols.push_back(molecules[i].getTotalMass());  
239 <      
239 >
240 >      zPos.push_back((*parameters)[searchResult].zPos);
241 >           kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
242 >      
243        molecules[i].getCOM(COM);
244      }
245      else
# Line 329 | Line 265 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
265      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
266    
267   #endif
268 <  
268 >
269 >  //get total number of unconstrained atoms
270 >  int nUnconsAtoms_local;
271 >  nUnconsAtoms_local = 0;
272 >  for(int i = 0; i < unconsMols.size(); i++)
273 >    nUnconsAtoms_local += unconsMols[i]->getNAtoms();
274 >    
275 > #ifndef IS_MPI
276 >  totNumOfUnconsAtoms = nUnconsAtoms_local;
277 > #else
278 >  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
279 > #endif
280 >
281 >  checkZConsState();
282 >
283 >  //  
284    fzOut = new ZConsWriter(zconsOutput.c_str());  
285    
286    if(!fzOut){
# Line 339 | Line 290 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
290      simError();
291    }
292    
342 –  fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
293   }
294  
295   template<typename T> ZConstraint<T>::~ZConstraint()
# Line 362 | Line 312 | template<typename T> void ZConstraint<T>::update()
312    
313    zconsMols.clear();
314    massOfZConsMols.clear();
315 <  refZ.clear();
315 >  zPos.clear();
316 >  kz.clear();
317    
318    unconsMols.clear();
319    massOfUnconsMols.clear();
# Line 376 | Line 327 | template<typename T> void ZConstraint<T>::update()
327      if(index > -1){
328      
329        zconsMols.push_back(&molecules[i]);      
330 +      zPos.push_back((*parameters)[index].zPos);
331 +        kz.push_back((*parameters)[index].kRatio * zForceConst);
332 +
333        massOfZConsMols.push_back(molecules[i].getTotalMass());  
334        
335        molecules[i].getCOM(COM);
382 –      refZ.push_back(allRefZ[index]);      
336      }
337      else
338      {
# Line 441 | Line 394 | template<typename T> int ZConstraint<T>::isZConstraint
394    index = mol->getGlobalIndex();
395    
396    low = 0;
397 <  high = indexOfAllZConsMols.size() - 1;
397 >  high = parameters->size() - 1;
398    
399    //Binary Search (we have sorted the array)  
400    while(low <= high){
401      mid = (low + high) /2;
402 <    if (indexOfAllZConsMols[mid] == index)
402 >    if ((*parameters)[mid].zconsIndex == index)
403        return mid;
404 <    else if (indexOfAllZConsMols[mid] > index )
404 >    else if ((*parameters)[mid].zconsIndex > index )
405         high = mid -1;
406      else    
407        low = mid + 1;
# Line 457 | Line 410 | template<typename T> int ZConstraint<T>::isZConstraint
410    return -1;
411   }
412  
413 < /** Function Name: integrateStep
414 < ** Parameter:
415 < **   int calcPot;
463 < **   int calcStress;
464 < ** Description:
465 < **  Advance One Step.
466 < ** Memo:
467 < **   The best way to implement z-constraint is to override integrateStep
468 < **   Overriding constrainB is not a good choice, since in integrateStep,
469 < **   constrainB is invoked by below line,
470 < **                  if(nConstrained) constrainB();
471 < **   For instance, we would like to apply z-constraint without bond contrain,
472 < **   In that case, if we override constrainB, Z-constrain method will never be executed;
413 > /**
414 > * Description:
415 > *  Reset the z coordinates
416   */
417 < template<typename T> void ZConstraint<T>::integrateStep( int calcPot, int calcStress )
475 < {
476 <  T::integrateStep( calcPot, calcStress );
477 <  resetZ();
417 > template<typename T> void ZConstraint<T>::integrate(){
418    
419 <  double currZConsTime = 0;
420 <  
421 <  //write out forces of z constraint
422 <  if( info->getTime() >= currZConsTime){  
423 <      fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
424 <  }    
419 >  //zero out the velocities of center of mass of unconstrained molecules
420 >  //and the velocities of center of mass of every single z-constrained molecueles
421 >  zeroOutVel();
422 >  
423 >  T::integrate();
424 >
425   }
426 +
427  
428 < /** Function Name: resetZ
429 < ** Description:
430 < **  Reset the z coordinates
428 > /**
429 > *
430 > *
431 > *
432 > *
433   */
434  
435 < template<typename T> void ZConstraint<T>::resetZ()
435 >
436 > template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
437 >
438 >  T::calcForce(calcPot, calcStress);
439 >
440 >  if (checkZConsState())
441 >    zeroOutVel();
442 >
443 >  //do zconstraint force;
444 >  if (haveFixedZMols())
445 >    this->doZconstraintForce();
446 >  
447 >  //use harmonical poteintial to move the molecules to the specified positions
448 >  if (haveMovingZMols())
449 >    //this->doHarmonic();
450 >  
451 >  fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
452 >      
453 > }
454 >
455 > template<typename T> double ZConstraint<T>::calcZSys()
456   {
457 <  double deltaZ;
458 <  double mzOfZCons;   //total sum of m*z of z-constrain molecules
459 <  double mzOfUncons; //total sum of m*z of unconstrain molecuels;
460 <  double totalMZOfZCons;
461 <  double totalMZOfUncons;
457 >  //calculate reference z coordinate for z-constraint molecules
458 >  double totalMass_local;
459 >  double totalMass;
460 >  double totalMZ_local;
461 >  double totalMZ;
462 >  double massOfUncons_local;
463 >  double massOfCurMol;
464    double COM[3];
465 +  
466 +  totalMass_local = 0;
467 +  totalMass = 0;
468 +  totalMZ_local = 0;
469 +  totalMZ = 0;
470 +  massOfUncons_local = 0;
471 +    
472 +  
473 +  for(int i = 0; i < nMols; i++){
474 +    massOfCurMol = molecules[i].getTotalMass();
475 +    molecules[i].getCOM(COM);
476 +    
477 +    totalMass_local += massOfCurMol;
478 +    totalMZ_local += massOfCurMol * COM[whichDirection];
479 +    
480 +    if(isZConstraintMol(&molecules[i]) == -1){
481 +    
482 +      massOfUncons_local += massOfCurMol;
483 +    }  
484 +    
485 +  }
486 +  
487 +  
488 + #ifdef IS_MPI  
489 +  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
490 +  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
491 +  MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
492 + #else
493 +  totalMass = totalMass_local;
494 +  totalMZ = totalMZ_local;
495 +  totalMassOfUncons = massOfUncons_local;
496 + #endif  
497 +
498    double zsys;
499 <  Atom** zconsAtoms;
499 >  zsys = totalMZ / totalMass;
500  
501 <  mzOfZCons = 0;
502 <  mzOfUncons  = 0;
501 >  return zsys;
502 > }
503 >
504 > /**
505 > *
506 > */
507 > template<typename T> void ZConstraint<T>::thermalize( void ){
508 >
509 >  T::thermalize();
510 >  zeroOutVel();
511 > }
512 >
513 > /**
514 > *
515 > *
516 > *
517 > */
518 >
519 > template<typename T> void ZConstraint<T>::zeroOutVel(){
520 >
521 >  Atom** fixedZAtoms;  
522 >  double COMvel[3];
523 >  double vel[3];
524    
525 +  //zero out the velocities of center of mass of fixed z-constrained molecules
526 +  
527    for(int i = 0; i < zconsMols.size(); i++){
528 <    mzOfZCons += massOfZConsMols[i] * refZ[i];    
528 >
529 >    if (states[i] == zcsFixed){
530 >
531 >        zconsMols[i]->getCOMvel(COMvel);      
532 >      fixedZAtoms = zconsMols[i]->getMyAtoms();
533 >          
534 >      for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
535 >        fixedZAtoms[j]->getVel(vel);
536 >          vel[whichDirection] -= COMvel[whichDirection];
537 >          fixedZAtoms[j]->setVel(vel);
538 >      }
539 >          
540 >    }
541 >        
542    }
543 +  
544 +  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
545 +  double MVzOfMovingMols_local;
546 +  double MVzOfMovingMols;
547 +  double totalMassOfMovingZMols_local;
548 +  double totalMassOfMovingZMols;
549 +      
550 +  MVzOfMovingMols_local = 0;
551 +  totalMassOfMovingZMols_local = 0;
552  
553 < #ifdef IS_MPI
554 <  MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
553 >  for(int i =0; i < unconsMols.size(); i++){
554 >    unconsMols[i]->getCOMvel(COMvel);
555 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
556 >  }
557 >
558 >  for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
559 >
560 >    if (states[i] == zcsMoving){
561 >      zconsMols[i]->getCOMvel(COMvel);
562 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
563 >      totalMassOfMovingZMols_local += massOfZConsMols[i];              
564 >    }
565 >                
566 >  }
567 >
568 > #ifndef IS_MPI
569 >  MVzOfMovingMols = MVzOfMovingMols_local;
570 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
571   #else
572 <  totalMZOfZCons = mzOfZCons;
572 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
573 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
574   #endif
575  
576 +  double vzOfMovingMols;
577 +  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
578 +
579 +  //modify the velocites of unconstrained molecules  
580 +  Atom** unconsAtoms;
581    for(int i = 0; i < unconsMols.size(); i++){
517 –    unconsMols[i]->getCOM(COM);
518 –    mzOfUncons += massOfUnconsMols[i] * COM[2];
519 –  }
582    
583 < #ifdef IS_MPI
584 <  MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
585 < #else
586 <  totalMZOfUncons = mzOfUncons;
587 < #endif  
583 >    unconsAtoms = unconsMols[i]->getMyAtoms();
584 >    for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
585 >      unconsAtoms[j]->getVel(vel);
586 >      vel[whichDirection] -= vzOfMovingMols;
587 >      unconsAtoms[j]->setVel(vel);
588 >    }
589    
590 <  zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons;
590 >  }  
591  
592 <  cout << "current time: " << info->getTime() <<endl;  
593 <  for(int i = 0; i < zconsMols.size(); i++){  
592 >  //modify the velocities of moving z-constrained molecuels
593 >  Atom** movingZAtoms;
594 >  for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
595 >
596 >    if (states[i] ==zcsMoving){
597    
598 <    zconsMols[i]->getCOM(COM);
598 >      movingZAtoms = zconsMols[i]->getMyAtoms();
599 >        for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
600 >        movingZAtoms[j]->getVel(vel);
601 >        vel[whichDirection] -= vzOfMovingMols;
602 >          movingZAtoms[j]->setVel(vel);
603 >        }
604 >          
605 >    }
606 >
607 >  }
608 >
609 > }
610 >
611 > template<typename T> void ZConstraint<T>::doZconstraintForce(){
612 >
613 >  Atom** zconsAtoms;
614 >  double totalFZ;
615 >  double totalFZ_local;
616 >  double COMvel[3];  
617 >  double COM[3];
618 >  double force[3];
619 >  double zsys;
620 >
621 >  int nMovingZMols_local;
622 >  int nMovingZMols;
623 >
624 >  //constrain the molecules which do not reach the specified positions  
625 >
626 >   zsys = calcZSys();
627 >   cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;  
628      
629 <    cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t";
630 <    deltaZ = zsys + refZ[i] - COM[2];
631 <    cout << "\tdistance: " << COM[2] +deltaZ - zsys;    
632 <    //update z coordinate    
633 <    zconsAtoms = zconsMols[i]->getMyAtoms();    
634 <    for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
635 <      zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ);  
636 <    }    
629 >  //Zero Out the force of z-contrained molecules    
630 >  totalFZ_local = 0;
631 >
632 >  //calculate the total z-contrained force of fixed z-contrained molecules
633 >  cout << "Fixed Molecules" << endl;
634 >  for(int i = 0; i < zconsMols.size(); i++){
635 >                
636 >    if (states[i] == zcsFixed){
637 >                
638 >      zconsMols[i]->getCOM(COM);
639 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
640 >
641 >      fz[i] = 0;      
642 >      for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
643 >        zconsAtoms[j]->getFrc(force);
644 >        fz[i] += force[whichDirection];      
645 >      }
646 >      totalFZ_local += fz[i];
647 >
648 >      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
649 >
650 >    }
651 >          
652 >  }
653 >
654 >  //calculate the number of atoms of moving z-constrained molecules
655 >  nMovingZMols_local = 0;
656 >  for(int i = 0; zconsMols.size(); i++){
657 >    if(states[i] == zcsMoving)
658 >        nMovingZMols_local += massOfZConsMols[i];
659 >  }
660 > #ifdef IS_MPI
661 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
662 >  MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
663 > #else
664 >  totalFZ = totalFZ_local;
665 >  nMovingZMols = nMovingZMols_local;
666 > #endif
667 >
668 >  force[0]= 0;
669 >  force[1]= 0;
670 >  force[2]= 0;
671 >  force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);
672 >
673 >  //modify the velocites of unconstrained molecules
674 >  for(int i = 0; i < unconsMols.size(); i++){
675 >    
676 >     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
677 >    
678 >     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)          
679 >       unconsAtoms[j]->addFrc(force);
680 >    
681 >  }      
682 >
683 > //modify the velocities of moving z-constrained molecules
684 >  for(int i = 0; i < zconsMols.size(); i++) {
685 >   if (states[i] == zcsMoving){
686 >                
687 >     Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
688 >
689 >     for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
690 >       movingZAtoms[j]->addFrc(force);
691 >     }
692 >  }
693 >
694 >  // apply negative to fixed z-constrained molecues;
695 >  force[0]= 0;
696 >  force[1]= 0;
697 >  force[2]= 0;
698 >
699 >  for(int i = 0; i < zconsMols.size(); i++){
700 >
701 >    if (states[i] == zcsFixed){  
702 >        
703 >      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
704 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
705      
706 <    //calculate z constrain force
707 <    fz[i] = massOfZConsMols[i]* deltaZ / dt2;
708 <    
709 <    cout << "\tforce: " << fz[i] << endl;
706 >      for(int j =0; j < nAtomOfCurZConsMol; j++) {
707 >        force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
708 >        zconsAtoms[j]->addFrc(force);
709 >      }
710 >                
711 >    }
712 >        
713 >  }
714 >
715 > }
716 >
717 > template<typename T> bool ZConstraint<T>::checkZConsState(){
718 >  double COM[3];
719 >  double diff;
720 >  
721 >  bool changed;
722 >  
723 >  changed = false;
724 >  
725 >  for(int i =0; i < zconsMols.size(); i++){
726 >
727 >    zconsMols[i]->getCOM(COM);
728 >    diff = fabs(COM[whichDirection] - zPos[i]);  
729 >    if (  diff <= zconsTol && states[i] == zcsMoving){
730 >      states[i] = zcsFixed;
731 >        changed = true;
732 >    }
733 >    else if ( diff > zconsTol && states[i] == zcsFixed){
734 >      states[i] = zcsMoving;
735 >        changed = true;  
736 >    }
737 >  
738    }
739  
740 <      
740 >  return changed;
741   }
742 +
743 + template<typename T> bool ZConstraint<T>::haveFixedZMols(){
744 +  for(int i = 0; i < zconsMols.size(); i++)
745 +    if (states[i] == zcsFixed)
746 +      return true;
747 +
748 +  return false;
749 + }
750 +
751 +
752 + /**
753 + *
754 + */
755 + template<typename T> bool ZConstraint<T>::haveMovingZMols(){
756 +  for(int i = 0; i < zconsMols.size(); i++)
757 +    if (states[i] == zcsMoving)
758 +      return true;
759 +
760 +  return false;
761 +  
762 + }
763 +
764 + /**
765 +  *
766 +  *
767 +  */
768 +
769 + template<typename T> void ZConstraint<T>::doHarmonic(){
770 +  double force[3];
771 +  double harmonicU;
772 +  double COM[3];
773 +  double diff;
774 +        
775 +  force[0] = 0;
776 +  force[1] = 0;
777 +  force[2] = 0;
778 +
779 +  cout << "Moving Molecules" << endl;  
780 +  for(int i = 0; i < zconsMols.size(); i++) {
781 +
782 +    if (states[i] == zcsMoving){
783 +      zconsMols[i]->getCOM(COM):
784 +      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
785 +                
786 +                diff = COM[whichDirection] -zPos[i];
787 +                
788 +      harmonicU = 0.5 * kz[i] * diff * diff;  
789 +                info->ltPot += harmonicU;
790 +
791 +                force[whichDirection] = - kz[i] * diff / zconsMols[i]->getNAtoms();
792 +                
793 +      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
794 +
795 +       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
796 +         movingZAtoms[j]->addFrc(force);
797 +    }
798 +
799 +  }
800 +
801 + }

Diff Legend

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