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

Comparing trunk/OOPSE/libmdtools/SimInfo.cpp (file contents):
Revision 621 by gezelter, Wed Jul 16 02:11:02 2003 UTC vs.
Revision 626 by mmeineke, Wed Jul 16 21:30:56 2003 UTC

# Line 36 | Line 36 | SimInfo::SimInfo(){
36    rCut = 0.0;
37    ecr = 0.0;
38    est = 0.0;
39 +  oldEcr = 0.0;
40 +  oldRcut = 0.0;
41  
42 +  haveOrigRcut = 0;
43 +  haveOrigEcr = 0;
44 +  boxIsInit = 0;
45 +  
46 +  
47 +
48    usePBC = 0;
49    useLJ = 0;
50    useSticky = 0;
# Line 75 | Line 83 | void SimInfo::setBoxM( double theBox[3][3] ){
83                           // [ 2 5 8 ]
84    double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
85  
86 +  
87 +  if( !boxIsInit ) boxIsInit = 1;
88  
89    for(i=0; i < 3; i++)
90      for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
91    
82 –  //  cerr
83 –  // << "setting Hmat ->\n"
84 –  // << "[ " << Hmat[0][0] << ", " << Hmat[0][1] << ", " << Hmat[0][2] << " ]\n"
85 –  // << "[ " << Hmat[1][0] << ", " << Hmat[1][1] << ", " << Hmat[1][2] << " ]\n"
86 –  // << "[ " << Hmat[2][0] << ", " << Hmat[2][1] << ", " << Hmat[2][2] << " ]\n";
87 –
92    calcBoxL();
93    calcHmatInv();
94  
# Line 97 | Line 101 | void SimInfo::setBoxM( double theBox[3][3] ){
101  
102    setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
103  
100 –  smallestBoxL = boxL[0];
101 –  if (boxL[1] < smallestBoxL) smallestBoxL = boxL[1];
102 –  if (boxL[2] > smallestBoxL) smallestBoxL = boxL[2];
103 –
104 –  maxCutoff = smallestBoxL / 2.0;
105 –
106 –  if (rList > maxCutoff) {
107 –    sprintf( painCave.errMsg,
108 –             "New Box size is forcing neighborlist radius down to %lf\n",
109 –             maxCutoff );
110 –    painCave.isFatal = 0;
111 –    simError();
112 –    rList = maxCutoff;
113 –
114 –    if (rCut > (rList - 1.0)) {
115 –      sprintf( painCave.errMsg,
116 –               "New Box size is forcing LJ cutoff radius down to %lf\n",
117 –               rList - 1.0 );
118 –      painCave.isFatal = 0;
119 –      simError();
120 –      rCut = rList - 1.0;
121 –    }
122 –
123 –    if( ecr > (rList - 1.0) ){
124 –      sprintf( painCave.errMsg,
125 –               "New Box size is forcing electrostaticCutoffRadius "
126 –               "down to %lf\n"
127 –               "electrostaticSkinThickness is now %lf\n",
128 –               rList - 1.0, 0.05*(rList-1.0) );
129 –      painCave.isFatal = 0;
130 –      simError();      
131 –      ecr = maxCutoff;
132 –      est = 0.05 * ecr;
133 –    }
134 –
135 –    // At least one of the radii changed, so we need a refresh:
136 –    refreshSim();
137 –  }    
104   }
105  
106  
# Line 310 | Line 276 | void SimInfo::calcBoxL( void ){
276    dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
277    dsq = dx*dx + dy*dy + dz*dz;
278    boxL[0] = sqrt( dsq );
279 +  maxCutoff = 0.5 * boxL[0];
280  
281    // boxLy
282    
283    dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
284    dsq = dx*dx + dy*dy + dz*dz;
285    boxL[1] = sqrt( dsq );
286 +  if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
287  
288    // boxLz
289    
290    dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
291    dsq = dx*dx + dy*dy + dz*dz;
292    boxL[2] = sqrt( dsq );
293 <  
293 >  if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
294 >
295   }
296  
297  
# Line 402 | Line 371 | void SimInfo::refreshSim(){
371    int isError;
372    int n_global;
373    int* excl;
405 –  
406 –  fInfo.rrf = 0.0;
407 –  fInfo.rt = 0.0;
408 –  fInfo.dielect = 0.0;
374  
375 <  fInfo.rlist = rList;
411 <  fInfo.rcut = rCut;
375 >  fInfo.dielect = 0.0;
376  
377    if( useDipole ){
414 –    fInfo.rrf = ecr;
415 –    fInfo.rt = ecr - est;
378      if( useReactionField )fInfo.dielect = dielectric;
379    }
380  
# Line 461 | Line 423 | void SimInfo::refreshSim(){
423  
424   }
425  
426 +
427 + void SimInfo::setRcut( double theRcut ){
428 +
429 +  if( !haveOrigRcut ){
430 +    haveOrigRcut = 1;
431 +    origRcut = theRcut;
432 +  }
433 +
434 +  rCut = theRcut;
435 +  checkCutOffs();
436 + }
437 +
438 + void SimInfo::setEcr( double theEcr ){
439 +
440 +  if( !haveOrigEcr ){
441 +    haveOrigEcr = 1;
442 +    origEcr = theEcr;
443 +  }
444 +
445 +  ecr = theEcr;
446 +  checkCutOffs();
447 + }
448 +
449 + void SimInfo::setEcr( double theEcr, double theEst ){
450 +
451 +  est = theEst;
452 +  setEcr( theEcr );
453 + }
454 +
455 +
456 + void SimInfo::checkCutOffs( void ){
457 +
458 +  int cutChanged = 0;
459 +
460 +  if( boxIsInit ){
461 +    
462 +    //we need to check cutOffs against the box
463 +    
464 +    if( maxCutoff > rCut ){
465 +      if( rCut < origRcut ){
466 +        rCut = origRcut;
467 +        if (rCut > maxCutoff) rCut = maxCutoff;
468 +        
469 +        sprintf( painCave.errMsg,
470 +                 "New Box size is setting the long range cutoff radius "
471 +                 "to %lf\n",
472 +                 rCut );
473 +        painCave.isFatal = 0;
474 +        simError();
475 +      }
476 +    }
477 +
478 +    if( maxCutoff > ecr ){
479 +      if( ecr < origEcr ){
480 +        rCut = origEcr;
481 +        if (ecr > maxCutoff) ecr = maxCutoff;
482 +        
483 +        sprintf( painCave.errMsg,
484 +                 "New Box size is setting the electrostaticCutoffRadius "
485 +                 "to %lf\n",
486 +                 ecr );
487 +        painCave.isFatal = 0;
488 +        simError();
489 +      }
490 +    }
491 +
492 +
493 +    if (rCut > maxCutoff) {
494 +      sprintf( painCave.errMsg,
495 +               "New Box size is setting the long range cutoff radius "
496 +               "to %lf\n",
497 +               maxCutoff );
498 +      painCave.isFatal = 0;
499 +      simError();
500 +      rCut = maxCutoff;
501 +    }
502 +
503 +    if( ecr > maxCutoff){
504 +      sprintf( painCave.errMsg,
505 +               "New Box size is setting the electrostaticCutoffRadius "
506 +               "to %lf\n",
507 +               maxCutoff  );
508 +      painCave.isFatal = 0;
509 +      simError();      
510 +      ecr = maxCutoff;
511 +    }
512 +
513 +    
514 +  }
515 +  
516 +
517 +  if( (oldEcr != ecr) || ( oldRcut != rCut ) ) cutChanged = 1;
518 +
519 +  // rlist is the 1.0 plus max( rcut, ecr )
520 +  
521 +  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
522 +
523 +  if( cutChanged ){
524 +    
525 +    notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
526 +  }
527 +
528 +  oldEcr = ecr;
529 +  oldRcut = rCut;
530 + }

Diff Legend

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