| 6 | 
  | 
 * redistribute this software in source and binary code form, provided | 
| 7 | 
  | 
 * that the following conditions are met: | 
| 8 | 
  | 
 * | 
| 9 | 
< | 
 * 1. Acknowledgement of the program authors must be made in any | 
| 10 | 
< | 
 *    publication of scientific results based in part on use of the | 
| 11 | 
< | 
 *    program.  An acceptable form of acknowledgement is citation of | 
| 12 | 
< | 
 *    the article in which the program was described (Matthew | 
| 13 | 
< | 
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher | 
| 14 | 
< | 
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented | 
| 15 | 
< | 
 *    Parallel Simulation Engine for Molecular Dynamics," | 
| 16 | 
< | 
 *    J. Comput. Chem. 26, pp. 252-271 (2005)) | 
| 17 | 
< | 
 * | 
| 18 | 
< | 
 * 2. Redistributions of source code must retain the above copyright | 
| 9 | 
> | 
 * 1. Redistributions of source code must retain the above copyright | 
| 10 | 
  | 
 *    notice, this list of conditions and the following disclaimer. | 
| 11 | 
  | 
 * | 
| 12 | 
< | 
 * 3. Redistributions in binary form must reproduce the above copyright | 
| 12 | 
> | 
 * 2. Redistributions in binary form must reproduce the above copyright | 
| 13 | 
  | 
 *    notice, this list of conditions and the following disclaimer in the | 
| 14 | 
  | 
 *    documentation and/or other materials provided with the | 
| 15 | 
  | 
 *    distribution. | 
| 28 | 
  | 
 * arising out of the use of or inability to use software, even if the | 
| 29 | 
  | 
 * University of Notre Dame has been advised of the possibility of | 
| 30 | 
  | 
 * such damages. | 
| 31 | 
+ | 
 * | 
| 32 | 
+ | 
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your | 
| 33 | 
+ | 
 * research, please cite the appropriate papers when you publish your | 
| 34 | 
+ | 
 * work.  Good starting points are: | 
| 35 | 
+ | 
 *                                                                       | 
| 36 | 
+ | 
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).              | 
| 37 | 
+ | 
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).           | 
| 38 | 
+ | 
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).           | 
| 39 | 
+ | 
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 | 
+ | 
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 | 
  | 
 */ | 
| 42 | 
  | 
  | 
| 43 | 
  | 
/** | 
| 44 | 
  | 
 * @file SimCreator.cpp | 
| 45 | 
  | 
 * @author tlin | 
| 46 | 
  | 
 * @date 11/03/2004 | 
| 46 | 
– | 
 * @time 13:51am | 
| 47 | 
  | 
 * @version 1.0 | 
| 48 | 
  | 
 */ | 
| 49 | 
  | 
 | 
| 50 | 
+ | 
#ifdef IS_MPI | 
| 51 | 
+ | 
#include "mpi.h" | 
| 52 | 
+ | 
#include "math/ParallelRandNumGen.hpp" | 
| 53 | 
+ | 
#endif | 
| 54 | 
+ | 
 | 
| 55 | 
+ | 
#include <exception> | 
| 56 | 
+ | 
#include <iostream> | 
| 57 | 
+ | 
#include <sstream> | 
| 58 | 
+ | 
#include <string> | 
| 59 | 
+ | 
 | 
| 60 | 
  | 
#include "brains/MoleculeCreator.hpp" | 
| 61 | 
  | 
#include "brains/SimCreator.hpp" | 
| 62 | 
  | 
#include "brains/SimSnapshotManager.hpp" | 
| 63 | 
  | 
#include "io/DumpReader.hpp" | 
| 64 | 
< | 
#include "io/parse_me.h" | 
| 55 | 
< | 
#include "UseTheForce/ForceFieldFactory.hpp" | 
| 64 | 
> | 
#include "brains/ForceField.hpp" | 
| 65 | 
  | 
#include "utils/simError.h" | 
| 66 | 
  | 
#include "utils/StringUtils.hpp" | 
| 67 | 
  | 
#include "math/SeqRandNumGen.hpp" | 
| 68 | 
< | 
#ifdef IS_MPI | 
| 69 | 
< | 
#include "io/mpiBASS.h" | 
| 70 | 
< | 
#include "math/ParallelRandNumGen.hpp" | 
| 71 | 
< | 
#endif | 
| 68 | 
> | 
#include "mdParser/MDLexer.hpp" | 
| 69 | 
> | 
#include "mdParser/MDParser.hpp" | 
| 70 | 
> | 
#include "mdParser/MDTreeParser.hpp" | 
| 71 | 
> | 
#include "mdParser/SimplePreprocessor.hpp" | 
| 72 | 
> | 
#include "antlr/ANTLRException.hpp" | 
| 73 | 
> | 
#include "antlr/TokenStreamRecognitionException.hpp" | 
| 74 | 
> | 
#include "antlr/TokenStreamIOException.hpp" | 
| 75 | 
> | 
#include "antlr/TokenStreamException.hpp" | 
| 76 | 
> | 
#include "antlr/RecognitionException.hpp" | 
| 77 | 
> | 
#include "antlr/CharStreamException.hpp" | 
| 78 | 
  | 
 | 
| 79 | 
< | 
namespace oopse { | 
| 79 | 
> | 
#include "antlr/MismatchedCharException.hpp" | 
| 80 | 
> | 
#include "antlr/MismatchedTokenException.hpp" | 
| 81 | 
> | 
#include "antlr/NoViableAltForCharException.hpp" | 
| 82 | 
> | 
#include "antlr/NoViableAltException.hpp" | 
| 83 | 
> | 
 | 
| 84 | 
> | 
#include "types/DirectionalAdapter.hpp" | 
| 85 | 
> | 
#include "types/MultipoleAdapter.hpp" | 
| 86 | 
> | 
#include "types/EAMAdapter.hpp" | 
| 87 | 
> | 
#include "types/SuttonChenAdapter.hpp" | 
| 88 | 
> | 
#include "types/PolarizableAdapter.hpp" | 
| 89 | 
> | 
#include "types/FixedChargeAdapter.hpp" | 
| 90 | 
> | 
#include "types/FluctuatingChargeAdapter.hpp" | 
| 91 | 
> | 
 | 
| 92 | 
> | 
 | 
| 93 | 
> | 
namespace OpenMD { | 
| 94 | 
  | 
   | 
| 95 | 
< | 
  void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps,  | 
| 96 | 
< | 
                             Globals* simParams){ | 
| 95 | 
> | 
  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){ | 
| 96 | 
> | 
    Globals* simParams = NULL; | 
| 97 | 
> | 
    try { | 
| 98 | 
> | 
 | 
| 99 | 
> | 
      // Create a preprocessor that preprocesses md file into an ostringstream | 
| 100 | 
> | 
      std::stringstream ppStream; | 
| 101 | 
> | 
#ifdef IS_MPI             | 
| 102 | 
> | 
      int streamSize; | 
| 103 | 
> | 
      const int masterNode = 0; | 
| 104 | 
> | 
 | 
| 105 | 
> | 
      if (worldRank == masterNode) { | 
| 106 | 
> | 
        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); | 
| 107 | 
> | 
#endif                  | 
| 108 | 
> | 
        SimplePreprocessor preprocessor; | 
| 109 | 
> | 
        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream); | 
| 110 | 
> | 
                 | 
| 111 | 
> | 
#ifdef IS_MPI             | 
| 112 | 
> | 
        //brocasting the stream size | 
| 113 | 
> | 
        streamSize = ppStream.str().size() +1; | 
| 114 | 
> | 
        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); | 
| 115 | 
> | 
        MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI::CHAR, masterNode); | 
| 116 | 
> | 
                            | 
| 117 | 
> | 
      } else { | 
| 118 | 
> | 
        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); | 
| 119 | 
> | 
 | 
| 120 | 
> | 
        //get stream size | 
| 121 | 
> | 
        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); | 
| 122 | 
> | 
 | 
| 123 | 
> | 
        char* buf = new char[streamSize]; | 
| 124 | 
> | 
        assert(buf); | 
| 125 | 
> | 
                 | 
| 126 | 
> | 
        //receive file content | 
| 127 | 
> | 
        MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode); | 
| 128 | 
> | 
                 | 
| 129 | 
> | 
        ppStream.str(buf); | 
| 130 | 
> | 
        delete [] buf; | 
| 131 | 
> | 
      } | 
| 132 | 
> | 
#endif             | 
| 133 | 
> | 
      // Create a scanner that reads from the input stream | 
| 134 | 
> | 
      MDLexer lexer(ppStream); | 
| 135 | 
> | 
      lexer.setFilename(filename); | 
| 136 | 
> | 
      lexer.initDeferredLineCount(); | 
| 137 | 
  | 
     | 
| 138 | 
< | 
#ifdef IS_MPI | 
| 138 | 
> | 
      // Create a parser that reads from the scanner | 
| 139 | 
> | 
      MDParser parser(lexer); | 
| 140 | 
> | 
      parser.setFilename(filename); | 
| 141 | 
> | 
 | 
| 142 | 
> | 
      // Create an observer that synchorizes file name change | 
| 143 | 
> | 
      FilenameObserver observer; | 
| 144 | 
> | 
      observer.setLexer(&lexer); | 
| 145 | 
> | 
      observer.setParser(&parser); | 
| 146 | 
> | 
      lexer.setObserver(&observer); | 
| 147 | 
  | 
     | 
| 148 | 
< | 
    if (worldRank == 0) { | 
| 149 | 
< | 
#endif // is_mpi | 
| 148 | 
> | 
      antlr::ASTFactory factory; | 
| 149 | 
> | 
      parser.initializeASTFactory(factory); | 
| 150 | 
> | 
      parser.setASTFactory(&factory); | 
| 151 | 
> | 
      parser.mdfile(); | 
| 152 | 
> | 
 | 
| 153 | 
> | 
      // Create a tree parser that reads information into Globals | 
| 154 | 
> | 
      MDTreeParser treeParser; | 
| 155 | 
> | 
      treeParser.initializeASTFactory(factory); | 
| 156 | 
> | 
      treeParser.setASTFactory(&factory); | 
| 157 | 
> | 
      simParams = treeParser.walkTree(parser.getAST()); | 
| 158 | 
> | 
    } | 
| 159 | 
> | 
 | 
| 160 | 
  | 
       | 
| 161 | 
< | 
      simParams->initalize(); | 
| 162 | 
< | 
      set_interface_stamps(stamps, simParams); | 
| 161 | 
> | 
    catch(antlr::MismatchedCharException& e) { | 
| 162 | 
> | 
      sprintf(painCave.errMsg,  | 
| 163 | 
> | 
              "parser exception: %s %s:%d:%d\n", | 
| 164 | 
> | 
              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 165 | 
> | 
      painCave.isFatal = 1; | 
| 166 | 
> | 
      simError();            | 
| 167 | 
> | 
    } | 
| 168 | 
> | 
    catch(antlr::MismatchedTokenException &e) { | 
| 169 | 
> | 
      sprintf(painCave.errMsg,  | 
| 170 | 
> | 
              "parser exception: %s %s:%d:%d\n", | 
| 171 | 
> | 
              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 172 | 
> | 
      painCave.isFatal = 1; | 
| 173 | 
> | 
      simError();    | 
| 174 | 
> | 
    } | 
| 175 | 
> | 
    catch(antlr::NoViableAltForCharException &e) { | 
| 176 | 
> | 
      sprintf(painCave.errMsg,  | 
| 177 | 
> | 
              "parser exception: %s %s:%d:%d\n", | 
| 178 | 
> | 
              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 179 | 
> | 
      painCave.isFatal = 1; | 
| 180 | 
> | 
      simError();    | 
| 181 | 
> | 
    } | 
| 182 | 
> | 
    catch(antlr::NoViableAltException &e) { | 
| 183 | 
> | 
      sprintf(painCave.errMsg,  | 
| 184 | 
> | 
              "parser exception: %s %s:%d:%d\n", | 
| 185 | 
> | 
              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 186 | 
> | 
      painCave.isFatal = 1; | 
| 187 | 
> | 
      simError();    | 
| 188 | 
> | 
    } | 
| 189 | 
  | 
       | 
| 190 | 
< | 
#ifdef IS_MPI | 
| 191 | 
< | 
       | 
| 192 | 
< | 
      mpiEventInit(); | 
| 193 | 
< | 
       | 
| 194 | 
< | 
#endif | 
| 195 | 
< | 
       | 
| 83 | 
< | 
      yacc_BASS(mdFileName.c_str()); | 
| 84 | 
< | 
       | 
| 85 | 
< | 
#ifdef IS_MPI | 
| 86 | 
< | 
       | 
| 87 | 
< | 
      throwMPIEvent(NULL); | 
| 88 | 
< | 
    } else { | 
| 89 | 
< | 
      set_interface_stamps(stamps, simParams); | 
| 90 | 
< | 
      mpiEventInit(); | 
| 91 | 
< | 
      MPIcheckPoint(); | 
| 92 | 
< | 
      mpiEventLoop(); | 
| 190 | 
> | 
    catch(antlr::TokenStreamRecognitionException& e) { | 
| 191 | 
> | 
      sprintf(painCave.errMsg,  | 
| 192 | 
> | 
              "parser exception: %s %s:%d:%d\n", | 
| 193 | 
> | 
              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 194 | 
> | 
      painCave.isFatal = 1; | 
| 195 | 
> | 
      simError();    | 
| 196 | 
  | 
    } | 
| 197 | 
< | 
     | 
| 198 | 
< | 
#endif | 
| 199 | 
< | 
     | 
| 197 | 
> | 
         | 
| 198 | 
> | 
    catch(antlr::TokenStreamIOException& e) { | 
| 199 | 
> | 
      sprintf(painCave.errMsg,  | 
| 200 | 
> | 
              "parser exception: %s\n", | 
| 201 | 
> | 
              e.getMessage().c_str()); | 
| 202 | 
> | 
      painCave.isFatal = 1; | 
| 203 | 
> | 
      simError(); | 
| 204 | 
> | 
    } | 
| 205 | 
> | 
         | 
| 206 | 
> | 
    catch(antlr::TokenStreamException& e) { | 
| 207 | 
> | 
      sprintf(painCave.errMsg,  | 
| 208 | 
> | 
              "parser exception: %s\n", | 
| 209 | 
> | 
              e.getMessage().c_str()); | 
| 210 | 
> | 
      painCave.isFatal = 1; | 
| 211 | 
> | 
      simError(); | 
| 212 | 
> | 
    }         | 
| 213 | 
> | 
    catch (antlr::RecognitionException& e) { | 
| 214 | 
> | 
      sprintf(painCave.errMsg,  | 
| 215 | 
> | 
              "parser exception: %s %s:%d:%d\n", | 
| 216 | 
> | 
              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 217 | 
> | 
      painCave.isFatal = 1; | 
| 218 | 
> | 
      simError();           | 
| 219 | 
> | 
    } | 
| 220 | 
> | 
    catch (antlr::CharStreamException& e) { | 
| 221 | 
> | 
      sprintf(painCave.errMsg,  | 
| 222 | 
> | 
              "parser exception: %s\n", | 
| 223 | 
> | 
              e.getMessage().c_str()); | 
| 224 | 
> | 
      painCave.isFatal = 1; | 
| 225 | 
> | 
      simError();         | 
| 226 | 
> | 
    } | 
| 227 | 
> | 
    catch (OpenMDException& e) { | 
| 228 | 
> | 
      sprintf(painCave.errMsg,  | 
| 229 | 
> | 
              "%s\n", | 
| 230 | 
> | 
              e.getMessage().c_str()); | 
| 231 | 
> | 
      painCave.isFatal = 1; | 
| 232 | 
> | 
      simError(); | 
| 233 | 
> | 
    } | 
| 234 | 
> | 
    catch (std::exception& e) { | 
| 235 | 
> | 
      sprintf(painCave.errMsg,  | 
| 236 | 
> | 
              "parser exception: %s\n", | 
| 237 | 
> | 
              e.what()); | 
| 238 | 
> | 
      painCave.isFatal = 1; | 
| 239 | 
> | 
      simError(); | 
| 240 | 
> | 
    } | 
| 241 | 
> | 
 | 
| 242 | 
> | 
    simParams->setMDfileVersion(mdFileVersion); | 
| 243 | 
> | 
    return simParams; | 
| 244 | 
  | 
  } | 
| 245 | 
  | 
   | 
| 246 | 
< | 
  SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) { | 
| 246 | 
> | 
  SimInfo*  SimCreator::createSim(const std::string & mdFileName,  | 
| 247 | 
> | 
                                  bool loadInitCoords) { | 
| 248 | 
  | 
     | 
| 249 | 
< | 
    MakeStamps * stamps = new MakeStamps(); | 
| 250 | 
< | 
     | 
| 251 | 
< | 
    Globals * simParams = new Globals(); | 
| 252 | 
< | 
     | 
| 249 | 
> | 
    const int bufferSize = 65535; | 
| 250 | 
> | 
    char buffer[bufferSize]; | 
| 251 | 
> | 
    int lineNo = 0; | 
| 252 | 
> | 
    std::string mdRawData; | 
| 253 | 
> | 
    int metaDataBlockStart = -1; | 
| 254 | 
> | 
    int metaDataBlockEnd = -1; | 
| 255 | 
> | 
    int i, j; | 
| 256 | 
> | 
    streamoff mdOffset; | 
| 257 | 
> | 
    int mdFileVersion; | 
| 258 | 
> | 
 | 
| 259 | 
> | 
    // Create a string for embedding the version information in the MetaData | 
| 260 | 
> | 
    std::string version; | 
| 261 | 
> | 
    version.assign("## Last run using OpenMD Version: "); | 
| 262 | 
> | 
    version.append(OPENMD_VERSION_MAJOR); | 
| 263 | 
> | 
    version.append("."); | 
| 264 | 
> | 
    version.append(OPENMD_VERSION_MINOR); | 
| 265 | 
> | 
 | 
| 266 | 
> | 
    std::string svnrev; | 
| 267 | 
> | 
    //convert a macro from compiler to a string in c++ | 
| 268 | 
> | 
    STR_DEFINE(svnrev, SVN_REV ); | 
| 269 | 
> | 
    version.append(" Revision: "); | 
| 270 | 
> | 
    // If there's no SVN revision, just call this the RELEASE revision. | 
| 271 | 
> | 
    if (!svnrev.empty()) { | 
| 272 | 
> | 
      version.append(svnrev); | 
| 273 | 
> | 
    } else { | 
| 274 | 
> | 
      version.append("RELEASE"); | 
| 275 | 
> | 
    } | 
| 276 | 
> | 
    | 
| 277 | 
> | 
#ifdef IS_MPI             | 
| 278 | 
> | 
    const int masterNode = 0; | 
| 279 | 
> | 
    if (worldRank == masterNode) { | 
| 280 | 
> | 
#endif  | 
| 281 | 
> | 
 | 
| 282 | 
> | 
      std::ifstream mdFile_; | 
| 283 | 
> | 
      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary); | 
| 284 | 
> | 
       | 
| 285 | 
> | 
      if (mdFile_.fail()) {  | 
| 286 | 
> | 
        sprintf(painCave.errMsg,  | 
| 287 | 
> | 
                "SimCreator: Cannot open file: %s\n",  | 
| 288 | 
> | 
                mdFileName.c_str());  | 
| 289 | 
> | 
        painCave.isFatal = 1;  | 
| 290 | 
> | 
        simError();  | 
| 291 | 
> | 
      }  | 
| 292 | 
> | 
 | 
| 293 | 
> | 
      mdFile_.getline(buffer, bufferSize); | 
| 294 | 
> | 
      ++lineNo; | 
| 295 | 
> | 
      std::string line = trimLeftCopy(buffer); | 
| 296 | 
> | 
      i = CaseInsensitiveFind(line, "<OpenMD"); | 
| 297 | 
> | 
      if (static_cast<size_t>(i) == string::npos) { | 
| 298 | 
> | 
        // try the older file strings to see if that works: | 
| 299 | 
> | 
        i = CaseInsensitiveFind(line, "<OOPSE"); | 
| 300 | 
> | 
      } | 
| 301 | 
> | 
       | 
| 302 | 
> | 
      if (static_cast<size_t>(i) == string::npos) { | 
| 303 | 
> | 
        // still no luck! | 
| 304 | 
> | 
        sprintf(painCave.errMsg,  | 
| 305 | 
> | 
                "SimCreator: File: %s is not a valid OpenMD file!\n",  | 
| 306 | 
> | 
                mdFileName.c_str());  | 
| 307 | 
> | 
        painCave.isFatal = 1;  | 
| 308 | 
> | 
        simError();  | 
| 309 | 
> | 
      } | 
| 310 | 
> | 
       | 
| 311 | 
> | 
      // found the correct opening string, now try to get the file | 
| 312 | 
> | 
      // format version number. | 
| 313 | 
> | 
 | 
| 314 | 
> | 
      StringTokenizer tokenizer(line, "=<> \t\n\r"); | 
| 315 | 
> | 
      std::string fileType = tokenizer.nextToken(); | 
| 316 | 
> | 
      toUpper(fileType); | 
| 317 | 
> | 
 | 
| 318 | 
> | 
      mdFileVersion = 0; | 
| 319 | 
> | 
 | 
| 320 | 
> | 
      if (fileType == "OPENMD") { | 
| 321 | 
> | 
        while (tokenizer.hasMoreTokens()) { | 
| 322 | 
> | 
          std::string token(tokenizer.nextToken()); | 
| 323 | 
> | 
          toUpper(token); | 
| 324 | 
> | 
          if (token == "VERSION") { | 
| 325 | 
> | 
            mdFileVersion = tokenizer.nextTokenAsInt(); | 
| 326 | 
> | 
            break; | 
| 327 | 
> | 
          } | 
| 328 | 
> | 
        } | 
| 329 | 
> | 
      } | 
| 330 | 
> | 
             | 
| 331 | 
> | 
      //scan through the input stream and find MetaData tag         | 
| 332 | 
> | 
      while(mdFile_.getline(buffer, bufferSize)) { | 
| 333 | 
> | 
        ++lineNo; | 
| 334 | 
> | 
         | 
| 335 | 
> | 
        std::string line = trimLeftCopy(buffer); | 
| 336 | 
> | 
        if (metaDataBlockStart == -1) { | 
| 337 | 
> | 
          i = CaseInsensitiveFind(line, "<MetaData>"); | 
| 338 | 
> | 
          if (i != string::npos) { | 
| 339 | 
> | 
            metaDataBlockStart = lineNo; | 
| 340 | 
> | 
            mdOffset = mdFile_.tellg(); | 
| 341 | 
> | 
          } | 
| 342 | 
> | 
        } else { | 
| 343 | 
> | 
          i = CaseInsensitiveFind(line, "</MetaData>"); | 
| 344 | 
> | 
          if (i != string::npos) { | 
| 345 | 
> | 
            metaDataBlockEnd = lineNo; | 
| 346 | 
> | 
          } | 
| 347 | 
> | 
        } | 
| 348 | 
> | 
      } | 
| 349 | 
> | 
 | 
| 350 | 
> | 
      if (metaDataBlockStart == -1) { | 
| 351 | 
> | 
        sprintf(painCave.errMsg,  | 
| 352 | 
> | 
                "SimCreator: File: %s did not contain a <MetaData> tag!\n",  | 
| 353 | 
> | 
                mdFileName.c_str());  | 
| 354 | 
> | 
        painCave.isFatal = 1;  | 
| 355 | 
> | 
        simError();  | 
| 356 | 
> | 
      } | 
| 357 | 
> | 
      if (metaDataBlockEnd == -1) { | 
| 358 | 
> | 
        sprintf(painCave.errMsg,  | 
| 359 | 
> | 
                "SimCreator: File: %s did not contain a closed MetaData block!\n",  | 
| 360 | 
> | 
                mdFileName.c_str());  | 
| 361 | 
> | 
        painCave.isFatal = 1;  | 
| 362 | 
> | 
        simError();  | 
| 363 | 
> | 
      } | 
| 364 | 
> | 
         | 
| 365 | 
> | 
      mdFile_.clear(); | 
| 366 | 
> | 
      mdFile_.seekg(0); | 
| 367 | 
> | 
      mdFile_.seekg(mdOffset); | 
| 368 | 
> | 
 | 
| 369 | 
> | 
      mdRawData.clear(); | 
| 370 | 
> | 
 | 
| 371 | 
> | 
      bool foundVersion = false; | 
| 372 | 
> | 
 | 
| 373 | 
> | 
      for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) { | 
| 374 | 
> | 
        mdFile_.getline(buffer, bufferSize); | 
| 375 | 
> | 
        std::string line = trimLeftCopy(buffer); | 
| 376 | 
> | 
        j = CaseInsensitiveFind(line, "## Last run using OpenMD Version"); | 
| 377 | 
> | 
        if (static_cast<size_t>(j) != string::npos) { | 
| 378 | 
> | 
          foundVersion = true; | 
| 379 | 
> | 
          mdRawData += version; | 
| 380 | 
> | 
        } else { | 
| 381 | 
> | 
          mdRawData += buffer; | 
| 382 | 
> | 
        } | 
| 383 | 
> | 
        mdRawData += "\n"; | 
| 384 | 
> | 
      } | 
| 385 | 
> | 
       | 
| 386 | 
> | 
      if (!foundVersion) mdRawData += version + "\n"; | 
| 387 | 
> | 
       | 
| 388 | 
> | 
      mdFile_.close(); | 
| 389 | 
> | 
 | 
| 390 | 
> | 
#ifdef IS_MPI | 
| 391 | 
> | 
    } | 
| 392 | 
> | 
#endif | 
| 393 | 
> | 
 | 
| 394 | 
> | 
    std::stringstream rawMetaDataStream(mdRawData); | 
| 395 | 
> | 
 | 
| 396 | 
  | 
    //parse meta-data file | 
| 397 | 
< | 
    parseFile(mdFileName, stamps, simParams); | 
| 397 | 
> | 
    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion, | 
| 398 | 
> | 
                                   metaDataBlockStart + 1); | 
| 399 | 
  | 
     | 
| 400 | 
  | 
    //create the force field | 
| 401 | 
< | 
    ForceField * ff = ForceFieldFactory::getInstance()->createForceField( | 
| 402 | 
< | 
                                                                         simParams->getForceField()); | 
| 111 | 
< | 
     | 
| 401 | 
> | 
    ForceField * ff = new ForceField(simParams->getForceField()); | 
| 402 | 
> | 
 | 
| 403 | 
  | 
    if (ff == NULL) { | 
| 404 | 
< | 
      sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n", | 
| 405 | 
< | 
              simParams->getForceField()); | 
| 404 | 
> | 
      sprintf(painCave.errMsg,  | 
| 405 | 
> | 
              "ForceField Factory can not create %s force field\n", | 
| 406 | 
> | 
              simParams->getForceField().c_str()); | 
| 407 | 
  | 
      painCave.isFatal = 1; | 
| 408 | 
  | 
      simError(); | 
| 409 | 
  | 
    } | 
| 432 | 
  | 
    }  | 
| 433 | 
  | 
     | 
| 434 | 
  | 
    ff->parse(forcefieldFileName); | 
| 143 | 
– | 
     | 
| 144 | 
– | 
    //extract the molecule stamps | 
| 145 | 
– | 
    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs; | 
| 146 | 
– | 
    compList(stamps, simParams, moleculeStampPairs); | 
| 147 | 
– | 
     | 
| 435 | 
  | 
    //create SimInfo | 
| 436 | 
< | 
    SimInfo * info = new SimInfo(stamps, moleculeStampPairs, ff, simParams); | 
| 436 | 
> | 
    SimInfo * info = new SimInfo(ff, simParams); | 
| 437 | 
> | 
 | 
| 438 | 
> | 
    info->setRawMetaData(mdRawData); | 
| 439 | 
  | 
      | 
| 440 | 
< | 
    //gather parameters (SimCreator only retrieves part of the parameters) | 
| 440 | 
> | 
    //gather parameters (SimCreator only retrieves part of the | 
| 441 | 
> | 
    //parameters) | 
| 442 | 
  | 
    gatherParameters(info, mdFileName); | 
| 443 | 
  | 
     | 
| 444 | 
  | 
    //divide the molecules and determine the global index of molecules | 
| 449 | 
  | 
    //create the molecules | 
| 450 | 
  | 
    createMolecules(info); | 
| 451 | 
  | 
     | 
| 452 | 
+ | 
    //find the storage layout | 
| 453 | 
+ | 
 | 
| 454 | 
+ | 
    int storageLayout = computeStorageLayout(info); | 
| 455 | 
+ | 
 | 
| 456 | 
+ | 
    //allocate memory for DataStorage(circular reference, need to | 
| 457 | 
+ | 
    //break it) | 
| 458 | 
+ | 
    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout)); | 
| 459 | 
  | 
     | 
| 460 | 
< | 
    //allocate memory for DataStorage(circular reference, need to break it) | 
| 461 | 
< | 
    info->setSnapshotManager(new SimSnapshotManager(info)); | 
| 462 | 
< | 
     | 
| 463 | 
< | 
    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the | 
| 464 | 
< | 
    //global index will never change again). Local indices of atoms and rigidbodies are already set by  | 
| 168 | 
< | 
    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.  | 
| 460 | 
> | 
    //set the global index of atoms, rigidbodies and cutoffgroups | 
| 461 | 
> | 
    //(only need to be set once, the global index will never change | 
| 462 | 
> | 
    //again). Local indices of atoms and rigidbodies are already set | 
| 463 | 
> | 
    //by MoleculeCreator class which actually delegates the | 
| 464 | 
> | 
    //responsibility to LocalIndexManager. | 
| 465 | 
  | 
    setGlobalIndex(info); | 
| 466 | 
  | 
     | 
| 467 | 
< | 
    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point | 
| 468 | 
< | 
    //atoms don't have the global index yet  (their global index are all initialized to -1). | 
| 469 | 
< | 
    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that | 
| 470 | 
< | 
    //we can determine the beginning global indices of atoms before they get created. | 
| 467 | 
> | 
    //Although addInteractionPairs is called inside SimInfo's addMolecule | 
| 468 | 
> | 
    //method, at that point atoms don't have the global index yet | 
| 469 | 
> | 
    //(their global index are all initialized to -1).  Therefore we | 
| 470 | 
> | 
    //have to call addInteractionPairs explicitly here. A way to work | 
| 471 | 
> | 
    //around is that we can determine the beginning global indices of | 
| 472 | 
> | 
    //atoms before they get created. | 
| 473 | 
  | 
    SimInfo::MoleculeIterator mi; | 
| 474 | 
  | 
    Molecule* mol; | 
| 475 | 
  | 
    for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 476 | 
< | 
      info->addExcludePairs(mol); | 
| 476 | 
> | 
      info->addInteractionPairs(mol); | 
| 477 | 
  | 
    } | 
| 478 | 
  | 
     | 
| 479 | 
  | 
    if (loadInitCoords) | 
| 480 | 
< | 
      loadCoordinates(info);     | 
| 183 | 
< | 
     | 
| 480 | 
> | 
      loadCoordinates(info, mdFileName);     | 
| 481 | 
  | 
    return info; | 
| 482 | 
  | 
  } | 
| 483 | 
  | 
   | 
| 484 | 
  | 
  void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) { | 
| 485 | 
  | 
     | 
| 486 | 
< | 
    //figure out the ouput file names | 
| 486 | 
> | 
    //figure out the output file names | 
| 487 | 
  | 
    std::string prefix; | 
| 488 | 
  | 
     | 
| 489 | 
  | 
#ifdef IS_MPI | 
| 512 | 
  | 
   | 
| 513 | 
  | 
#ifdef IS_MPI | 
| 514 | 
  | 
  void SimCreator::divideMolecules(SimInfo *info) { | 
| 515 | 
< | 
    double numerator; | 
| 219 | 
< | 
    double denominator; | 
| 220 | 
< | 
    double precast; | 
| 221 | 
< | 
    double x; | 
| 222 | 
< | 
    double y; | 
| 223 | 
< | 
    double a; | 
| 224 | 
< | 
    int old_atoms; | 
| 225 | 
< | 
    int add_atoms; | 
| 226 | 
< | 
    int new_atoms; | 
| 227 | 
< | 
    int nTarget; | 
| 228 | 
< | 
    int done; | 
| 229 | 
< | 
    int i; | 
| 230 | 
< | 
    int j; | 
| 231 | 
< | 
    int loops; | 
| 232 | 
< | 
    int which_proc; | 
| 515 | 
> | 
    RealType a; | 
| 516 | 
  | 
    int nProcessors; | 
| 517 | 
  | 
    std::vector<int> atomsPerProc; | 
| 518 | 
  | 
    int nGlobalMols = info->getNGlobalMolecules(); | 
| 519 | 
< | 
    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition: | 
| 519 | 
> | 
    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an | 
| 520 | 
> | 
                                                    // error | 
| 521 | 
> | 
                                                    // condition: | 
| 522 | 
  | 
     | 
| 523 | 
< | 
    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors); | 
| 523 | 
> | 
    nProcessors = MPI::COMM_WORLD.Get_size(); | 
| 524 | 
  | 
     | 
| 525 | 
  | 
    if (nProcessors > nGlobalMols) { | 
| 526 | 
  | 
      sprintf(painCave.errMsg, | 
| 529 | 
  | 
              "\tthe number of molecules.  This will not result in a \n" | 
| 530 | 
  | 
              "\tusable division of atoms for force decomposition.\n" | 
| 531 | 
  | 
              "\tEither try a smaller number of processors, or run the\n" | 
| 532 | 
< | 
              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols); | 
| 532 | 
> | 
              "\tsingle-processor version of OpenMD.\n", nProcessors,  | 
| 533 | 
> | 
              nGlobalMols); | 
| 534 | 
  | 
       | 
| 535 | 
  | 
      painCave.isFatal = 1; | 
| 536 | 
  | 
      simError(); | 
| 537 | 
  | 
    } | 
| 538 | 
  | 
     | 
| 253 | 
– | 
    int seedValue; | 
| 539 | 
  | 
    Globals * simParams = info->getSimParams(); | 
| 540 | 
< | 
    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator | 
| 540 | 
> | 
    SeqRandNumGen* myRandom; //divide labor does not need Parallel | 
| 541 | 
> | 
                             //random number generator | 
| 542 | 
  | 
    if (simParams->haveSeed()) { | 
| 543 | 
< | 
      seedValue = simParams->getSeed(); | 
| 543 | 
> | 
      int seedValue = simParams->getSeed(); | 
| 544 | 
  | 
      myRandom = new SeqRandNumGen(seedValue); | 
| 545 | 
  | 
    }else { | 
| 546 | 
  | 
      myRandom = new SeqRandNumGen(); | 
| 553 | 
  | 
    atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0); | 
| 554 | 
  | 
     | 
| 555 | 
  | 
    if (worldRank == 0) { | 
| 556 | 
< | 
      numerator = info->getNGlobalAtoms(); | 
| 557 | 
< | 
      denominator = nProcessors; | 
| 558 | 
< | 
      precast = numerator / denominator; | 
| 559 | 
< | 
      nTarget = (int)(precast + 0.5); | 
| 556 | 
> | 
      RealType numerator = info->getNGlobalAtoms(); | 
| 557 | 
> | 
      RealType denominator = nProcessors; | 
| 558 | 
> | 
      RealType precast = numerator / denominator; | 
| 559 | 
> | 
      int nTarget = (int)(precast + 0.5); | 
| 560 | 
  | 
       | 
| 561 | 
< | 
      for(i = 0; i < nGlobalMols; i++) { | 
| 562 | 
< | 
        done = 0; | 
| 563 | 
< | 
        loops = 0; | 
| 561 | 
> | 
      for(int i = 0; i < nGlobalMols; i++) { | 
| 562 | 
> | 
 | 
| 563 | 
> | 
        int done = 0; | 
| 564 | 
> | 
        int loops = 0; | 
| 565 | 
  | 
         | 
| 566 | 
  | 
        while (!done) { | 
| 567 | 
  | 
          loops++; | 
| 568 | 
  | 
           | 
| 569 | 
  | 
          // Pick a processor at random | 
| 570 | 
  | 
           | 
| 571 | 
< | 
          which_proc = (int) (myRandom->rand() * nProcessors); | 
| 571 | 
> | 
          int which_proc = (int) (myRandom->rand() * nProcessors); | 
| 572 | 
  | 
           | 
| 573 | 
  | 
          //get the molecule stamp first | 
| 574 | 
  | 
          int stampId = info->getMoleculeStampId(i); | 
| 575 | 
  | 
          MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId); | 
| 576 | 
  | 
           | 
| 577 | 
  | 
          // How many atoms does this processor have so far? | 
| 578 | 
< | 
          old_atoms = atomsPerProc[which_proc]; | 
| 579 | 
< | 
          add_atoms = moleculeStamp->getNAtoms(); | 
| 580 | 
< | 
          new_atoms = old_atoms + add_atoms; | 
| 578 | 
> | 
          int old_atoms = atomsPerProc[which_proc]; | 
| 579 | 
> | 
          int add_atoms = moleculeStamp->getNAtoms(); | 
| 580 | 
> | 
          int new_atoms = old_atoms + add_atoms; | 
| 581 | 
  | 
           | 
| 582 | 
  | 
          // If we've been through this loop too many times, we need | 
| 583 | 
  | 
          // to just give up and assign the molecule to this processor | 
| 584 | 
  | 
          // and be done with it.  | 
| 585 | 
  | 
           | 
| 586 | 
  | 
          if (loops > 100) { | 
| 587 | 
+ | 
 | 
| 588 | 
  | 
            sprintf(painCave.errMsg, | 
| 589 | 
< | 
                    "I've tried 100 times to assign molecule %d to a " | 
| 590 | 
< | 
                    " processor, but can't find a good spot.\n" | 
| 591 | 
< | 
                    "I'm assigning it at random to processor %d.\n", | 
| 589 | 
> | 
                    "There have been 100 attempts to assign molecule %d to an\n" | 
| 590 | 
> | 
                    "\tunderworked processor, but there's no good place to\n" | 
| 591 | 
> | 
                    "\tleave it.  OpenMD is assigning it at random to processor %d.\n", | 
| 592 | 
  | 
                    i, which_proc); | 
| 593 | 
< | 
             | 
| 593 | 
> | 
            | 
| 594 | 
  | 
            painCave.isFatal = 0; | 
| 595 | 
+ | 
            painCave.severity = OPENMD_INFO; | 
| 596 | 
  | 
            simError(); | 
| 597 | 
  | 
             | 
| 598 | 
  | 
            molToProcMap[i] = which_proc; | 
| 621 | 
  | 
          //           Pacc(x) = exp(- a * x) | 
| 622 | 
  | 
          // where a = penalty / (average atoms per molecule) | 
| 623 | 
  | 
           | 
| 624 | 
< | 
          x = (double)(new_atoms - nTarget); | 
| 625 | 
< | 
          y = myRandom->rand(); | 
| 624 | 
> | 
          RealType x = (RealType)(new_atoms - nTarget); | 
| 625 | 
> | 
          RealType y = myRandom->rand(); | 
| 626 | 
  | 
           | 
| 627 | 
  | 
          if (y < exp(- a * x)) { | 
| 628 | 
  | 
            molToProcMap[i] = which_proc; | 
| 637 | 
  | 
      } | 
| 638 | 
  | 
       | 
| 639 | 
  | 
      delete myRandom; | 
| 640 | 
< | 
       | 
| 640 | 
> | 
 | 
| 641 | 
  | 
      // Spray out this nonsense to all other processors: | 
| 642 | 
< | 
       | 
| 354 | 
< | 
      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); | 
| 642 | 
> | 
      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0); | 
| 643 | 
  | 
    } else { | 
| 644 | 
  | 
       | 
| 645 | 
  | 
      // Listen to your marching orders from processor 0: | 
| 646 | 
< | 
       | 
| 647 | 
< | 
      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); | 
| 646 | 
> | 
      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0); | 
| 647 | 
> | 
 | 
| 648 | 
  | 
    } | 
| 649 | 
  | 
     | 
| 650 | 
  | 
    info->setMolToProcMap(molToProcMap); | 
| 651 | 
  | 
    sprintf(checkPointMsg, | 
| 652 | 
  | 
            "Successfully divided the molecules among the processors.\n"); | 
| 653 | 
< | 
    MPIcheckPoint(); | 
| 653 | 
> | 
    errorCheckPoint(); | 
| 654 | 
  | 
  } | 
| 655 | 
  | 
   | 
| 656 | 
  | 
#endif | 
| 667 | 
  | 
#endif | 
| 668 | 
  | 
         | 
| 669 | 
  | 
        stampId = info->getMoleculeStampId(i); | 
| 670 | 
< | 
        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId), | 
| 671 | 
< | 
                                                   stampId, i, info->getLocalIndexManager()); | 
| 670 | 
> | 
        Molecule * mol = molCreator.createMolecule(info->getForceField(),  | 
| 671 | 
> | 
                                                   info->getMoleculeStamp(stampId), | 
| 672 | 
> | 
                                                   stampId, i,  | 
| 673 | 
> | 
                                                   info->getLocalIndexManager()); | 
| 674 | 
  | 
         | 
| 675 | 
  | 
        info->addMolecule(mol); | 
| 676 | 
  | 
         | 
| 682 | 
  | 
       | 
| 683 | 
  | 
    } //end for(int i=0)    | 
| 684 | 
  | 
  } | 
| 395 | 
– | 
   | 
| 396 | 
– | 
  void SimCreator::compList(MakeStamps *stamps, Globals* simParams, | 
| 397 | 
– | 
                            std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) { | 
| 398 | 
– | 
    int i; | 
| 399 | 
– | 
    char * id; | 
| 400 | 
– | 
    MoleculeStamp * currentStamp; | 
| 401 | 
– | 
    Component** the_components = simParams->getComponents(); | 
| 402 | 
– | 
    int n_components = simParams->getNComponents(); | 
| 685 | 
  | 
     | 
| 686 | 
< | 
    if (!simParams->haveNMol()) { | 
| 405 | 
< | 
      // we don't have the total number of molecules, so we assume it is | 
| 406 | 
< | 
      // given in each component | 
| 407 | 
< | 
       | 
| 408 | 
< | 
      for(i = 0; i < n_components; i++) { | 
| 409 | 
< | 
        if (!the_components[i]->haveNMol()) { | 
| 410 | 
< | 
          // we have a problem | 
| 411 | 
< | 
          sprintf(painCave.errMsg, | 
| 412 | 
< | 
                  "SimCreator Error. No global NMol or component NMol given.\n" | 
| 413 | 
< | 
                  "\tCannot calculate the number of atoms.\n"); | 
| 414 | 
< | 
           | 
| 415 | 
< | 
          painCave.isFatal = 1; | 
| 416 | 
< | 
          simError(); | 
| 417 | 
< | 
        } | 
| 418 | 
< | 
         | 
| 419 | 
< | 
        id = the_components[i]->getType(); | 
| 686 | 
> | 
  int SimCreator::computeStorageLayout(SimInfo* info) { | 
| 687 | 
  | 
 | 
| 688 | 
< | 
        currentStamp = stamps->getMolStamp(id); | 
| 689 | 
< | 
        if (currentStamp == NULL) { | 
| 690 | 
< | 
          sprintf(painCave.errMsg, | 
| 691 | 
< | 
                  "SimCreator error: Component \"%s\" was not found in the " | 
| 692 | 
< | 
                  "list of declared molecules\n", id); | 
| 693 | 
< | 
           | 
| 694 | 
< | 
          painCave.isFatal = 1; | 
| 695 | 
< | 
          simError(); | 
| 696 | 
< | 
        } | 
| 697 | 
< | 
         | 
| 698 | 
< | 
        moleculeStampPairs.push_back( | 
| 699 | 
< | 
                                     std::make_pair(currentStamp, the_components[i]->getNMol())); | 
| 700 | 
< | 
      } //end for (i = 0; i < n_components; i++) | 
| 701 | 
< | 
    } else { | 
| 702 | 
< | 
      sprintf(painCave.errMsg, "SimSetup error.\n" | 
| 703 | 
< | 
              "\tSorry, the ability to specify total" | 
| 704 | 
< | 
              " nMols and then give molfractions in the components\n" | 
| 705 | 
< | 
              "\tis not currently supported." | 
| 706 | 
< | 
              " Please give nMol in the components.\n"); | 
| 707 | 
< | 
       | 
| 708 | 
< | 
      painCave.isFatal = 1; | 
| 709 | 
< | 
      simError(); | 
| 688 | 
> | 
    Globals* simParams = info->getSimParams(); | 
| 689 | 
> | 
    int nRigidBodies = info->getNGlobalRigidBodies(); | 
| 690 | 
> | 
    set<AtomType*> atomTypes = info->getSimulatedAtomTypes(); | 
| 691 | 
> | 
    set<AtomType*>::iterator i; | 
| 692 | 
> | 
    bool hasDirectionalAtoms = false; | 
| 693 | 
> | 
    bool hasFixedCharge = false; | 
| 694 | 
> | 
    bool hasDipoles = false;     | 
| 695 | 
> | 
    bool hasQuadrupoles = false;     | 
| 696 | 
> | 
    bool hasPolarizable = false;     | 
| 697 | 
> | 
    bool hasFluctuatingCharge = false;     | 
| 698 | 
> | 
    bool hasMetallic = false; | 
| 699 | 
> | 
    int storageLayout = 0; | 
| 700 | 
> | 
    storageLayout |= DataStorage::dslPosition; | 
| 701 | 
> | 
    storageLayout |= DataStorage::dslVelocity; | 
| 702 | 
> | 
    storageLayout |= DataStorage::dslForce; | 
| 703 | 
> | 
 | 
| 704 | 
> | 
    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 705 | 
> | 
 | 
| 706 | 
> | 
      DirectionalAdapter da = DirectionalAdapter( (*i) ); | 
| 707 | 
> | 
      MultipoleAdapter ma = MultipoleAdapter( (*i) ); | 
| 708 | 
> | 
      EAMAdapter ea = EAMAdapter( (*i) ); | 
| 709 | 
> | 
      SuttonChenAdapter sca = SuttonChenAdapter( (*i) ); | 
| 710 | 
> | 
      PolarizableAdapter pa = PolarizableAdapter( (*i) ); | 
| 711 | 
> | 
      FixedChargeAdapter fca = FixedChargeAdapter( (*i) ); | 
| 712 | 
> | 
      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) ); | 
| 713 | 
> | 
 | 
| 714 | 
> | 
      if (da.isDirectional()){ | 
| 715 | 
> | 
        hasDirectionalAtoms = true; | 
| 716 | 
> | 
      } | 
| 717 | 
> | 
      if (ma.isDipole()){ | 
| 718 | 
> | 
        hasDipoles = true; | 
| 719 | 
> | 
      } | 
| 720 | 
> | 
      if (ma.isQuadrupole()){ | 
| 721 | 
> | 
        hasQuadrupoles = true; | 
| 722 | 
> | 
      } | 
| 723 | 
> | 
      if (ea.isEAM() || sca.isSuttonChen()){ | 
| 724 | 
> | 
        hasMetallic = true; | 
| 725 | 
> | 
      } | 
| 726 | 
> | 
      if ( fca.isFixedCharge() ){ | 
| 727 | 
> | 
        hasFixedCharge = true; | 
| 728 | 
> | 
      } | 
| 729 | 
> | 
      if ( fqa.isFluctuatingCharge() ){ | 
| 730 | 
> | 
        hasFluctuatingCharge = true; | 
| 731 | 
> | 
      } | 
| 732 | 
> | 
      if ( pa.isPolarizable() ){ | 
| 733 | 
> | 
        hasPolarizable = true; | 
| 734 | 
> | 
      } | 
| 735 | 
  | 
    } | 
| 736 | 
  | 
     | 
| 737 | 
< | 
#ifdef IS_MPI | 
| 737 | 
> | 
    if (nRigidBodies > 0 || hasDirectionalAtoms) { | 
| 738 | 
> | 
      storageLayout |= DataStorage::dslAmat; | 
| 739 | 
> | 
      if(storageLayout & DataStorage::dslVelocity) { | 
| 740 | 
> | 
        storageLayout |= DataStorage::dslAngularMomentum; | 
| 741 | 
> | 
      } | 
| 742 | 
> | 
      if (storageLayout & DataStorage::dslForce) { | 
| 743 | 
> | 
        storageLayout |= DataStorage::dslTorque; | 
| 744 | 
> | 
      } | 
| 745 | 
> | 
    } | 
| 746 | 
> | 
    if (hasDipoles) { | 
| 747 | 
> | 
      storageLayout |= DataStorage::dslDipole; | 
| 748 | 
> | 
    } | 
| 749 | 
> | 
    if (hasQuadrupoles) { | 
| 750 | 
> | 
      storageLayout |= DataStorage::dslQuadrupole; | 
| 751 | 
> | 
    } | 
| 752 | 
> | 
    if (hasFixedCharge || hasFluctuatingCharge) { | 
| 753 | 
> | 
      storageLayout |= DataStorage::dslSkippedCharge; | 
| 754 | 
> | 
    } | 
| 755 | 
> | 
    if (hasMetallic) { | 
| 756 | 
> | 
      storageLayout |= DataStorage::dslDensity; | 
| 757 | 
> | 
      storageLayout |= DataStorage::dslFunctional; | 
| 758 | 
> | 
      storageLayout |= DataStorage::dslFunctionalDerivative; | 
| 759 | 
> | 
    } | 
| 760 | 
> | 
    if (hasPolarizable) { | 
| 761 | 
> | 
      storageLayout |= DataStorage::dslElectricField; | 
| 762 | 
> | 
    } | 
| 763 | 
> | 
    if (hasFluctuatingCharge){ | 
| 764 | 
> | 
      storageLayout |= DataStorage::dslFlucQPosition; | 
| 765 | 
> | 
      if(storageLayout & DataStorage::dslVelocity) { | 
| 766 | 
> | 
        storageLayout |= DataStorage::dslFlucQVelocity; | 
| 767 | 
> | 
      } | 
| 768 | 
> | 
      if (storageLayout & DataStorage::dslForce) { | 
| 769 | 
> | 
        storageLayout |= DataStorage::dslFlucQForce; | 
| 770 | 
> | 
      } | 
| 771 | 
> | 
    } | 
| 772 | 
  | 
     | 
| 773 | 
< | 
    strcpy(checkPointMsg, "Component stamps successfully extracted\n"); | 
| 774 | 
< | 
    MPIcheckPoint(); | 
| 775 | 
< | 
     | 
| 776 | 
< | 
#endif // is_mpi | 
| 777 | 
< | 
     | 
| 773 | 
> | 
    // if the user has asked for them, make sure we've got the memory for the | 
| 774 | 
> | 
    // objects defined. | 
| 775 | 
> | 
 | 
| 776 | 
> | 
    if (simParams->getOutputParticlePotential()) { | 
| 777 | 
> | 
      storageLayout |= DataStorage::dslParticlePot; | 
| 778 | 
> | 
    } | 
| 779 | 
> | 
 | 
| 780 | 
> | 
    if (simParams->havePrintHeatFlux()) { | 
| 781 | 
> | 
      if (simParams->getPrintHeatFlux()) { | 
| 782 | 
> | 
        storageLayout |= DataStorage::dslParticlePot; | 
| 783 | 
> | 
      } | 
| 784 | 
> | 
    } | 
| 785 | 
> | 
 | 
| 786 | 
> | 
    if (simParams->getOutputElectricField() | simParams->haveElectricField()) { | 
| 787 | 
> | 
      storageLayout |= DataStorage::dslElectricField; | 
| 788 | 
> | 
    } | 
| 789 | 
> | 
 | 
| 790 | 
> | 
    if (simParams->getOutputFluctuatingCharges()) { | 
| 791 | 
> | 
      storageLayout |= DataStorage::dslFlucQPosition; | 
| 792 | 
> | 
      storageLayout |= DataStorage::dslFlucQVelocity; | 
| 793 | 
> | 
      storageLayout |= DataStorage::dslFlucQForce; | 
| 794 | 
> | 
    } | 
| 795 | 
> | 
 | 
| 796 | 
> | 
    info->setStorageLayout(storageLayout); | 
| 797 | 
> | 
 | 
| 798 | 
> | 
    return storageLayout; | 
| 799 | 
  | 
  } | 
| 800 | 
< | 
   | 
| 800 | 
> | 
 | 
| 801 | 
  | 
  void SimCreator::setGlobalIndex(SimInfo *info) { | 
| 802 | 
  | 
    SimInfo::MoleculeIterator mi; | 
| 803 | 
  | 
    Molecule::AtomIterator ai; | 
| 804 | 
  | 
    Molecule::RigidBodyIterator ri; | 
| 805 | 
  | 
    Molecule::CutoffGroupIterator ci; | 
| 806 | 
< | 
    Molecule * mol; | 
| 807 | 
< | 
    Atom * atom; | 
| 808 | 
< | 
    RigidBody * rb; | 
| 809 | 
< | 
    CutoffGroup * cg; | 
| 806 | 
> | 
    Molecule::BondIterator boi; | 
| 807 | 
> | 
    Molecule::BendIterator bei; | 
| 808 | 
> | 
    Molecule::TorsionIterator ti; | 
| 809 | 
> | 
    Molecule::InversionIterator ii; | 
| 810 | 
> | 
    Molecule::IntegrableObjectIterator  ioi; | 
| 811 | 
> | 
    Molecule* mol; | 
| 812 | 
> | 
    Atom* atom; | 
| 813 | 
> | 
    RigidBody* rb; | 
| 814 | 
> | 
    CutoffGroup* cg; | 
| 815 | 
> | 
    Bond* bond; | 
| 816 | 
> | 
    Bend* bend; | 
| 817 | 
> | 
    Torsion* torsion; | 
| 818 | 
> | 
    Inversion* inversion; | 
| 819 | 
  | 
    int beginAtomIndex; | 
| 820 | 
  | 
    int beginRigidBodyIndex; | 
| 821 | 
  | 
    int beginCutoffGroupIndex; | 
| 822 | 
+ | 
    int beginBondIndex; | 
| 823 | 
+ | 
    int beginBendIndex; | 
| 824 | 
+ | 
    int beginTorsionIndex; | 
| 825 | 
+ | 
    int beginInversionIndex; | 
| 826 | 
  | 
    int nGlobalAtoms = info->getNGlobalAtoms(); | 
| 827 | 
< | 
     | 
| 468 | 
< | 
#ifndef IS_MPI | 
| 827 | 
> | 
    int nGlobalRigidBodies = info->getNGlobalRigidBodies(); | 
| 828 | 
  | 
     | 
| 829 | 
  | 
    beginAtomIndex = 0; | 
| 830 | 
< | 
    beginRigidBodyIndex = 0; | 
| 830 | 
> | 
    // The rigid body indices begin immediately after the atom indices: | 
| 831 | 
> | 
    beginRigidBodyIndex = info->getNGlobalAtoms(); | 
| 832 | 
  | 
    beginCutoffGroupIndex = 0; | 
| 833 | 
< | 
     | 
| 834 | 
< | 
#else | 
| 835 | 
< | 
     | 
| 836 | 
< | 
    int nproc; | 
| 837 | 
< | 
    int myNode; | 
| 838 | 
< | 
     | 
| 479 | 
< | 
    myNode = worldRank; | 
| 480 | 
< | 
    MPI_Comm_size(MPI_COMM_WORLD, &nproc); | 
| 481 | 
< | 
     | 
| 482 | 
< | 
    std::vector < int > tmpAtomsInProc(nproc, 0); | 
| 483 | 
< | 
    std::vector < int > tmpRigidBodiesInProc(nproc, 0); | 
| 484 | 
< | 
    std::vector < int > tmpCutoffGroupsInProc(nproc, 0); | 
| 485 | 
< | 
    std::vector < int > NumAtomsInProc(nproc, 0); | 
| 486 | 
< | 
    std::vector < int > NumRigidBodiesInProc(nproc, 0); | 
| 487 | 
< | 
    std::vector < int > NumCutoffGroupsInProc(nproc, 0); | 
| 488 | 
< | 
     | 
| 489 | 
< | 
    tmpAtomsInProc[myNode] = info->getNAtoms(); | 
| 490 | 
< | 
    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies(); | 
| 491 | 
< | 
    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups(); | 
| 492 | 
< | 
     | 
| 493 | 
< | 
    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups | 
| 494 | 
< | 
    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT, | 
| 495 | 
< | 
                  MPI_SUM, MPI_COMM_WORLD); | 
| 496 | 
< | 
    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc, | 
| 497 | 
< | 
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 498 | 
< | 
    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc, | 
| 499 | 
< | 
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 500 | 
< | 
     | 
| 501 | 
< | 
    beginAtomIndex = 0; | 
| 502 | 
< | 
    beginRigidBodyIndex = 0; | 
| 503 | 
< | 
    beginCutoffGroupIndex = 0; | 
| 504 | 
< | 
     | 
| 505 | 
< | 
    for(int i = 0; i < myNode; i++) { | 
| 506 | 
< | 
      beginAtomIndex += NumAtomsInProc[i]; | 
| 507 | 
< | 
      beginRigidBodyIndex += NumRigidBodiesInProc[i]; | 
| 508 | 
< | 
      beginCutoffGroupIndex += NumCutoffGroupsInProc[i]; | 
| 509 | 
< | 
    } | 
| 510 | 
< | 
     | 
| 511 | 
< | 
#endif | 
| 512 | 
< | 
     | 
| 513 | 
< | 
    //rigidbody's index begins right after atom's | 
| 514 | 
< | 
    beginRigidBodyIndex += info->getNGlobalAtoms(); | 
| 515 | 
< | 
     | 
| 516 | 
< | 
    for(mol = info->beginMolecule(mi); mol != NULL; | 
| 517 | 
< | 
        mol = info->nextMolecule(mi)) { | 
| 833 | 
> | 
    beginBondIndex = 0; | 
| 834 | 
> | 
    beginBendIndex = 0; | 
| 835 | 
> | 
    beginTorsionIndex = 0; | 
| 836 | 
> | 
    beginInversionIndex = 0; | 
| 837 | 
> | 
    | 
| 838 | 
> | 
    for(int i = 0; i < info->getNGlobalMolecules(); i++) { | 
| 839 | 
  | 
       | 
| 840 | 
< | 
      //local index(index in DataStorge) of atom is important | 
| 841 | 
< | 
      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 842 | 
< | 
        atom->setGlobalIndex(beginAtomIndex++); | 
| 840 | 
> | 
#ifdef IS_MPI       | 
| 841 | 
> | 
      if (info->getMolToProc(i) == worldRank) { | 
| 842 | 
> | 
#endif         | 
| 843 | 
> | 
        // stuff to do if I own this molecule | 
| 844 | 
> | 
        mol = info->getMoleculeByGlobalIndex(i); | 
| 845 | 
> | 
 | 
| 846 | 
> | 
        // The local index(index in DataStorge) of the atom is important: | 
| 847 | 
> | 
        for(atom = mol->beginAtom(ai); atom != NULL;  | 
| 848 | 
> | 
            atom = mol->nextAtom(ai)) { | 
| 849 | 
> | 
          atom->setGlobalIndex(beginAtomIndex++); | 
| 850 | 
> | 
        } | 
| 851 | 
> | 
         | 
| 852 | 
> | 
        for(rb = mol->beginRigidBody(ri); rb != NULL; | 
| 853 | 
> | 
            rb = mol->nextRigidBody(ri)) { | 
| 854 | 
> | 
          rb->setGlobalIndex(beginRigidBodyIndex++); | 
| 855 | 
> | 
        } | 
| 856 | 
> | 
         | 
| 857 | 
> | 
        // The local index of other objects only depends on the order | 
| 858 | 
> | 
        // of traversal: | 
| 859 | 
> | 
        for(cg = mol->beginCutoffGroup(ci); cg != NULL; | 
| 860 | 
> | 
            cg = mol->nextCutoffGroup(ci)) { | 
| 861 | 
> | 
          cg->setGlobalIndex(beginCutoffGroupIndex++); | 
| 862 | 
> | 
        }         | 
| 863 | 
> | 
        for(bond = mol->beginBond(boi); bond != NULL; | 
| 864 | 
> | 
            bond = mol->nextBond(boi)) { | 
| 865 | 
> | 
          bond->setGlobalIndex(beginBondIndex++); | 
| 866 | 
> | 
        }         | 
| 867 | 
> | 
        for(bend = mol->beginBend(bei); bend != NULL; | 
| 868 | 
> | 
            bend = mol->nextBend(bei)) { | 
| 869 | 
> | 
          bend->setGlobalIndex(beginBendIndex++); | 
| 870 | 
> | 
        }         | 
| 871 | 
> | 
        for(torsion = mol->beginTorsion(ti); torsion != NULL; | 
| 872 | 
> | 
            torsion = mol->nextTorsion(ti)) { | 
| 873 | 
> | 
          torsion->setGlobalIndex(beginTorsionIndex++); | 
| 874 | 
> | 
        }         | 
| 875 | 
> | 
        for(inversion = mol->beginInversion(ii); inversion != NULL; | 
| 876 | 
> | 
            inversion = mol->nextInversion(ii)) { | 
| 877 | 
> | 
          inversion->setGlobalIndex(beginInversionIndex++); | 
| 878 | 
> | 
        }         | 
| 879 | 
> | 
         | 
| 880 | 
> | 
#ifdef IS_MPI         | 
| 881 | 
> | 
      }  else { | 
| 882 | 
> | 
 | 
| 883 | 
> | 
        // stuff to do if I don't own this molecule | 
| 884 | 
> | 
         | 
| 885 | 
> | 
        int stampId = info->getMoleculeStampId(i); | 
| 886 | 
> | 
        MoleculeStamp* stamp = info->getMoleculeStamp(stampId); | 
| 887 | 
> | 
 | 
| 888 | 
> | 
        beginAtomIndex += stamp->getNAtoms(); | 
| 889 | 
> | 
        beginRigidBodyIndex += stamp->getNRigidBodies(); | 
| 890 | 
> | 
        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms(); | 
| 891 | 
> | 
        beginBondIndex += stamp->getNBonds(); | 
| 892 | 
> | 
        beginBendIndex += stamp->getNBends(); | 
| 893 | 
> | 
        beginTorsionIndex += stamp->getNTorsions(); | 
| 894 | 
> | 
        beginInversionIndex += stamp->getNInversions(); | 
| 895 | 
  | 
      } | 
| 896 | 
< | 
       | 
| 897 | 
< | 
      for(rb = mol->beginRigidBody(ri); rb != NULL; | 
| 898 | 
< | 
          rb = mol->nextRigidBody(ri)) { | 
| 899 | 
< | 
        rb->setGlobalIndex(beginRigidBodyIndex++); | 
| 527 | 
< | 
      } | 
| 528 | 
< | 
       | 
| 529 | 
< | 
      //local index of cutoff group is trivial, it only depends on the order of travesing | 
| 530 | 
< | 
      for(cg = mol->beginCutoffGroup(ci); cg != NULL; | 
| 531 | 
< | 
          cg = mol->nextCutoffGroup(ci)) { | 
| 532 | 
< | 
        cg->setGlobalIndex(beginCutoffGroupIndex++); | 
| 533 | 
< | 
      } | 
| 534 | 
< | 
    } | 
| 535 | 
< | 
     | 
| 896 | 
> | 
#endif           | 
| 897 | 
> | 
 | 
| 898 | 
> | 
    } //end for(int i=0)   | 
| 899 | 
> | 
 | 
| 900 | 
  | 
    //fill globalGroupMembership | 
| 901 | 
  | 
    std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0); | 
| 902 | 
< | 
    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {         | 
| 903 | 
< | 
      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { | 
| 904 | 
< | 
         | 
| 902 | 
> | 
    for(mol = info->beginMolecule(mi); mol != NULL;  | 
| 903 | 
> | 
        mol = info->nextMolecule(mi)) {         | 
| 904 | 
> | 
      for (cg = mol->beginCutoffGroup(ci); cg != NULL;  | 
| 905 | 
> | 
           cg = mol->nextCutoffGroup(ci)) {         | 
| 906 | 
  | 
        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | 
| 907 | 
  | 
          globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex(); | 
| 908 | 
  | 
        } | 
| 909 | 
  | 
         | 
| 910 | 
  | 
      }        | 
| 911 | 
  | 
    } | 
| 912 | 
< | 
     | 
| 912 | 
> | 
    | 
| 913 | 
  | 
#ifdef IS_MPI     | 
| 914 | 
  | 
    // Since the globalGroupMembership has been zero filled and we've only | 
| 915 | 
  | 
    // poked values into the atoms we know, we can do an Allreduce | 
| 916 | 
  | 
    // to get the full globalGroupMembership array (We think). | 
| 917 | 
  | 
    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 | 
| 918 | 
  | 
    // docs said we could. | 
| 919 | 
< | 
    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0); | 
| 920 | 
< | 
    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms, | 
| 921 | 
< | 
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 919 | 
> | 
    std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0); | 
| 920 | 
> | 
    MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],  | 
| 921 | 
> | 
                              &tmpGroupMembership[0], nGlobalAtoms, | 
| 922 | 
> | 
                              MPI::INT, MPI::SUM); | 
| 923 | 
  | 
    info->setGlobalGroupMembership(tmpGroupMembership); | 
| 924 | 
  | 
#else | 
| 925 | 
  | 
    info->setGlobalGroupMembership(globalGroupMembership); | 
| 926 | 
  | 
#endif | 
| 927 | 
  | 
     | 
| 928 | 
  | 
    //fill molMembership | 
| 929 | 
< | 
    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0); | 
| 929 | 
> | 
    std::vector<int> globalMolMembership(info->getNGlobalAtoms() +  | 
| 930 | 
> | 
                                         info->getNGlobalRigidBodies(), 0); | 
| 931 | 
  | 
     | 
| 932 | 
< | 
    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 933 | 
< | 
       | 
| 932 | 
> | 
    for(mol = info->beginMolecule(mi); mol != NULL;  | 
| 933 | 
> | 
        mol = info->nextMolecule(mi)) { | 
| 934 | 
  | 
      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 935 | 
  | 
        globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex(); | 
| 936 | 
  | 
      } | 
| 937 | 
+ | 
      for (rb = mol->beginRigidBody(ri); rb != NULL;  | 
| 938 | 
+ | 
           rb = mol->nextRigidBody(ri)) { | 
| 939 | 
+ | 
        globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex(); | 
| 940 | 
+ | 
      } | 
| 941 | 
  | 
    } | 
| 942 | 
  | 
     | 
| 943 | 
  | 
#ifdef IS_MPI | 
| 944 | 
< | 
    std::vector<int> tmpMolMembership(nGlobalAtoms, 0); | 
| 944 | 
> | 
    std::vector<int> tmpMolMembership(info->getNGlobalAtoms() +  | 
| 945 | 
> | 
                                      info->getNGlobalRigidBodies(), 0); | 
| 946 | 
> | 
    MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],  | 
| 947 | 
> | 
                              nGlobalAtoms + nGlobalRigidBodies, | 
| 948 | 
> | 
                              MPI::INT, MPI::SUM); | 
| 949 | 
  | 
     | 
| 575 | 
– | 
    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms, | 
| 576 | 
– | 
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 577 | 
– | 
     | 
| 950 | 
  | 
    info->setGlobalMolMembership(tmpMolMembership); | 
| 951 | 
  | 
#else | 
| 952 | 
  | 
    info->setGlobalMolMembership(globalMolMembership); | 
| 953 | 
  | 
#endif | 
| 954 | 
+ | 
 | 
| 955 | 
+ | 
    // nIOPerMol holds the number of integrable objects per molecule | 
| 956 | 
+ | 
    // here the molecules are listed by their global indices. | 
| 957 | 
+ | 
 | 
| 958 | 
+ | 
    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0); | 
| 959 | 
+ | 
    for (mol = info->beginMolecule(mi); mol != NULL; | 
| 960 | 
+ | 
         mol = info->nextMolecule(mi)) { | 
| 961 | 
+ | 
      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();        | 
| 962 | 
+ | 
    } | 
| 963 | 
  | 
     | 
| 964 | 
< | 
  } | 
| 965 | 
< | 
   | 
| 966 | 
< | 
  void SimCreator::loadCoordinates(SimInfo* info) { | 
| 967 | 
< | 
    Globals* simParams; | 
| 968 | 
< | 
    simParams = info->getSimParams(); | 
| 964 | 
> | 
#ifdef IS_MPI | 
| 965 | 
> | 
    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0); | 
| 966 | 
> | 
    MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],  | 
| 967 | 
> | 
                              info->getNGlobalMolecules(), MPI::INT, MPI::SUM); | 
| 968 | 
> | 
#else | 
| 969 | 
> | 
    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol; | 
| 970 | 
> | 
#endif     | 
| 971 | 
> | 
 | 
| 972 | 
> | 
    std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules()); | 
| 973 | 
  | 
     | 
| 974 | 
< | 
    if (!simParams->haveInitialConfig()) { | 
| 975 | 
< | 
      sprintf(painCave.errMsg, | 
| 976 | 
< | 
              "Cannot intialize a simulation without an initial configuration file.\n"); | 
| 977 | 
< | 
      painCave.isFatal = 1;; | 
| 593 | 
< | 
      simError(); | 
| 974 | 
> | 
    int startingIndex = 0; | 
| 975 | 
> | 
    for (int i = 0; i < info->getNGlobalMolecules(); i++) { | 
| 976 | 
> | 
      startingIOIndexForMol[i] = startingIndex; | 
| 977 | 
> | 
      startingIndex += numIntegrableObjectsPerMol[i]; | 
| 978 | 
  | 
    } | 
| 979 | 
  | 
     | 
| 980 | 
< | 
    DumpReader reader(info, simParams->getInitialConfig()); | 
| 980 | 
> | 
    std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); | 
| 981 | 
> | 
    for (mol = info->beginMolecule(mi); mol != NULL;  | 
| 982 | 
> | 
         mol = info->nextMolecule(mi)) { | 
| 983 | 
> | 
      int myGlobalIndex = mol->getGlobalIndex(); | 
| 984 | 
> | 
      int globalIO = startingIOIndexForMol[myGlobalIndex]; | 
| 985 | 
> | 
      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL; | 
| 986 | 
> | 
           sd = mol->nextIntegrableObject(ioi)) { | 
| 987 | 
> | 
        sd->setGlobalIntegrableObjectIndex(globalIO); | 
| 988 | 
> | 
        IOIndexToIntegrableObject[globalIO] = sd; | 
| 989 | 
> | 
        globalIO++; | 
| 990 | 
> | 
      } | 
| 991 | 
> | 
    } | 
| 992 | 
> | 
       | 
| 993 | 
> | 
    info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject); | 
| 994 | 
> | 
     | 
| 995 | 
> | 
  } | 
| 996 | 
> | 
   | 
| 997 | 
> | 
  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { | 
| 998 | 
> | 
     | 
| 999 | 
> | 
    DumpReader reader(info, mdFileName); | 
| 1000 | 
  | 
    int nframes = reader.getNFrames(); | 
| 1001 | 
  | 
     | 
| 1002 | 
  | 
    if (nframes > 0) { | 
| 1005 | 
  | 
      //invalid initial coordinate file | 
| 1006 | 
  | 
      sprintf(painCave.errMsg,  | 
| 1007 | 
  | 
              "Initial configuration file %s should at least contain one frame\n", | 
| 1008 | 
< | 
              simParams->getInitialConfig()); | 
| 1008 | 
> | 
              mdFileName.c_str()); | 
| 1009 | 
  | 
      painCave.isFatal = 1; | 
| 1010 | 
  | 
      simError(); | 
| 1011 | 
  | 
    } | 
| 609 | 
– | 
     | 
| 1012 | 
  | 
    //copy the current snapshot to previous snapshot | 
| 1013 | 
  | 
    info->getSnapshotManager()->advance(); | 
| 1014 | 
  | 
  } | 
| 1015 | 
  | 
   | 
| 1016 | 
< | 
} //end namespace oopse | 
| 1016 | 
> | 
} //end namespace OpenMD | 
| 1017 | 
  | 
 | 
| 1018 | 
  | 
 |