| 41 |
|
ForceField* ff, Globals* globals) : |
| 42 |
|
forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0), |
| 43 |
|
nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0), |
| 44 |
< |
nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL) { |
| 44 |
> |
nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL), |
| 45 |
> |
fortranInitialized_(false) { |
| 46 |
|
|
| 47 |
|
std::vector<std::pair<MoleculeStamp*, int> >::iterator i; |
| 48 |
|
int nCutoffAtoms; // number of atoms belong to cutoff groups |
| 160 |
|
|
| 161 |
|
Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { |
| 162 |
|
i = molecules_.begin(); |
| 163 |
< |
return i == molecules_.end() ? NULL : *i; |
| 163 |
> |
return i == molecules_.end() ? NULL : i->second; |
| 164 |
|
} |
| 165 |
|
|
| 166 |
|
Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { |
| 167 |
|
++i; |
| 168 |
< |
return i == molecules_.end() ? NULL : *i; |
| 168 |
> |
return i == molecules_.end() ? NULL : i->second; |
| 169 |
|
} |
| 170 |
|
|
| 171 |
|
|
| 363 |
|
|
| 364 |
|
void SimInfo::update() { |
| 365 |
|
|
| 366 |
+ |
setupSimType(); |
| 367 |
+ |
|
| 368 |
|
#ifdef IS_MPI |
| 369 |
|
setupFortranParallel(); |
| 370 |
|
#endif |
| 371 |
|
|
| 372 |
|
setupFortranSim(); |
| 373 |
|
|
| 374 |
+ |
setupCutoff(); |
| 375 |
+ |
|
| 376 |
+ |
//notify fortran whether reaction field is used or not. It is deprecated now |
| 377 |
+ |
//int isError = 0; |
| 378 |
+ |
//initFortranFF( &useReactionField, &isError ); |
| 379 |
+ |
|
| 380 |
+ |
//if(isError){ |
| 381 |
+ |
// sprintf( painCave.errMsg, |
| 382 |
+ |
// "SimCreator::initFortran() error: There was an error initializing the forceField in fortran.\n" ); |
| 383 |
+ |
// painCave.isFatal = 1; |
| 384 |
+ |
// simError(); |
| 385 |
+ |
//} |
| 386 |
+ |
|
| 387 |
|
calcNdf(); |
| 388 |
|
calcNdfRaw(); |
| 389 |
|
calcNdfTrans(); |
| 390 |
+ |
|
| 391 |
+ |
fortranInitialized_ = true; |
| 392 |
|
} |
| 393 |
|
|
| 394 |
|
std::set<AtomType*> SimInfo::getUniqueAtomTypes() { |
| 562 |
|
} |
| 563 |
|
} |
| 564 |
|
|
| 565 |
+ |
//fill molMembershipArray |
| 566 |
+ |
//molMembershipArray is filled by SimCreator |
| 567 |
+ |
std::vector<int> molMembershipArray(nGlobalAtoms_); |
| 568 |
+ |
for (int i = 0; i < nGlobalAtoms_; i++) { |
| 569 |
+ |
molMembershipArray.push_back(globalMolMembership_[i] + 1); |
| 570 |
+ |
} |
| 571 |
+ |
|
| 572 |
|
//setup fortran simulation |
| 573 |
|
//gloalExcludes and molMembershipArray should go away (They are never used) |
| 574 |
|
//why the hell fortran need to know molecule? |
| 654 |
|
|
| 655 |
|
#endif |
| 656 |
|
|
| 657 |
+ |
double SimInfo::calcMaxCutoffRadius() { |
| 658 |
+ |
|
| 659 |
+ |
|
| 660 |
+ |
std::vector<AtomType*> atomTypes; |
| 661 |
+ |
std::vector<AtomType*>::iterator i; |
| 662 |
+ |
std::vector<double> cutoffRadius; |
| 663 |
+ |
|
| 664 |
+ |
//get the unique atom types |
| 665 |
+ |
atomTypes = getUniqueAtomTypes(); |
| 666 |
+ |
|
| 667 |
+ |
//query the max cutoff radius among these atom types |
| 668 |
+ |
for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
| 669 |
+ |
cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); |
| 670 |
+ |
} |
| 671 |
+ |
|
| 672 |
+ |
double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end()); |
| 673 |
+ |
#ifdef IS_MPI |
| 674 |
+ |
//pick the max cutoff radius among the processors |
| 675 |
+ |
#endif |
| 676 |
+ |
|
| 677 |
+ |
return maxCutoffRadius; |
| 678 |
+ |
} |
| 679 |
+ |
|
| 680 |
+ |
void SimInfo::setupCutoff() { |
| 681 |
+ |
double rcut_; //cutoff radius |
| 682 |
+ |
double rsw_; //switching radius |
| 683 |
+ |
|
| 684 |
+ |
if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
| 685 |
+ |
|
| 686 |
+ |
if (!globals_->haveRcut()){ |
| 687 |
+ |
sprintf(painCave.errMsg, |
| 688 |
+ |
"SimCreator Warning: No value was set for the cutoffRadius.\n" |
| 689 |
+ |
"\tOOPSE will use a default value of 15.0 angstroms" |
| 690 |
+ |
"\tfor the cutoffRadius.\n"); |
| 691 |
+ |
painCave.isFatal = 0; |
| 692 |
+ |
simError(); |
| 693 |
+ |
rcut_ = 15.0; |
| 694 |
+ |
} else{ |
| 695 |
+ |
rcut_ = globals_->getRcut(); |
| 696 |
+ |
} |
| 697 |
+ |
|
| 698 |
+ |
if (!globals_->haveRsw()){ |
| 699 |
+ |
sprintf(painCave.errMsg, |
| 700 |
+ |
"SimCreator Warning: No value was set for switchingRadius.\n" |
| 701 |
+ |
"\tOOPSE will use a default value of\n" |
| 702 |
+ |
"\t0.95 * cutoffRadius for the switchingRadius\n"); |
| 703 |
+ |
painCave.isFatal = 0; |
| 704 |
+ |
simError(); |
| 705 |
+ |
rsw_ = 0.95 * rcut_; |
| 706 |
+ |
} else{ |
| 707 |
+ |
rsw_ = globals_->getRsw(); |
| 708 |
+ |
} |
| 709 |
+ |
|
| 710 |
+ |
} else { |
| 711 |
+ |
// if charge, dipole or reaction field is not used and the cutofff radius is not specified in |
| 712 |
+ |
//meta-data file, the maximum cutoff radius calculated from forcefiled will be used |
| 713 |
+ |
|
| 714 |
+ |
if (globals_->haveRcut()) { |
| 715 |
+ |
rcut_ = globals_->getRcut(); |
| 716 |
+ |
} else { |
| 717 |
+ |
//set cutoff radius to the maximum cutoff radius based on atom types in the whole system |
| 718 |
+ |
rcut_ = calcMaxCutoffRadius(); |
| 719 |
+ |
} |
| 720 |
+ |
|
| 721 |
+ |
if (globals_->haveRsw()) { |
| 722 |
+ |
rsw_ = globals_->getRsw() |
| 723 |
+ |
} else { |
| 724 |
+ |
rsw_ = rcut_; |
| 725 |
+ |
} |
| 726 |
+ |
|
| 727 |
+ |
} |
| 728 |
+ |
|
| 729 |
+ |
double rnblist = rcut_ + 1; // skin of neighbor list |
| 730 |
+ |
|
| 731 |
+ |
//Pass these cutoff radius etc. to fortran. This function should be called once and only once |
| 732 |
+ |
notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); |
| 733 |
+ |
} |
| 734 |
+ |
|
| 735 |
|
void SimInfo::addProperty(GenericData* genData) { |
| 736 |
|
properties_.addProperty(genData); |
| 737 |
|
} |
