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 1093 by tim, Wed Mar 17 14:22:59 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   #endif
360  
361 +  
362    //get total masss of unconstraint molecules
363    double totalMassOfUncons_local;
364    totalMassOfUncons_local = 0;
365  
366 <  for (int i = 0; i < (int) (unconsMols.size()); i++)
366 >  for (size_t i = 0; i < unconsMols.size(); i++)
367      totalMassOfUncons_local += unconsMols[i]->getTotalMass();
368  
369   #ifndef IS_MPI
# Line 366 | Line 390 | template<typename T> ZConstraint<T>::~ZConstraint(){
390   }
391  
392   template<typename T> ZConstraint<T>::~ZConstraint(){
369 –  if (fz){
370 –    delete[] fz;
371 –  }
393  
373 –  if (curZPos){
374 –    delete[] curZPos;
375 –  }
376 –
377 –  if (indexOfZConsMols){
378 –    delete[] indexOfZConsMols;
379 –  }
380 –
394    if (fzOut){
395      delete fzOut;
396    }
# Line 401 | Line 414 | template<typename T> void ZConstraint<T>::update(){
414    massOfZConsMols.clear();
415    zPos.clear();
416    kz.clear();
417 +  cantPos.clear();
418 +  cantVel.clear();
419  
420    unconsMols.clear();
421    massOfUnconsMols.clear();
# Line 414 | Line 429 | template<typename T> void ZConstraint<T>::update(){
429        zconsMols.push_back(&molecules[i]);      
430        zPos.push_back((*parameters)[index].zPos);
431        kz.push_back((*parameters)[index].kRatio * zForceConst);
432 <      massOfZConsMols.push_back(molecules[i].getTotalMass());  
432 >      massOfZConsMols.push_back(molecules[i].getTotalMass());
433 >      
434 >      if(usingSMD)
435 >        cantVel.push_back((*parameters)[index].cantVel);
436  
419 –      molecules[i].getCOM(COM);
437      }
438      else{
439        unconsMols.push_back(&molecules[i]);
# Line 424 | Line 441 | template<typename T> void ZConstraint<T>::update(){
441      }
442    }
443  
444 <
445 <  //The reason to declare fz and indexOfZconsMols as pointer to array is
446 <  // that we want to make the MPI communication simple
447 <  if (fz){
448 <    delete[] fz;
449 <  }
433 <
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;
444 >  fz.resize(zconsMols.size());
445 >  curZPos.resize(zconsMols.size());
446 >  indexOfZConsMols.resize(zconsMols.size());  
447 >
448 >  for (size_t i = 0; i < zconsMols.size(); i++){
449 >    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
450    }
451 <
451 >    
452    //determine the states of z-constraint molecules
453    for (int i = 0; i < (int) (zconsMols.size()); i++){
454 +
455      zconsMols[i]->getCOM(COM);
456 +    
457      if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
458        states.push_back(zcsFixed);
459  
# Line 476 | Line 466 | template<typename T> void ZConstraint<T>::update(){
466        if (hasZConsGap)
467          endFixTime.push_back(INFINITE_TIME);
468      }
469 +
470 +    if(usingSMD)
471 +      cantPos.push_back(COM[whichDirection]);        
472    }
473    //
474    forcePolicy->update();
# Line 580 | Line 573 | template<typename T> void ZConstraint<T>::calcForce(in
573      this->doZconstraintForce();
574    }
575  
576 <  //use harmonical poteintial to move the molecules to the specified positions
576 >  //use external force to move the molecules to the specified positions
577    if (haveMovingZMols()){
578 <    this->doHarmonic();
578 >    if (usingSMD)
579 >      this->doHarmonic(cantPos);
580 >    else
581 >      this->doHarmonic(zPos);      
582    }
583  
584    //write out forces and current positions of z-constraint molecules
# Line 602 | Line 598 | template<typename T> void ZConstraint<T>::calcForce(in
598          }
599        }
600      }
601 <    fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz,
602 <                   curZPos, &zPos[0]);
601 >    fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
602 >                   &curZPos[0], &zPos[0]);
603      curZconsTime += zconsTime;
604    }
605  
# Line 887 | Line 883 | template<typename T> void ZConstraint<T>::doHarmonic()
883    *
884    */
885  
886 < template<typename T> void ZConstraint<T>::doHarmonic(){
886 > template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
887    double force[3];
888    double harmonicU;
889    double harmonicF;
# Line 908 | Line 904 | template<typename T> void ZConstraint<T>::doHarmonic()
904        //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
905        //     << "\tcurrent zpos: " << COM[whichDirection] << endl;
906  
907 <      diff = COM[whichDirection] - zPos[i];
907 >      diff = COM[whichDirection] - resPos[i];
908  
909        harmonicU = 0.5 * kz[i] * diff * diff;  
910        info->lrPot += harmonicU;
# Line 1249 | Line 1245 | template<typename T> void ZConstraint<T>::updateZPos()
1245  
1246   template<typename T> void ZConstraint<T>::updateZPos(){
1247    double curTime;
1248 <
1248 >  double COM[3];
1249 >  
1250    curTime = info->getTime();
1251  
1252    for (size_t i = 0; i < zconsMols.size(); i++){
1253 +
1254      if (states[i] == zcsFixed && curTime >= endFixTime[i]){
1255        zPos[i] += zconsGap;
1256 +
1257 +      if (usingSMD){
1258 +        zconsMols[i]->getCOM(COM);
1259 +        cantPos[i] = COM[whichDirection];
1260 +      }
1261 +      
1262      }
1263 +    
1264    }
1265 +  
1266   }
1267 +
1268 + template<typename T> void ZConstraint<T>::updateCantPos(){
1269 +  double curTime;
1270 +  double dt;
1271 +
1272 +  curTime = info->getTime();
1273 +  dt = info->dt;
1274 +
1275 +  for (size_t i = 0; i < zconsMols.size(); i++){
1276 +    if (states[i] == zcsMoving){
1277 +      cantPos[i] += cantVel[i] * dt;
1278 +    }
1279 +  }
1280 +
1281 + }

Diff Legend

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