| 513 |
|
} |
| 514 |
|
} |
| 515 |
|
|
| 516 |
+ |
//fill globalGroupMembership |
| 517 |
|
std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0); |
| 518 |
|
for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
| 519 |
|
for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
| 524 |
|
|
| 525 |
|
} |
| 526 |
|
} |
| 527 |
< |
|
| 527 |
> |
|
| 528 |
> |
#ifdef IS_MPI |
| 529 |
|
// Since the globalGroupMembership has been zero filled and we've only |
| 530 |
|
// poked values into the atoms we know, we can do an Allreduce |
| 531 |
|
// to get the full globalGroupMembership array (We think). |
| 536 |
|
info->getGlobalGroupMembershipPointer(), |
| 537 |
|
info->getNGlobalAtoms(), |
| 538 |
|
MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
| 539 |
+ |
#else |
| 540 |
+ |
std::copy(globalGroupMembership.begin(), globalGroupMembership.end(), |
| 541 |
+ |
info->getGlobalGroupMembershipPointer()); |
| 542 |
+ |
#endif |
| 543 |
|
|
| 544 |
|
//fill molMembership |
| 545 |
|
std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0); |
| 550 |
|
globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex(); |
| 551 |
|
} |
| 552 |
|
|
| 553 |
+ |
#ifdef IS_MPI |
| 554 |
|
MPI_Allreduce(&globalMolMembership[0], |
| 555 |
|
info->getGlobalMolMembershipPointer(), |
| 556 |
|
info->getNGlobalAtoms(), |
| 557 |
|
MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
| 558 |
< |
|
| 558 |
> |
#else |
| 559 |
> |
std::copy(globalMolMembership.begin(), globalMolMembership.end(), |
| 560 |
> |
info->getGlobalMolMembershipPointer()); |
| 561 |
> |
#endif |
| 562 |
> |
|
| 563 |
|
} |
| 564 |
|
|
| 565 |
|
|
| 566 |
|
void SimCreator::initFortran(SimInfo* info) { |
| 567 |
+ |
|
| 568 |
+ |
//setup fortran simulation (Actually, in parallel version it will initialze the parallel part of fortran first) |
| 569 |
|
info->update(); |
| 570 |
|
|
| 571 |
< |
std::vector<AtomType*> atomTypes; |
| 572 |
< |
std::vector<AtomType*>::iterator i; |
| 573 |
< |
std::vector<double> cutoffRadius; |
| 574 |
< |
ForceField* ff; |
| 575 |
< |
|
| 563 |
< |
//get the unique atom types |
| 564 |
< |
atomTypes = info->getUniqueAtomTypes(); |
| 565 |
< |
|
| 566 |
< |
ff = info->getForceField(); |
| 571 |
> |
//notify fortran whether reaction field is used or not |
| 572 |
> |
//deprecated |
| 573 |
> |
int isError; |
| 574 |
> |
isError = 0; |
| 575 |
> |
initFortranFF( &useReactionField, &isError ); |
| 576 |
|
|
| 577 |
< |
//query the max cutoff radius among these atom types |
| 578 |
< |
for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
| 579 |
< |
cutoffRadius.push_back(ff->(*i)); |
| 577 |
> |
if(isError){ |
| 578 |
> |
sprintf( painCave.errMsg, |
| 579 |
> |
"SimCreator::initFortran() error: There was an error initializing the forceField in fortran.\n" ); |
| 580 |
> |
painCave.isFatal = 1; |
| 581 |
> |
simError(); |
| 582 |
|
} |
| 583 |
+ |
|
| 584 |
+ |
//figure out the cutoff radius and pass it to fortran |
| 585 |
+ |
setCutoff(info); |
| 586 |
+ |
} |
| 587 |
|
|
| 588 |
< |
double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end()); |
| 588 |
> |
void SimCreator::setCutoff(SimInfo* info) { |
| 589 |
> |
Globals* globals = info->getGlobals(); |
| 590 |
> |
double rcut; //cutoff radius |
| 591 |
> |
double rsw; //switching radius |
| 592 |
|
|
| 593 |
< |
#ifdef IS_MPI |
| 594 |
< |
//pick the max cutoff radius among the processors |
| 595 |
< |
#endif |
| 593 |
> |
simtype* st = info->getSimType(); |
| 594 |
> |
|
| 595 |
> |
if (st->SIM_uses_Charges | st->SIM_uses_Dipoles | st->SIM_uses_RF) { |
| 596 |
> |
|
| 597 |
> |
if (!globals->haveRcut()){ |
| 598 |
> |
sprintf(painCave.errMsg, |
| 599 |
> |
"SimCreator Warning: No value was set for the cutoffRadius.\n" |
| 600 |
> |
"\tOOPSE will use a default value of 15.0 angstroms" |
| 601 |
> |
"\tfor the cutoffRadius.\n"); |
| 602 |
> |
painCave.isFatal = 0; |
| 603 |
> |
simError(); |
| 604 |
> |
rcut = 15.0; |
| 605 |
> |
} else{ |
| 606 |
> |
rcut = globals->getRcut(); |
| 607 |
> |
} |
| 608 |
|
|
| 609 |
+ |
if (!globals->haveRsw()){ |
| 610 |
+ |
sprintf(painCave.errMsg, |
| 611 |
+ |
"SimCreator Warning: No value was set for switchingRadius.\n" |
| 612 |
+ |
"\tOOPSE will use a default value of\n" |
| 613 |
+ |
"\t0.95 * cutoffRadius for the switchingRadius\n"); |
| 614 |
+ |
painCave.isFatal = 0; |
| 615 |
+ |
simError(); |
| 616 |
+ |
rsw = 0.95 * rcut; |
| 617 |
+ |
} else{ |
| 618 |
+ |
rsw = globals->getRsw(); |
| 619 |
+ |
} |
| 620 |
+ |
|
| 621 |
+ |
} else { |
| 622 |
+ |
// if charge, dipole or reaction field is not used and the cutofff radius is not specified in |
| 623 |
+ |
//meta-data file, the maximum cutoff radius calculated from forcefiled will be used |
| 624 |
+ |
|
| 625 |
+ |
if (globals->haveRcut()) { |
| 626 |
+ |
rcut = globals->getRcut(); |
| 627 |
+ |
} else { |
| 628 |
+ |
//set cutoff radius to the maximum cutoff radius based on atom types in the whole system |
| 629 |
+ |
rcut = info->calcMaxCutoffRadius(); |
| 630 |
+ |
} |
| 631 |
+ |
|
| 632 |
+ |
if (globals->haveRsw()) { |
| 633 |
+ |
rsw = globals->getRsw() |
| 634 |
+ |
} else { |
| 635 |
+ |
rsw = rcut; |
| 636 |
+ |
} |
| 637 |
|
|
| 638 |
+ |
} |
| 639 |
|
|
| 640 |
+ |
//store them into Siminfo |
| 641 |
+ |
info->setRcut(rcut); |
| 642 |
+ |
info->setRsw(rsw); |
| 643 |
+ |
|
| 644 |
+ |
double rnblist = rcut + 1; // skin of neighbor list |
| 645 |
+ |
|
| 646 |
+ |
//Pass these cutoff radius etc. to fortran. This function should be called once and only once |
| 647 |
+ |
notifyFortranCutoffs(&rcut, &rsw, &rnblist); |
| 648 |
|
} |
| 649 |
|
|
| 650 |
|
void SimCreator::loadCoordinates(SimInfo* info) { |
| 664 |
|
if (nframes > 0) { |
| 665 |
|
reader.readFrame(info, nframes - 1); |
| 666 |
|
} else { |
| 667 |
< |
//invalid initial coordinate file |
| 667 |
> |
//invalid initial coordinate file |
| 668 |
|
sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n", |
| 669 |
|
globals->getInitialConfig()); |
| 670 |
|
painCave.isFatal = 1; |
