| 14 | 
  | 
 | 
| 15 | 
  | 
#include "MatVec3.h" | 
| 16 | 
  | 
 | 
| 17 | 
+ | 
#include "ConstraintManager.hpp" | 
| 18 | 
+ | 
 | 
| 19 | 
  | 
#ifdef IS_MPI | 
| 20 | 
  | 
#include "mpiSimulation.hpp" | 
| 21 | 
  | 
#endif | 
| 44 | 
  | 
  thermalTime = 0.0; | 
| 45 | 
  | 
  currentTime = 0.0; | 
| 46 | 
  | 
  rCut = 0.0; | 
| 47 | 
< | 
  ecr = 0.0; | 
| 46 | 
< | 
  est = 0.0; | 
| 47 | 
> | 
  rSw = 0.0; | 
| 48 | 
  | 
 | 
| 49 | 
  | 
  haveRcut = 0; | 
| 50 | 
< | 
  haveEcr = 0; | 
| 50 | 
> | 
  haveRsw = 0; | 
| 51 | 
  | 
  boxIsInit = 0; | 
| 52 | 
  | 
   | 
| 53 | 
  | 
  resetTime = 1e99; | 
| 64 | 
  | 
  useReactionField = 0; | 
| 65 | 
  | 
  useGB = 0; | 
| 66 | 
  | 
  useEAM = 0; | 
| 67 | 
+ | 
  useSolidThermInt = 0; | 
| 68 | 
+ | 
  useLiquidThermInt = 0; | 
| 69 | 
  | 
 | 
| 70 | 
+ | 
  haveCutoffGroups = false; | 
| 71 | 
+ | 
 | 
| 72 | 
  | 
  excludes = Exclude::Instance(); | 
| 73 | 
  | 
 | 
| 74 | 
  | 
  myConfiguration = new SimState(); | 
| 76 | 
  | 
  has_minimizer = false; | 
| 77 | 
  | 
  the_minimizer =NULL; | 
| 78 | 
  | 
 | 
| 79 | 
+ | 
  ngroup = 0; | 
| 80 | 
+ | 
 | 
| 81 | 
+ | 
  consMan = NULL; | 
| 82 | 
+ | 
   | 
| 83 | 
  | 
  wrapMeSimInfo( this ); | 
| 84 | 
  | 
} | 
| 85 | 
  | 
 | 
| 92 | 
  | 
   | 
| 93 | 
  | 
  for(i = properties.begin(); i != properties.end(); i++) | 
| 94 | 
  | 
    delete (*i).second; | 
| 95 | 
< | 
     | 
| 95 | 
> | 
 | 
| 96 | 
> | 
  if (!consMan) | 
| 97 | 
> | 
    delete consMan;   | 
| 98 | 
  | 
} | 
| 99 | 
  | 
 | 
| 100 | 
  | 
void SimInfo::setBox(double newBox[3]) { | 
| 197 | 
  | 
 | 
| 198 | 
  | 
  if( oldOrtho != orthoRhombic ){ | 
| 199 | 
  | 
     | 
| 200 | 
< | 
    if( orthoRhombic ){ | 
| 200 | 
> | 
    if( orthoRhombic ) { | 
| 201 | 
  | 
      sprintf( painCave.errMsg, | 
| 202 | 
  | 
               "OOPSE is switching from the default Non-Orthorhombic\n" | 
| 203 | 
  | 
               "\tto the faster Orthorhombic periodic boundary computations.\n" | 
| 205 | 
  | 
               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n" | 
| 206 | 
  | 
               "\tvariable ( currently set to %G ) smaller.\n", | 
| 207 | 
  | 
               orthoTolerance); | 
| 208 | 
+ | 
      painCave.severity = OOPSE_INFO; | 
| 209 | 
  | 
      simError(); | 
| 210 | 
  | 
    } | 
| 211 | 
  | 
    else { | 
| 217 | 
  | 
               "\tthe Orthorhombic computations, make the orthoBoxTolerance\n" | 
| 218 | 
  | 
               "\tvariable ( currently set to %G ) larger.\n", | 
| 219 | 
  | 
               orthoTolerance); | 
| 220 | 
+ | 
      painCave.severity = OOPSE_WARNING; | 
| 221 | 
  | 
      simError(); | 
| 222 | 
  | 
    } | 
| 223 | 
  | 
  } | 
| 362 | 
  | 
 | 
| 363 | 
  | 
  ndf = ndf - 3 - nZconstraints; | 
| 364 | 
  | 
 | 
| 352 | 
– | 
  std::cerr << "ndf = " << ndf; | 
| 353 | 
– | 
 | 
| 365 | 
  | 
  return ndf; | 
| 366 | 
  | 
} | 
| 367 | 
  | 
 | 
| 452 | 
  | 
 | 
| 453 | 
  | 
  n_exclude = excludes->getSize(); | 
| 454 | 
  | 
  excl = excludes->getFortranArray(); | 
| 455 | 
< | 
 | 
| 455 | 
> | 
   | 
| 456 | 
  | 
#ifdef IS_MPI | 
| 457 | 
< | 
  n_global = mpiSim->getTotAtoms(); | 
| 457 | 
> | 
  n_global = mpiSim->getNAtomsGlobal(); | 
| 458 | 
  | 
#else | 
| 459 | 
  | 
  n_global = n_atoms; | 
| 460 | 
  | 
#endif | 
| 461 | 
< | 
 | 
| 461 | 
> | 
   | 
| 462 | 
  | 
  isError = 0; | 
| 463 | 
< | 
 | 
| 463 | 
> | 
   | 
| 464 | 
> | 
  getFortranGroupArrays(this, FglobalGroupMembership, mfact); | 
| 465 | 
> | 
  //it may not be a good idea to pass the address of first element in vector | 
| 466 | 
> | 
  //since c++ standard does not require vector to be stored continuously in meomory | 
| 467 | 
> | 
  //Most of the compilers will organize the memory of vector continuously | 
| 468 | 
  | 
  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,  | 
| 469 | 
< | 
                  &nGlobalExcludes, globalExcludes, molMembershipArray,  | 
| 470 | 
< | 
                  &isError ); | 
| 469 | 
> | 
                  &nGlobalExcludes, globalExcludes, molMembershipArray,  | 
| 470 | 
> | 
                  &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError);  | 
| 471 | 
  | 
 | 
| 472 | 
  | 
  if( isError ){ | 
| 473 | 
< | 
 | 
| 474 | 
< | 
    sprintf( painCave.errMsg, | 
| 475 | 
< | 
             "There was an error setting the simulation information in fortran.\n" ); | 
| 473 | 
> | 
     | 
| 474 | 
> | 
    sprintf( painCave.errMsg, | 
| 475 | 
> | 
             "There was an error setting the simulation information in fortran.\n" ); | 
| 476 | 
  | 
    painCave.isFatal = 1; | 
| 477 | 
+ | 
    painCave.severity = OOPSE_ERROR; | 
| 478 | 
  | 
    simError(); | 
| 479 | 
  | 
  } | 
| 480 | 
< | 
 | 
| 480 | 
> | 
   | 
| 481 | 
  | 
#ifdef IS_MPI | 
| 482 | 
  | 
  sprintf( checkPointMsg, | 
| 483 | 
  | 
           "succesfully sent the simulation information to fortran.\n"); | 
| 484 | 
  | 
  MPIcheckPoint(); | 
| 485 | 
  | 
#endif // is_mpi | 
| 486 | 
< | 
 | 
| 486 | 
> | 
   | 
| 487 | 
  | 
  this->ndf = this->getNDF(); | 
| 488 | 
  | 
  this->ndfRaw = this->getNDFraw(); | 
| 489 | 
  | 
  this->ndfTrans = this->getNDFtranslational(); | 
| 490 | 
  | 
} | 
| 491 | 
  | 
 | 
| 492 | 
  | 
void SimInfo::setDefaultRcut( double theRcut ){ | 
| 493 | 
< | 
 | 
| 493 | 
> | 
   | 
| 494 | 
  | 
  haveRcut = 1; | 
| 495 | 
  | 
  rCut = theRcut; | 
| 496 | 
< | 
 | 
| 481 | 
< | 
  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; | 
| 482 | 
< | 
 | 
| 483 | 
< | 
  notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); | 
| 484 | 
< | 
} | 
| 485 | 
< | 
 | 
| 486 | 
< | 
void SimInfo::setDefaultEcr( double theEcr ){ | 
| 487 | 
< | 
 | 
| 488 | 
< | 
  haveEcr = 1; | 
| 489 | 
< | 
  ecr = theEcr; | 
| 496 | 
> | 
  rList = rCut + 1.0;  | 
| 497 | 
  | 
   | 
| 498 | 
< | 
  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; | 
| 492 | 
< | 
 | 
| 493 | 
< | 
  notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); | 
| 498 | 
> | 
  notifyFortranCutOffs( &rCut, &rSw, &rList ); | 
| 499 | 
  | 
} | 
| 500 | 
  | 
 | 
| 501 | 
< | 
void SimInfo::setDefaultEcr( double theEcr, double theEst ){ | 
| 501 | 
> | 
void SimInfo::setDefaultRcut( double theRcut, double theRsw ){ | 
| 502 | 
  | 
 | 
| 503 | 
< | 
  est = theEst; | 
| 504 | 
< | 
  setDefaultEcr( theEcr ); | 
| 503 | 
> | 
  rSw = theRsw; | 
| 504 | 
> | 
  setDefaultRcut( theRcut ); | 
| 505 | 
  | 
} | 
| 506 | 
  | 
 | 
| 507 | 
  | 
 | 
| 513 | 
  | 
     | 
| 514 | 
  | 
    if( rCut > maxCutoff ){ | 
| 515 | 
  | 
      sprintf( painCave.errMsg, | 
| 516 | 
< | 
               "LJrcut is too large for the current periodic box.\n" | 
| 517 | 
< | 
               "\tCurrent Value of LJrcut = %G at time %G\n " | 
| 516 | 
> | 
               "cutoffRadius is too large for the current periodic box.\n" | 
| 517 | 
> | 
               "\tCurrent Value of cutoffRadius = %G at time %G\n " | 
| 518 | 
  | 
               "\tThis is larger than half of at least one of the\n" | 
| 519 | 
  | 
               "\tperiodic box vectors.  Right now, the Box matrix is:\n" | 
| 520 | 
< | 
               "\n, %G" | 
| 520 | 
> | 
               "\n" | 
| 521 | 
  | 
               "\t[ %G %G %G ]\n" | 
| 522 | 
  | 
               "\t[ %G %G %G ]\n" | 
| 523 | 
  | 
               "\t[ %G %G %G ]\n", | 
| 524 | 
< | 
               rCut, currentTime, maxCutoff, | 
| 524 | 
> | 
               rCut, currentTime, | 
| 525 | 
  | 
               Hmat[0][0], Hmat[0][1], Hmat[0][2], | 
| 526 | 
  | 
               Hmat[1][0], Hmat[1][1], Hmat[1][2], | 
| 527 | 
  | 
               Hmat[2][0], Hmat[2][1], Hmat[2][2]); | 
| 528 | 
+ | 
      painCave.severity = OOPSE_ERROR; | 
| 529 | 
  | 
      painCave.isFatal = 1; | 
| 530 | 
  | 
      simError(); | 
| 531 | 
< | 
    } | 
| 526 | 
< | 
     | 
| 527 | 
< | 
    if( haveEcr ){ | 
| 528 | 
< | 
      if( ecr > maxCutoff ){ | 
| 529 | 
< | 
        sprintf( painCave.errMsg, | 
| 530 | 
< | 
                 "electrostaticCutoffRadius is too large for the current\n" | 
| 531 | 
< | 
                 "\tperiodic box.\n\n" | 
| 532 | 
< | 
                 "\tCurrent Value of ECR = %G at time %G\n " | 
| 533 | 
< | 
                 "\tThis is larger than half of at least one of the\n" | 
| 534 | 
< | 
                 "\tperiodic box vectors.  Right now, the Box matrix is:\n" | 
| 535 | 
< | 
                 "\n" | 
| 536 | 
< | 
                 "\t[ %G %G %G ]\n" | 
| 537 | 
< | 
                 "\t[ %G %G %G ]\n" | 
| 538 | 
< | 
                 "\t[ %G %G %G ]\n", | 
| 539 | 
< | 
                 ecr, currentTime, | 
| 540 | 
< | 
                 Hmat[0][0], Hmat[0][1], Hmat[0][2], | 
| 541 | 
< | 
                 Hmat[1][0], Hmat[1][1], Hmat[1][2], | 
| 542 | 
< | 
                 Hmat[2][0], Hmat[2][1], Hmat[2][2]); | 
| 543 | 
< | 
        painCave.isFatal = 1; | 
| 544 | 
< | 
        simError(); | 
| 545 | 
< | 
      } | 
| 546 | 
< | 
    } | 
| 531 | 
> | 
    }     | 
| 532 | 
  | 
  } else { | 
| 533 | 
  | 
    // initialize this stuff before using it, OK? | 
| 534 | 
  | 
    sprintf( painCave.errMsg, | 
| 535 | 
  | 
             "Trying to check cutoffs without a box.\n" | 
| 536 | 
  | 
             "\tOOPSE should have better programmers than that.\n" ); | 
| 537 | 
+ | 
    painCave.severity = OOPSE_ERROR; | 
| 538 | 
  | 
    painCave.isFatal = 1; | 
| 539 | 
  | 
    simError();       | 
| 540 | 
  | 
  } | 
| 577 | 
  | 
    return NULL;   | 
| 578 | 
  | 
} | 
| 579 | 
  | 
 | 
| 580 | 
+ | 
 | 
| 581 | 
+ | 
void SimInfo::getFortranGroupArrays(SimInfo* info,  | 
| 582 | 
+ | 
                                    vector<int>& FglobalGroupMembership, | 
| 583 | 
+ | 
                                    vector<double>& mfact){ | 
| 584 | 
+ | 
   | 
| 585 | 
+ | 
  Molecule* myMols; | 
| 586 | 
+ | 
  Atom** myAtoms; | 
| 587 | 
+ | 
  int numAtom; | 
| 588 | 
+ | 
  double mtot; | 
| 589 | 
+ | 
  int numMol; | 
| 590 | 
+ | 
  int numCutoffGroups; | 
| 591 | 
+ | 
  CutoffGroup* myCutoffGroup; | 
| 592 | 
+ | 
  vector<CutoffGroup*>::iterator iterCutoff; | 
| 593 | 
+ | 
  Atom* cutoffAtom; | 
| 594 | 
+ | 
  vector<Atom*>::iterator iterAtom; | 
| 595 | 
+ | 
  int atomIndex; | 
| 596 | 
+ | 
  double totalMass; | 
| 597 | 
+ | 
   | 
| 598 | 
+ | 
  mfact.clear(); | 
| 599 | 
+ | 
  FglobalGroupMembership.clear(); | 
| 600 | 
+ | 
   | 
| 601 | 
+ | 
 | 
| 602 | 
+ | 
  // Fix the silly fortran indexing problem | 
| 603 | 
+ | 
#ifdef IS_MPI | 
| 604 | 
+ | 
  numAtom = mpiSim->getNAtomsGlobal(); | 
| 605 | 
+ | 
#else | 
| 606 | 
+ | 
  numAtom = n_atoms; | 
| 607 | 
+ | 
#endif | 
| 608 | 
+ | 
  for (int i = 0; i < numAtom; i++)  | 
| 609 | 
+ | 
    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1); | 
| 610 | 
+ | 
   | 
| 611 | 
+ | 
 | 
| 612 | 
+ | 
  myMols = info->molecules; | 
| 613 | 
+ | 
  numMol = info->n_mol; | 
| 614 | 
+ | 
  for(int i  = 0; i < numMol; i++){ | 
| 615 | 
+ | 
    numCutoffGroups = myMols[i].getNCutoffGroups(); | 
| 616 | 
+ | 
    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff);  | 
| 617 | 
+ | 
        myCutoffGroup != NULL;  | 
| 618 | 
+ | 
        myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){ | 
| 619 | 
+ | 
 | 
| 620 | 
+ | 
      totalMass = myCutoffGroup->getMass(); | 
| 621 | 
+ | 
       | 
| 622 | 
+ | 
      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom);  | 
| 623 | 
+ | 
          cutoffAtom != NULL;  | 
| 624 | 
+ | 
          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ | 
| 625 | 
+ | 
        mfact.push_back(cutoffAtom->getMass()/totalMass); | 
| 626 | 
+ | 
      }   | 
| 627 | 
+ | 
    } | 
| 628 | 
+ | 
  } | 
| 629 | 
+ | 
 | 
| 630 | 
+ | 
} |