| 64 |  | #include "selection/SelectionManager.hpp" | 
| 65 |  | #include "io/ForceFieldOptions.hpp" | 
| 66 |  | #include "UseTheForce/ForceField.hpp" | 
| 67 | + | #include "nonbonded/InteractionManager.hpp" | 
| 68 |  |  | 
| 69 |  |  | 
| 70 |  | #ifdef IS_MPI | 
| 72 |  | #include "UseTheForce/DarkSide/simParallel_interface.h" | 
| 73 |  | #endif | 
| 74 |  |  | 
| 75 | + | using namespace std; | 
| 76 |  | namespace OpenMD { | 
| 75 | – | std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { | 
| 76 | – | std::map<int, std::set<int> >::iterator i = container.find(index); | 
| 77 | – | std::set<int> result; | 
| 78 | – | if (i != container.end()) { | 
| 79 | – | result = i->second; | 
| 80 | – | } | 
| 81 | – |  | 
| 82 | – | return result; | 
| 83 | – | } | 
| 77 |  |  | 
| 78 |  | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | 
| 79 |  | forceField_(ff), simParams_(simParams), | 
| 83 |  | nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0), | 
| 84 |  | nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), | 
| 85 |  | nConstraints_(0), sman_(NULL), fortranInitialized_(false), | 
| 86 | < | calcBoxDipole_(false), useAtomicVirial_(true) { | 
| 87 | < |  | 
| 88 | < |  | 
| 89 | < | MoleculeStamp* molStamp; | 
| 90 | < | int nMolWithSameStamp; | 
| 91 | < | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | 
| 92 | < | int nGroups = 0;      //total cutoff groups defined in meta-data file | 
| 93 | < | CutoffGroupStamp* cgStamp; | 
| 94 | < | RigidBodyStamp* rbStamp; | 
| 95 | < | int nRigidAtoms = 0; | 
| 96 | < |  | 
| 97 | < | std::vector<Component*> components = simParams->getComponents(); | 
| 86 | > | calcBoxDipole_(false), useAtomicVirial_(true) { | 
| 87 | > |  | 
| 88 | > | MoleculeStamp* molStamp; | 
| 89 | > | int nMolWithSameStamp; | 
| 90 | > | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | 
| 91 | > | int nGroups = 0;       //total cutoff groups defined in meta-data file | 
| 92 | > | CutoffGroupStamp* cgStamp; | 
| 93 | > | RigidBodyStamp* rbStamp; | 
| 94 | > | int nRigidAtoms = 0; | 
| 95 | > |  | 
| 96 | > | vector<Component*> components = simParams->getComponents(); | 
| 97 | > |  | 
| 98 | > | for (vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { | 
| 99 | > | molStamp = (*i)->getMoleculeStamp(); | 
| 100 | > | nMolWithSameStamp = (*i)->getNMol(); | 
| 101 |  |  | 
| 102 | < | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { | 
| 103 | < | molStamp = (*i)->getMoleculeStamp(); | 
| 104 | < | nMolWithSameStamp = (*i)->getNMol(); | 
| 105 | < |  | 
| 106 | < | addMoleculeStamp(molStamp, nMolWithSameStamp); | 
| 107 | < |  | 
| 108 | < | //calculate atoms in molecules | 
| 109 | < | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | 
| 110 | < |  | 
| 111 | < | //calculate atoms in cutoff groups | 
| 112 | < | int nAtomsInGroups = 0; | 
| 113 | < | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | 
| 118 | < |  | 
| 119 | < | for (int j=0; j < nCutoffGroupsInStamp; j++) { | 
| 120 | < | cgStamp = molStamp->getCutoffGroupStamp(j); | 
| 121 | < | nAtomsInGroups += cgStamp->getNMembers(); | 
| 122 | < | } | 
| 123 | < |  | 
| 124 | < | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; | 
| 125 | < |  | 
| 126 | < | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; | 
| 127 | < |  | 
| 128 | < | //calculate atoms in rigid bodies | 
| 129 | < | int nAtomsInRigidBodies = 0; | 
| 130 | < | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | 
| 131 | < |  | 
| 132 | < | for (int j=0; j < nRigidBodiesInStamp; j++) { | 
| 133 | < | rbStamp = molStamp->getRigidBodyStamp(j); | 
| 134 | < | nAtomsInRigidBodies += rbStamp->getNMembers(); | 
| 135 | < | } | 
| 136 | < |  | 
| 137 | < | nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; | 
| 138 | < | nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; | 
| 139 | < |  | 
| 102 | > | addMoleculeStamp(molStamp, nMolWithSameStamp); | 
| 103 | > |  | 
| 104 | > | //calculate atoms in molecules | 
| 105 | > | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | 
| 106 | > |  | 
| 107 | > | //calculate atoms in cutoff groups | 
| 108 | > | int nAtomsInGroups = 0; | 
| 109 | > | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | 
| 110 | > |  | 
| 111 | > | for (int j=0; j < nCutoffGroupsInStamp; j++) { | 
| 112 | > | cgStamp = molStamp->getCutoffGroupStamp(j); | 
| 113 | > | nAtomsInGroups += cgStamp->getNMembers(); | 
| 114 |  | } | 
| 115 | < |  | 
| 116 | < | //every free atom (atom does not belong to cutoff groups) is a cutoff | 
| 117 | < | //group therefore the total number of cutoff groups in the system is | 
| 118 | < | //equal to the total number of atoms minus number of atoms belong to | 
| 119 | < | //cutoff group defined in meta-data file plus the number of cutoff | 
| 120 | < | //groups defined in meta-data file | 
| 121 | < | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | 
| 122 | < |  | 
| 123 | < | //every free atom (atom does not belong to rigid bodies) is an | 
| 124 | < | //integrable object therefore the total number of integrable objects | 
| 125 | < | //in the system is equal to the total number of atoms minus number of | 
| 126 | < | //atoms belong to rigid body defined in meta-data file plus the number | 
| 127 | < | //of rigid bodies defined in meta-data file | 
| 128 | < | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms | 
| 129 | < | + nGlobalRigidBodies_; | 
| 130 | < |  | 
| 131 | < | nGlobalMols_ = molStampIds_.size(); | 
| 132 | < | molToProcMap_.resize(nGlobalMols_); | 
| 133 | < | } | 
| 134 | < |  | 
| 115 | > |  | 
| 116 | > | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; | 
| 117 | > |  | 
| 118 | > | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; | 
| 119 | > |  | 
| 120 | > | //calculate atoms in rigid bodies | 
| 121 | > | int nAtomsInRigidBodies = 0; | 
| 122 | > | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | 
| 123 | > |  | 
| 124 | > | for (int j=0; j < nRigidBodiesInStamp; j++) { | 
| 125 | > | rbStamp = molStamp->getRigidBodyStamp(j); | 
| 126 | > | nAtomsInRigidBodies += rbStamp->getNMembers(); | 
| 127 | > | } | 
| 128 | > |  | 
| 129 | > | nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; | 
| 130 | > | nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; | 
| 131 | > |  | 
| 132 | > | } | 
| 133 | > |  | 
| 134 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff | 
| 135 | > | //group therefore the total number of cutoff groups in the system is | 
| 136 | > | //equal to the total number of atoms minus number of atoms belong to | 
| 137 | > | //cutoff group defined in meta-data file plus the number of cutoff | 
| 138 | > | //groups defined in meta-data file | 
| 139 | > | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | 
| 140 | > |  | 
| 141 | > | //every free atom (atom does not belong to rigid bodies) is an | 
| 142 | > | //integrable object therefore the total number of integrable objects | 
| 143 | > | //in the system is equal to the total number of atoms minus number of | 
| 144 | > | //atoms belong to rigid body defined in meta-data file plus the number | 
| 145 | > | //of rigid bodies defined in meta-data file | 
| 146 | > | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms | 
| 147 | > | + nGlobalRigidBodies_; | 
| 148 | > |  | 
| 149 | > | nGlobalMols_ = molStampIds_.size(); | 
| 150 | > | molToProcMap_.resize(nGlobalMols_); | 
| 151 | > | } | 
| 152 | > |  | 
| 153 |  | SimInfo::~SimInfo() { | 
| 154 | < | std::map<int, Molecule*>::iterator i; | 
| 154 | > | map<int, Molecule*>::iterator i; | 
| 155 |  | for (i = molecules_.begin(); i != molecules_.end(); ++i) { | 
| 156 |  | delete i->second; | 
| 157 |  | } | 
| 162 |  | delete forceField_; | 
| 163 |  | } | 
| 164 |  |  | 
| 173 | – | int SimInfo::getNGlobalConstraints() { | 
| 174 | – | int nGlobalConstraints; | 
| 175 | – | #ifdef IS_MPI | 
| 176 | – | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, | 
| 177 | – | MPI_COMM_WORLD); | 
| 178 | – | #else | 
| 179 | – | nGlobalConstraints =  nConstraints_; | 
| 180 | – | #endif | 
| 181 | – | return nGlobalConstraints; | 
| 182 | – | } | 
| 165 |  |  | 
| 166 |  | bool SimInfo::addMolecule(Molecule* mol) { | 
| 167 |  | MoleculeIterator i; | 
| 168 | < |  | 
| 168 | > |  | 
| 169 |  | i = molecules_.find(mol->getGlobalIndex()); | 
| 170 |  | if (i == molecules_.end() ) { | 
| 171 | < |  | 
| 172 | < | molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); | 
| 173 | < |  | 
| 171 | > |  | 
| 172 | > | molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); | 
| 173 | > |  | 
| 174 |  | nAtoms_ += mol->getNAtoms(); | 
| 175 |  | nBonds_ += mol->getNBonds(); | 
| 176 |  | nBends_ += mol->getNBends(); | 
| 180 |  | nIntegrableObjects_ += mol->getNIntegrableObjects(); | 
| 181 |  | nCutoffGroups_ += mol->getNCutoffGroups(); | 
| 182 |  | nConstraints_ += mol->getNConstraintPairs(); | 
| 183 | < |  | 
| 183 | > |  | 
| 184 |  | addInteractionPairs(mol); | 
| 185 | < |  | 
| 185 | > |  | 
| 186 |  | return true; | 
| 187 |  | } else { | 
| 188 |  | return false; | 
| 189 |  | } | 
| 190 |  | } | 
| 191 | < |  | 
| 191 | > |  | 
| 192 |  | bool SimInfo::removeMolecule(Molecule* mol) { | 
| 193 |  | MoleculeIterator i; | 
| 194 |  | i = molecules_.find(mol->getGlobalIndex()); | 
| 216 |  | } else { | 
| 217 |  | return false; | 
| 218 |  | } | 
| 237 | – |  | 
| 238 | – |  | 
| 219 |  | } | 
| 220 |  |  | 
| 221 |  |  | 
| 233 |  | void SimInfo::calcNdf() { | 
| 234 |  | int ndf_local; | 
| 235 |  | MoleculeIterator i; | 
| 236 | < | std::vector<StuntDouble*>::iterator j; | 
| 236 | > | vector<StuntDouble*>::iterator j; | 
| 237 |  | Molecule* mol; | 
| 238 |  | StuntDouble* integrableObject; | 
| 239 |  |  | 
| 284 |  | int ndfRaw_local; | 
| 285 |  |  | 
| 286 |  | MoleculeIterator i; | 
| 287 | < | std::vector<StuntDouble*>::iterator j; | 
| 287 | > | vector<StuntDouble*>::iterator j; | 
| 288 |  | Molecule* mol; | 
| 289 |  | StuntDouble* integrableObject; | 
| 290 |  |  | 
| 333 |  |  | 
| 334 |  | void SimInfo::addInteractionPairs(Molecule* mol) { | 
| 335 |  | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | 
| 336 | < | std::vector<Bond*>::iterator bondIter; | 
| 337 | < | std::vector<Bend*>::iterator bendIter; | 
| 338 | < | std::vector<Torsion*>::iterator torsionIter; | 
| 339 | < | std::vector<Inversion*>::iterator inversionIter; | 
| 336 | > | vector<Bond*>::iterator bondIter; | 
| 337 | > | vector<Bend*>::iterator bendIter; | 
| 338 | > | vector<Torsion*>::iterator torsionIter; | 
| 339 | > | vector<Inversion*>::iterator inversionIter; | 
| 340 |  | Bond* bond; | 
| 341 |  | Bend* bend; | 
| 342 |  | Torsion* torsion; | 
| 354 |  | // always be excluded.  These are done at the bottom of this | 
| 355 |  | // function. | 
| 356 |  |  | 
| 357 | < | std::map<int, std::set<int> > atomGroups; | 
| 357 | > | map<int, set<int> > atomGroups; | 
| 358 |  | Molecule::RigidBodyIterator rbIter; | 
| 359 |  | RigidBody* rb; | 
| 360 |  | Molecule::IntegrableObjectIterator ii; | 
| 366 |  |  | 
| 367 |  | if (integrableObject->isRigidBody()) { | 
| 368 |  | rb = static_cast<RigidBody*>(integrableObject); | 
| 369 | < | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 370 | < | std::set<int> rigidAtoms; | 
| 369 | > | vector<Atom*> atoms = rb->getAtoms(); | 
| 370 | > | set<int> rigidAtoms; | 
| 371 |  | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | 
| 372 |  | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 373 |  | } | 
| 374 |  | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | 
| 375 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 375 | > | atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 376 |  | } | 
| 377 |  | } else { | 
| 378 | < | std::set<int> oneAtomSet; | 
| 378 | > | set<int> oneAtomSet; | 
| 379 |  | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 380 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 380 | > | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 381 |  | } | 
| 382 |  | } | 
| 383 |  |  | 
| 480 |  |  | 
| 481 |  | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 482 |  | rb = mol->nextRigidBody(rbIter)) { | 
| 483 | < | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 483 | > | vector<Atom*> atoms = rb->getAtoms(); | 
| 484 |  | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | 
| 485 |  | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | 
| 486 |  | a = atoms[i]->getGlobalIndex(); | 
| 494 |  |  | 
| 495 |  | void SimInfo::removeInteractionPairs(Molecule* mol) { | 
| 496 |  | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | 
| 497 | < | std::vector<Bond*>::iterator bondIter; | 
| 498 | < | std::vector<Bend*>::iterator bendIter; | 
| 499 | < | std::vector<Torsion*>::iterator torsionIter; | 
| 500 | < | std::vector<Inversion*>::iterator inversionIter; | 
| 497 | > | vector<Bond*>::iterator bondIter; | 
| 498 | > | vector<Bend*>::iterator bendIter; | 
| 499 | > | vector<Torsion*>::iterator torsionIter; | 
| 500 | > | vector<Inversion*>::iterator inversionIter; | 
| 501 |  | Bond* bond; | 
| 502 |  | Bend* bend; | 
| 503 |  | Torsion* torsion; | 
| 507 |  | int c; | 
| 508 |  | int d; | 
| 509 |  |  | 
| 510 | < | std::map<int, std::set<int> > atomGroups; | 
| 510 | > | map<int, set<int> > atomGroups; | 
| 511 |  | Molecule::RigidBodyIterator rbIter; | 
| 512 |  | RigidBody* rb; | 
| 513 |  | Molecule::IntegrableObjectIterator ii; | 
| 519 |  |  | 
| 520 |  | if (integrableObject->isRigidBody()) { | 
| 521 |  | rb = static_cast<RigidBody*>(integrableObject); | 
| 522 | < | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 523 | < | std::set<int> rigidAtoms; | 
| 522 | > | vector<Atom*> atoms = rb->getAtoms(); | 
| 523 | > | set<int> rigidAtoms; | 
| 524 |  | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | 
| 525 |  | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 526 |  | } | 
| 527 |  | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | 
| 528 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 528 | > | atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 529 |  | } | 
| 530 |  | } else { | 
| 531 | < | std::set<int> oneAtomSet; | 
| 531 | > | set<int> oneAtomSet; | 
| 532 |  | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 533 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 533 | > | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 534 |  | } | 
| 535 |  | } | 
| 536 |  |  | 
| 633 |  |  | 
| 634 |  | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 635 |  | rb = mol->nextRigidBody(rbIter)) { | 
| 636 | < | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 636 | > | vector<Atom*> atoms = rb->getAtoms(); | 
| 637 |  | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | 
| 638 |  | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | 
| 639 |  | a = atoms[i]->getGlobalIndex(); | 
| 659 |  | void SimInfo::update() { | 
| 660 |  |  | 
| 661 |  | setupSimType(); | 
| 662 | + | setupCutoffRadius(); | 
| 663 | + | setupSwitchingRadius(); | 
| 664 | + | setupCutoffMethod(); | 
| 665 | + | setupSkinThickness(); | 
| 666 | + | setupSwitchingFunction(); | 
| 667 | + | setupAccumulateBoxDipole(); | 
| 668 |  |  | 
| 669 |  | #ifdef IS_MPI | 
| 670 |  | setupFortranParallel(); | 
| 671 |  | #endif | 
| 686 | – |  | 
| 672 |  | setupFortranSim(); | 
| 673 | + | fortranInitialized_ = true; | 
| 674 |  |  | 
| 689 | – | //setup fortran force field | 
| 690 | – | /** @deprecate */ | 
| 691 | – | int isError = 0; | 
| 692 | – |  | 
| 693 | – | setupCutoff(); | 
| 694 | – |  | 
| 695 | – | setupElectrostaticSummationMethod( isError ); | 
| 696 | – | setupSwitchingFunction(); | 
| 697 | – | setupAccumulateBoxDipole(); | 
| 698 | – |  | 
| 699 | – | if(isError){ | 
| 700 | – | sprintf( painCave.errMsg, | 
| 701 | – | "ForceField error: There was an error initializing the forceField in fortran.\n" ); | 
| 702 | – | painCave.isFatal = 1; | 
| 703 | – | simError(); | 
| 704 | – | } | 
| 705 | – |  | 
| 675 |  | calcNdf(); | 
| 676 |  | calcNdfRaw(); | 
| 677 |  | calcNdfTrans(); | 
| 709 | – |  | 
| 710 | – | fortranInitialized_ = true; | 
| 678 |  | } | 
| 679 | < |  | 
| 680 | < | std::set<AtomType*> SimInfo::getUniqueAtomTypes() { | 
| 679 | > |  | 
| 680 | > | set<AtomType*> SimInfo::getSimulatedAtomTypes() { | 
| 681 |  | SimInfo::MoleculeIterator mi; | 
| 682 |  | Molecule* mol; | 
| 683 |  | Molecule::AtomIterator ai; | 
| 684 |  | Atom* atom; | 
| 685 | < | std::set<AtomType*> atomTypes; | 
| 686 | < |  | 
| 685 | > | set<AtomType*> atomTypes; | 
| 686 | > |  | 
| 687 |  | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | 
| 688 | < |  | 
| 688 | > |  | 
| 689 |  | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 690 |  | atomTypes.insert(atom->getAtomType()); | 
| 691 |  | } | 
| 692 | < |  | 
| 692 | > |  | 
| 693 |  | } | 
| 694 | < |  | 
| 694 | > |  | 
| 695 |  | return atomTypes; | 
| 696 |  | } | 
| 697 |  |  | 
| 698 | < | void SimInfo::setupSimType() { | 
| 699 | < | std::set<AtomType*>::iterator i; | 
| 700 | < | std::set<AtomType*> atomTypes; | 
| 701 | < | atomTypes = getUniqueAtomTypes(); | 
| 702 | < |  | 
| 703 | < | int useLennardJones = 0; | 
| 704 | < | int useElectrostatic = 0; | 
| 705 | < | int useEAM = 0; | 
| 706 | < | int useSC = 0; | 
| 707 | < | int useCharge = 0; | 
| 708 | < | int useDirectional = 0; | 
| 709 | < | int useDipole = 0; | 
| 710 | < | int useGayBerne = 0; | 
| 711 | < | int useSticky = 0; | 
| 712 | < | int useStickyPower = 0; | 
| 713 | < | int useShape = 0; | 
| 714 | < | int useFLARB = 0; //it is not in AtomType yet | 
| 715 | < | int useDirectionalAtom = 0; | 
| 716 | < | int useElectrostatics = 0; | 
| 717 | < | //usePBC and useRF are from simParams | 
| 718 | < | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | 
| 719 | < | int useRF; | 
| 720 | < | int useSF; | 
| 721 | < | int useSP; | 
| 722 | < | int useBoxDipole; | 
| 723 | < |  | 
| 724 | < | std::string myMethod; | 
| 725 | < |  | 
| 726 | < | // set the useRF logical | 
| 727 | < | useRF = 0; | 
| 728 | < | useSF = 0; | 
| 729 | < | useSP = 0; | 
| 730 | < | useBoxDipole = 0; | 
| 731 | < |  | 
| 732 | < | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 733 | < | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 734 | < | toUpper(myMethod); | 
| 735 | < | if (myMethod == "REACTION_FIELD"){ | 
| 736 | < | useRF = 1; | 
| 770 | < | } else if (myMethod == "SHIFTED_FORCE"){ | 
| 771 | < | useSF = 1; | 
| 772 | < | } else if (myMethod == "SHIFTED_POTENTIAL"){ | 
| 773 | < | useSP = 1; | 
| 774 | < | } | 
| 698 | > | /** | 
| 699 | > | * setupCutoffRadius | 
| 700 | > | * | 
| 701 | > | *  If the cutoffRadius was explicitly set, use that value. | 
| 702 | > | *  If the cutoffRadius was not explicitly set: | 
| 703 | > | *      Are there electrostatic atoms?  Use 12.0 Angstroms. | 
| 704 | > | *      No electrostatic atoms?  Poll the atom types present in the | 
| 705 | > | *      simulation for suggested cutoff values (e.g. 2.5 * sigma). | 
| 706 | > | *      Use the maximum suggested value that was found. | 
| 707 | > | */ | 
| 708 | > | void SimInfo::setupCutoffRadius() { | 
| 709 | > |  | 
| 710 | > | if (simParams_->haveCutoffRadius()) { | 
| 711 | > | cutoffRadius_ = simParams_->getCutoffRadius(); | 
| 712 | > | } else { | 
| 713 | > | if (usesElectrostaticAtoms_) { | 
| 714 | > | sprintf(painCave.errMsg, | 
| 715 | > | "SimInfo Warning: No value was set for the cutoffRadius.\n" | 
| 716 | > | "\tOpenMD will use a default value of 12.0 angstroms" | 
| 717 | > | "\tfor the cutoffRadius.\n"); | 
| 718 | > | painCave.isFatal = 0; | 
| 719 | > | simError(); | 
| 720 | > | cutoffRadius_ = 12.0; | 
| 721 | > | } else { | 
| 722 | > | RealType thisCut; | 
| 723 | > | set<AtomType*>::iterator i; | 
| 724 | > | set<AtomType*> atomTypes; | 
| 725 | > | atomTypes = getSimulatedAtomTypes(); | 
| 726 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 727 | > | thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i)); | 
| 728 | > | cutoffRadius_ = max(thisCut, cutoffRadius_); | 
| 729 | > | } | 
| 730 | > | sprintf(painCave.errMsg, | 
| 731 | > | "SimInfo Warning: No value was set for the cutoffRadius.\n" | 
| 732 | > | "\tOpenMD will use %lf angstroms.\n", | 
| 733 | > | cutoffRadius_); | 
| 734 | > | painCave.isFatal = 0; | 
| 735 | > | simError(); | 
| 736 | > | } | 
| 737 |  | } | 
| 738 | + |  | 
| 739 | + | InteractionManager::Instance()->setCutoffRadius(cutoffRadius_); | 
| 740 | + | } | 
| 741 | + |  | 
| 742 | + | /** | 
| 743 | + | * setupSwitchingRadius | 
| 744 | + | * | 
| 745 | + | *  If the switchingRadius was explicitly set, use that value (but check it) | 
| 746 | + | *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ | 
| 747 | + | */ | 
| 748 | + | void SimInfo::setupSwitchingRadius() { | 
| 749 |  |  | 
| 750 | < | if (simParams_->haveAccumulateBoxDipole()) | 
| 751 | < | if (simParams_->getAccumulateBoxDipole()) | 
| 752 | < | useBoxDipole = 1; | 
| 750 | > | if (simParams_->haveSwitchingRadius()) { | 
| 751 | > | switchingRadius_ = simParams_->getSwitchingRadius(); | 
| 752 | > | if (switchingRadius_ > cutoffRadius_) { | 
| 753 | > | sprintf(painCave.errMsg, | 
| 754 | > | "SimInfo Error: switchingRadius (%f) is larger than cutoffRadius(%f)\n", | 
| 755 | > | switchingRadius_, cutoffRadius_); | 
| 756 | > | painCave.isFatal = 1; | 
| 757 | > | simError(); | 
| 758 |  |  | 
| 759 | + | } | 
| 760 | + | } else { | 
| 761 | + | switchingRadius_ = 0.85 * cutoffRadius_; | 
| 762 | + | sprintf(painCave.errMsg, | 
| 763 | + | "SimInfo Warning: No value was set for the switchingRadius.\n" | 
| 764 | + | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" | 
| 765 | + | "\tswitchingRadius = %f. for this simulation\n", switchingRadius_); | 
| 766 | + | painCave.isFatal = 0; | 
| 767 | + | simError(); | 
| 768 | + | } | 
| 769 | + | InteractionManager::Instance()->setSwitchingRadius(switchingRadius_); | 
| 770 | + | } | 
| 771 | + |  | 
| 772 | + | /** | 
| 773 | + | * setupSkinThickness | 
| 774 | + | * | 
| 775 | + | *  If the skinThickness was explicitly set, use that value (but check it) | 
| 776 | + | *  If the skinThickness was not explicitly set: use 1.0 angstroms | 
| 777 | + | */ | 
| 778 | + | void SimInfo::setupSkinThickness() { | 
| 779 | + | if (simParams_->haveSkinThickness()) { | 
| 780 | + | skinThickness_ = simParams_->getSkinThickness(); | 
| 781 | + | } else { | 
| 782 | + | skinThickness_ = 1.0; | 
| 783 | + | sprintf(painCave.errMsg, | 
| 784 | + | "SimInfo Warning: No value was set for the skinThickness.\n" | 
| 785 | + | "\tOpenMD will use a default value of %f Angstroms\n" | 
| 786 | + | "\tfor this simulation\n", skinThickness_); | 
| 787 | + | painCave.isFatal = 0; | 
| 788 | + | simError(); | 
| 789 | + | } | 
| 790 | + | } | 
| 791 | + |  | 
| 792 | + | void SimInfo::setupSimType() { | 
| 793 | + | set<AtomType*>::iterator i; | 
| 794 | + | set<AtomType*> atomTypes; | 
| 795 | + | atomTypes = getSimulatedAtomTypes(); | 
| 796 | + |  | 
| 797 |  | useAtomicVirial_ = simParams_->getUseAtomicVirial(); | 
| 798 |  |  | 
| 799 | + | int usesElectrostatic = 0; | 
| 800 | + | int usesMetallic = 0; | 
| 801 | + | int usesDirectional = 0; | 
| 802 |  | //loop over all of the atom types | 
| 803 |  | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 804 | < | useLennardJones |= (*i)->isLennardJones(); | 
| 805 | < | useElectrostatic |= (*i)->isElectrostatic(); | 
| 806 | < | useEAM |= (*i)->isEAM(); | 
| 788 | < | useSC |= (*i)->isSC(); | 
| 789 | < | useCharge |= (*i)->isCharge(); | 
| 790 | < | useDirectional |= (*i)->isDirectional(); | 
| 791 | < | useDipole |= (*i)->isDipole(); | 
| 792 | < | useGayBerne |= (*i)->isGayBerne(); | 
| 793 | < | useSticky |= (*i)->isSticky(); | 
| 794 | < | useStickyPower |= (*i)->isStickyPower(); | 
| 795 | < | useShape |= (*i)->isShape(); | 
| 804 | > | usesElectrostatic |= (*i)->isElectrostatic(); | 
| 805 | > | usesMetallic |= (*i)->isMetal(); | 
| 806 | > | usesDirectional |= (*i)->isDirectional(); | 
| 807 |  | } | 
| 808 |  |  | 
| 798 | – | if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { | 
| 799 | – | useDirectionalAtom = 1; | 
| 800 | – | } | 
| 801 | – |  | 
| 802 | – | if (useCharge || useDipole) { | 
| 803 | – | useElectrostatics = 1; | 
| 804 | – | } | 
| 805 | – |  | 
| 809 |  | #ifdef IS_MPI | 
| 810 |  | int temp; | 
| 811 | + | temp = usesDirectional; | 
| 812 | + | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 813 |  |  | 
| 814 | < | temp = usePBC; | 
| 815 | < | MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 814 | > | temp = usesMetallic; | 
| 815 | > | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 816 |  |  | 
| 817 | < | temp = useDirectionalAtom; | 
| 818 | < | MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 814 | < |  | 
| 815 | < | temp = useLennardJones; | 
| 816 | < | MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 817 | < |  | 
| 818 | < | temp = useElectrostatics; | 
| 819 | < | MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 820 | < |  | 
| 821 | < | temp = useCharge; | 
| 822 | < | MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 823 | < |  | 
| 824 | < | temp = useDipole; | 
| 825 | < | MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 826 | < |  | 
| 827 | < | temp = useSticky; | 
| 828 | < | MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 829 | < |  | 
| 830 | < | temp = useStickyPower; | 
| 831 | < | MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 832 | < |  | 
| 833 | < | temp = useGayBerne; | 
| 834 | < | MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 835 | < |  | 
| 836 | < | temp = useEAM; | 
| 837 | < | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 838 | < |  | 
| 839 | < | temp = useSC; | 
| 840 | < | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 841 | < |  | 
| 842 | < | temp = useShape; | 
| 843 | < | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 844 | < |  | 
| 845 | < | temp = useFLARB; | 
| 846 | < | MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 847 | < |  | 
| 848 | < | temp = useRF; | 
| 849 | < | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 850 | < |  | 
| 851 | < | temp = useSF; | 
| 852 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 853 | < |  | 
| 854 | < | temp = useSP; | 
| 855 | < | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 856 | < |  | 
| 857 | < | temp = useBoxDipole; | 
| 858 | < | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 859 | < |  | 
| 860 | < | temp = useAtomicVirial_; | 
| 861 | < | MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 862 | < |  | 
| 817 | > | temp = usesElectrostatic; | 
| 818 | > | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 819 |  | #endif | 
| 820 | < | fInfo_.SIM_uses_PBC = usePBC; | 
| 821 | < | fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; | 
| 822 | < | fInfo_.SIM_uses_LennardJones = useLennardJones; | 
| 823 | < | fInfo_.SIM_uses_Electrostatics = useElectrostatics; | 
| 824 | < | fInfo_.SIM_uses_Charges = useCharge; | 
| 825 | < | fInfo_.SIM_uses_Dipoles = useDipole; | 
| 870 | < | fInfo_.SIM_uses_Sticky = useSticky; | 
| 871 | < | fInfo_.SIM_uses_StickyPower = useStickyPower; | 
| 872 | < | fInfo_.SIM_uses_GayBerne = useGayBerne; | 
| 873 | < | fInfo_.SIM_uses_EAM = useEAM; | 
| 874 | < | fInfo_.SIM_uses_SC = useSC; | 
| 875 | < | fInfo_.SIM_uses_Shapes = useShape; | 
| 876 | < | fInfo_.SIM_uses_FLARB = useFLARB; | 
| 877 | < | fInfo_.SIM_uses_RF = useRF; | 
| 878 | < | fInfo_.SIM_uses_SF = useSF; | 
| 879 | < | fInfo_.SIM_uses_SP = useSP; | 
| 880 | < | fInfo_.SIM_uses_BoxDipole = useBoxDipole; | 
| 881 | < | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; | 
| 820 | > | fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_; | 
| 821 | > | fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; | 
| 822 | > | fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; | 
| 823 | > | fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; | 
| 824 | > | fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; | 
| 825 | > | fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; | 
| 826 |  | } | 
| 827 |  |  | 
| 828 |  | void SimInfo::setupFortranSim() { | 
| 829 |  | int isError; | 
| 830 |  | int nExclude, nOneTwo, nOneThree, nOneFour; | 
| 831 | < | std::vector<int> fortranGlobalGroupMembership; | 
| 831 | > | vector<int> fortranGlobalGroupMembership; | 
| 832 |  |  | 
| 833 | + | notifyFortranSkinThickness(&skinThickness_); | 
| 834 | + |  | 
| 835 | + | int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; | 
| 836 | + | int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; | 
| 837 | + | notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf); | 
| 838 | + |  | 
| 839 |  | isError = 0; | 
| 840 |  |  | 
| 841 |  | //globalGroupMembership_ is filled by SimCreator | 
| 844 |  | } | 
| 845 |  |  | 
| 846 |  | //calculate mass ratio of cutoff group | 
| 847 | < | std::vector<RealType> mfact; | 
| 847 | > | vector<RealType> mfact; | 
| 848 |  | SimInfo::MoleculeIterator mi; | 
| 849 |  | Molecule* mol; | 
| 850 |  | Molecule::CutoffGroupIterator ci; | 
| 871 |  | } | 
| 872 |  |  | 
| 873 |  | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) | 
| 874 | < | std::vector<int> identArray; | 
| 874 | > | vector<int> identArray; | 
| 875 |  |  | 
| 876 |  | //to avoid memory reallocation, reserve enough space identArray | 
| 877 |  | identArray.reserve(getNAtoms()); | 
| 884 |  |  | 
| 885 |  | //fill molMembershipArray | 
| 886 |  | //molMembershipArray is filled by SimCreator | 
| 887 | < | std::vector<int> molMembershipArray(nGlobalAtoms_); | 
| 887 | > | vector<int> molMembershipArray(nGlobalAtoms_); | 
| 888 |  | for (int i = 0; i < nGlobalAtoms_; i++) { | 
| 889 |  | molMembershipArray[i] = globalMolMembership_[i] + 1; | 
| 890 |  | } | 
| 937 |  | void SimInfo::setupFortranParallel() { | 
| 938 |  | #ifdef IS_MPI | 
| 939 |  | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex | 
| 940 | < | std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0); | 
| 941 | < | std::vector<int> localToGlobalCutoffGroupIndex; | 
| 940 | > | vector<int> localToGlobalAtomIndex(getNAtoms(), 0); | 
| 941 | > | vector<int> localToGlobalCutoffGroupIndex; | 
| 942 |  | SimInfo::MoleculeIterator mi; | 
| 943 |  | Molecule::AtomIterator ai; | 
| 944 |  | Molecule::CutoffGroupIterator ci; | 
| 988 |  | errorCheckPoint(); | 
| 989 |  |  | 
| 990 |  | #endif | 
| 1041 | – | } | 
| 1042 | – |  | 
| 1043 | – | void SimInfo::setupCutoff() { | 
| 1044 | – |  | 
| 1045 | – | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); | 
| 1046 | – |  | 
| 1047 | – | // Check the cutoff policy | 
| 1048 | – | int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default | 
| 1049 | – |  | 
| 1050 | – | // Set LJ shifting bools to false | 
| 1051 | – | ljsp_ = 0; | 
| 1052 | – | ljsf_ = 0; | 
| 1053 | – |  | 
| 1054 | – | std::string myPolicy; | 
| 1055 | – | if (forceFieldOptions_.haveCutoffPolicy()){ | 
| 1056 | – | myPolicy = forceFieldOptions_.getCutoffPolicy(); | 
| 1057 | – | }else if (simParams_->haveCutoffPolicy()) { | 
| 1058 | – | myPolicy = simParams_->getCutoffPolicy(); | 
| 1059 | – | } | 
| 1060 | – |  | 
| 1061 | – | if (!myPolicy.empty()){ | 
| 1062 | – | toUpper(myPolicy); | 
| 1063 | – | if (myPolicy == "MIX") { | 
| 1064 | – | cp = MIX_CUTOFF_POLICY; | 
| 1065 | – | } else { | 
| 1066 | – | if (myPolicy == "MAX") { | 
| 1067 | – | cp = MAX_CUTOFF_POLICY; | 
| 1068 | – | } else { | 
| 1069 | – | if (myPolicy == "TRADITIONAL") { | 
| 1070 | – | cp = TRADITIONAL_CUTOFF_POLICY; | 
| 1071 | – | } else { | 
| 1072 | – | // throw error | 
| 1073 | – | sprintf( painCave.errMsg, | 
| 1074 | – | "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); | 
| 1075 | – | painCave.isFatal = 1; | 
| 1076 | – | simError(); | 
| 1077 | – | } | 
| 1078 | – | } | 
| 1079 | – | } | 
| 1080 | – | } | 
| 1081 | – | notifyFortranCutoffPolicy(&cp); | 
| 1082 | – |  | 
| 1083 | – | // Check the Skin Thickness for neighborlists | 
| 1084 | – | RealType skin; | 
| 1085 | – | if (simParams_->haveSkinThickness()) { | 
| 1086 | – | skin = simParams_->getSkinThickness(); | 
| 1087 | – | notifyFortranSkinThickness(&skin); | 
| 1088 | – | } | 
| 1089 | – |  | 
| 1090 | – | // Check if the cutoff was set explicitly: | 
| 1091 | – | if (simParams_->haveCutoffRadius()) { | 
| 1092 | – | rcut_ = simParams_->getCutoffRadius(); | 
| 1093 | – | if (simParams_->haveSwitchingRadius()) { | 
| 1094 | – | rsw_  = simParams_->getSwitchingRadius(); | 
| 1095 | – | } else { | 
| 1096 | – | if (fInfo_.SIM_uses_Charges | | 
| 1097 | – | fInfo_.SIM_uses_Dipoles | | 
| 1098 | – | fInfo_.SIM_uses_RF) { | 
| 1099 | – |  | 
| 1100 | – | rsw_ = 0.85 * rcut_; | 
| 1101 | – | sprintf(painCave.errMsg, | 
| 1102 | – | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 1103 | – | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" | 
| 1104 | – | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 1105 | – | painCave.isFatal = 0; | 
| 1106 | – | simError(); | 
| 1107 | – | } else { | 
| 1108 | – | rsw_ = rcut_; | 
| 1109 | – | sprintf(painCave.errMsg, | 
| 1110 | – | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 1111 | – | "\tOpenMD will use the same value as the cutoffRadius.\n" | 
| 1112 | – | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 1113 | – | painCave.isFatal = 0; | 
| 1114 | – | simError(); | 
| 1115 | – | } | 
| 1116 | – | } | 
| 1117 | – |  | 
| 1118 | – | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1119 | – | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1120 | – | toUpper(myMethod); | 
| 1121 | – |  | 
| 1122 | – | if (myMethod == "SHIFTED_POTENTIAL") { | 
| 1123 | – | ljsp_ = 1; | 
| 1124 | – | } else if (myMethod == "SHIFTED_FORCE") { | 
| 1125 | – | ljsf_ = 1; | 
| 1126 | – | } | 
| 1127 | – | } | 
| 1128 | – |  | 
| 1129 | – | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); | 
| 1130 | – |  | 
| 1131 | – | } else { | 
| 1132 | – |  | 
| 1133 | – | // For electrostatic atoms, we'll assume a large safe value: | 
| 1134 | – | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 1135 | – | sprintf(painCave.errMsg, | 
| 1136 | – | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 1137 | – | "\tOpenMD will use a default value of 15.0 angstroms" | 
| 1138 | – | "\tfor the cutoffRadius.\n"); | 
| 1139 | – | painCave.isFatal = 0; | 
| 1140 | – | simError(); | 
| 1141 | – | rcut_ = 15.0; | 
| 1142 | – |  | 
| 1143 | – | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1144 | – | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1145 | – | toUpper(myMethod); | 
| 1146 | – |  | 
| 1147 | – | // For the time being, we're tethering the LJ shifted behavior to the | 
| 1148 | – | // electrostaticSummationMethod keyword options | 
| 1149 | – | if (myMethod == "SHIFTED_POTENTIAL") { | 
| 1150 | – | ljsp_ = 1; | 
| 1151 | – | } else if (myMethod == "SHIFTED_FORCE") { | 
| 1152 | – | ljsf_ = 1; | 
| 1153 | – | } | 
| 1154 | – | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { | 
| 1155 | – | if (simParams_->haveSwitchingRadius()){ | 
| 1156 | – | sprintf(painCave.errMsg, | 
| 1157 | – | "SimInfo Warning: A value was set for the switchingRadius\n" | 
| 1158 | – | "\teven though the electrostaticSummationMethod was\n" | 
| 1159 | – | "\tset to %s\n", myMethod.c_str()); | 
| 1160 | – | painCave.isFatal = 1; | 
| 1161 | – | simError(); | 
| 1162 | – | } | 
| 1163 | – | } | 
| 1164 | – | } | 
| 1165 | – |  | 
| 1166 | – | if (simParams_->haveSwitchingRadius()){ | 
| 1167 | – | rsw_ = simParams_->getSwitchingRadius(); | 
| 1168 | – | } else { | 
| 1169 | – | sprintf(painCave.errMsg, | 
| 1170 | – | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 1171 | – | "\tOpenMD will use a default value of\n" | 
| 1172 | – | "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 1173 | – | painCave.isFatal = 0; | 
| 1174 | – | simError(); | 
| 1175 | – | rsw_ = 0.85 * rcut_; | 
| 1176 | – | } | 
| 1177 | – |  | 
| 1178 | – | Electrostatic::setElectrostaticCutoffRadius(rcut_, rsw_); | 
| 1179 | – | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); | 
| 1180 | – |  | 
| 1181 | – | } else { | 
| 1182 | – | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so | 
| 1183 | – | // We'll punt and let fortran figure out the cutoffs later. | 
| 1184 | – |  | 
| 1185 | – | notifyFortranYouAreOnYourOwn(); | 
| 1186 | – |  | 
| 1187 | – | } | 
| 1188 | – | } | 
| 991 |  | } | 
| 1190 | – |  | 
| 1191 | – | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | 
| 1192 | – |  | 
| 1193 | – | int errorOut; | 
| 1194 | – | ElectrostaticSummationMethod esm = NONE; | 
| 1195 | – | ElectrostaticScreeningMethod sm = UNDAMPED; | 
| 1196 | – | RealType alphaVal; | 
| 1197 | – | RealType dielectric; | 
| 1198 | – |  | 
| 1199 | – | errorOut = isError; | 
| 1200 | – |  | 
| 1201 | – | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1202 | – | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1203 | – | toUpper(myMethod); | 
| 1204 | – | if (myMethod == "NONE") { | 
| 1205 | – | esm = NONE; | 
| 1206 | – | } else { | 
| 1207 | – | if (myMethod == "SWITCHING_FUNCTION") { | 
| 1208 | – | esm = SWITCHING_FUNCTION; | 
| 1209 | – | } else { | 
| 1210 | – | if (myMethod == "SHIFTED_POTENTIAL") { | 
| 1211 | – | esm = SHIFTED_POTENTIAL; | 
| 1212 | – | } else { | 
| 1213 | – | if (myMethod == "SHIFTED_FORCE") { | 
| 1214 | – | esm = SHIFTED_FORCE; | 
| 1215 | – | } else { | 
| 1216 | – | if (myMethod == "REACTION_FIELD") { | 
| 1217 | – | esm = REACTION_FIELD; | 
| 1218 | – | dielectric = simParams_->getDielectric(); | 
| 1219 | – | if (!simParams_->haveDielectric()) { | 
| 1220 | – | // throw warning | 
| 1221 | – | sprintf( painCave.errMsg, | 
| 1222 | – | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" | 
| 1223 | – | "\tA default value of %f will be used for the dielectric.\n", dielectric); | 
| 1224 | – | painCave.isFatal = 0; | 
| 1225 | – | simError(); | 
| 1226 | – | } | 
| 1227 | – | } else { | 
| 1228 | – | // throw error | 
| 1229 | – | sprintf( painCave.errMsg, | 
| 1230 | – | "SimInfo error: Unknown electrostaticSummationMethod.\n" | 
| 1231 | – | "\t(Input file specified %s .)\n" | 
| 1232 | – | "\telectrostaticSummationMethod must be one of: \"none\",\n" | 
| 1233 | – | "\t\"shifted_potential\", \"shifted_force\", or \n" | 
| 1234 | – | "\t\"reaction_field\".\n", myMethod.c_str() ); | 
| 1235 | – | painCave.isFatal = 1; | 
| 1236 | – | simError(); | 
| 1237 | – | } | 
| 1238 | – | } | 
| 1239 | – | } | 
| 1240 | – | } | 
| 1241 | – | } | 
| 1242 | – | } | 
| 1243 | – |  | 
| 1244 | – | if (simParams_->haveElectrostaticScreeningMethod()) { | 
| 1245 | – | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); | 
| 1246 | – | toUpper(myScreen); | 
| 1247 | – | if (myScreen == "UNDAMPED") { | 
| 1248 | – | sm = UNDAMPED; | 
| 1249 | – | } else { | 
| 1250 | – | if (myScreen == "DAMPED") { | 
| 1251 | – | sm = DAMPED; | 
| 1252 | – | if (!simParams_->haveDampingAlpha()) { | 
| 1253 | – | // first set a cutoff dependent alpha value | 
| 1254 | – | // we assume alpha depends linearly with rcut from 0 to 20.5 ang | 
| 1255 | – | alphaVal = 0.5125 - rcut_* 0.025; | 
| 1256 | – | // for values rcut > 20.5, alpha is zero | 
| 1257 | – | if (alphaVal < 0) alphaVal = 0; | 
| 1258 | – |  | 
| 1259 | – | // throw warning | 
| 1260 | – | sprintf( painCave.errMsg, | 
| 1261 | – | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | 
| 1262 | – | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); | 
| 1263 | – | painCave.isFatal = 0; | 
| 1264 | – | simError(); | 
| 1265 | – | } else { | 
| 1266 | – | alphaVal = simParams_->getDampingAlpha(); | 
| 1267 | – | } | 
| 1268 | – |  | 
| 1269 | – | } else { | 
| 1270 | – | // throw error | 
| 1271 | – | sprintf( painCave.errMsg, | 
| 1272 | – | "SimInfo error: Unknown electrostaticScreeningMethod.\n" | 
| 1273 | – | "\t(Input file specified %s .)\n" | 
| 1274 | – | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" | 
| 1275 | – | "or \"damped\".\n", myScreen.c_str() ); | 
| 1276 | – | painCave.isFatal = 1; | 
| 1277 | – | simError(); | 
| 1278 | – | } | 
| 1279 | – | } | 
| 1280 | – | } | 
| 1281 | – |  | 
| 992 |  |  | 
| 1283 | – | Electrostatic::setElectrostaticSummationMethod( esm ); | 
| 1284 | – | Electrostatic::setElectrostaticScreeningMethod( sm ); | 
| 1285 | – | Electrostatic::setDampingAlpha( alphaVal ); | 
| 1286 | – | Electrostatic::setReactionFieldDielectric( dielectric ); | 
| 1287 | – | initFortranFF( &errorOut ); | 
| 1288 | – | } | 
| 993 |  |  | 
| 994 |  | void SimInfo::setupSwitchingFunction() { | 
| 995 |  | int ft = CUBIC; | 
| 996 | < |  | 
| 996 | > |  | 
| 997 |  | if (simParams_->haveSwitchingFunctionType()) { | 
| 998 | < | std::string funcType = simParams_->getSwitchingFunctionType(); | 
| 998 | > | string funcType = simParams_->getSwitchingFunctionType(); | 
| 999 |  | toUpper(funcType); | 
| 1000 |  | if (funcType == "CUBIC") { | 
| 1001 |  | ft = CUBIC; | 
| 1022 |  | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | 
| 1023 |  | if ( simParams_->haveAccumulateBoxDipole() ) | 
| 1024 |  | if ( simParams_->getAccumulateBoxDipole() ) { | 
| 1321 | – | setAccumulateBoxDipole(); | 
| 1025 |  | calcBoxDipole_ = true; | 
| 1026 |  | } | 
| 1027 |  |  | 
| 1031 |  | properties_.addProperty(genData); | 
| 1032 |  | } | 
| 1033 |  |  | 
| 1034 | < | void SimInfo::removeProperty(const std::string& propName) { | 
| 1034 | > | void SimInfo::removeProperty(const string& propName) { | 
| 1035 |  | properties_.removeProperty(propName); | 
| 1036 |  | } | 
| 1037 |  |  | 
| 1039 |  | properties_.clearProperties(); | 
| 1040 |  | } | 
| 1041 |  |  | 
| 1042 | < | std::vector<std::string> SimInfo::getPropertyNames() { | 
| 1042 | > | vector<string> SimInfo::getPropertyNames() { | 
| 1043 |  | return properties_.getPropertyNames(); | 
| 1044 |  | } | 
| 1045 |  |  | 
| 1046 | < | std::vector<GenericData*> SimInfo::getProperties() { | 
| 1046 | > | vector<GenericData*> SimInfo::getProperties() { | 
| 1047 |  | return properties_.getProperties(); | 
| 1048 |  | } | 
| 1049 |  |  | 
| 1050 | < | GenericData* SimInfo::getPropertyByName(const std::string& propName) { | 
| 1050 | > | GenericData* SimInfo::getPropertyByName(const string& propName) { | 
| 1051 |  | return properties_.getPropertyByName(propName); | 
| 1052 |  | } | 
| 1053 |  |  | 
| 1130 |  |  | 
| 1131 |  | } | 
| 1132 |  |  | 
| 1133 | < | std::ostream& operator <<(std::ostream& o, SimInfo& info) { | 
| 1133 | > | ostream& operator <<(ostream& o, SimInfo& info) { | 
| 1134 |  |  | 
| 1135 |  | return o; | 
| 1136 |  | } | 
| 1280 |  | return IOIndexToIntegrableObject.at(index); | 
| 1281 |  | } | 
| 1282 |  |  | 
| 1283 | < | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { | 
| 1283 | > | void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) { | 
| 1284 |  | IOIndexToIntegrableObject= v; | 
| 1285 |  | } | 
| 1286 |  |  | 
| 1322 |  | return; | 
| 1323 |  | } | 
| 1324 |  | /* | 
| 1325 | < | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { | 
| 1325 | > | void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) { | 
| 1326 |  | assert( v.size() == nAtoms_ + nRigidBodies_); | 
| 1327 |  | sdByGlobalIndex_ = v; | 
| 1328 |  | } | 
| 1332 |  | return sdByGlobalIndex_.at(index); | 
| 1333 |  | } | 
| 1334 |  | */ | 
| 1335 | + | int SimInfo::getNGlobalConstraints() { | 
| 1336 | + | int nGlobalConstraints; | 
| 1337 | + | #ifdef IS_MPI | 
| 1338 | + | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, | 
| 1339 | + | MPI_COMM_WORLD); | 
| 1340 | + | #else | 
| 1341 | + | nGlobalConstraints =  nConstraints_; | 
| 1342 | + | #endif | 
| 1343 | + | return nGlobalConstraints; | 
| 1344 | + | } | 
| 1345 | + |  | 
| 1346 |  | }//end namespace OpenMD | 
| 1347 |  |  |