| 1 | tim | 1710 | /* | 
| 2 |  |  | * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project | 
| 3 |  |  | * | 
| 4 |  |  | * Contact: oopse@oopse.org | 
| 5 |  |  | * | 
| 6 |  |  | * This program is free software; you can redistribute it and/or | 
| 7 |  |  | * modify it under the terms of the GNU Lesser General Public License | 
| 8 |  |  | * as published by the Free Software Foundation; either version 2.1 | 
| 9 |  |  | * of the License, or (at your option) any later version. | 
| 10 |  |  | * All we ask is that proper credit is given for our work, which includes | 
| 11 |  |  | * - but is not limited to - adding the above copyright notice to the beginning | 
| 12 |  |  | * of your source code files, and to any copyright notice that you may distribute | 
| 13 |  |  | * with programs based on this work. | 
| 14 |  |  | * | 
| 15 |  |  | * This program is distributed in the hope that it will be useful, | 
| 16 |  |  | * but WITHOUT ANY WARRANTY; without even the implied warranty of | 
| 17 |  |  | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the | 
| 18 |  |  | * GNU Lesser General Public License for more details. | 
| 19 |  |  | * | 
| 20 |  |  | * You should have received a copy of the GNU Lesser General Public License | 
| 21 |  |  | * along with this program; if not, write to the Free Software | 
| 22 |  |  | * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA. | 
| 23 |  |  | * | 
| 24 |  |  | */ | 
| 25 | gezelter | 1490 |  | 
| 26 | tim | 1710 | /** | 
| 27 |  |  | * @file SimInfo.cpp | 
| 28 |  |  | * @author    tlin | 
| 29 |  |  | * @date  11/02/2004 | 
| 30 |  |  | * @version 1.0 | 
| 31 |  |  | */ | 
| 32 | gezelter | 1490 |  | 
| 33 | tim | 1710 | #include <algorithm> | 
| 34 |  |  |  | 
| 35 | tim | 1492 | #include "brains/SimInfo.hpp" | 
| 36 | tim | 1710 | #include "utils/MemoryUtils.hpp" | 
| 37 | gezelter | 1490 |  | 
| 38 | tim | 1710 | namespace oopse { | 
| 39 | gezelter | 1490 |  | 
| 40 | tim | 1733 | SimInfo::SimInfo(const std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, | 
| 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) { | 
| 45 | gezelter | 1490 |  | 
| 46 | tim | 1733 | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; | 
| 47 |  |  | int nCutoffAtoms; // number of atoms belong to cutoff groups | 
| 48 |  |  | int ngroups;          //total cutoff groups defined in meta-data file | 
| 49 |  |  | MoleculeStamp* molStamp; | 
| 50 |  |  | int nMolWithSameStamp; | 
| 51 |  |  | CutoffGroupStamp* cgStamp; | 
| 52 |  |  | int nAtomsIngroups; | 
| 53 |  |  | int nCutoffGroupsInStamp; | 
| 54 |  |  |  | 
| 55 |  |  | nGlobalAtoms_ =  0; | 
| 56 |  |  | ngroups = 0; | 
| 57 |  |  |  | 
| 58 |  |  | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { | 
| 59 |  |  | molStamp = i->first; | 
| 60 |  |  | nMolWithSameStamp = i->second; | 
| 61 |  |  |  | 
| 62 |  |  | addMoleculeStamp(molStamp, nMolWithSameStamp); | 
| 63 |  |  |  | 
| 64 |  |  | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | 
| 65 |  |  |  | 
| 66 |  |  | nAtomsIngroups = 0; | 
| 67 |  |  | nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | 
| 68 |  |  |  | 
| 69 |  |  | for (int j=0; j < nCutoffGroupsInStamp; j++) { | 
| 70 |  |  | cgStamp = molStamp->getCutoffGroup(j); | 
| 71 |  |  | nAtomsIngroups += cgStamp->getNMembers(); | 
| 72 |  |  | } | 
| 73 |  |  |  | 
| 74 |  |  | ngroups += *nMolWithSameStamp; | 
| 75 |  |  | nCutoffAtoms += nAtomsIngroups * nMolWithSameStamp; | 
| 76 |  |  | } | 
| 77 |  |  |  | 
| 78 |  |  | //every free atom (atom does not belong to cutoff groups) is a cutoff group | 
| 79 |  |  | //therefore the total number of cutoff groups in the system is equal to | 
| 80 |  |  | //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data | 
| 81 |  |  | //file plus the number of cutoff groups defined in meta-data file | 
| 82 |  |  | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + ngroups; | 
| 83 |  |  |  | 
| 84 |  |  | //initialize globalGroupMembership_, every element of this array will be 0 | 
| 85 |  |  | globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0); | 
| 86 |  |  |  | 
| 87 |  |  | nGlobalMols_ = molStampIds_.size(); | 
| 88 |  |  |  | 
| 89 |  |  | #ifdef IS_MPI | 
| 90 |  |  | molToProcMap_.resize(nGlobalMols_); | 
| 91 |  |  | #endif | 
| 92 |  |  |  | 
| 93 | gezelter | 1490 | } | 
| 94 |  |  |  | 
| 95 | tim | 1710 | SimInfo::~SimInfo() { | 
| 96 | tim | 1726 | //MemoryUtils::deleteVectorOfPointer(molecules_); | 
| 97 | tim | 1733 |  | 
| 98 |  |  | MemoryUtils::deleteVectorOfPointer(moleculeStamps_); | 
| 99 |  |  |  | 
| 100 | tim | 1710 | delete sman_; | 
| 101 | tim | 1733 | delete globals_; | 
| 102 |  |  | delete forceField_; | 
| 103 | gezelter | 1490 |  | 
| 104 |  |  | } | 
| 105 |  |  |  | 
| 106 |  |  |  | 
| 107 | tim | 1710 | bool SimInfo::addMolecule(Molecule* mol) { | 
| 108 | tim | 1726 | MoleculeIterator i; | 
| 109 | tim | 1733 |  | 
| 110 |  |  | i = molecules_.find(mol->getGlobalIndex()); | 
| 111 | tim | 1710 | if (i != molecules_.end() ) { | 
| 112 | gezelter | 1490 |  | 
| 113 | tim | 1726 | molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); | 
| 114 |  |  |  | 
| 115 | tim | 1710 | nAtoms_ += mol->getNAtoms(); | 
| 116 |  |  | nBonds_ += mol->getNBonds(); | 
| 117 |  |  | nBends_ += mol->getNBends(); | 
| 118 |  |  | nTorsions_ += mol->getNTorsions(); | 
| 119 |  |  | nRigidBodies_ += mol->getNRigidBodies(); | 
| 120 |  |  | nIntegrableObjects_ += mol->getNIntegrableObjects(); | 
| 121 |  |  | nCutoffGroups_ += mol->getNCutoffGroups(); | 
| 122 |  |  | nConstraints_ += mol->getNConstraints(); | 
| 123 | gezelter | 1490 |  | 
| 124 | tim | 1710 | return true; | 
| 125 |  |  | } else { | 
| 126 |  |  | return false; | 
| 127 | gezelter | 1490 | } | 
| 128 |  |  | } | 
| 129 |  |  |  | 
| 130 | tim | 1710 | bool SimInfo::removeMolecule(Molecule* mol) { | 
| 131 | tim | 1726 | MoleculeIterator i; | 
| 132 | tim | 1733 | i = molecules_.find(mol->getGlobalIndex()); | 
| 133 | gezelter | 1490 |  | 
| 134 | tim | 1710 | if (i != molecules_.end() ) { | 
| 135 | tim | 1726 |  | 
| 136 | tim | 1733 | assert(mol == i->second); | 
| 137 |  |  |  | 
| 138 | tim | 1710 | nAtoms_ -= mol->getNAtoms(); | 
| 139 |  |  | nBonds_ -= mol->getNBonds(); | 
| 140 |  |  | nBends_ -= mol->getNBends(); | 
| 141 |  |  | nTorsions_ -= mol->getNTorsions(); | 
| 142 |  |  | nRigidBodies_ -= mol->getNRigidBodies(); | 
| 143 |  |  | nIntegrableObjects_ -= mol->getNIntegrableObjects(); | 
| 144 |  |  | nCutoffGroups_ -= mol->getNCutoffGroups(); | 
| 145 |  |  | nConstraints_ -= mol->getNConstraints(); | 
| 146 | gezelter | 1490 |  | 
| 147 | tim | 1726 | molecules_.erase(mol->getGlobalIndex()); | 
| 148 |  |  |  | 
| 149 |  |  | delete mol; | 
| 150 |  |  |  | 
| 151 | tim | 1710 | return true; | 
| 152 |  |  | } else { | 
| 153 |  |  | return false; | 
| 154 | gezelter | 1490 | } | 
| 155 |  |  |  | 
| 156 |  |  |  | 
| 157 | tim | 1710 | } | 
| 158 | gezelter | 1490 |  | 
| 159 | tim | 1710 |  | 
| 160 | tim | 1726 | Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { | 
| 161 | tim | 1710 | i = molecules_.begin(); | 
| 162 |  |  | return i == molecules_.end() ? NULL : *i; | 
| 163 |  |  | } | 
| 164 | gezelter | 1490 |  | 
| 165 | tim | 1726 | Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { | 
| 166 | tim | 1710 | ++i; | 
| 167 |  |  | return i == molecules_.end() ? NULL : *i; | 
| 168 | gezelter | 1490 | } | 
| 169 |  |  |  | 
| 170 |  |  |  | 
| 171 | tim | 1722 | void SimInfo::calcNdf() { | 
| 172 | tim | 1710 | int ndf_local; | 
| 173 | tim | 1726 | MoleculeIterator i; | 
| 174 | tim | 1710 | std::vector<StuntDouble*>::iterator j; | 
| 175 |  |  | Molecule* mol; | 
| 176 |  |  | StuntDouble* integrableObject; | 
| 177 | gezelter | 1490 |  | 
| 178 | tim | 1710 | ndf_local = 0; | 
| 179 | gezelter | 1490 |  | 
| 180 | tim | 1710 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 181 |  |  | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; | 
| 182 |  |  | integrableObject = mol->nextIntegrableObject(j)) { | 
| 183 | gezelter | 1490 |  | 
| 184 | tim | 1710 | ndf_local += 3; | 
| 185 | gezelter | 1490 |  | 
| 186 | tim | 1710 | if (integrableObject->isDirectional()) { | 
| 187 |  |  | if (integrableObject->isLinear()) { | 
| 188 |  |  | ndf_local += 2; | 
| 189 |  |  | } else { | 
| 190 |  |  | ndf_local += 3; | 
| 191 |  |  | } | 
| 192 |  |  | } | 
| 193 |  |  |  | 
| 194 |  |  | }//end for (integrableObject) | 
| 195 |  |  | }// end for (mol) | 
| 196 | gezelter | 1490 |  | 
| 197 | tim | 1710 | // n_constraints is local, so subtract them on each processor | 
| 198 |  |  | ndf_local -= nConstraints_; | 
| 199 | gezelter | 1490 |  | 
| 200 |  |  | #ifdef IS_MPI | 
| 201 | tim | 1710 | MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 202 | gezelter | 1490 | #else | 
| 203 | tim | 1710 | ndf_ = ndf_local; | 
| 204 | gezelter | 1490 | #endif | 
| 205 |  |  |  | 
| 206 | tim | 1733 | // nZconstraints_ is global, as are the 3 COM translations for the | 
| 207 | tim | 1710 | // entire system: | 
| 208 | tim | 1733 | ndf_ = ndf_ - 3 - nZconstraints_; | 
| 209 | gezelter | 1490 |  | 
| 210 |  |  | } | 
| 211 |  |  |  | 
| 212 | tim | 1722 | void SimInfo::calcNdfRaw() { | 
| 213 | tim | 1710 | int ndfRaw_local; | 
| 214 | gezelter | 1490 |  | 
| 215 | tim | 1726 | MoleculeIterator i; | 
| 216 | tim | 1710 | std::vector<StuntDouble*>::iterator j; | 
| 217 |  |  | Molecule* mol; | 
| 218 |  |  | StuntDouble* integrableObject; | 
| 219 | gezelter | 1490 |  | 
| 220 | tim | 1710 | // Raw degrees of freedom that we have to set | 
| 221 |  |  | ndfRaw_local = 0; | 
| 222 | gezelter | 1490 |  | 
| 223 | tim | 1710 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 224 |  |  | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; | 
| 225 |  |  | integrableObject = mol->nextIntegrableObject(j)) { | 
| 226 | gezelter | 1490 |  | 
| 227 | tim | 1710 | ndfRaw_local += 3; | 
| 228 | gezelter | 1490 |  | 
| 229 | tim | 1710 | if (integrableObject->isDirectional()) { | 
| 230 |  |  | if (integrableObject->isLinear()) { | 
| 231 |  |  | ndfRaw_local += 2; | 
| 232 |  |  | } else { | 
| 233 |  |  | ndfRaw_local += 3; | 
| 234 |  |  | } | 
| 235 |  |  | } | 
| 236 |  |  |  | 
| 237 |  |  | } | 
| 238 |  |  | } | 
| 239 |  |  |  | 
| 240 | gezelter | 1490 | #ifdef IS_MPI | 
| 241 | tim | 1710 | MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 242 | gezelter | 1490 | #else | 
| 243 | tim | 1710 | ndfRaw_ = ndfRaw_local; | 
| 244 | gezelter | 1490 | #endif | 
| 245 |  |  | } | 
| 246 |  |  |  | 
| 247 | tim | 1722 | void SimInfo::calcNdfTrans() { | 
| 248 | tim | 1710 | int ndfTrans_local; | 
| 249 | gezelter | 1490 |  | 
| 250 | tim | 1710 | ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_; | 
| 251 | gezelter | 1490 |  | 
| 252 |  |  |  | 
| 253 |  |  | #ifdef IS_MPI | 
| 254 | tim | 1710 | MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 255 | gezelter | 1490 | #else | 
| 256 | tim | 1710 | ndfTrans_ = ndfTrans_local; | 
| 257 | gezelter | 1490 | #endif | 
| 258 |  |  |  | 
| 259 | tim | 1733 | ndfTrans_ = ndfTrans_ - 3 - nZconstraints_; | 
| 260 | gezelter | 1490 |  | 
| 261 |  |  | } | 
| 262 |  |  |  | 
| 263 | tim | 1719 | void SimInfo::addExcludePairs(Molecule* mol) { | 
| 264 |  |  | std::vector<Bond*>::iterator bondIter; | 
| 265 |  |  | std::vector<Bend*>::iterator bendIter; | 
| 266 |  |  | std::vector<Torsion*>::iterator torsionIter; | 
| 267 |  |  | Bond* bond; | 
| 268 |  |  | Bend* bend; | 
| 269 |  |  | Torsion* torsion; | 
| 270 |  |  | int a; | 
| 271 |  |  | int b; | 
| 272 |  |  | int c; | 
| 273 |  |  | int d; | 
| 274 |  |  |  | 
| 275 |  |  | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 276 |  |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 277 |  |  | b = bond->getAtomB()->getGlobalIndex(); | 
| 278 |  |  | exclude_.addPair(a, b); | 
| 279 |  |  | } | 
| 280 | gezelter | 1490 |  | 
| 281 | tim | 1719 | for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { | 
| 282 |  |  | a = bend->getAtomA()->getGlobalIndex(); | 
| 283 |  |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 284 |  |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 285 |  |  |  | 
| 286 |  |  | exclude_.addPair(a, b); | 
| 287 |  |  | exclude_.addPair(a, c); | 
| 288 |  |  | exclude_.addPair(b, c); | 
| 289 |  |  | } | 
| 290 |  |  |  | 
| 291 |  |  | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) { | 
| 292 |  |  | a = torsion->getAtomA()->getGlobalIndex(); | 
| 293 |  |  | b = torsion->getAtomB()->getGlobalIndex(); | 
| 294 |  |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 295 |  |  | d = torsion->getAtomD()->getGlobalIndex(); | 
| 296 |  |  |  | 
| 297 |  |  | exclude_.addPair(a, b); | 
| 298 |  |  | exclude_.addPair(a, c); | 
| 299 |  |  | exclude_.addPair(a, d); | 
| 300 |  |  | exclude_.addPair(b, c); | 
| 301 |  |  | exclude_.addPair(b, d); | 
| 302 |  |  | exclude_.addPair(c, d); | 
| 303 |  |  | } | 
| 304 |  |  |  | 
| 305 |  |  |  | 
| 306 |  |  | } | 
| 307 |  |  |  | 
| 308 |  |  | void SimInfo::removeExcludePairs(Molecule* mol) { | 
| 309 |  |  | std::vector<Bond*>::iterator bondIter; | 
| 310 |  |  | std::vector<Bend*>::iterator bendIter; | 
| 311 |  |  | std::vector<Torsion*>::iterator torsionIter; | 
| 312 |  |  | Bond* bond; | 
| 313 |  |  | Bend* bend; | 
| 314 |  |  | Torsion* torsion; | 
| 315 |  |  | int a; | 
| 316 |  |  | int b; | 
| 317 |  |  | int c; | 
| 318 |  |  | int d; | 
| 319 |  |  |  | 
| 320 |  |  | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 321 |  |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 322 |  |  | b = bond->getAtomB()->getGlobalIndex(); | 
| 323 |  |  | exclude_.removePair(a, b); | 
| 324 |  |  | } | 
| 325 |  |  |  | 
| 326 |  |  | for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { | 
| 327 |  |  | a = bend->getAtomA()->getGlobalIndex(); | 
| 328 |  |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 329 |  |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 330 |  |  |  | 
| 331 |  |  | exclude_.removePair(a, b); | 
| 332 |  |  | exclude_.removePair(a, c); | 
| 333 |  |  | exclude_.removePair(b, c); | 
| 334 |  |  | } | 
| 335 |  |  |  | 
| 336 |  |  | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) { | 
| 337 |  |  | a = torsion->getAtomA()->getGlobalIndex(); | 
| 338 |  |  | b = torsion->getAtomB()->getGlobalIndex(); | 
| 339 |  |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 340 |  |  | d = torsion->getAtomD()->getGlobalIndex(); | 
| 341 |  |  |  | 
| 342 |  |  | exclude_.removePair(a, b); | 
| 343 |  |  | exclude_.removePair(a, c); | 
| 344 |  |  | exclude_.removePair(a, d); | 
| 345 |  |  | exclude_.removePair(b, c); | 
| 346 |  |  | exclude_.removePair(b, d); | 
| 347 |  |  | exclude_.removePair(c, d); | 
| 348 |  |  | } | 
| 349 |  |  |  | 
| 350 |  |  | } | 
| 351 |  |  |  | 
| 352 |  |  |  | 
| 353 | tim | 1725 | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { | 
| 354 |  |  | int curStampId; | 
| 355 |  |  |  | 
| 356 |  |  | //index from 0 | 
| 357 |  |  | curStampId = molStampIds_.size(); | 
| 358 |  |  |  | 
| 359 |  |  | moleculeStamps_.push_back(molStamp); | 
| 360 |  |  | molStampIds_.insert(molStampIds_.end(), nmol, curStampId) | 
| 361 |  |  | } | 
| 362 |  |  |  | 
| 363 | tim | 1727 | void SimInfo::update() { | 
| 364 |  |  |  | 
| 365 | tim | 1733 | #ifdef IS_MPI | 
| 366 |  |  | setupFortranParallel(); | 
| 367 |  |  | #endif | 
| 368 | tim | 1727 |  | 
| 369 | tim | 1733 | setupFortranSim(); | 
| 370 |  |  |  | 
| 371 |  |  | calcNdf(); | 
| 372 |  |  | calcNdfRaw(); | 
| 373 |  |  | calcNdfTrans(); | 
| 374 |  |  | } | 
| 375 |  |  |  | 
| 376 |  |  | std::set<AtomType*> SimInfo::getUniqueAtomTypes() { | 
| 377 |  |  | typename SimInfo::MoleculeIterator mi; | 
| 378 |  |  | Molecule* mol; | 
| 379 |  |  | typename Molecule::AtomIterator ai; | 
| 380 |  |  | Atom* atom; | 
| 381 |  |  | std::set<AtomType*> atomTypes; | 
| 382 |  |  |  | 
| 383 |  |  | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | 
| 384 |  |  |  | 
| 385 |  |  | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 386 |  |  | atomTypes.insert(atom->getAtomType()); | 
| 387 |  |  | } | 
| 388 |  |  |  | 
| 389 |  |  | } | 
| 390 |  |  |  | 
| 391 |  |  | return atomTypes; | 
| 392 |  |  | } | 
| 393 |  |  |  | 
| 394 |  |  | void SimInfo::setupSimType() { | 
| 395 |  |  | std::set<AtomType*>::iterator i; | 
| 396 |  |  | std::set<AtomType*> atomTypes; | 
| 397 |  |  | atomTypes = getUniqueAtomTypes(); | 
| 398 | tim | 1727 |  | 
| 399 | tim | 1733 | int useLennardJones = 0; | 
| 400 |  |  | int useElectrostatic = 0; | 
| 401 |  |  | int useEAM = 0; | 
| 402 |  |  | int useCharge = 0; | 
| 403 |  |  | int useDirectional = 0; | 
| 404 |  |  | int useDipole = 0; | 
| 405 |  |  | int useGayBerne = 0; | 
| 406 |  |  | int useSticky = 0; | 
| 407 |  |  | int useShape = 0; | 
| 408 |  |  | int useFLARB = 0; //it is not in AtomType yet | 
| 409 |  |  | int useDirectionalAtom = 0; | 
| 410 |  |  | int useElectrostatics = 0; | 
| 411 |  |  | //usePBC and useRF are from globals | 
| 412 |  |  | bool usePBC = globals_->getPBC(); | 
| 413 |  |  | bool useRF = globals_->getUseRF(); | 
| 414 |  |  |  | 
| 415 |  |  | //loop over all of the atom types | 
| 416 |  |  | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 417 |  |  | useLennardJones |= i->isLennardJones(); | 
| 418 |  |  | useElectrostatic |= i->isElectrostatic(); | 
| 419 |  |  | useEAM |= i->isEAM(); | 
| 420 |  |  | useCharge |= i->isCharge(); | 
| 421 |  |  | useDirectional |= i->isDirectional(); | 
| 422 |  |  | useDipole |= i->isDipole(); | 
| 423 |  |  | useGayBerne |= i->isGayBerne(); | 
| 424 |  |  | useSticky |= i->isSticky(); | 
| 425 |  |  | useShape |= i->isShape(); | 
| 426 |  |  | } | 
| 427 |  |  |  | 
| 428 |  |  | if (useSticky || useDipole || useGayBerne || useShape) { | 
| 429 |  |  | useDirectionalAtom = 1; | 
| 430 |  |  | } | 
| 431 |  |  |  | 
| 432 |  |  | if (useCharge || useDipole) { | 
| 433 |  |  | useElectrostatics = 1; | 
| 434 |  |  | } | 
| 435 |  |  |  | 
| 436 |  |  | #ifdef IS_MPI | 
| 437 |  |  | int temp; | 
| 438 |  |  |  | 
| 439 |  |  | temp = usePBC; | 
| 440 |  |  | MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 441 |  |  |  | 
| 442 |  |  | temp = useDirectionalAtom; | 
| 443 |  |  | MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 444 |  |  |  | 
| 445 |  |  | temp = useLennardJones; | 
| 446 |  |  | MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 447 |  |  |  | 
| 448 |  |  | temp = useElectrostatics; | 
| 449 |  |  | MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 450 |  |  |  | 
| 451 |  |  | temp = useCharge; | 
| 452 |  |  | MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 453 |  |  |  | 
| 454 |  |  | temp = useDipole; | 
| 455 |  |  | MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 456 |  |  |  | 
| 457 |  |  | temp = useSticky; | 
| 458 |  |  | MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 459 |  |  |  | 
| 460 |  |  | temp = useGayBerne; | 
| 461 |  |  | MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 462 |  |  |  | 
| 463 |  |  | temp = useEAM; | 
| 464 |  |  | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 465 |  |  |  | 
| 466 |  |  | temp = useShape; | 
| 467 |  |  | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 468 |  |  |  | 
| 469 |  |  | temp = useFLARB; | 
| 470 |  |  | MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 471 |  |  |  | 
| 472 |  |  | temp = useRF; | 
| 473 |  |  | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 474 |  |  |  | 
| 475 |  |  | #endif | 
| 476 |  |  |  | 
| 477 |  |  | fInfo_.SIM_uses_PBC = usePBC; | 
| 478 |  |  | fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; | 
| 479 |  |  | fInfo_.SIM_uses_LennardJones = useLennardJones; | 
| 480 |  |  | fInfo_.SIM_uses_Electrostatics = useElectrostatics; | 
| 481 |  |  | fInfo_.SIM_uses_Charges = useCharge; | 
| 482 |  |  | fInfo_.SIM_uses_Dipoles = useDipole; | 
| 483 |  |  | fInfo_.SIM_uses_Sticky = useSticky; | 
| 484 |  |  | fInfo_.SIM_uses_GayBerne = useGayBerne; | 
| 485 |  |  | fInfo_.SIM_uses_EAM = useEAM; | 
| 486 |  |  | fInfo_.SIM_uses_Shapes = useShape; | 
| 487 |  |  | fInfo_.SIM_uses_FLARB = useFLARB; | 
| 488 |  |  | fInfo_.SIM_uses_RF = useRF; | 
| 489 |  |  |  | 
| 490 |  |  | if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { | 
| 491 |  |  | fInfo_.dielect = dielectric; | 
| 492 |  |  | } else { | 
| 493 |  |  | fInfo_.dielect = 0.0; | 
| 494 |  |  | } | 
| 495 |  |  |  | 
| 496 |  |  | } | 
| 497 |  |  |  | 
| 498 |  |  | void SimInfo::setupFortranSim() { | 
| 499 |  |  | int isError; | 
| 500 |  |  | int nExclude; | 
| 501 |  |  | std::vector<int> fortranGlobalGroupMembership; | 
| 502 |  |  |  | 
| 503 |  |  | nExclude = exclude_.getSize(); | 
| 504 |  |  | isError = 0; | 
| 505 |  |  |  | 
| 506 |  |  | //globalGroupMembership_ is filled by SimCreator | 
| 507 |  |  | for (int i = 0; i < nGlobalAtoms_; i++) { | 
| 508 |  |  | fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); | 
| 509 |  |  | } | 
| 510 |  |  |  | 
| 511 |  |  | //calculate mass ratio of cutoff group | 
| 512 |  |  | std::vector<double> mfact; | 
| 513 |  |  | typename SimInfo::MoleculeIterator mi; | 
| 514 |  |  | Molecule* mol; | 
| 515 |  |  | typename Molecule::CutoffGroupIterator ci; | 
| 516 |  |  | CutoffGroup* cg; | 
| 517 |  |  | typename Molecule::AtomIterator ai; | 
| 518 |  |  | Atom* atom; | 
| 519 |  |  | double totalMass; | 
| 520 |  |  |  | 
| 521 |  |  | //to avoid memory reallocation, reserve enough space for mfact | 
| 522 |  |  | mfact.reserve(getNCutoffGroups()); | 
| 523 |  |  |  | 
| 524 |  |  | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | 
| 525 |  |  | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { | 
| 526 |  |  |  | 
| 527 |  |  | totalMass = cg->getMass(); | 
| 528 |  |  | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | 
| 529 |  |  | mfact.push_back(atom->getMass()/totalMass); | 
| 530 |  |  | } | 
| 531 |  |  |  | 
| 532 |  |  | } | 
| 533 |  |  | } | 
| 534 |  |  |  | 
| 535 |  |  | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) | 
| 536 |  |  | std::vector<int> identArray; | 
| 537 |  |  |  | 
| 538 |  |  | //to avoid memory reallocation, reserve enough space identArray | 
| 539 |  |  | identArray.reserve(getNAtoms()); | 
| 540 |  |  |  | 
| 541 |  |  | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | 
| 542 |  |  | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 543 |  |  | identArray.push_back(atom->getIdent()); | 
| 544 |  |  | } | 
| 545 |  |  | } | 
| 546 |  |  |  | 
| 547 |  |  | //setup fortran simulation | 
| 548 |  |  | //gloalExcludes and molMembershipArray should go away (They are never used) | 
| 549 |  |  | //why the hell fortran need to know molecule? | 
| 550 |  |  | //OOPSE = Object-Obfuscated Parallel Simulation Engine | 
| 551 |  |  |  | 
| 552 |  |  | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, exclude_->getExcludeList(), | 
| 553 |  |  | &nGlobalExcludes, globalExcludes, molMembershipArray, | 
| 554 |  |  | &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError); | 
| 555 |  |  |  | 
| 556 |  |  | if( isError ){ | 
| 557 |  |  |  | 
| 558 |  |  | sprintf( painCave.errMsg, | 
| 559 |  |  | "There was an error setting the simulation information in fortran.\n" ); | 
| 560 |  |  | painCave.isFatal = 1; | 
| 561 |  |  | painCave.severity = OOPSE_ERROR; | 
| 562 |  |  | simError(); | 
| 563 |  |  | } | 
| 564 |  |  |  | 
| 565 |  |  | #ifdef IS_MPI | 
| 566 |  |  | sprintf( checkPointMsg, | 
| 567 |  |  | "succesfully sent the simulation information to fortran.\n"); | 
| 568 |  |  | MPIcheckPoint(); | 
| 569 |  |  | #endif // is_mpi | 
| 570 |  |  | } | 
| 571 |  |  |  | 
| 572 |  |  |  | 
| 573 |  |  | #ifdef IS_MPI | 
| 574 |  |  | void SimInfo::setupFortranParallel() { | 
| 575 |  |  |  | 
| 576 | tim | 1727 | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex | 
| 577 |  |  | std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0); | 
| 578 |  |  | std::vector<int> localToGlobalCutoffGroupIndex; | 
| 579 |  |  | typename SimInfo::MoleculeIterator mi; | 
| 580 |  |  | typename Molecule::AtomIterator ai; | 
| 581 |  |  | typename Molecule::CutoffGroupIterator ci; | 
| 582 |  |  | Molecule* mol; | 
| 583 |  |  | Atom* atom; | 
| 584 |  |  | CutoffGroup* cg; | 
| 585 |  |  | mpiSimData parallelData; | 
| 586 |  |  | int isError; | 
| 587 |  |  |  | 
| 588 |  |  | for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) { | 
| 589 |  |  |  | 
| 590 |  |  | //local index(index in DataStorge) of atom is important | 
| 591 |  |  | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 592 |  |  | localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; | 
| 593 |  |  | } | 
| 594 |  |  |  | 
| 595 |  |  | //local index of cutoff group is trivial, it only depends on the order of travesing | 
| 596 |  |  | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { | 
| 597 |  |  | localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); | 
| 598 |  |  | } | 
| 599 |  |  |  | 
| 600 |  |  | } | 
| 601 |  |  |  | 
| 602 | tim | 1733 | //fill up mpiSimData struct | 
| 603 |  |  | parallelData.nMolGlobal = getNGlobalMolecules(); | 
| 604 |  |  | parallelData.nMolLocal = getNMolecules(); | 
| 605 |  |  | parallelData.nAtomsGlobal = getNGlobalAtoms(); | 
| 606 |  |  | parallelData.nAtomsLocal = getNAtoms(); | 
| 607 |  |  | parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); | 
| 608 |  |  | parallelData.nGroupsLocal = getNCutoffGroups(); | 
| 609 | tim | 1727 | parallelData.myNode = worldRank; | 
| 610 |  |  | MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors)); | 
| 611 | tim | 1733 |  | 
| 612 |  |  | //pass mpiSimData struct and index arrays to fortran | 
| 613 |  |  | setFsimParallel(parallelData, &(parallelData->nAtomsLocal), | 
| 614 | tim | 1727 | &localToGlobalAtomIndex[0],  &(parallelData->nGroupsLocal), | 
| 615 |  |  | &localToGlobalCutoffGroupIndex[0], &isError); | 
| 616 |  |  |  | 
| 617 |  |  | if (isError) { | 
| 618 |  |  | sprintf(painCave.errMsg, | 
| 619 |  |  | "mpiRefresh errror: fortran didn't like something we gave it.\n"); | 
| 620 |  |  | painCave.isFatal = 1; | 
| 621 |  |  | simError(); | 
| 622 |  |  | } | 
| 623 |  |  |  | 
| 624 |  |  | sprintf(checkPointMsg, " mpiRefresh successful.\n"); | 
| 625 |  |  | MPIcheckPoint(); | 
| 626 |  |  |  | 
| 627 |  |  |  | 
| 628 |  |  | } | 
| 629 |  |  |  | 
| 630 | tim | 1733 | #endif | 
| 631 |  |  |  | 
| 632 |  |  | void SimInfo::addProperty(GenericData* genData) { | 
| 633 |  |  | properties_.addProperty(genData); | 
| 634 |  |  | } | 
| 635 |  |  |  | 
| 636 |  |  | void SimInfo::removeProperty(const std::string& propName) { | 
| 637 |  |  | properties_.removeProperty(propName); | 
| 638 |  |  | } | 
| 639 |  |  |  | 
| 640 |  |  | void SimInfo::clearProperties() { | 
| 641 |  |  | properties_.clearProperties(); | 
| 642 |  |  | } | 
| 643 |  |  |  | 
| 644 |  |  | std::vector<std::string> SimInfo::getPropertyNames() { | 
| 645 |  |  | return properties_.getPropertyNames(); | 
| 646 |  |  | } | 
| 647 |  |  |  | 
| 648 |  |  | std::vector<GenericData*> SimInfo::getProperties() { | 
| 649 |  |  | return properties_.getProperties(); | 
| 650 |  |  | } | 
| 651 |  |  |  | 
| 652 |  |  | GenericData* SimInfo::getPropertyByName(const std::string& propName) { | 
| 653 |  |  | return properties_.getPropertyByName(propName); | 
| 654 |  |  | } | 
| 655 |  |  |  | 
| 656 |  |  |  | 
| 657 | tim | 1725 | std::ostream& operator <<(ostream& o, SimInfo& info) { | 
| 658 |  |  |  | 
| 659 |  |  | return o; | 
| 660 |  |  | } | 
| 661 |  |  |  | 
| 662 | tim | 1710 | }//end namespace oopse |