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 1091 by tim, Tue Mar 16 19:22:56 2004 UTC vs.
Revision 1141 by tim, Wed Apr 28 23:09:32 2004 UTC

# Line 5 | Line 5 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
5   const double INFINITE_TIME = 10e30;
6   template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo,
7                                                   ForceFields* the_ff): T(theInfo, the_ff),
8 –                                                                       indexOfZConsMols(NULL),
9 –                                                                       fz(NULL),
10 –                                                                       curZPos(NULL),
8                                                                         fzOut(NULL),
9                                                                         curZconsTime(0),
10                                                                         forcePolicy(NULL),
11 +                                                                       usingSMD(false),
12                                                                         hasZConsGap(false){
13    //get properties from SimInfo
14    GenericData* data;
# Line 21 | Line 19 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
19    DoubleData* fixtime;
20    StringData* policy;
21    StringData* filename;
22 +  IntData* smdFlag;
23    double COM[3];
24  
25    //by default, the direction of constraint is z
# Line 150 | Line 149 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
149  
150    if (data){
151      gap = dynamic_cast<DoubleData*>(data);
153 –  }
152  
153 <  if (!gap){
154 <    sprintf(painCave.errMsg,
155 <            "ZConstraint error: Can not get property from SimInfo\n");
156 <    painCave.isFatal = 1;
157 <    simError();
153 >    if (!gap){
154 >      sprintf(painCave.errMsg,
155 >              "ZConstraint error: Can not get property from SimInfo\n");
156 >      painCave.isFatal = 1;
157 >      simError();
158 >    }
159 >    else{
160 >      this->hasZConsGap = true;
161 >      this->zconsGap = gap->getData();
162 >    }
163    }
161 –  else{
162 –    this->hasZConsGap = true;
163 –    this->zconsGap = gap->getData();
164 –  }
164  
165 +
166 +
167    data = info->getProperty(ZCONSFIXTIME_ID);
168  
169    if (data){
170      fixtime = dynamic_cast<DoubleData*>(data);
171 +    if (!fixtime){
172 +      sprintf(painCave.errMsg,
173 +              "ZConstraint error: Can not get zconsFixTime from SimInfo\n");
174 +      painCave.isFatal = 1;
175 +      simError();
176 +    }
177 +    else{
178 +      this->zconsFixTime = fixtime->getData();
179 +    }
180    }
181 <
182 <  if (!fixtime){
183 <    sprintf(painCave.errMsg,
184 <            "ZConstraint error: Can not get property from SimInfo\n");
185 <    painCave.isFatal = 1;
176 <    simError();
181 >  else if(hasZConsGap){
182 >      sprintf(painCave.errMsg,
183 >              "ZConstraint error: must set fixtime if already set zconsGap\n");
184 >      painCave.isFatal = 1;
185 >      simError();
186    }
187 <  else{
188 <    this->zconsFixTime = fixtime->getData();
187 >
188 >
189 >
190 >  data = info->getProperty(ZCONSUSINGSMD_ID);
191 >
192 >  if (data){
193 >    smdFlag = dynamic_cast<IntData*>(data);
194 >
195 >    if (!smdFlag){
196 >      sprintf(painCave.errMsg,
197 >              "ZConstraint error: Can not get property from SimInfo\n");
198 >      painCave.isFatal = 1;
199 >      simError();
200 >    }
201 >    else{
202 >      this->usingSMD= smdFlag->getData() ? true : false;
203 >    }
204 >
205    }
206  
207  
# Line 292 | Line 317 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
317        massOfZConsMols.push_back(molecules[i].getTotalMass());  
318  
319        zPos.push_back((*parameters)[searchResult].zPos);
295 –      //       cout << "index: "<< (*parameters)[searchResult].zconsIndex
296 –      //              <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
297 –
320        kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
321 <      molecules[i].getCOM(COM);
321 >      
322 >      if(usingSMD)
323 >        cantVel.push_back((*parameters)[searchResult].cantVel);
324 >
325      }
326      else{
327        unconsMols.push_back(&molecules[i]);
# Line 304 | Line 329 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
329      }
330    }
331  
332 <  fz = new double[zconsMols.size()];
333 <  curZPos = new double[zconsMols.size()];
334 <  indexOfZConsMols = new int [zconsMols.size()];
332 >  fz.resize(zconsMols.size());
333 >  curZPos.resize(zconsMols.size());
334 >  indexOfZConsMols.resize(zconsMols.size());  
335  
311 –  if (!fz || !curZPos || !indexOfZConsMols){
312 –    sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
313 –    painCave.isFatal = 1;
314 –    simError();
315 –  }
316 –
336    //determine the states of z-constraint molecules
337 <  for (int i = 0; i < (int) (zconsMols.size()); i++){
337 >  for (size_t i = 0; i < zconsMols.size(); i++){
338      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
339  
340      zconsMols[i]->getCOM(COM);
341 +    
342      if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
343        states.push_back(zcsFixed);
344  
# Line 331 | Line 351 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
351        if (hasZConsGap)
352          endFixTime.push_back(INFINITE_TIME);
353      }
354 +
355 +    if(usingSMD)
356 +      cantPos.push_back(COM[whichDirection]);    
357    }
358  
359 +  if(usingSMD)
360 +    prevCantPos = cantPos;
361   #endif
362  
363 +  
364    //get total masss of unconstraint molecules
365    double totalMassOfUncons_local;
366    totalMassOfUncons_local = 0;
367  
368 <  for (int i = 0; i < (int) (unconsMols.size()); i++)
368 >  for (size_t i = 0; i < unconsMols.size(); i++)
369      totalMassOfUncons_local += unconsMols[i]->getTotalMass();
370  
371   #ifndef IS_MPI
# Line 366 | Line 392 | template<typename T> ZConstraint<T>::~ZConstraint(){
392   }
393  
394   template<typename T> ZConstraint<T>::~ZConstraint(){
369 –  if (fz){
370 –    delete[] fz;
371 –  }
395  
373 –  if (curZPos){
374 –    delete[] curZPos;
375 –  }
376 –
377 –  if (indexOfZConsMols){
378 –    delete[] indexOfZConsMols;
379 –  }
380 –
396    if (fzOut){
397      delete fzOut;
398    }
# Line 401 | Line 416 | template<typename T> void ZConstraint<T>::update(){
416    massOfZConsMols.clear();
417    zPos.clear();
418    kz.clear();
419 +  cantPos.clear();
420 +  cantVel.clear();
421  
422    unconsMols.clear();
423    massOfUnconsMols.clear();
# Line 414 | Line 431 | template<typename T> void ZConstraint<T>::update(){
431        zconsMols.push_back(&molecules[i]);      
432        zPos.push_back((*parameters)[index].zPos);
433        kz.push_back((*parameters)[index].kRatio * zForceConst);
434 <      massOfZConsMols.push_back(molecules[i].getTotalMass());  
434 >      massOfZConsMols.push_back(molecules[i].getTotalMass());
435 >      
436 >      if(usingSMD)
437 >        cantVel.push_back((*parameters)[index].cantVel);
438  
419 –      molecules[i].getCOM(COM);
439      }
440      else{
441        unconsMols.push_back(&molecules[i]);
# Line 424 | Line 443 | template<typename T> void ZConstraint<T>::update(){
443      }
444    }
445  
446 <
447 <  //The reason to declare fz and indexOfZconsMols as pointer to array is
448 <  // that we want to make the MPI communication simple
449 <  if (fz){
450 <    delete[] fz;
446 >  fz.resize(zconsMols.size());
447 >  curZPos.resize(zconsMols.size());
448 >  indexOfZConsMols.resize(zconsMols.size());  
449 >
450 >  for (size_t i = 0; i < zconsMols.size(); i++){
451 >    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
452    }
453 <
434 <  if (curZPos){
435 <    delete[] curZPos;
436 <  }
437 <
438 <  if (indexOfZConsMols){
439 <    delete[] indexOfZConsMols;
440 <  }
441 <
442 <  if (zconsMols.size() > 0){
443 <    fz = new double[zconsMols.size()];
444 <    curZPos = new double[zconsMols.size()];
445 <    indexOfZConsMols = new int[zconsMols.size()];
446 <
447 <    if (!fz || !curZPos || !indexOfZConsMols){
448 <      sprintf(painCave.errMsg,
449 <              "Memory allocation failure in class Zconstraint\n");
450 <      painCave.isFatal = 1;
451 <      simError();
452 <    }
453 <
454 <    for (int i = 0; i < (int) (zconsMols.size()); i++){
455 <      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
456 <    }
457 <  }
458 <  else{
459 <    fz = NULL;
460 <    curZPos = NULL;
461 <    indexOfZConsMols = NULL;
462 <  }
463 <
453 >    
454    //determine the states of z-constraint molecules
455    for (int i = 0; i < (int) (zconsMols.size()); i++){
456 +
457      zconsMols[i]->getCOM(COM);
458 +    
459      if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
460        states.push_back(zcsFixed);
461  
# Line 476 | Line 468 | template<typename T> void ZConstraint<T>::update(){
468        if (hasZConsGap)
469          endFixTime.push_back(INFINITE_TIME);
470      }
471 +
472 +    if(usingSMD)
473 +      cantPos.push_back(COM[whichDirection]);        
474    }
475 +
476 +  if(usingSMD)
477 +  prevCantPos = cantPos;
478 +
479    //
480    forcePolicy->update();
481   }
# Line 580 | Line 579 | template<typename T> void ZConstraint<T>::calcForce(in
579      this->doZconstraintForce();
580    }
581  
582 <  //use harmonical poteintial to move the molecules to the specified positions
582 >  //use external force to move the molecules to the specified positions
583    if (haveMovingZMols()){
584 <    this->doHarmonic();
584 >    if (usingSMD)
585 >      this->doHarmonic(cantPos);
586 >    else
587 >      this->doHarmonic(zPos);      
588    }
589  
590    //write out forces and current positions of z-constraint molecules
# Line 602 | Line 604 | template<typename T> void ZConstraint<T>::calcForce(in
604          }
605        }
606      }
607 <    fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz,
608 <                   curZPos, &zPos[0]);
607 >    fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
608 >                   &curZPos[0], &zPos[0]);
609      curZconsTime += zconsTime;
610    }
611  
# Line 887 | Line 889 | template<typename T> void ZConstraint<T>::doZconstrain
889    *
890    */
891  
892 < template<typename T> void ZConstraint<T>::doHarmonic(){
892 > template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
893    double force[3];
894    double harmonicU;
895    double harmonicF;
# Line 908 | Line 910 | template<typename T> void ZConstraint<T>::doHarmonic()
910        //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
911        //     << "\tcurrent zpos: " << COM[whichDirection] << endl;
912  
913 <      diff = COM[whichDirection] - zPos[i];
913 >      diff = COM[whichDirection] - resPos[i];
914  
915        harmonicU = 0.5 * kz[i] * diff * diff;  
916        info->lrPot += harmonicU;
# Line 978 | Line 980 | template<typename T> bool ZConstraint<T>::checkZConsSt
980      if (diff <= zconsTol && states[i] == zcsMoving){
981        states[i] = zcsFixed;
982        changed_local = 1;
983 +
984 +      if(usingSMD)
985 +        prevCantPos = cantPos;
986  
987        if (hasZConsGap)
988          endFixTime[i] = info->getTime() + zconsFixTime;
# Line 986 | Line 991 | template<typename T> bool ZConstraint<T>::checkZConsSt
991        states[i] = zcsMoving;
992        changed_local = 1;  
993  
994 +      if(usingSMD)
995 +         cantPos = prevCantPos;
996 +      
997        if (hasZConsGap)
998          endFixTime[i] = INFINITE_TIME;
999      }
# Line 1249 | Line 1257 | template<typename T> void ZConstraint<T>::updateZPos()
1257  
1258   template<typename T> void ZConstraint<T>::updateZPos(){
1259    double curTime;
1260 <
1260 >  double COM[3];
1261 >  
1262    curTime = info->getTime();
1263  
1264    for (size_t i = 0; i < zconsMols.size(); i++){
1265 +
1266      if (states[i] == zcsFixed && curTime >= endFixTime[i]){
1267        zPos[i] += zconsGap;
1268 +
1269 +      if (usingSMD){
1270 +        zconsMols[i]->getCOM(COM);
1271 +        cantPos[i] = COM[whichDirection];
1272 +      }
1273 +      
1274      }
1275 +    
1276    }
1277 +  
1278   }
1279 +
1280 + template<typename T> void ZConstraint<T>::updateCantPos(){
1281 +  double curTime;
1282 +  double dt;
1283 +
1284 +  curTime = info->getTime();
1285 +  dt = info->dt;
1286 +
1287 +  for (size_t i = 0; i < zconsMols.size(); i++){
1288 +    if (states[i] == zcsMoving){
1289 +      cantPos[i] += cantVel[i] * dt;
1290 +    }
1291 +  }
1292 +
1293 + }

Diff Legend

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