ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/OOPSE/libmdtools/SimInfo.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/SimInfo.cpp (file contents):
Revision 618 by mmeineke, Tue Jul 15 21:34:56 2003 UTC vs.
Revision 670 by mmeineke, Thu Aug 7 21:47:18 2003 UTC

# Line 33 | Line 33 | SimInfo::SimInfo(){
33    the_integrator = NULL;
34    setTemp = 0;
35    thermalTime = 0.0;
36 +  currentTime = 0.0;
37    rCut = 0.0;
38    ecr = 0.0;
39 +  est = 0.0;
40 +  oldEcr = 0.0;
41 +  oldRcut = 0.0;
42  
43 +  haveOrigRcut = 0;
44 +  haveOrigEcr = 0;
45 +  boxIsInit = 0;
46 +  
47 +  
48 +
49    usePBC = 0;
50    useLJ = 0;
51    useSticky = 0;
# Line 44 | Line 54 | SimInfo::SimInfo(){
54    useGB = 0;
55    useEAM = 0;
56  
57 +  myConfiguration = new SimState();
58 +
59    wrapMeSimInfo( this );
60   }
61  
62 +
63 + SimInfo::~SimInfo(){
64 +
65 +  delete myConfiguration;
66 +
67 +  map<string, GenericData*>::iterator i;
68 +  
69 +  for(i = properties.begin(); i != properties.end(); i++)
70 +    delete (*i).second;
71 +    
72 + }
73 +
74   void SimInfo::setBox(double newBox[3]) {
75    
76    int i, j;
# Line 74 | Line 98 | void SimInfo::setBoxM( double theBox[3][3] ){
98                           // [ 2 5 8 ]
99    double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
100  
101 +  
102 +  if( !boxIsInit ) boxIsInit = 1;
103  
104    for(i=0; i < 3; i++)
105      for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
106    
81 –  //  cerr
82 –  // << "setting Hmat ->\n"
83 –  // << "[ " << Hmat[0][0] << ", " << Hmat[0][1] << ", " << Hmat[0][2] << " ]\n"
84 –  // << "[ " << Hmat[1][0] << ", " << Hmat[1][1] << ", " << Hmat[1][2] << " ]\n"
85 –  // << "[ " << Hmat[2][0] << ", " << Hmat[2][1] << ", " << Hmat[2][2] << " ]\n";
86 –
107    calcBoxL();
108    calcHmatInv();
109  
# Line 96 | Line 116 | void SimInfo::setBoxM( double theBox[3][3] ){
116  
117    setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
118  
99 –  smallestBoxL = boxLx;
100 –  if (boxLy < smallestBoxL) smallestBoxL = boxLy;
101 –  if (boxLz < smallestBoxL) smallestBoxL = boxLz;
102 –
103 –  maxCutoff = smallestBoxL / 2.0;
104 –
105 –  if (rList > maxCutoff) {
106 –    sprintf( painCave.errMsg,
107 –             "New Box size is forcing neighborlist radius down to %lf\n",
108 –             maxCutoff );
109 –    painCave.isFatal = 0;
110 –    simError();
111 –
112 –    rList = maxCutoff;
113 –
114 –    sprintf( painCave.errMsg,
115 –             "New Box size is forcing cutoff radius down to %lf\n",
116 –             maxCutoff - 1.0 );
117 –    painCave.isFatal = 0;
118 –    simError();
119 –
120 –    rCut = rList - 1.0;
121 –
122 –    // list radius changed so we have to refresh the simulation structure.
123 –    refreshSim();
124 –  }
125 –
126 –  if( ecr > maxCutoff ){
127 –
128 –    sprintf( painCave.errMsg,
129 –             "New Box size is forcing electrostatic cutoff radius "
130 –             "down to %lf\n",
131 –             maxCutoff );
132 –    painCave.isFatal = 0;
133 –    simError();
134 –
135 –    ecr = maxCutoff;
136 –    est = 0.05 * ecr;
137 –
138 –    refreshSim();
139 –  }
140 –    
119   }
120  
121  
# Line 312 | Line 290 | void SimInfo::calcBoxL( void ){
290    
291    dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
292    dsq = dx*dx + dy*dy + dz*dz;
293 <  boxLx = sqrt( dsq );
293 >  boxL[0] = sqrt( dsq );
294 >  maxCutoff = 0.5 * boxL[0];
295  
296    // boxLy
297    
298    dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
299    dsq = dx*dx + dy*dy + dz*dz;
300 <  boxLy = sqrt( dsq );
300 >  boxL[1] = sqrt( dsq );
301 >  if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
302  
303    // boxLz
304    
305    dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
306    dsq = dx*dx + dy*dy + dz*dz;
307 <  boxLz = sqrt( dsq );
307 >  boxL[2] = sqrt( dsq );
308 >  if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
309    
310 +  checkCutOffs();
311 +
312   }
313  
314  
# Line 405 | Line 388 | void SimInfo::refreshSim(){
388    int isError;
389    int n_global;
390    int* excl;
391 <  
409 <  fInfo.rrf = 0.0;
410 <  fInfo.rt = 0.0;
391 >
392    fInfo.dielect = 0.0;
393  
413 –  fInfo.rlist = rList;
414 –  fInfo.rcut = rCut;
415 –
394    if( useDipole ){
417 –    fInfo.rrf = ecr;
418 –    fInfo.rt = ecr - est;
395      if( useReactionField )fInfo.dielect = dielectric;
396    }
397  
# Line 461 | Line 437 | void SimInfo::refreshSim(){
437  
438    this->ndf = this->getNDF();
439    this->ndfRaw = this->getNDFraw();
440 +
441 + }
442 +
443 +
444 + void SimInfo::setRcut( double theRcut ){
445 +
446 +  if( !haveOrigRcut ){
447 +    haveOrigRcut = 1;
448 +    origRcut = theRcut;
449 +  }
450 +
451 +  rCut = theRcut;
452 +  checkCutOffs();
453 + }
454 +
455 + void SimInfo::setEcr( double theEcr ){
456 +
457 +  if( !haveOrigEcr ){
458 +    haveOrigEcr = 1;
459 +    origEcr = theEcr;
460 +  }
461 +
462 +  ecr = theEcr;
463 +  checkCutOffs();
464 + }
465 +
466 + void SimInfo::setEcr( double theEcr, double theEst ){
467 +
468 +  est = theEst;
469 +  setEcr( theEcr );
470 + }
471 +
472  
473 + void SimInfo::checkCutOffs( void ){
474 +
475 +  int cutChanged = 0;
476 +
477 +
478 +
479 +  if( boxIsInit ){
480 +    
481 +    //we need to check cutOffs against the box
482 +  
483 +    if(( maxCutoff > rCut )&&(usePBC)){
484 +      if( rCut < origRcut ){
485 +        rCut = origRcut;
486 +        if (rCut > maxCutoff) rCut = maxCutoff;
487 +        
488 +        sprintf( painCave.errMsg,
489 +                 "New Box size is setting the long range cutoff radius "
490 +                 "to %lf\n",
491 +                 rCut );
492 +        painCave.isFatal = 0;
493 +        simError();
494 +      }
495 +    }
496 +
497 +    if( maxCutoff > ecr ){
498 +      if( ecr < origEcr ){
499 +        rCut = origEcr;
500 +        if (ecr > maxCutoff) ecr = maxCutoff;
501 +        
502 +        sprintf( painCave.errMsg,
503 +                 "New Box size is setting the electrostaticCutoffRadius "
504 +                 "to %lf\n",
505 +                 ecr );
506 +        painCave.isFatal = 0;
507 +        simError();
508 +      }
509 +    }
510 +
511 +
512 +    if ((rCut > maxCutoff)&&(usePBC)) {
513 +      sprintf( painCave.errMsg,
514 +               "New Box size is setting the long range cutoff radius "
515 +               "to %lf\n",
516 +               maxCutoff );
517 +      painCave.isFatal = 0;
518 +      simError();
519 +      rCut = maxCutoff;
520 +    }
521 +
522 +    if( ecr > maxCutoff){
523 +      sprintf( painCave.errMsg,
524 +               "New Box size is setting the electrostaticCutoffRadius "
525 +               "to %lf\n",
526 +               maxCutoff  );
527 +      painCave.isFatal = 0;
528 +      simError();      
529 +      ecr = maxCutoff;
530 +    }
531 +
532 +    
533 +  }
534 +  
535 +
536 +  if( (oldEcr != ecr) || ( oldRcut != rCut ) ) cutChanged = 1;
537 +
538 +  // rlist is the 1.0 plus max( rcut, ecr )
539 +  
540 +  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
541 +
542 +  if( cutChanged ){
543 +    
544 +    notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
545 +  }
546 +
547 +  oldEcr = ecr;
548 +  oldRcut = rCut;
549   }
550  
551 + void SimInfo::addProperty(GenericData* prop){
552 +
553 +  map<string, GenericData*>::iterator result;
554 +  result = properties.find(prop->getID());
555 +  
556 +  //we can't simply use  properties[prop->getID()] = prop,
557 +  //it will cause memory leak if we already contain a propery which has the same name of prop
558 +  
559 +  if(result != properties.end()){
560 +    
561 +    delete (*result).second;
562 +    (*result).second = prop;
563 +      
564 +  }
565 +  else{
566 +
567 +    properties[prop->getID()] = prop;
568 +
569 +  }
570 +    
571 + }
572 +
573 + GenericData* SimInfo::getProperty(const string& propName){
574 +
575 +  map<string, GenericData*>::iterator result;
576 +  
577 +  //string lowerCaseName = ();
578 +  
579 +  result = properties.find(propName);
580 +  
581 +  if(result != properties.end())
582 +    return (*result).second;  
583 +  else  
584 +    return NULL;  
585 + }
586 +
587 + vector<GenericData*> SimInfo::getProperties(){
588 +
589 +  vector<GenericData*> result;
590 +  map<string, GenericData*>::iterator i;
591 +  
592 +  for(i = properties.begin(); i != properties.end(); i++)
593 +    result.push_back((*i).second);
594 +    
595 +  return result;
596 + }
597 +
598 +

Diff Legend

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