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 696 by tim, Thu Aug 14 16:16:39 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 <  }
141 <
266 <  //Brocast the allRefZ to slave nodes;
267 <  double* allRefZBuf;
268 <  int nZConsMols;
269 <  nZConsMols = indexOfAllZConsMols.size();
270 <  
271 <  allRefZBuf = new double[nZConsMols];
272 <  
273 <  if(worldRank == 0){
139 >    else{
140 >      
141 >      parameters = zConsParaData->getData();
142  
143 <    for(int i = 0; i < nZConsMols; i++)
144 <      allRefZBuf[i] = allRefZ[i];
145 <  }    
146 <  
147 <    MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
148 <  
149 <  if(worldRank != 0){
150 <    
151 <    for(int i = 0; i < nZConsMols; i++)
152 <      allRefZ.push_back(allRefZBuf[i]);  
153 <  }
154 <  
155 <  delete[] allRefZBuf;
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
457 < */
454 > /**
455 > *
456 > *
457 > *
458 > *
459 > */
460 > template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
461 >  double zsys;
462  
463 < template<typename T> void ZConstraint<T>::resetZ()
463 >  T::calcForce(calcPot, calcStress);
464 >
465 >  if (checkZConsState())
466 >    zeroOutVel();
467 >  
468 >  zsys = calcZSys();
469 >  cout << "---------------------------------------------------------------------" <<endl;
470 >  cout << "current time: " << info->getTime() << endl;
471 >  cout << "center of mass at z: " << zsys << endl;      
472 >  //cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
473 >  cout << "before calcForce, the COMVel of system is " << calcSysCOMVel() <<endl;
474 >
475 >  //cout <<      "before doZConstraintForce, totalForce is " << calcTotalForce() << endl;
476 >
477 >  //do zconstraint force;
478 >  if (haveFixedZMols())
479 >    this->doZconstraintForce();
480 >    
481 >  //use harmonical poteintial to move the molecules to the specified positions
482 >  if (haveMovingZMols())
483 >    this->doHarmonic();
484 >
485 >  //cout <<      "after doHarmonic, totalForce is " << calcTotalForce() << endl;
486 >
487 >  fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
488 >  //cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
489 >  cout << "after calcForce, the COMVel of system is " << calcSysCOMVel() <<endl;
490 > }
491 >
492 > template<typename T> double ZConstraint<T>::calcZSys()
493   {
494 <  double deltaZ;
495 <  double mzOfZCons;   //total sum of m*z of z-constrain molecules
496 <  double mzOfUncons; //total sum of m*z of unconstrain molecuels;
497 <  double totalMZOfZCons;
498 <  double totalMZOfUncons;
494 >  //calculate reference z coordinate for z-constraint molecules
495 >  double totalMass_local;
496 >  double totalMass;
497 >  double totalMZ_local;
498 >  double totalMZ;
499 >  double massOfCurMol;
500    double COM[3];
501 +        
502 +  totalMass_local = 0;
503 +  totalMZ_local = 0;
504 +  
505 +  for(int i = 0; i < nMols; i++){
506 +    massOfCurMol = molecules[i].getTotalMass();
507 +    molecules[i].getCOM(COM);
508 +    
509 +    totalMass_local += massOfCurMol;
510 +    totalMZ_local += massOfCurMol * COM[whichDirection];
511 +
512 +  }
513 +
514 +  
515 + #ifdef IS_MPI  
516 +  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
517 +  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
518 + #else
519 +  totalMass = totalMass_local;
520 +  totalMZ = totalMZ_local;
521 + #endif  
522 +
523    double zsys;
524 <  Atom** zconsAtoms;
524 >  zsys = totalMZ / totalMass;
525  
526 <  mzOfZCons = 0;
527 <  mzOfUncons  = 0;
526 >  return zsys;
527 > }
528 >
529 > /**
530 > *
531 > */
532 > template<typename T> void ZConstraint<T>::thermalize( void ){
533 >
534 >  T::thermalize();
535 >  zeroOutVel();
536 > }
537 >
538 > /**
539 > *
540 > *
541 > *
542 > */
543 >
544 > template<typename T> void ZConstraint<T>::zeroOutVel(){
545 >
546 >  Atom** fixedZAtoms;  
547 >  double COMvel[3];
548 >  double vel[3];
549 >
550 >  //zero out the velocities of center of mass of fixed z-constrained molecules
551    
552    for(int i = 0; i < zconsMols.size(); i++){
553 <    mzOfZCons += massOfZConsMols[i] * refZ[i];    
553 >
554 >    if (states[i] == zcsFixed){
555 >
556 >           zconsMols[i]->getCOMvel(COMvel);      
557 >                //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
558 >
559 >      fixedZAtoms = zconsMols[i]->getMyAtoms();
560 >          
561 >      for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
562 >        fixedZAtoms[j]->getVel(vel);
563 >             vel[whichDirection] -= COMvel[whichDirection];
564 >             fixedZAtoms[j]->setVel(vel);
565 >      }
566 >
567 >                zconsMols[i]->getCOMvel(COMvel);
568 >                //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
569 >    }
570 >        
571    }
572  
573 < #ifdef IS_MPI
574 <  MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
573 >        //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;      
574 >        cout << "before resetting the COMVel of sytem is " << calcSysCOMVel() << endl;  
575 >                  
576 >  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
577 >  double MVzOfMovingMols_local;
578 >  double MVzOfMovingMols;
579 >  double totalMassOfMovingZMols_local;
580 >  double totalMassOfMovingZMols;
581 >      
582 >  MVzOfMovingMols_local = 0;
583 >  totalMassOfMovingZMols_local = 0;
584 >
585 >  for(int i =0; i < unconsMols.size(); i++){
586 >    unconsMols[i]->getCOMvel(COMvel);
587 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
588 >  }
589 >
590 >  for(int i = 0; i < zconsMols.size(); i++){
591 >    if (states[i] == zcsMoving){
592 >      zconsMols[i]->getCOMvel(COMvel);
593 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
594 >      totalMassOfMovingZMols_local += massOfZConsMols[i];              
595 >    }
596 >                
597 >  }
598 >
599 > #ifndef IS_MPI
600 >  MVzOfMovingMols = MVzOfMovingMols_local;
601 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
602   #else
603 <  totalMZOfZCons = mzOfZCons;
603 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
604 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
605   #endif
606  
607 +  double vzOfMovingMols;
608 +  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
609 +
610 +  //modify the velocites of unconstrained molecules  
611 +  Atom** unconsAtoms;
612    for(int i = 0; i < unconsMols.size(); i++){
517 –    unconsMols[i]->getCOM(COM);
518 –    mzOfUncons += massOfUnconsMols[i] * COM[2];
519 –  }
613    
614 < #ifdef IS_MPI
615 <  MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
616 < #else
617 <  totalMZOfUncons = mzOfUncons;
618 < #endif  
614 >    unconsAtoms = unconsMols[i]->getMyAtoms();
615 >    for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
616 >      unconsAtoms[j]->getVel(vel);
617 >      vel[whichDirection] -= vzOfMovingMols;
618 >      unconsAtoms[j]->setVel(vel);
619 >    }
620    
621 <  zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons;
621 >  }  
622  
623 <  cout << "current time: " << info->getTime() <<endl;  
624 <  for(int i = 0; i < zconsMols.size(); i++){  
623 >  //modify the velocities of moving z-constrained molecuels
624 >  Atom** movingZAtoms;
625 >  for(int i = 0; i < zconsMols.size(); i++){
626 >
627 >    if (states[i] ==zcsMoving){
628    
629 <    zconsMols[i]->getCOM(COM);
629 >      movingZAtoms = zconsMols[i]->getMyAtoms();
630 >           for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
631 >        movingZAtoms[j]->getVel(vel);
632 >        vel[whichDirection] -= vzOfMovingMols;
633 >             movingZAtoms[j]->setVel(vel);
634 >          }
635 >          
636 >   }
637 >
638 > }
639 >
640 >        cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() <<endl;
641 >
642 > }
643 >
644 > template<typename T> void ZConstraint<T>::doZconstraintForce(){
645 >
646 >  Atom** zconsAtoms;
647 >  double totalFZ;
648 >  double totalFZ_local;
649 >  double COMvel[3];  
650 >  double COM[3];
651 >  double force[3];
652 >
653 >
654 >
655 >  //constrain the molecules which do not reach the specified positions  
656      
657 <    cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t";
658 <    deltaZ = zsys + refZ[i] - COM[2];
659 <    cout << "\tdistance: " << COM[2] +deltaZ - zsys;    
660 <    //update z coordinate    
661 <    zconsAtoms = zconsMols[i]->getMyAtoms();    
662 <    for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
663 <      zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ);  
664 <    }    
657 >  //Zero Out the force of z-contrained molecules    
658 >  totalFZ_local = 0;
659 >
660 >  //calculate the total z-contrained force of fixed z-contrained molecules
661 >  cout << "Fixed Molecules" << endl;
662 >  for(int i = 0; i < zconsMols.size(); i++){
663 >                
664 >    if (states[i] == zcsFixed){
665 >                
666 >      zconsMols[i]->getCOM(COM);
667 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
668 >
669 >      fz[i] = 0;      
670 >      for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
671 >        zconsAtoms[j]->getFrc(force);
672 >        fz[i] += force[whichDirection];      
673 >      }
674 >      totalFZ_local += fz[i];
675 >
676 >      cout << "index: " << indexOfZConsMols[i]
677 >                                <<"\tcurrent zpos: " << COM[whichDirection]
678 >                                << "\tcurrent fz: " <<fz[i] << endl;
679 >
680 >    }
681 >          
682 >  }
683 >        
684 >  // apply negative to fixed z-constrained molecues;
685 >  force[0]= 0;
686 >  force[1]= 0;
687 >  force[2]= 0;
688 >
689 >  for(int i = 0; i < zconsMols.size(); i++){
690 >
691 >    if (states[i] == zcsFixed){  
692 >        
693 >      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
694 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
695      
696 <    //calculate z constrain force
697 <    fz[i] = massOfZConsMols[i]* deltaZ / dt2;
698 <    
699 <    cout << "\tforce: " << fz[i] << endl;
696 >      for(int j =0; j < nAtomOfCurZConsMol; j++) {
697 >        force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
698 >        zconsAtoms[j]->addFrc(force);
699 >      }
700 >                
701 >    }
702 >        
703 >  }
704 >
705 >  cout << "after zero out z-constraint force on fixed z-constraint molecuels "
706 >                  << "total force is " << calcTotalForce() << endl;
707 >  //calculate the number of atoms of moving z-constrained molecules
708 >  int nMovingZAtoms_local;
709 >  int nMovingZAtoms;
710 >        
711 >  nMovingZAtoms_local = 0;
712 >  for(int i = 0; i < zconsMols.size(); i++)
713 >    if(states[i] == zcsMoving)
714 >           nMovingZAtoms_local += zconsMols[i]->getNAtoms();
715 >  
716 > #ifdef IS_MPI
717 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
718 >  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
719 > #else
720 >  totalFZ = totalFZ_local;
721 >  nMovingZAtoms = nMovingZAtoms_local;
722 > #endif
723 >
724 >  force[0]= 0;
725 >  force[1]= 0;
726 >  force[2]= 0;
727 >  force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
728 >
729 >  //modify the forces of unconstrained molecules
730 >  int accessCount = 0;
731 >  for(int i = 0; i < unconsMols.size(); i++){
732 >    
733 >     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
734 >    
735 >     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)          
736 >       unconsAtoms[j]->addFrc(force);
737 >    
738 >  }      
739 >
740 > //modify the forces of moving z-constrained molecules
741 >  for(int i = 0; i < zconsMols.size(); i++) {
742 >    if (states[i] == zcsMoving){
743 >                
744 >      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
745 >
746 >      for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)
747 >        movingZAtoms[j]->addFrc(force);
748 >    }
749    }
750 +        
751 +  cout << "after substracting z-constraint force from moving molecuels "
752 +                  << "total force is " << calcTotalForce()  << endl;
753  
754 + }
755 +
756 + template<typename T> bool ZConstraint<T>::checkZConsState(){
757 +  double COM[3];
758 +  double diff;
759 +  
760 +  bool changed;
761 +  
762 +  changed = false;
763 +  
764 +  for(int i =0; i < zconsMols.size(); i++){
765 +
766 +    zconsMols[i]->getCOM(COM);
767 +    diff = fabs(COM[whichDirection] - zPos[i]);  
768 +    if (  diff <= zconsTol && states[i] == zcsMoving){
769 +      states[i] = zcsFixed;
770 +        changed = true;
771 +    }
772 +    else if ( diff > zconsTol && states[i] == zcsFixed){
773 +      states[i] = zcsMoving;
774 +        changed = true;  
775 +    }
776 +  
777 +  }
778 +
779 +  return changed;
780 + }
781 +
782 + template<typename T> bool ZConstraint<T>::haveFixedZMols(){
783 +  for(int i = 0; i < zconsMols.size(); i++)
784 +    if (states[i] == zcsFixed)
785 +      return true;
786 +
787 +  return false;
788 + }
789 +
790 +
791 + /**
792 + *
793 + */
794 + template<typename T> bool ZConstraint<T>::haveMovingZMols(){
795 +  for(int i = 0; i < zconsMols.size(); i++)
796 +    if (states[i] == zcsMoving)
797 +      return true;
798 +
799 +  return false;
800 +  
801 + }
802 +
803 + /**
804 +  *
805 +  *
806 +  */
807 +
808 + template<typename T> void ZConstraint<T>::doHarmonic(){
809 +  double force[3];
810 +  double harmonicU;
811 +  double harmonicF;
812 +  double COM[3];
813 +  double diff;
814 +  double totalFZ;
815 +        
816 +  force[0] = 0;
817 +  force[1] = 0;
818 +  force[2] = 0;
819 +
820 +  totalFZ = 0;
821 +
822 +  cout << "Moving Molecules" << endl;  
823 +  for(int i = 0; i < zconsMols.size(); i++) {
824 +
825 +    if (states[i] == zcsMoving){
826 +      zconsMols[i]->getCOM(COM);
827 +      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
828 +                
829 +                diff = COM[whichDirection] -zPos[i];
830 +                
831 +      harmonicU = 0.5 * kz[i] * diff * diff;  
832 +                info->lrPot += harmonicU;
833 +
834 +      harmonicF =  - kz[i] * diff;
835 +      totalFZ += harmonicF;
836 +
837 +       //
838 +                force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
839 +                
840 +      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
841 +
842 +       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
843 +         movingZAtoms[j]->addFrc(force);
844 +    }
845 +
846 +  }
847 +        
848 +  force[0]= 0;
849 +  force[1]= 0;
850 +  force[2]= 0;
851 +  force[whichDirection] = -totalFZ /totNumOfUnconsAtoms;
852 +
853 +  //modify the forces of unconstrained molecules
854 +  for(int i = 0; i < unconsMols.size(); i++){
855 +    
856 +     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
857 +    
858 +     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)          
859 +       unconsAtoms[j]->addFrc(force);    
860 +  }  
861 +
862 + }
863 +
864 + template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
865 + {
866 +  double MVzOfMovingMols_local;
867 +  double MVzOfMovingMols;
868 +  double totalMassOfMovingZMols_local;
869 +  double totalMassOfMovingZMols;
870 +  double COMvel[3];
871        
872 +  MVzOfMovingMols_local = 0;
873 +  totalMassOfMovingZMols_local = 0;
874 +
875 +  for(int i =0; i < unconsMols.size(); i++){
876 +    unconsMols[i]->getCOMvel(COMvel);
877 +    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
878 +  }
879 +
880 +  for(int i = 0; i < zconsMols.size(); i++){
881 +
882 +    if (states[i] == zcsMoving){
883 +      zconsMols[i]->getCOMvel(COMvel);
884 +      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
885 +      totalMassOfMovingZMols_local += massOfZConsMols[i];              
886 +    }
887 +                
888 +  }
889 +
890 + #ifndef IS_MPI
891 +  MVzOfMovingMols = MVzOfMovingMols_local;
892 +  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
893 + #else
894 +  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
895 +  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
896 + #endif
897 +
898 +  double vzOfMovingMols;
899 +  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
900 +
901 +  return vzOfMovingMols;
902   }
903 +
904 +
905 + template<typename T> double ZConstraint<T>::calcSysCOMVel()
906 + {
907 +  double COMvel[3];
908 +  double tempMVz = 0;
909 +                
910 +  for(int i =0 ; i < nMols; i++){
911 +    molecules[i].getCOMvel(COMvel);
912 +         tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection];
913 +  }
914 +
915 +  double massOfZCons_local;
916 +  double massOfZCons;
917 +
918 +  massOfZCons_local = 0;
919 +        
920 +  for(int i = 0; i < massOfZConsMols.size(); i++){
921 +    massOfZCons_local += massOfZConsMols[i];
922 +  }
923 + #ifndef IS_MPI
924 +  massOfZCons = massOfZCons_local;
925 + #else
926 +  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
927 + #endif
928 +
929 +  return tempMVz /(totalMassOfUncons + massOfZCons);
930 + }
931 +
932 + template<typename T> double ZConstraint<T>::calcTotalForce(){
933 +
934 +  double force[3];  
935 +  double totalForce_local;
936 +  double totalForce;
937 +
938 +  totalForce_local = 0;
939 +
940 +  for(int i = 0; i < nAtoms; i++){
941 +    atoms[i]->getFrc(force);
942 +    totalForce_local += force[whichDirection];
943 +  }
944 +
945 + #ifndef IS_MPI
946 +  totalForce = totalForce_local;
947 + #else
948 +  MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
949 + #endif
950 +
951 +  return totalForce;
952 +
953 + }

Diff Legend

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