# | 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 47 | Line 57 | SimInfo::SimInfo(){ | |
57 | wrapMeSimInfo( this ); | |
58 | } | |
59 | ||
60 | + | SimInfo::~SimInfo(){ |
61 | + | |
62 | + | map<string, GenericData*>::iterator i; |
63 | + | |
64 | + | for(i = properties.begin(); i != properties.end(); i++) |
65 | + | delete (*i).second; |
66 | + | |
67 | + | |
68 | + | } |
69 | + | |
70 | void SimInfo::setBox(double newBox[3]) { | |
71 | ||
72 | int i, j; | |
# | Line 74 | Line 94 | void SimInfo::setBoxM( double theBox[3][3] ){ | |
94 | // [ 2 5 8 ] | |
95 | double FortranHmatInv[9]; // the inverted Hmat (for Fortran); | |
96 | ||
97 | + | |
98 | + | if( !boxIsInit ) boxIsInit = 1; |
99 | ||
100 | for(i=0; i < 3; i++) | |
101 | for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j]; | |
102 | ||
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 | – | |
103 | calcBoxL(); | |
104 | calcHmatInv(); | |
105 | ||
# | Line 96 | Line 112 | void SimInfo::setBoxM( double theBox[3][3] ){ | |
112 | ||
113 | setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic); | |
114 | ||
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 | – | |
115 | } | |
116 | ||
117 | ||
# | Line 312 | Line 286 | void SimInfo::calcBoxL( void ){ | |
286 | ||
287 | dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0]; | |
288 | dsq = dx*dx + dy*dy + dz*dz; | |
289 | < | boxLx = sqrt( dsq ); |
289 | > | boxL[0] = sqrt( dsq ); |
290 | > | maxCutoff = 0.5 * boxL[0]; |
291 | ||
292 | // boxLy | |
293 | ||
294 | dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1]; | |
295 | dsq = dx*dx + dy*dy + dz*dz; | |
296 | < | boxLy = sqrt( dsq ); |
296 | > | boxL[1] = sqrt( dsq ); |
297 | > | if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1]; |
298 | ||
299 | // boxLz | |
300 | ||
301 | dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2]; | |
302 | dsq = dx*dx + dy*dy + dz*dz; | |
303 | < | boxLz = sqrt( dsq ); |
304 | < | |
303 | > | boxL[2] = sqrt( dsq ); |
304 | > | if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2]; |
305 | > | |
306 | } | |
307 | ||
308 | ||
# | Line 405 | Line 382 | void SimInfo::refreshSim(){ | |
382 | int isError; | |
383 | int n_global; | |
384 | int* excl; | |
385 | < | |
409 | < | fInfo.rrf = 0.0; |
410 | < | fInfo.rt = 0.0; |
385 | > | |
386 | fInfo.dielect = 0.0; | |
387 | ||
413 | – | fInfo.rlist = rList; |
414 | – | fInfo.rcut = rCut; |
415 | – | |
388 | if( useDipole ){ | |
417 | – | fInfo.rrf = ecr; |
418 | – | fInfo.rt = ecr - est; |
389 | if( useReactionField )fInfo.dielect = dielectric; | |
390 | } | |
391 | ||
# | Line 461 | Line 431 | void SimInfo::refreshSim(){ | |
431 | ||
432 | this->ndf = this->getNDF(); | |
433 | this->ndfRaw = this->getNDFraw(); | |
434 | + | |
435 | + | } |
436 | + | |
437 | + | |
438 | + | void SimInfo::setRcut( double theRcut ){ |
439 | + | |
440 | + | if( !haveOrigRcut ){ |
441 | + | haveOrigRcut = 1; |
442 | + | origRcut = theRcut; |
443 | + | } |
444 | + | |
445 | + | rCut = theRcut; |
446 | + | checkCutOffs(); |
447 | + | } |
448 | + | |
449 | + | void SimInfo::setEcr( double theEcr ){ |
450 | + | |
451 | + | if( !haveOrigEcr ){ |
452 | + | haveOrigEcr = 1; |
453 | + | origEcr = theEcr; |
454 | + | } |
455 | + | |
456 | + | ecr = theEcr; |
457 | + | checkCutOffs(); |
458 | + | } |
459 | + | |
460 | + | void SimInfo::setEcr( double theEcr, double theEst ){ |
461 | + | |
462 | + | est = theEst; |
463 | + | setEcr( theEcr ); |
464 | + | } |
465 | + | |
466 | + | |
467 | + | void SimInfo::checkCutOffs( void ){ |
468 | ||
469 | + | int cutChanged = 0; |
470 | + | |
471 | + | if( boxIsInit ){ |
472 | + | |
473 | + | //we need to check cutOffs against the box |
474 | + | |
475 | + | if( maxCutoff > rCut ){ |
476 | + | if( rCut < origRcut ){ |
477 | + | rCut = origRcut; |
478 | + | if (rCut > maxCutoff) rCut = maxCutoff; |
479 | + | |
480 | + | sprintf( painCave.errMsg, |
481 | + | "New Box size is setting the long range cutoff radius " |
482 | + | "to %lf\n", |
483 | + | rCut ); |
484 | + | painCave.isFatal = 0; |
485 | + | simError(); |
486 | + | } |
487 | + | } |
488 | + | |
489 | + | if( maxCutoff > ecr ){ |
490 | + | if( ecr < origEcr ){ |
491 | + | rCut = origEcr; |
492 | + | if (ecr > maxCutoff) ecr = maxCutoff; |
493 | + | |
494 | + | sprintf( painCave.errMsg, |
495 | + | "New Box size is setting the electrostaticCutoffRadius " |
496 | + | "to %lf\n", |
497 | + | ecr ); |
498 | + | painCave.isFatal = 0; |
499 | + | simError(); |
500 | + | } |
501 | + | } |
502 | + | |
503 | + | |
504 | + | if (rCut > maxCutoff) { |
505 | + | sprintf( painCave.errMsg, |
506 | + | "New Box size is setting the long range cutoff radius " |
507 | + | "to %lf\n", |
508 | + | maxCutoff ); |
509 | + | painCave.isFatal = 0; |
510 | + | simError(); |
511 | + | rCut = maxCutoff; |
512 | + | } |
513 | + | |
514 | + | if( ecr > maxCutoff){ |
515 | + | sprintf( painCave.errMsg, |
516 | + | "New Box size is setting the electrostaticCutoffRadius " |
517 | + | "to %lf\n", |
518 | + | maxCutoff ); |
519 | + | painCave.isFatal = 0; |
520 | + | simError(); |
521 | + | ecr = maxCutoff; |
522 | + | } |
523 | + | |
524 | + | |
525 | + | } |
526 | + | |
527 | + | |
528 | + | if( (oldEcr != ecr) || ( oldRcut != rCut ) ) cutChanged = 1; |
529 | + | |
530 | + | // rlist is the 1.0 plus max( rcut, ecr ) |
531 | + | |
532 | + | ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; |
533 | + | |
534 | + | if( cutChanged ){ |
535 | + | |
536 | + | notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); |
537 | + | } |
538 | + | |
539 | + | oldEcr = ecr; |
540 | + | oldRcut = rCut; |
541 | } | |
542 | ||
543 | + | void SimInfo::addProperty(GenericData* prop){ |
544 | + | |
545 | + | map<string, GenericData*>::iterator result; |
546 | + | result = properties.find(prop->getID()); |
547 | + | |
548 | + | //we can't simply use properties[prop->getID()] = prop, |
549 | + | //it will cause memory leak if we already contain a propery which has the same name of prop |
550 | + | |
551 | + | if(result != properties.end()){ |
552 | + | |
553 | + | delete (*result).second; |
554 | + | (*result).second = prop; |
555 | + | |
556 | + | } |
557 | + | else{ |
558 | + | |
559 | + | properties[prop->getID()] = prop; |
560 | + | |
561 | + | } |
562 | + | |
563 | + | } |
564 | + | |
565 | + | GenericData* SimInfo::getProperty(const string& propName){ |
566 | + | |
567 | + | map<string, GenericData*>::iterator result; |
568 | + | |
569 | + | //string lowerCaseName = (); |
570 | + | |
571 | + | result = properties.find(propName); |
572 | + | |
573 | + | if(result != properties.end()) |
574 | + | return (*result).second; |
575 | + | else |
576 | + | return NULL; |
577 | + | } |
578 | + | |
579 | + | vector<GenericData*> SimInfo::getProperties(){ |
580 | + | |
581 | + | vector<GenericData*> result; |
582 | + | map<string, GenericData*>::iterator i; |
583 | + | |
584 | + | for(i = properties.begin(); i != properties.end(); i++) |
585 | + | result.push_back((*i).second); |
586 | + | |
587 | + | return result; |
588 | + | } |
589 | + | |
590 | + |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |