ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/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 693 by tim, Wed Aug 13 19:21:53 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;
135 <  totalMass = 0;
136 <  totalMZ_local = 0;
137 <  totalMZ = 0;
138 <  massOfUncons_local = 0;
139 <    
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 <    
153 <  }
154 <  
155 <  
156 < #ifdef IS_MPI  
157 <  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  
165 <
166 <  double zsys;
167 <  zsys = totalMZ / totalMass;
168 <
169 < #ifndef IS_MPI  
170 <  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);  
175 <    }
176 <    
177 <  }
178 < #else
179 <
180 <  int whichNode;
181 <  enum CommType { RequestMolZPos, EndOfRequest} status;
182 <  //int status;
183 <  double zpos;
184 <  int localIndex;
185 <  MPI_Status ierr;
186 <  int tag = 0;
187 <  
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++){
92 >  if(!data) {
93        
94 <      whichNode = MolToProcMap[indexOfAllZConsMols[i]];
95 <      globalIndexOfCurMol = indexOfAllZConsMols[i];
96 <      
97 <      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 <    
94 >    sprintf( painCave.errMsg,
95 >               "ZConstraint error: can not get tolerance \n");
96 >    painCave.isFatal = 1;
97 >    simError();      
98    }
99    else{
100    
101 <    int whichMol;
102 <    bool done = false;
101 >    tolerance = dynamic_cast<DoubleData*>(data);
102 >    
103 >    if(!tolerance){
104  
105 <    while (!done){  
106 <      
107 <      MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
108 <    
109 <      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 <      
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 +        
117 +  //retrieve index of z-constraint molecules
118 +  data = info->getProperty(ZCONSPARADATA_ID);
119 +  if(!data) {
120  
121 <  //Brocast the allRefZ to slave nodes;
122 <  double* allRefZBuf;
123 <  int nZConsMols;
124 <  nZConsMols = indexOfAllZConsMols.size();
125 <  
271 <  allRefZBuf = new double[nZConsMols];
272 <  
273 <  if(worldRank == 0){
274 <
275 <    for(int i = 0; i < nZConsMols; i++)
276 <      allRefZBuf[i] = allRefZ[i];
277 <  }    
278 <  
279 <    MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
280 <  
281 <  if(worldRank != 0){
282 <    
283 <    for(int i = 0; i < nZConsMols; i++)
284 <      allRefZ.push_back(allRefZBuf[i]);  
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 <  delete[] allRefZBuf;
129 >    zConsParaData = dynamic_cast<ZConsParaData*>(data);
130 >    
131 >    if(!zConsParaData){
132 >
133 >      sprintf( painCave.errMsg,
134 >                 "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
135 >      painCave.isFatal = 1;
136 >      simError();  
137 >    
138 >    }
139 >    else{
140 >      
141 >      parameters = zConsParaData->getData();
142 >
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 >      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() - 1].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[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
183 >                 molecules[j].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 j = 0; j < nMols; j++)
201 >               if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
202 >                 molecules[j].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 325 | Line 261 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
261      simError();
262    }
263  
264 <  for(int i = 0; i < zconsMols.size(); i++)
264 >  //determine the states of z-constraint molecules
265 >  for(int i = 0; i < zconsMols.size(); i++){
266      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
267 +
268 +         zconsMols[i]->getCOM(COM);
269 +    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
270 +                states.push_back(zcsFixed);
271 +         else
272 +                states.push_back(zcsMoving);
273 +  }
274    
275   #endif
276 +
277 +  //get total masss of unconstraint molecules
278 +  double totalMassOfUncons_local;
279 +  totalMassOfUncons_local = 0;
280    
281 +  for(int i = 0; i < unconsMols.size(); i++)
282 +    totalMassOfUncons_local += unconsMols[i]->getTotalMass();
283 +    
284 + #ifndef IS_MPI
285 +  totalMassOfUncons = totalMassOfUncons_local;
286 + #else
287 +  MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
288 + #endif
289 +
290 +
291 +  //get total number of unconstrained atoms
292 +  int nUnconsAtoms_local;
293 +  nUnconsAtoms_local = 0;
294 +  for(int i = 0; i < unconsMols.size(); i++)
295 +    nUnconsAtoms_local += unconsMols[i]->getNAtoms();
296 +    
297 + #ifndef IS_MPI
298 +  totNumOfUnconsAtoms = nUnconsAtoms_local;
299 + #else
300 +  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
301 + #endif  
302 +
303 +  // creat zconsWriter  
304    fzOut = new ZConsWriter(zconsOutput.c_str());  
305    
306    if(!fzOut){
# Line 339 | Line 310 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
310      simError();
311    }
312    
342 –  fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
313   }
314  
315   template<typename T> ZConstraint<T>::~ZConstraint()
# Line 362 | Line 332 | template<typename T> void ZConstraint<T>::update()
332    
333    zconsMols.clear();
334    massOfZConsMols.clear();
335 <  refZ.clear();
335 >  zPos.clear();
336 >  kz.clear();
337    
338    unconsMols.clear();
339    massOfUnconsMols.clear();
# Line 376 | Line 347 | template<typename T> void ZConstraint<T>::update()
347      if(index > -1){
348      
349        zconsMols.push_back(&molecules[i]);      
350 +      zPos.push_back((*parameters)[index].zPos);
351 +        kz.push_back((*parameters)[index].kRatio * zForceConst);
352 +                        
353        massOfZConsMols.push_back(molecules[i].getTotalMass());  
354        
355        molecules[i].getCOM(COM);
382 –      refZ.push_back(allRefZ[index]);      
356      }
357      else
358      {
# Line 389 | Line 362 | template<typename T> void ZConstraint<T>::update()
362  
363      }
364    }
365 +
366 +  //determine the states of z-constraint molecules
367 +  for(int i = 0; i < zconsMols.size(); i++){
368 +           zconsMols[i]->getCOM(COM);
369 +      if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
370 +                  states.push_back(zcsFixed);
371 +           else
372 +                  states.push_back(zcsMoving);
373 +  }
374 +
375      
376    //The reason to declare fz and indexOfZconsMols as pointer to array is
377    // that we want to make the MPI communication simple
# Line 441 | Line 424 | template<typename T> int ZConstraint<T>::isZConstraint
424    index = mol->getGlobalIndex();
425    
426    low = 0;
427 <  high = indexOfAllZConsMols.size() - 1;
427 >  high = parameters->size() - 1;
428    
429    //Binary Search (we have sorted the array)  
430    while(low <= high){
431      mid = (low + high) /2;
432 <    if (indexOfAllZConsMols[mid] == index)
432 >    if ((*parameters)[mid].zconsIndex == index)
433        return mid;
434 <    else if (indexOfAllZConsMols[mid] > index )
434 >    else if ((*parameters)[mid].zconsIndex > index )
435         high = mid -1;
436      else    
437        low = mid + 1;
# Line 457 | Line 440 | template<typename T> int ZConstraint<T>::isZConstraint
440    return -1;
441   }
442  
443 < /** Function Name: integrateStep
461 < ** Parameter:
462 < **   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;
473 < */
474 < template<typename T> void ZConstraint<T>::integrateStep( int calcPot, int calcStress )
475 < {
476 <  T::integrateStep( calcPot, calcStress );
477 <  resetZ();
443 > template<typename T> void ZConstraint<T>::integrate(){
444    
445 <  double currZConsTime = 0;
446 <  
447 <  //write out forces of z constraint
448 <  if( info->getTime() >= currZConsTime){  
449 <      fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
450 <  }    
445 >  //zero out the velocities of center of mass of unconstrained molecules
446 >  //and the velocities of center of mass of every single z-constrained molecueles
447 >  zeroOutVel();
448 >  
449 >  T::integrate();
450 >
451   }
452 +
453  
454 < /** Function Name: resetZ
455 < ** Description:
456 < **  Reset the z coordinates
454 > /**
455 > *
456 > *
457 > *
458 > *
459   */
460  
461 < template<typename T> void ZConstraint<T>::resetZ()
461 >
462 > template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
463 >  double zsys;
464 >
465 >  T::calcForce(calcPot, calcStress);
466 >
467 >  if (checkZConsState())
468 >  zeroOutVel();
469 >  
470 >  zsys = calcZSys();
471 >  cout << "---------------------------------------------------------------------" <<endl;
472 >  cout << "current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;      
473 >  cout << "before calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
474 >        
475 >
476 >  //do zconstraint force;
477 >  if (haveFixedZMols())
478 >    this->doZconstraintForce();
479 >
480 >
481 >      
482 >  //use harmonical poteintial to move the molecules to the specified positions
483 >  if (haveMovingZMols())
484 >    this->doHarmonic();
485 >  
486 >  fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
487 >  cout << "after calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
488 > }
489 >
490 > template<typename T> double ZConstraint<T>::calcZSys()
491   {
492 <  double deltaZ;
493 <  double mzOfZCons;   //total sum of m*z of z-constrain molecules
494 <  double mzOfUncons; //total sum of m*z of unconstrain molecuels;
495 <  double totalMZOfZCons;
496 <  double totalMZOfUncons;
492 >  //calculate reference z coordinate for z-constraint molecules
493 >  double totalMass_local;
494 >  double totalMass;
495 >  double totalMZ_local;
496 >  double totalMZ;
497 >  double massOfUncons_local;
498 >  double massOfCurMol;
499    double COM[3];
500 +  
501 +  totalMass_local = 0;
502 +  totalMass = 0;
503 +  totalMZ_local = 0;
504 +  totalMZ = 0;
505 +  massOfUncons_local = 0;
506 +    
507 +  
508 +  for(int i = 0; i < nMols; i++){
509 +    massOfCurMol = molecules[i].getTotalMass();
510 +    molecules[i].getCOM(COM);
511 +    
512 +    totalMass_local += massOfCurMol;
513 +    totalMZ_local += massOfCurMol * COM[whichDirection];
514 +    
515 +    if(isZConstraintMol(&molecules[i]) == -1){
516 +    
517 +      massOfUncons_local += massOfCurMol;
518 +    }  
519 +    
520 +  }
521 +  
522 +  
523 + #ifdef IS_MPI  
524 +  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
525 +  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
526 +  MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
527 + #else
528 +  totalMass = totalMass_local;
529 +  totalMZ = totalMZ_local;
530 +  totalMassOfUncons = massOfUncons_local;
531 + #endif  
532 +
533    double zsys;
534 <  Atom** zconsAtoms;
534 >  zsys = totalMZ / totalMass;
535  
536 <  mzOfZCons = 0;
537 <  mzOfUncons  = 0;
536 >  return zsys;
537 > }
538 >
539 > /**
540 > *
541 > */
542 > template<typename T> void ZConstraint<T>::thermalize( void ){
543 >
544 >  T::thermalize();
545 >  zeroOutVel();
546 > }
547 >
548 > /**
549 > *
550 > *
551 > *
552 > */
553 >
554 > template<typename T> void ZConstraint<T>::zeroOutVel(){
555 >
556 >  Atom** fixedZAtoms;  
557 >  double COMvel[3];
558 >  double vel[3];
559 >
560 >  double tempMass = 0;
561 >  double tempMVz =0;
562 >  double tempVz = 0;
563 >  for(int i = 0; i < nMols; i++){
564 >    molecules[i].getCOMvel(COMvel);
565 >    tempMass +=molecules[i].getTotalMass();  
566 >         tempMVz += molecules[i].getTotalMass() * COMvel[whichDirection];
567 >         tempVz += COMvel[whichDirection];
568 >  }
569 >  cout << "initial velocity of center of mass is " << tempMVz / tempMass << endl;
570 >
571 >  //zero out the velocities of center of mass of fixed z-constrained molecules
572    
573    for(int i = 0; i < zconsMols.size(); i++){
574 <    mzOfZCons += massOfZConsMols[i] * refZ[i];    
574 >
575 >    if (states[i] == zcsFixed){
576 >
577 >           zconsMols[i]->getCOMvel(COMvel);      
578 >                cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
579 >
580 >      fixedZAtoms = zconsMols[i]->getMyAtoms();
581 >          
582 >      for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
583 >        fixedZAtoms[j]->getVel(vel);
584 >             vel[whichDirection] -= COMvel[whichDirection];
585 >             fixedZAtoms[j]->setVel(vel);
586 >      }
587 >
588 >                zconsMols[i]->getCOMvel(COMvel);
589 >                cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
590 >    }
591 >        
592    }
593  
594 < #ifdef IS_MPI
595 <  MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
594 >        cout << "before resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;    
595 >                  
596 >  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
597 >  double MVzOfMovingMols_local;
598 >  double MVzOfMovingMols;
599 >  double totalMassOfMovingZMols_local;
600 >  double totalMassOfMovingZMols;
601 >      
602 >  MVzOfMovingMols_local = 0;
603 >  totalMassOfMovingZMols_local = 0;
604 >
605 >  for(int i =0; i < unconsMols.size(); i++){
606 >    unconsMols[i]->getCOMvel(COMvel);
607 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
608 >  }
609 >
610 >  for(int i = 0; i < zconsMols.size(); i++){
611 >    if (states[i] == zcsMoving){
612 >      zconsMols[i]->getCOMvel(COMvel);
613 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
614 >      totalMassOfMovingZMols_local += massOfZConsMols[i];              
615 >    }
616 >                
617 >  }
618 >
619 > #ifndef IS_MPI
620 >  MVzOfMovingMols = MVzOfMovingMols_local;
621 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
622   #else
623 <  totalMZOfZCons = mzOfZCons;
623 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
624 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
625   #endif
626  
627 +  double vzOfMovingMols;
628 +  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
629 +
630 +  //modify the velocites of unconstrained molecules  
631 +  Atom** unconsAtoms;
632    for(int i = 0; i < unconsMols.size(); i++){
517 –    unconsMols[i]->getCOM(COM);
518 –    mzOfUncons += massOfUnconsMols[i] * COM[2];
519 –  }
633    
634 < #ifdef IS_MPI
635 <  MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
636 < #else
637 <  totalMZOfUncons = mzOfUncons;
638 < #endif  
634 >    unconsAtoms = unconsMols[i]->getMyAtoms();
635 >    for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
636 >      unconsAtoms[j]->getVel(vel);
637 >      vel[whichDirection] -= vzOfMovingMols;
638 >      unconsAtoms[j]->setVel(vel);
639 >    }
640    
641 <  zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons;
641 >  }  
642  
643 <  cout << "current time: " << info->getTime() <<endl;  
644 <  for(int i = 0; i < zconsMols.size(); i++){  
643 >  //modify the velocities of moving z-constrained molecuels
644 >  Atom** movingZAtoms;
645 >  for(int i = 0; i < zconsMols.size(); i++){
646 >
647 >    if (states[i] ==zcsMoving){
648    
649 <    zconsMols[i]->getCOM(COM);
649 >      movingZAtoms = zconsMols[i]->getMyAtoms();
650 >           for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
651 >        movingZAtoms[j]->getVel(vel);
652 >        vel[whichDirection] -= vzOfMovingMols;
653 >             movingZAtoms[j]->setVel(vel);
654 >          }
655 >          
656 >   }
657 >
658 > }
659 >
660 >        cout << "after resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
661 >
662 > }
663 >
664 > template<typename T> void ZConstraint<T>::doZconstraintForce(){
665 >
666 >  Atom** zconsAtoms;
667 >  double totalFZ;
668 >  double totalFZ_local;
669 >  double COMvel[3];  
670 >  double COM[3];
671 >  double force[3];
672 >
673 >  int nMovingZMols_local;
674 >  int nMovingZMols;
675 >
676 >  //constrain the molecules which do not reach the specified positions  
677      
678 <    cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t";
679 <    deltaZ = zsys + refZ[i] - COM[2];
680 <    cout << "\tdistance: " << COM[2] +deltaZ - zsys;    
681 <    //update z coordinate    
682 <    zconsAtoms = zconsMols[i]->getMyAtoms();    
683 <    for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
684 <      zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ);  
685 <    }    
678 >  //Zero Out the force of z-contrained molecules    
679 >  totalFZ_local = 0;
680 >
681 >  //calculate the total z-contrained force of fixed z-contrained molecules
682 >  cout << "Fixed Molecules" << endl;
683 >  for(int i = 0; i < zconsMols.size(); i++){
684 >                
685 >    if (states[i] == zcsFixed){
686 >                
687 >      zconsMols[i]->getCOM(COM);
688 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
689 >
690 >      fz[i] = 0;      
691 >      for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
692 >        zconsAtoms[j]->getFrc(force);
693 >        fz[i] += force[whichDirection];      
694 >      }
695 >      totalFZ_local += fz[i];
696 >
697 >      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
698 >
699 >    }
700 >          
701 >  }
702 >
703 >  //calculate the number of atoms of moving z-constrained molecules
704 >  nMovingZMols_local = 0;
705 >  for(int i = 0; i < zconsMols.size(); i++)
706 >    if(states[i] == zcsMoving)
707 >           nMovingZMols_local += massOfZConsMols[i];
708 >  
709 > #ifdef IS_MPI
710 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
711 >  MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
712 > #else
713 >  totalFZ = totalFZ_local;
714 >  nMovingZMols = nMovingZMols_local;
715 > #endif
716 >
717 >  force[0]= 0;
718 >  force[1]= 0;
719 >  force[2]= 0;
720 >  force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);
721 >
722 >  //modify the forces of unconstrained molecules
723 >  for(int i = 0; i < unconsMols.size(); i++){
724 >    
725 >     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
726 >    
727 >     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)          
728 >       unconsAtoms[j]->addFrc(force);
729 >    
730 >  }      
731 >
732 > //modify the forces of moving z-constrained molecules
733 >  for(int i = 0; i < zconsMols.size(); i++) {
734 >   if (states[i] == zcsMoving){
735 >                
736 >     Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
737 >
738 >     for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
739 >       movingZAtoms[j]->addFrc(force);
740 >     }
741 >  }
742 >
743 >  // apply negative to fixed z-constrained molecues;
744 >  force[0]= 0;
745 >  force[1]= 0;
746 >  force[2]= 0;
747 >
748 >  for(int i = 0; i < zconsMols.size(); i++){
749 >
750 >    if (states[i] == zcsFixed){  
751 >        
752 >      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
753 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
754      
755 <    //calculate z constrain force
756 <    fz[i] = massOfZConsMols[i]* deltaZ / dt2;
757 <    
758 <    cout << "\tforce: " << fz[i] << endl;
755 >      for(int j =0; j < nAtomOfCurZConsMol; j++) {
756 >        force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
757 >        zconsAtoms[j]->addFrc(force);
758 >      }
759 >                
760 >    }
761 >        
762 >  }
763 >
764 > }
765 >
766 > template<typename T> bool ZConstraint<T>::checkZConsState(){
767 >  double COM[3];
768 >  double diff;
769 >  
770 >  bool changed;
771 >  
772 >  changed = false;
773 >  
774 >  for(int i =0; i < zconsMols.size(); i++){
775 >
776 >    zconsMols[i]->getCOM(COM);
777 >    diff = fabs(COM[whichDirection] - zPos[i]);  
778 >    if (  diff <= zconsTol && states[i] == zcsMoving){
779 >      states[i] = zcsFixed;
780 >        changed = true;
781 >    }
782 >    else if ( diff > zconsTol && states[i] == zcsFixed){
783 >      states[i] = zcsMoving;
784 >        changed = true;  
785 >    }
786 >  
787    }
788  
789 +  return changed;
790 + }
791 +
792 + template<typename T> bool ZConstraint<T>::haveFixedZMols(){
793 +  for(int i = 0; i < zconsMols.size(); i++)
794 +    if (states[i] == zcsFixed)
795 +      return true;
796 +
797 +  return false;
798 + }
799 +
800 +
801 + /**
802 + *
803 + */
804 + template<typename T> bool ZConstraint<T>::haveMovingZMols(){
805 +  for(int i = 0; i < zconsMols.size(); i++)
806 +    if (states[i] == zcsMoving)
807 +      return true;
808 +
809 +  return false;
810 +  
811 + }
812 +
813 + /**
814 +  *
815 +  *
816 +  */
817 +
818 + template<typename T> void ZConstraint<T>::doHarmonic(){
819 +  double force[3];
820 +  double harmonicU;
821 +  double harmonicF;
822 +  double COM[3];
823 +  double diff;
824 +  double totalFZ;
825 +        
826 +  force[0] = 0;
827 +  force[1] = 0;
828 +  force[2] = 0;
829 +
830 +  totalFZ = 0;
831 +
832 +  cout << "Moving Molecules" << endl;  
833 +  for(int i = 0; i < zconsMols.size(); i++) {
834 +
835 +    if (states[i] == zcsMoving){
836 +      zconsMols[i]->getCOM(COM);
837 +      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
838 +                
839 +                diff = COM[whichDirection] -zPos[i];
840 +                
841 +      harmonicU = 0.5 * kz[i] * diff * diff;  
842 +                info->lrPot += harmonicU;
843 +
844 +      harmonicF =  - kz[i] * diff / zconsMols[i]->getNAtoms();
845 +                force[whichDirection] = harmonicF;
846 +      totalFZ += harmonicF;
847 +                
848 +      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
849 +
850 +       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
851 +         movingZAtoms[j]->addFrc(force);
852 +    }
853 +
854 +  }
855 +
856 +  force[0]= 0;
857 +  force[1]= 0;
858 +  force[2]= 0;
859 +  force[whichDirection] = -totalFZ /totNumOfUnconsAtoms;
860 +
861 +  //modify the forces of unconstrained molecules
862 +  for(int i = 0; i < unconsMols.size(); i++){
863 +    
864 +     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
865 +    
866 +     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)          
867 +       unconsAtoms[j]->addFrc(force);    
868 +  }  
869 +
870 + }
871 +
872 + template<typename T> double ZConstraint<T>::calcCOMVel()
873 + {
874 +  double MVzOfMovingMols_local;
875 +  double MVzOfMovingMols;
876 +  double totalMassOfMovingZMols_local;
877 +  double totalMassOfMovingZMols;
878 +  double COMvel[3];
879        
880 +  MVzOfMovingMols_local = 0;
881 +  totalMassOfMovingZMols_local = 0;
882 +
883 +  for(int i =0; i < unconsMols.size(); i++){
884 +    unconsMols[i]->getCOMvel(COMvel);
885 +    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
886 +  }
887 +
888 +  for(int i = 0; i < zconsMols.size(); i++){
889 +
890 +    if (states[i] == zcsMoving){
891 +      zconsMols[i]->getCOMvel(COMvel);
892 +      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
893 +      totalMassOfMovingZMols_local += massOfZConsMols[i];              
894 +    }
895 +                
896 +  }
897 +
898 + #ifndef IS_MPI
899 +  MVzOfMovingMols = MVzOfMovingMols_local;
900 +  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
901 + #else
902 +  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
903 +  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
904 + #endif
905 +
906 +  double vzOfMovingMols;
907 +  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
908 +
909 +  return vzOfMovingMols;
910   }
911 +
912 +
913 + template<typename T> double ZConstraint<T>::calcCOMVel2()
914 + {
915 +  double COMvel[3];
916 +  double tempMVz = 0;
917 +  int index;
918 +                
919 +  for(int i =0 ; i < nMols; i++){
920 +         index = isZConstraintMol(&molecules[i]);
921 +    if( index == -1){
922 +       molecules[i].getCOMvel(COMvel);
923 +            tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection];
924 +    }
925 +         else if(states[i] == zcsMoving){
926 +       zconsMols[index]->getCOMvel(COMvel);
927 +            tempMVz += massOfZConsMols[index]*COMvel[whichDirection];    
928 +         }
929 +  }
930 +        
931 +  return tempMVz /totalMassOfUncons;
932 + }

Diff Legend

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