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 619 by mmeineke, Tue Jul 15 22:22:41 2003 UTC vs.
Revision 767 by tim, Tue Sep 16 20:02:11 2003 UTC

# Line 26 | Line 26 | SimInfo::SimInfo(){
26   SimInfo::SimInfo(){
27    excludes = NULL;
28    n_constraints = 0;
29 +  nZconstraints = 0;
30    n_oriented = 0;
31    n_dipoles = 0;
32    ndf = 0;
33    ndfRaw = 0;
34 +  nZconstraints = 0;
35    the_integrator = NULL;
36    setTemp = 0;
37    thermalTime = 0.0;
38 +  currentTime = 0.0;
39    rCut = 0.0;
40 +  origRcut = -1.0;
41    ecr = 0.0;
42 +  origEcr = -1.0;
43    est = 0.0;
44 +  oldEcr = 0.0;
45 +  oldRcut = 0.0;
46  
47 +  haveOrigRcut = 0;
48 +  haveOrigEcr = 0;
49 +  boxIsInit = 0;
50 +  
51 +  
52 +
53    usePBC = 0;
54    useLJ = 0;
55    useSticky = 0;
# Line 45 | Line 58 | SimInfo::SimInfo(){
58    useGB = 0;
59    useEAM = 0;
60  
61 +  myConfiguration = new SimState();
62 +
63    wrapMeSimInfo( this );
64   }
65  
66 +
67 + SimInfo::~SimInfo(){
68 +
69 +  delete myConfiguration;
70 +
71 +  map<string, GenericData*>::iterator i;
72 +  
73 +  for(i = properties.begin(); i != properties.end(); i++)
74 +    delete (*i).second;
75 +    
76 + }
77 +
78   void SimInfo::setBox(double newBox[3]) {
79    
80    int i, j;
# Line 75 | Line 102 | void SimInfo::setBoxM( double theBox[3][3] ){
102                           // [ 2 5 8 ]
103    double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
104  
105 +  
106 +  if( !boxIsInit ) boxIsInit = 1;
107  
108    for(i=0; i < 3; i++)
109      for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
110    
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 –
111    calcBoxL();
112    calcHmatInv();
113  
# Line 97 | Line 120 | void SimInfo::setBoxM( double theBox[3][3] ){
120  
121    setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
122  
100 –  smallestBoxL = boxLx;
101 –  if (boxLy < smallestBoxL) smallestBoxL = boxLy;
102 –  if (boxLz < smallestBoxL) smallestBoxL = boxLz;
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 –
113 –    rList = maxCutoff;
114 –
115 –    sprintf( painCave.errMsg,
116 –             "New Box size is forcing cutoff radius down to %lf\n",
117 –             maxCutoff - 1.0 );
118 –    painCave.isFatal = 0;
119 –    simError();
120 –
121 –    rCut = rList - 1.0;
122 –
123 –    // list radius changed so we have to refresh the simulation structure.
124 –    refreshSim();
125 –  }
126 –
127 –  if( ecr > maxCutoff ){
128 –
129 –    sprintf( painCave.errMsg,
130 –             "New Box size is forcing electrostatic cutoff radius "
131 –             "down to %lf\n",
132 –             maxCutoff );
133 –    painCave.isFatal = 0;
134 –    simError();
135 –
136 –    ecr = maxCutoff;
137 –    est = 0.05 * ecr;
138 –
139 –    refreshSim();
140 –  }
141 –    
123   }
124  
125  
# Line 313 | Line 294 | void SimInfo::calcBoxL( void ){
294    
295    dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
296    dsq = dx*dx + dy*dy + dz*dz;
297 <  boxLx = sqrt( dsq );
297 >  boxL[0] = sqrt( dsq );
298 >  maxCutoff = 0.5 * boxL[0];
299  
300    // boxLy
301    
302    dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
303    dsq = dx*dx + dy*dy + dz*dz;
304 <  boxLy = sqrt( dsq );
304 >  boxL[1] = sqrt( dsq );
305 >  if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
306  
307    // boxLz
308    
309    dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
310    dsq = dx*dx + dy*dy + dz*dz;
311 <  boxLz = sqrt( dsq );
311 >  boxL[2] = sqrt( dsq );
312 >  if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
313    
314 +  checkCutOffs();
315 +
316   }
317  
318  
# Line 380 | Line 366 | int SimInfo::getNDF(){
366    ndf = ndf_local;
367   #endif
368  
369 <  ndf = ndf - 3;
369 >  ndf = ndf - 3 - nZconstraints;
370  
371    return ndf;
372   }
# Line 399 | Line 385 | int SimInfo::getNDFraw() {
385  
386    return ndfRaw;
387   }
388 <
388 >
389 > int SimInfo::getNDFtranslational() {
390 >  int ndfTrans_local, ndfTrans;
391 >
392 >  ndfTrans_local = 3 * n_atoms - n_constraints;
393 >
394 > #ifdef IS_MPI
395 >  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
396 > #else
397 >  ndfTrans = ndfTrans_local;
398 > #endif
399 >
400 >  ndfTrans = ndfTrans - 3 - nZconstraints;
401 >
402 >  return ndfTrans;
403 > }
404 >
405   void SimInfo::refreshSim(){
406  
407    simtype fInfo;
408    int isError;
409    int n_global;
410    int* excl;
411 <  
410 <  fInfo.rrf = 0.0;
411 <  fInfo.rt = 0.0;
411 >
412    fInfo.dielect = 0.0;
413  
414 –  fInfo.rlist = rList;
415 –  fInfo.rcut = rCut;
416 –
414    if( useDipole ){
418 –    fInfo.rrf = ecr;
419 –    fInfo.rt = ecr - est;
415      if( useReactionField )fInfo.dielect = dielectric;
416    }
417  
# Line 462 | Line 457 | void SimInfo::refreshSim(){
457  
458    this->ndf = this->getNDF();
459    this->ndfRaw = this->getNDFraw();
460 +  this->ndfTrans = this->getNDFtranslational();
461 + }
462 +
463 +
464 + void SimInfo::setRcut( double theRcut ){
465 +
466 +  if( !haveOrigRcut ){
467 +    haveOrigRcut = 1;
468 +    origRcut = theRcut;
469 +  }
470 +
471 +  rCut = theRcut;
472 +  checkCutOffs();
473 + }
474 +
475 + void SimInfo::setEcr( double theEcr ){
476 +
477 +  if( !haveOrigEcr ){
478 +    haveOrigEcr = 1;
479 +    origEcr = theEcr;
480 +  }
481 +
482 +  ecr = theEcr;
483 +  checkCutOffs();
484 + }
485 +
486 + void SimInfo::setEcr( double theEcr, double theEst ){
487  
488 +  est = theEst;
489 +  setEcr( theEcr );
490   }
491  
492 +
493 + void SimInfo::checkCutOffs( void ){
494 +
495 +  int cutChanged = 0;
496 +
497 +
498 +
499 +  if( boxIsInit ){
500 +    
501 +    //we need to check cutOffs against the box
502 +  
503 +    if(( maxCutoff > rCut )&&(usePBC)){
504 +      if( rCut < origRcut ){
505 +        rCut = origRcut;
506 +        if (rCut > maxCutoff) rCut = maxCutoff;
507 +        
508 +        sprintf( painCave.errMsg,
509 +                 "New Box size is setting the long range cutoff radius "
510 +                 "to %lf\n",
511 +                 rCut );
512 +        painCave.isFatal = 0;
513 +        simError();
514 +      }
515 +    }
516 +
517 +    if( maxCutoff > ecr ){
518 +      if( ecr < origEcr ){
519 +        ecr = origEcr;
520 +        if (ecr > maxCutoff) ecr = maxCutoff;
521 +        
522 +        sprintf( painCave.errMsg,
523 +                 "New Box size is setting the electrostaticCutoffRadius "
524 +                 "to %lf\n",
525 +                 ecr );
526 +        painCave.isFatal = 0;
527 +        simError();
528 +      }
529 +    }
530 +
531 +
532 +    if ((rCut > maxCutoff)&&(usePBC)) {
533 +      sprintf( painCave.errMsg,
534 +               "New Box size is setting the long range cutoff radius "
535 +               "to %lf\n",
536 +               maxCutoff );
537 +      painCave.isFatal = 0;
538 +      simError();
539 +      rCut = maxCutoff;
540 +    }
541 +
542 +    if( ecr > maxCutoff){
543 +      sprintf( painCave.errMsg,
544 +               "New Box size is setting the electrostaticCutoffRadius "
545 +               "to %lf\n",
546 +               maxCutoff  );
547 +      painCave.isFatal = 0;
548 +      simError();      
549 +      ecr = maxCutoff;
550 +    }
551 +
552 +
553 +    if( (oldEcr != ecr) || ( oldRcut != rCut ) ) cutChanged = 1;
554 +
555 +    // rlist is the 1.0 plus max( rcut, ecr )
556 +    
557 +    ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
558 +    
559 +    if( cutChanged ){
560 +      
561 +      notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
562 +    }
563 +    
564 +    oldEcr = ecr;
565 +    oldRcut = rCut;
566 +
567 +  } else {
568 +    // initialize this stuff before using it, OK?
569 +      sprintf( painCave.errMsg,
570 +               "Trying to check cutoffs without a box. Be smarter.\n" );
571 +      painCave.isFatal = 1;
572 +      simError();      
573 +  }
574 +
575 + }
576 +
577 + void SimInfo::addProperty(GenericData* prop){
578 +
579 +  map<string, GenericData*>::iterator result;
580 +  result = properties.find(prop->getID());
581 +  
582 +  //we can't simply use  properties[prop->getID()] = prop,
583 +  //it will cause memory leak if we already contain a propery which has the same name of prop
584 +  
585 +  if(result != properties.end()){
586 +    
587 +    delete (*result).second;
588 +    (*result).second = prop;
589 +      
590 +  }
591 +  else{
592 +
593 +    properties[prop->getID()] = prop;
594 +
595 +  }
596 +    
597 + }
598 +
599 + GenericData* SimInfo::getProperty(const string& propName){
600 +
601 +  map<string, GenericData*>::iterator result;
602 +  
603 +  //string lowerCaseName = ();
604 +  
605 +  result = properties.find(propName);
606 +  
607 +  if(result != properties.end())
608 +    return (*result).second;  
609 +  else  
610 +    return NULL;  
611 + }
612 +
613 + vector<GenericData*> SimInfo::getProperties(){
614 +
615 +  vector<GenericData*> result;
616 +  map<string, GenericData*>::iterator i;
617 +  
618 +  for(i = properties.begin(); i != properties.end(); i++)
619 +    result.push_back((*i).second);
620 +    
621 +  return result;
622 + }
623 +
624 + double SimInfo::matTrace3(double m[3][3]){
625 +  double trace;
626 +  trace = m[0][0] + m[1][1] + m[2][2];
627 +
628 +  return trace;
629 + }

Diff Legend

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