| 1 | chrisfen | 2101 | /* | 
| 2 |  |  | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. | 
| 3 |  |  | * | 
| 4 |  |  | * The University of Notre Dame grants you ("Licensee") a | 
| 5 |  |  | * non-exclusive, royalty free, license to use, modify and | 
| 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 | 
| 19 |  |  | *    notice, this list of conditions and the following disclaimer. | 
| 20 |  |  | * | 
| 21 |  |  | * 3. Redistributions in binary form must reproduce the above copyright | 
| 22 |  |  | *    notice, this list of conditions and the following disclaimer in the | 
| 23 |  |  | *    documentation and/or other materials provided with the | 
| 24 |  |  | *    distribution. | 
| 25 |  |  | * | 
| 26 |  |  | * This software is provided "AS IS," without a warranty of any | 
| 27 |  |  | * kind. All express or implied conditions, representations and | 
| 28 |  |  | * warranties, including any implied warranty of merchantability, | 
| 29 |  |  | * fitness for a particular purpose or non-infringement, are hereby | 
| 30 |  |  | * excluded.  The University of Notre Dame and its licensors shall not | 
| 31 |  |  | * be liable for any damages suffered by licensee as a result of | 
| 32 |  |  | * using, modifying or distributing the software or its | 
| 33 |  |  | * derivatives. In no event will the University of Notre Dame or its | 
| 34 |  |  | * licensors be liable for any lost revenue, profit or data, or for | 
| 35 |  |  | * direct, indirect, special, consequential, incidental or punitive | 
| 36 |  |  | * damages, however caused and regardless of the theory of liability, | 
| 37 |  |  | * arising out of the use of or inability to use software, even if the | 
| 38 |  |  | * University of Notre Dame has been advised of the possibility of | 
| 39 |  |  | * such damages. | 
| 40 |  |  | */ | 
| 41 |  |  |  | 
| 42 |  |  | #define _LARGEFILE_SOURCE64 | 
| 43 |  |  | #define _FILE_OFFSET_BITS 64 | 
| 44 |  |  |  | 
| 45 |  |  | #include <sys/types.h> | 
| 46 |  |  | #include <sys/stat.h> | 
| 47 |  |  |  | 
| 48 |  |  | #include <iostream> | 
| 49 |  |  | #include <math.h> | 
| 50 |  |  |  | 
| 51 |  |  | #include <stdio.h> | 
| 52 |  |  | #include <stdlib.h> | 
| 53 |  |  | #include <string.h> | 
| 54 |  |  |  | 
| 55 |  |  | #include "primitives/Molecule.hpp" | 
| 56 |  |  | #include "utils/MemoryUtils.hpp" | 
| 57 |  |  | #include "utils/StringTokenizer.hpp" | 
| 58 |  |  | #include "io/RestReader.hpp" | 
| 59 |  |  | #include "utils/simError.h" | 
| 60 |  |  |  | 
| 61 |  |  | #ifdef IS_MPI | 
| 62 |  |  | #include <mpi.h> | 
| 63 |  |  | #define TAKE_THIS_TAG_CHAR 0 | 
| 64 |  |  | #define TAKE_THIS_TAG_INT 1 | 
| 65 |  |  | #define TAKE_THIS_TAG_DOUBLE 2 | 
| 66 |  |  | #endif // is_mpi | 
| 67 |  |  |  | 
| 68 |  |  | namespace oopse { | 
| 69 |  |  |  | 
| 70 |  |  | RestReader::RestReader( SimInfo* info ) : info_(info){ | 
| 71 | chrisfen | 2115 |  | 
| 72 | chrisfen | 2101 | idealName = "idealCrystal.in"; | 
| 73 | chrisfen | 2903 |  | 
| 74 | chrisfen | 2101 | #ifdef IS_MPI | 
| 75 |  |  | if (worldRank == 0) { | 
| 76 |  |  | #endif | 
| 77 | chrisfen | 2903 |  | 
| 78 |  |  | inIdealFile = new std::ifstream(idealName.c_str()); | 
| 79 |  |  |  | 
| 80 |  |  | if(inIdealFile->fail()){ | 
| 81 | chrisfen | 2101 | sprintf(painCave.errMsg, | 
| 82 | chrisfen | 2903 | "RestReader: Cannot open file: %s\n", | 
| 83 |  |  | idealName.c_str()); | 
| 84 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 85 |  |  | simError(); | 
| 86 |  |  | } | 
| 87 |  |  |  | 
| 88 |  |  | #ifdef IS_MPI | 
| 89 |  |  | } | 
| 90 |  |  | strcpy( checkPointMsg, | 
| 91 |  |  | "File \"idealCrystal.in\" opened successfully for reading." ); | 
| 92 |  |  | MPIcheckPoint(); | 
| 93 |  |  | #endif | 
| 94 | chrisfen | 2115 |  | 
| 95 | chrisfen | 2101 | return; | 
| 96 |  |  | } | 
| 97 |  |  |  | 
| 98 |  |  | RestReader :: ~RestReader( ){ | 
| 99 |  |  | #ifdef IS_MPI | 
| 100 |  |  | if (worldRank == 0) { | 
| 101 |  |  | #endif | 
| 102 | chrisfen | 2903 |  | 
| 103 |  |  | delete inIdealFile; | 
| 104 |  |  | delete inAngFile; | 
| 105 |  |  |  | 
| 106 | chrisfen | 2101 | #ifdef IS_MPI | 
| 107 |  |  | } | 
| 108 | chrisfen | 2903 | strcpy( checkPointMsg, | 
| 109 |  |  | "File idealCrystal.in (and .zang0 if present) closed successfully." ); | 
| 110 | chrisfen | 2101 | MPIcheckPoint(); | 
| 111 |  |  | #endif | 
| 112 |  |  |  | 
| 113 |  |  | return; | 
| 114 |  |  | } | 
| 115 |  |  |  | 
| 116 |  |  |  | 
| 117 |  |  | void RestReader :: readIdealCrystal(){ | 
| 118 | chrisfen | 2868 |  | 
| 119 | chrisfen | 2101 | #ifdef IS_MPI | 
| 120 | chrisfen | 2903 | int which_node; | 
| 121 |  |  | int i, j; | 
| 122 | chrisfen | 2101 | #endif //is_mpi | 
| 123 |  |  |  | 
| 124 |  |  | const int BUFFERSIZE = 2000; // size of the read buffer | 
| 125 |  |  | int nTotObjs; // the number of atoms | 
| 126 |  |  | char read_buffer[BUFFERSIZE]; //the line buffer for reading | 
| 127 |  |  |  | 
| 128 |  |  | char *parseErr; | 
| 129 |  |  |  | 
| 130 |  |  | std::vector<StuntDouble*> integrableObjects; | 
| 131 |  |  |  | 
| 132 |  |  | Molecule* mol; | 
| 133 |  |  | StuntDouble* integrableObject; | 
| 134 |  |  | SimInfo::MoleculeIterator mi; | 
| 135 |  |  | Molecule::IntegrableObjectIterator ii; | 
| 136 |  |  |  | 
| 137 |  |  | #ifndef IS_MPI | 
| 138 |  |  |  | 
| 139 | chrisfen | 2903 | inIdealFile->getline(read_buffer, sizeof(read_buffer)); | 
| 140 |  |  |  | 
| 141 |  |  | if( inIdealFile->eof() ){ | 
| 142 | chrisfen | 2101 | sprintf( painCave.errMsg, | 
| 143 |  |  | "RestraintReader error: error reading 1st line of \"%s\"\n", | 
| 144 | chrisfen | 2903 | idealName.c_str() ); | 
| 145 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 146 |  |  | simError(); | 
| 147 |  |  | } | 
| 148 |  |  |  | 
| 149 |  |  | nTotObjs = atoi( read_buffer ); | 
| 150 |  |  |  | 
| 151 |  |  | if( nTotObjs != info_->getNGlobalIntegrableObjects() ){ | 
| 152 |  |  | sprintf( painCave.errMsg, | 
| 153 |  |  | "RestraintReader error. %s nIntegrable, %d, " | 
| 154 |  |  | "does not match the meta-data file's nIntegrable, %d.\n", | 
| 155 | chrisfen | 2903 | idealName.c_str(), | 
| 156 |  |  | nTotObjs, | 
| 157 | chrisfen | 2101 | info_->getNGlobalIntegrableObjects()); | 
| 158 |  |  | painCave.isFatal = 1; | 
| 159 |  |  | simError(); | 
| 160 |  |  | } | 
| 161 |  |  |  | 
| 162 |  |  | // skip over the comment line | 
| 163 | chrisfen | 2903 | inIdealFile->getline(read_buffer, sizeof(read_buffer)); | 
| 164 |  |  |  | 
| 165 |  |  | if( inIdealFile->eof() ){ | 
| 166 | chrisfen | 2101 | sprintf( painCave.errMsg, | 
| 167 | chrisfen | 2903 | "error in reading commment in %s\n", | 
| 168 |  |  | idealName.c_str()); | 
| 169 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 170 |  |  | simError(); | 
| 171 |  |  | } | 
| 172 |  |  |  | 
| 173 |  |  | // parse the ideal crystal lines | 
| 174 |  |  | /* | 
| 175 |  |  | * Note: we assume that there is a one-to-one correspondence between | 
| 176 |  |  | * integrable objects and lines in the idealCrystal.in file.  Thermodynamic | 
| 177 |  |  | * integration is only supported for simple rigid bodies. | 
| 178 |  |  | */ | 
| 179 |  |  | for (mol = info_->beginMolecule(mi); mol != NULL; | 
| 180 |  |  | mol = info_->nextMolecule(mi)) { | 
| 181 |  |  |  | 
| 182 |  |  | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 183 |  |  | integrableObject != NULL; | 
| 184 |  |  | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 185 | chrisfen | 2903 |  | 
| 186 |  |  | inIdealFile->getline(read_buffer, sizeof(read_buffer)); | 
| 187 |  |  |  | 
| 188 |  |  | if( inIdealFile->eof() ){ | 
| 189 | chrisfen | 2101 | sprintf(painCave.errMsg, | 
| 190 |  |  | "RestReader Error: error in reading file %s\n" | 
| 191 |  |  | "natoms  = %d; index = %d\n" | 
| 192 |  |  | "error reading the line from the file.\n", | 
| 193 | chrisfen | 2903 | idealName.c_str(), nTotObjs, | 
| 194 | chrisfen | 2101 | integrableObject->getGlobalIndex() ); | 
| 195 |  |  | painCave.isFatal = 1; | 
| 196 |  |  | simError(); | 
| 197 |  |  | } | 
| 198 |  |  |  | 
| 199 |  |  | parseErr = parseIdealLine( read_buffer, integrableObject); | 
| 200 |  |  | if( parseErr != NULL ){ | 
| 201 |  |  | strcpy( painCave.errMsg, parseErr ); | 
| 202 |  |  | painCave.isFatal = 1; | 
| 203 |  |  | simError(); | 
| 204 |  |  | } | 
| 205 |  |  | } | 
| 206 |  |  | } | 
| 207 |  |  |  | 
| 208 |  |  | // MPI Section of code.......... | 
| 209 |  |  | #else //IS_MPI | 
| 210 |  |  |  | 
| 211 |  |  | // first thing first, suspend fatalities. | 
| 212 |  |  | painCave.isEventLoop = 1; | 
| 213 |  |  |  | 
| 214 |  |  | int masterNode = 0; | 
| 215 |  |  |  | 
| 216 |  |  | MPI_Status istatus; | 
| 217 |  |  | int nCurObj; | 
| 218 |  |  | int nitems; | 
| 219 | chrisfen | 2107 | int haveError; | 
| 220 |  |  |  | 
| 221 |  |  | nTotObjs = info_->getNGlobalIntegrableObjects(); | 
| 222 | chrisfen | 2101 | haveError = 0; | 
| 223 | chrisfen | 2115 |  | 
| 224 | chrisfen | 2101 | if (worldRank == masterNode) { | 
| 225 | chrisfen | 2903 | inIdealFile->getline(read_buffer, sizeof(read_buffer)); | 
| 226 |  |  |  | 
| 227 |  |  | if( inIdealFile->eof() ){ | 
| 228 | chrisfen | 2101 | sprintf( painCave.errMsg, | 
| 229 | chrisfen | 2903 | "Error reading 1st line of %s \n ",idealName.c_str()); | 
| 230 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 231 |  |  | simError(); | 
| 232 |  |  | } | 
| 233 |  |  |  | 
| 234 |  |  | nitems = atoi( read_buffer ); | 
| 235 |  |  |  | 
| 236 |  |  | // Check to see that the number of integrable objects in the | 
| 237 |  |  | // intial configuration file is the same as derived from the | 
| 238 |  |  | // meta-data file. | 
| 239 |  |  | if( nTotObjs != nitems){ | 
| 240 |  |  | sprintf( painCave.errMsg, | 
| 241 |  |  | "RestraintReader Error. %s nIntegrable, %d, " | 
| 242 |  |  | "does not match the meta-data file's nIntegrable, %d.\n", | 
| 243 | chrisfen | 2903 | idealName.c_str(), nTotObjs, | 
| 244 | chrisfen | 2101 | info_->getNGlobalIntegrableObjects()); | 
| 245 |  |  | painCave.isFatal = 1; | 
| 246 |  |  | simError(); | 
| 247 |  |  | } | 
| 248 |  |  |  | 
| 249 |  |  | // skip over the comment line | 
| 250 | chrisfen | 2903 | inIdealFile->getline(read_buffer, sizeof(read_buffer)); | 
| 251 |  |  |  | 
| 252 |  |  | if( inIdealFile->eof() ){ | 
| 253 | chrisfen | 2101 | sprintf( painCave.errMsg, | 
| 254 | chrisfen | 2903 | "error in reading commment in %s\n", idealName.c_str()); | 
| 255 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 256 |  |  | simError(); | 
| 257 |  |  | } | 
| 258 | chrisfen | 2115 |  | 
| 259 | chrisfen | 2868 | MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD); | 
| 260 |  |  |  | 
| 261 | chrisfen | 2101 | for (i=0 ; i < info_->getNGlobalMolecules(); i++) { | 
| 262 | chrisfen | 2903 | which_node = info_->getMolToProc(i); | 
| 263 | chrisfen | 2101 |  | 
| 264 |  |  | if(which_node == masterNode){ | 
| 265 |  |  | //molecules belong to master node | 
| 266 |  |  |  | 
| 267 | chrisfen | 2107 | mol = info_->getMoleculeByGlobalIndex(i); | 
| 268 | chrisfen | 2101 |  | 
| 269 | chrisfen | 2107 | if(mol == NULL) { | 
| 270 |  |  | sprintf(painCave.errMsg, | 
| 271 | gezelter | 2204 | "RestReader Error: Molecule not found on node %d!\n", | 
| 272 |  |  | worldRank); | 
| 273 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 274 |  |  | simError(); | 
| 275 |  |  | } | 
| 276 |  |  |  | 
| 277 |  |  | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 278 |  |  | integrableObject != NULL; | 
| 279 |  |  | integrableObject = mol->nextIntegrableObject(ii)){ | 
| 280 |  |  |  | 
| 281 | chrisfen | 2903 | inIdealFile->getline(read_buffer, sizeof(read_buffer)); | 
| 282 | chrisfen | 2101 |  | 
| 283 | chrisfen | 2903 | if( inIdealFile->eof() ){ | 
| 284 | chrisfen | 2101 | sprintf(painCave.errMsg, | 
| 285 |  |  | "RestReader Error: error in reading file %s\n" | 
| 286 |  |  | "natoms  = %d; index = %d\n" | 
| 287 |  |  | "error reading the line from the file.\n", | 
| 288 | chrisfen | 2903 | idealName.c_str(), nTotObjs, i ); | 
| 289 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 290 |  |  | simError(); | 
| 291 |  |  | } | 
| 292 | chrisfen | 2115 |  | 
| 293 |  |  | parseIdealLine(read_buffer, integrableObject); | 
| 294 |  |  |  | 
| 295 | chrisfen | 2101 | } | 
| 296 | chrisfen | 2868 |  | 
| 297 | chrisfen | 2101 | } else { | 
| 298 |  |  | //molecule belongs to slave nodes | 
| 299 |  |  |  | 
| 300 |  |  | MPI_Recv(&nCurObj, 1, MPI_INT, which_node, | 
| 301 |  |  | TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus); | 
| 302 |  |  |  | 
| 303 | chrisfen | 2868 | for(j = 0; j < nCurObj; j++){ | 
| 304 | chrisfen | 2903 | inIdealFile->getline(read_buffer, sizeof(read_buffer)); | 
| 305 |  |  |  | 
| 306 |  |  | if( inIdealFile->eof() ){ | 
| 307 | chrisfen | 2101 | sprintf(painCave.errMsg, | 
| 308 |  |  | "RestReader Error: error in reading file %s\n" | 
| 309 |  |  | "natoms  = %d; index = %d\n" | 
| 310 |  |  | "error reading the line from the file.\n", | 
| 311 | chrisfen | 2903 | idealName.c_str(), nTotObjs, i ); | 
| 312 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 313 |  |  | simError(); | 
| 314 |  |  | } | 
| 315 |  |  |  | 
| 316 | chrisfen | 2107 | MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node, | 
| 317 | chrisfen | 2101 | TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD); | 
| 318 |  |  | } | 
| 319 |  |  | } | 
| 320 |  |  | } | 
| 321 |  |  | } else { | 
| 322 | chrisfen | 2868 | //actions taken at slave nodes | 
| 323 |  |  | MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD); | 
| 324 |  |  |  | 
| 325 | chrisfen | 2101 | for (i=0 ; i < info_->getNGlobalMolecules(); i++) { | 
| 326 | chrisfen | 2868 | int which_node = info_->getMolToProc(i); | 
| 327 |  |  |  | 
| 328 | chrisfen | 2101 | if(which_node == worldRank){ | 
| 329 |  |  | //molecule with global index i belongs to this processor | 
| 330 |  |  |  | 
| 331 | chrisfen | 2107 | mol = info_->getMoleculeByGlobalIndex(i); | 
| 332 | chrisfen | 2101 |  | 
| 333 | chrisfen | 2107 | if(mol == NULL) { | 
| 334 | chrisfen | 2101 | sprintf(painCave.errMsg, | 
| 335 |  |  | "RestReader Error: molecule not found on node %d\n", | 
| 336 |  |  | worldRank); | 
| 337 |  |  | painCave.isFatal = 1; | 
| 338 |  |  | simError(); | 
| 339 |  |  | } | 
| 340 |  |  |  | 
| 341 | chrisfen | 2107 | nCurObj = mol->getNIntegrableObjects(); | 
| 342 | chrisfen | 2101 |  | 
| 343 |  |  | MPI_Send(&nCurObj, 1, MPI_INT, masterNode, | 
| 344 |  |  | TAKE_THIS_TAG_INT, MPI_COMM_WORLD); | 
| 345 |  |  |  | 
| 346 |  |  | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 347 |  |  | integrableObject != NULL; | 
| 348 |  |  | integrableObject = mol->nextIntegrableObject(ii)){ | 
| 349 |  |  |  | 
| 350 |  |  | MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, | 
| 351 |  |  | TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus); | 
| 352 |  |  |  | 
| 353 |  |  | parseErr = parseIdealLine(read_buffer, integrableObject); | 
| 354 |  |  |  | 
| 355 |  |  | if( parseErr != NULL ){ | 
| 356 |  |  | strcpy( painCave.errMsg, parseErr ); | 
| 357 |  |  | simError(); | 
| 358 |  |  | } | 
| 359 |  |  | } | 
| 360 |  |  | } | 
| 361 |  |  | } | 
| 362 |  |  | } | 
| 363 |  |  | #endif | 
| 364 |  |  | } | 
| 365 |  |  |  | 
| 366 |  |  | char* RestReader::parseIdealLine(char* readLine, StuntDouble* sd){ | 
| 367 |  |  |  | 
| 368 | tim | 2759 | RealType pos[3];        // position place holders | 
| 369 |  |  | RealType q[4];          // the quaternions | 
| 370 |  |  | RealType RfromQ[3][3];  // the rotation matrix | 
| 371 |  |  | RealType normalize;     // to normalize the reference unit vector | 
| 372 |  |  | RealType uX, uY, uZ;    // reference unit vector place holders | 
| 373 |  |  | RealType uselessToken; | 
| 374 | chrisfen | 2101 | StringTokenizer tokenizer(readLine); | 
| 375 |  |  | int nTokens; | 
| 376 |  |  |  | 
| 377 |  |  | nTokens = tokenizer.countTokens(); | 
| 378 | chrisfen | 2115 |  | 
| 379 | chrisfen | 2101 | if (nTokens < 14) { | 
| 380 |  |  | sprintf(painCave.errMsg, | 
| 381 |  |  | "RestReader Error: Not enough Tokens.\n"); | 
| 382 |  |  | painCave.isFatal = 1; | 
| 383 |  |  | simError(); | 
| 384 |  |  | } | 
| 385 |  |  |  | 
| 386 |  |  | std::string name = tokenizer.nextToken(); | 
| 387 | chrisfen | 2115 |  | 
| 388 | chrisfen | 2101 | if (name != sd->getType()) { | 
| 389 |  |  |  | 
| 390 |  |  | sprintf(painCave.errMsg, | 
| 391 |  |  | "RestReader Error: Atom type [%s] in %s does not " | 
| 392 |  |  | "match Atom Type [%s] in .md file.\n", | 
| 393 | chrisfen | 2903 | name.c_str(), idealName.c_str(), | 
| 394 | chrisfen | 2101 | sd->getType().c_str()); | 
| 395 |  |  | painCave.isFatal = 1; | 
| 396 |  |  | simError(); | 
| 397 |  |  | } | 
| 398 |  |  |  | 
| 399 |  |  | pos[0] = tokenizer.nextTokenAsDouble(); | 
| 400 |  |  | pos[1] = tokenizer.nextTokenAsDouble(); | 
| 401 |  |  | pos[2] = tokenizer.nextTokenAsDouble(); | 
| 402 | chrisfen | 2115 |  | 
| 403 | chrisfen | 2101 | // store the positions in the stuntdouble as generic data doubles | 
| 404 |  |  | DoubleGenericData* refPosX = new DoubleGenericData(); | 
| 405 |  |  | refPosX->setID("refPosX"); | 
| 406 |  |  | refPosX->setData(pos[0]); | 
| 407 |  |  | sd->addProperty(refPosX); | 
| 408 | chrisfen | 2115 |  | 
| 409 | chrisfen | 2101 | DoubleGenericData* refPosY = new DoubleGenericData(); | 
| 410 |  |  | refPosY->setID("refPosY"); | 
| 411 |  |  | refPosY->setData(pos[1]); | 
| 412 |  |  | sd->addProperty(refPosY); | 
| 413 |  |  |  | 
| 414 |  |  | DoubleGenericData* refPosZ = new DoubleGenericData(); | 
| 415 |  |  | refPosZ->setID("refPosZ"); | 
| 416 |  |  | refPosZ->setData(pos[2]); | 
| 417 |  |  | sd->addProperty(refPosZ); | 
| 418 | chrisfen | 2115 |  | 
| 419 | chrisfen | 2101 | // we don't need the velocities | 
| 420 |  |  | uselessToken = tokenizer.nextTokenAsDouble(); | 
| 421 |  |  | uselessToken = tokenizer.nextTokenAsDouble(); | 
| 422 |  |  | uselessToken = tokenizer.nextTokenAsDouble(); | 
| 423 |  |  |  | 
| 424 |  |  | if (sd->isDirectional()) { | 
| 425 |  |  |  | 
| 426 |  |  | q[0] = tokenizer.nextTokenAsDouble(); | 
| 427 |  |  | q[1] = tokenizer.nextTokenAsDouble(); | 
| 428 |  |  | q[2] = tokenizer.nextTokenAsDouble(); | 
| 429 |  |  | q[3] = tokenizer.nextTokenAsDouble(); | 
| 430 |  |  |  | 
| 431 |  |  | // now build the rotation matrix and find the unit vectors | 
| 432 |  |  | RfromQ[0][0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3]; | 
| 433 |  |  | RfromQ[0][1] = 2*(q[1]*q[2] + q[0]*q[3]); | 
| 434 |  |  | RfromQ[0][2] = 2*(q[1]*q[3] - q[0]*q[2]); | 
| 435 |  |  | RfromQ[1][0] = 2*(q[1]*q[2] - q[0]*q[3]); | 
| 436 |  |  | RfromQ[1][1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3]; | 
| 437 |  |  | RfromQ[1][2] = 2*(q[2]*q[3] + q[0]*q[1]); | 
| 438 |  |  | RfromQ[2][0] = 2*(q[1]*q[3] + q[0]*q[2]); | 
| 439 |  |  | RfromQ[2][1] = 2*(q[2]*q[3] - q[0]*q[1]); | 
| 440 |  |  | RfromQ[2][2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3]; | 
| 441 |  |  |  | 
| 442 |  |  | normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1] | 
| 443 |  |  | + RfromQ[2][2]*RfromQ[2][2]); | 
| 444 |  |  | uX = RfromQ[2][0]/normalize; | 
| 445 |  |  | uY = RfromQ[2][1]/normalize; | 
| 446 |  |  | uZ = RfromQ[2][2]/normalize; | 
| 447 |  |  |  | 
| 448 |  |  | // store reference unit vectors as generic data in the stuntdouble | 
| 449 |  |  | DoubleGenericData* refVectorX = new DoubleGenericData(); | 
| 450 |  |  | refVectorX->setID("refVectorX"); | 
| 451 |  |  | refVectorX->setData(uX); | 
| 452 |  |  | sd->addProperty(refVectorX); | 
| 453 |  |  |  | 
| 454 |  |  | DoubleGenericData* refVectorY = new DoubleGenericData(); | 
| 455 |  |  | refVectorY->setID("refVectorY"); | 
| 456 |  |  | refVectorY->setData(uY); | 
| 457 |  |  | sd->addProperty(refVectorY); | 
| 458 |  |  |  | 
| 459 |  |  | DoubleGenericData* refVectorZ = new DoubleGenericData(); | 
| 460 |  |  | refVectorZ->setID("refVectorZ"); | 
| 461 |  |  | refVectorZ->setData(uZ); | 
| 462 |  |  | sd->addProperty(refVectorZ); | 
| 463 |  |  | } | 
| 464 |  |  |  | 
| 465 |  |  | // we don't need the angular velocities, so let's exit the line | 
| 466 |  |  | return NULL; | 
| 467 |  |  | } | 
| 468 |  |  |  | 
| 469 |  |  | void RestReader::readZangle(){ | 
| 470 |  |  |  | 
| 471 |  |  | int i; | 
| 472 |  |  | int isPresent; | 
| 473 |  |  |  | 
| 474 |  |  | Molecule* mol; | 
| 475 |  |  | StuntDouble* integrableObject; | 
| 476 |  |  | SimInfo::MoleculeIterator mi; | 
| 477 |  |  | Molecule::IntegrableObjectIterator ii; | 
| 478 |  |  |  | 
| 479 |  |  | #ifdef IS_MPI | 
| 480 | chrisfen | 2903 | int which_node; | 
| 481 | chrisfen | 2101 | MPI_Status istatus; | 
| 482 |  |  | #endif //is_mpi | 
| 483 |  |  |  | 
| 484 |  |  | const int BUFFERSIZE = 2000; // size of the read buffer | 
| 485 | chrisfen | 2903 | unsigned int nTotObjs; // the number of atoms | 
| 486 | chrisfen | 2101 | char read_buffer[BUFFERSIZE]; //the line buffer for reading | 
| 487 |  |  |  | 
| 488 |  |  | std::vector<StuntDouble*> vecParticles; | 
| 489 | tim | 2759 | std::vector<RealType> tempZangs; | 
| 490 | chrisfen | 2101 |  | 
| 491 | chrisfen | 2903 | angFile = info_->getRestFileName(); | 
| 492 | chrisfen | 2101 |  | 
| 493 | chrisfen | 2903 | angFile += "0"; | 
| 494 | chrisfen | 2101 |  | 
| 495 |  |  | // open the omega value file for reading | 
| 496 |  |  | #ifdef IS_MPI | 
| 497 |  |  | if (worldRank == 0) { | 
| 498 |  |  | #endif | 
| 499 |  |  | isPresent = 1; | 
| 500 | chrisfen | 2903 |  | 
| 501 |  |  | inAngFile = new std::ifstream(angFile.c_str()); | 
| 502 |  |  |  | 
| 503 |  |  | if(inAngFile->fail()){ | 
| 504 | chrisfen | 2101 | sprintf(painCave.errMsg, | 
| 505 |  |  | "Restraints Warning: %s file is not present\n" | 
| 506 |  |  | "\tAll omega values will be initialized to zero. If the\n" | 
| 507 |  |  | "\tsimulation is starting from the idealCrystal.in reference\n" | 
| 508 |  |  | "\tconfiguration, this is the desired action. If this is not\n" | 
| 509 |  |  | "\tthe case, the energy calculations will be incorrect.\n", | 
| 510 | chrisfen | 2903 | angFile.c_str()); | 
| 511 | chrisfen | 2101 | painCave.severity = OOPSE_WARNING; | 
| 512 |  |  | painCave.isFatal = 0; | 
| 513 |  |  | simError(); | 
| 514 |  |  | isPresent = 0; | 
| 515 |  |  | } | 
| 516 |  |  |  | 
| 517 |  |  | #ifdef IS_MPI | 
| 518 |  |  | // let the other nodes know the status of the file search | 
| 519 |  |  | MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD); | 
| 520 |  |  | #endif // is_mpi | 
| 521 |  |  |  | 
| 522 |  |  | if (!isPresent) { | 
| 523 |  |  | zeroZangle(); | 
| 524 |  |  | return; | 
| 525 |  |  | } | 
| 526 |  |  |  | 
| 527 |  |  | #ifdef IS_MPI | 
| 528 | chrisfen | 2903 | if (!isPresent) { | 
| 529 |  |  | // master node zeroes out its zAngles if .zang0 isn't present | 
| 530 |  |  | zeroZangle(); | 
| 531 |  |  | return; | 
| 532 |  |  | } | 
| 533 |  |  |  | 
| 534 | chrisfen | 2101 | } | 
| 535 |  |  |  | 
| 536 |  |  | // listen to node 0 to see if we should exit this function | 
| 537 |  |  | if (worldRank != 0) { | 
| 538 |  |  | MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD); | 
| 539 |  |  | if (!isPresent) { | 
| 540 |  |  | zeroZangle(); | 
| 541 |  |  | return; | 
| 542 |  |  | } | 
| 543 |  |  | } | 
| 544 |  |  |  | 
| 545 |  |  | strcpy( checkPointMsg, "zAngle file opened successfully for reading." ); | 
| 546 |  |  | MPIcheckPoint(); | 
| 547 |  |  | #endif | 
| 548 |  |  |  | 
| 549 |  |  | #ifndef IS_MPI | 
| 550 |  |  |  | 
| 551 | chrisfen | 2903 | // read the first line and die if there is a failure | 
| 552 |  |  | inAngFile->getline(read_buffer, sizeof(read_buffer)); | 
| 553 |  |  |  | 
| 554 |  |  | if( inAngFile->eof() ){ | 
| 555 | chrisfen | 2101 | sprintf( painCave.errMsg, | 
| 556 |  |  | "RestraintReader error: error reading 1st line of \"%s\"\n", | 
| 557 | chrisfen | 2903 | angFile.c_str() ); | 
| 558 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 559 |  |  | simError(); | 
| 560 |  |  | } | 
| 561 | chrisfen | 2903 |  | 
| 562 |  |  | // read the file and load the values into a vector | 
| 563 |  |  | inAngFile->getline(read_buffer, sizeof(read_buffer)); | 
| 564 | chrisfen | 2101 |  | 
| 565 | chrisfen | 2903 | while ( !inAngFile->eof() ) { | 
| 566 | chrisfen | 2101 | // check for and ignore blank lines | 
| 567 |  |  | if ( read_buffer != NULL ) | 
| 568 |  |  | tempZangs.push_back( atof(read_buffer) ); | 
| 569 | chrisfen | 2903 |  | 
| 570 |  |  | inAngFile->getline(read_buffer, sizeof(read_buffer)); | 
| 571 | chrisfen | 2101 | } | 
| 572 |  |  |  | 
| 573 |  |  | nTotObjs = info_->getNGlobalIntegrableObjects(); | 
| 574 |  |  |  | 
| 575 |  |  | if( nTotObjs != tempZangs.size() ){ | 
| 576 |  |  | sprintf( painCave.errMsg, | 
| 577 |  |  | "RestraintReader zAngle reading error. %s nIntegrable, %d, " | 
| 578 | chrisfen | 2903 | "does not match the meta-data file's nIntegrable, %i.\n", | 
| 579 |  |  | angFile.c_str(), | 
| 580 |  |  | tempZangs.size(), | 
| 581 |  |  | nTotObjs ); | 
| 582 | chrisfen | 2101 | painCave.isFatal = 1; | 
| 583 |  |  | simError(); | 
| 584 |  |  | } | 
| 585 |  |  |  | 
| 586 |  |  | // load the zAngles into the integrable objects | 
| 587 |  |  | i = 0; | 
| 588 |  |  | for (mol = info_->beginMolecule(mi); mol != NULL; | 
| 589 |  |  | mol = info_->nextMolecule(mi)) { | 
| 590 |  |  |  | 
| 591 |  |  | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 592 |  |  | integrableObject != NULL; | 
| 593 |  |  | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 594 |  |  |  | 
| 595 |  |  | integrableObject->setZangle(tempZangs[i]); | 
| 596 |  |  | i++; | 
| 597 |  |  | } | 
| 598 |  |  | } | 
| 599 |  |  |  | 
| 600 |  |  | // MPI Section of code.......... | 
| 601 |  |  | #else //IS_MPI | 
| 602 |  |  |  | 
| 603 |  |  | // first thing first, suspend fatalities. | 
| 604 |  |  | painCave.isEventLoop = 1; | 
| 605 | chrisfen | 2107 |  | 
| 606 |  |  | int masterNode = 0; | 
| 607 | chrisfen | 2903 |  | 
| 608 | chrisfen | 2107 | int haveError; | 
| 609 | chrisfen | 2868 | int intObjIndex; | 
| 610 |  |  | int intObjIndexTransfer; | 
| 611 | chrisfen | 2107 |  | 
| 612 | chrisfen | 2903 | int j; | 
| 613 | chrisfen | 2101 | int nCurObj; | 
| 614 | chrisfen | 2903 | RealType angleTransfer; | 
| 615 | chrisfen | 2101 |  | 
| 616 | chrisfen | 2107 | nTotObjs = info_->getNGlobalIntegrableObjects(); | 
| 617 | chrisfen | 2101 | haveError = 0; | 
| 618 | chrisfen | 2107 |  | 
| 619 |  |  | if (worldRank == masterNode) { | 
| 620 | chrisfen | 2101 |  | 
| 621 | chrisfen | 2903 | inAngFile->getline(read_buffer, sizeof(read_buffer)); | 
| 622 |  |  |  | 
| 623 |  |  | if( inAngFile->eof() ){ | 
| 624 | chrisfen | 2101 | sprintf( painCave.errMsg, | 
| 625 | chrisfen | 2903 | "Error reading 1st line of %s \n ",angFile.c_str()); | 
| 626 | chrisfen | 2101 | haveError = 1; | 
| 627 |  |  | simError(); | 
| 628 |  |  | } | 
| 629 |  |  |  | 
| 630 | chrisfen | 2903 | // let the master node read the file and load the temporary angle vector | 
| 631 |  |  | inAngFile->getline(read_buffer, sizeof(read_buffer)); | 
| 632 |  |  |  | 
| 633 |  |  | while ( !inAngFile->eof() ) { | 
| 634 | chrisfen | 2101 | // check for and ignore blank lines | 
| 635 |  |  | if ( read_buffer != NULL ) | 
| 636 |  |  | tempZangs.push_back( atof(read_buffer) ); | 
| 637 | chrisfen | 2903 |  | 
| 638 |  |  | inAngFile->getline(read_buffer, sizeof(read_buffer)); | 
| 639 | chrisfen | 2101 | } | 
| 640 | chrisfen | 2868 |  | 
| 641 | chrisfen | 2101 | // Check to see that the number of integrable objects in the | 
| 642 |  |  | // intial configuration file is the same as derived from the | 
| 643 |  |  | // meta-data file. | 
| 644 |  |  | if( nTotObjs != tempZangs.size() ){ | 
| 645 |  |  | sprintf( painCave.errMsg, | 
| 646 |  |  | "RestraintReader zAngle reading Error. %s nIntegrable, %d, " | 
| 647 |  |  | "does not match the meta-data file's nIntegrable, %d.\n", | 
| 648 | chrisfen | 2903 | angFile.c_str(), | 
| 649 |  |  | tempZangs.size(), | 
| 650 |  |  | nTotObjs); | 
| 651 | chrisfen | 2101 | haveError= 1; | 
| 652 |  |  | simError(); | 
| 653 |  |  | } | 
| 654 |  |  |  | 
| 655 | chrisfen | 2107 | // At this point, node 0 has a tempZangs vector completed, and | 
| 656 |  |  | // everyone else has nada | 
| 657 | chrisfen | 2101 |  | 
| 658 | chrisfen | 2107 | for (i=0 ; i < info_->getNGlobalMolecules(); i++) { | 
| 659 |  |  | // Get the Node number which has this atom | 
| 660 |  |  | which_node = info_->getMolToProc(i); | 
| 661 |  |  |  | 
| 662 | chrisfen | 2814 | if (which_node == masterNode) { | 
| 663 | chrisfen | 2107 | mol = info_->getMoleculeByGlobalIndex(i); | 
| 664 | chrisfen | 2814 |  | 
| 665 | chrisfen | 2107 | if(mol == NULL) { | 
| 666 |  |  | strcpy(painCave.errMsg, "Molecule not found on node 0!"); | 
| 667 |  |  | haveError = 1; | 
| 668 |  |  | simError(); | 
| 669 |  |  | } | 
| 670 | chrisfen | 2814 |  | 
| 671 | chrisfen | 2107 | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 672 |  |  | integrableObject != NULL; | 
| 673 |  |  | integrableObject = mol->nextIntegrableObject(ii)){ | 
| 674 | chrisfen | 2992 | intObjIndex = integrableObject->getGlobalIntegrableObjectIndex(); | 
| 675 | chrisfen | 2868 | integrableObject->setZangle(tempZangs[intObjIndex]); | 
| 676 | chrisfen | 2107 | } | 
| 677 |  |  |  | 
| 678 |  |  | } else { | 
| 679 |  |  | // I am MASTER OF THE UNIVERSE, but I don't own this molecule | 
| 680 | chrisfen | 2868 | // listen for the number of integrableObjects in the molecule | 
| 681 | chrisfen | 2107 | MPI_Recv(&nCurObj, 1, MPI_INT, which_node, | 
| 682 |  |  | TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus); | 
| 683 |  |  |  | 
| 684 | chrisfen | 2868 | for(j=0; j < nCurObj; j++){ | 
| 685 |  |  | // listen for which integrableObject we need to send the value for | 
| 686 |  |  | MPI_Recv(&intObjIndexTransfer, 1, MPI_INT, which_node, | 
| 687 |  |  | TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus); | 
| 688 |  |  | angleTransfer = tempZangs[intObjIndexTransfer]; | 
| 689 |  |  | // send the value to the node so it can initialize the object | 
| 690 | tim | 2759 | MPI_Send(&angleTransfer, 1, MPI_REALTYPE, which_node, | 
| 691 | chrisfen | 2107 | TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD); | 
| 692 |  |  | } | 
| 693 |  |  | } | 
| 694 |  |  | } | 
| 695 |  |  | } else { | 
| 696 |  |  | // I am SLAVE TO THE MASTER | 
| 697 |  |  | for (i=0 ; i < info_->getNGlobalMolecules(); i++) { | 
| 698 | chrisfen | 2903 | which_node = info_->getMolToProc(i); | 
| 699 | chrisfen | 2107 |  | 
| 700 |  |  | if (which_node == worldRank) { | 
| 701 |  |  |  | 
| 702 |  |  | // BUT I OWN THIS MOLECULE!!! | 
| 703 |  |  |  | 
| 704 |  |  | mol = info_->getMoleculeByGlobalIndex(i); | 
| 705 |  |  |  | 
| 706 |  |  | if(mol == NULL) { | 
| 707 |  |  | sprintf(painCave.errMsg, | 
| 708 |  |  | "RestReader Error: molecule not found on node %d\n", | 
| 709 |  |  | worldRank); | 
| 710 |  |  | painCave.isFatal = 1; | 
| 711 | chrisfen | 2101 | simError(); | 
| 712 |  |  | } | 
| 713 | chrisfen | 2107 |  | 
| 714 |  |  | nCurObj = mol->getNIntegrableObjects(); | 
| 715 | chrisfen | 2868 | // send the number of integrableObjects in the molecule | 
| 716 | chrisfen | 2107 | MPI_Send(&nCurObj, 1, MPI_INT, 0, | 
| 717 |  |  | TAKE_THIS_TAG_INT, MPI_COMM_WORLD); | 
| 718 |  |  |  | 
| 719 |  |  | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 720 |  |  | integrableObject != NULL; | 
| 721 |  |  | integrableObject = mol->nextIntegrableObject(ii)){ | 
| 722 | chrisfen | 2992 | intObjIndexTransfer = integrableObject->getGlobalIntegrableObjectIndex(); | 
| 723 | chrisfen | 2868 | // send the global index of the integrableObject | 
| 724 |  |  | MPI_Send(&intObjIndexTransfer, 1, MPI_INT, 0, | 
| 725 |  |  | TAKE_THIS_TAG_INT, MPI_COMM_WORLD); | 
| 726 |  |  | // listen for the value we want to set locally | 
| 727 | tim | 2759 | MPI_Recv(&angleTransfer, 1, MPI_REALTYPE, 0, | 
| 728 | chrisfen | 2107 | TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus); | 
| 729 |  |  |  | 
| 730 |  |  | integrableObject->setZangle(angleTransfer); | 
| 731 |  |  | } | 
| 732 |  |  | } | 
| 733 | chrisfen | 2101 | } | 
| 734 | chrisfen | 2107 | } | 
| 735 | chrisfen | 2101 | #endif | 
| 736 |  |  | } | 
| 737 |  |  |  | 
| 738 |  |  | void RestReader :: zeroZangle(){ | 
| 739 |  |  |  | 
| 740 |  |  | Molecule* mol; | 
| 741 |  |  | StuntDouble* integrableObject; | 
| 742 |  |  | SimInfo::MoleculeIterator mi; | 
| 743 |  |  | Molecule::IntegrableObjectIterator ii; | 
| 744 | chrisfen | 2992 |  | 
| 745 | chrisfen | 2101 | #ifndef IS_MPI | 
| 746 |  |  | // set all zAngles to 0.0 | 
| 747 |  |  | for (mol = info_->beginMolecule(mi); mol != NULL; | 
| 748 |  |  | mol = info_->nextMolecule(mi)) | 
| 749 |  |  |  | 
| 750 |  |  | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 751 |  |  | integrableObject != NULL; | 
| 752 |  |  | integrableObject = mol->nextIntegrableObject(ii)) | 
| 753 |  |  | integrableObject->setZangle( 0.0 ); | 
| 754 |  |  |  | 
| 755 |  |  |  | 
| 756 |  |  | // MPI Section of code.......... | 
| 757 |  |  | #else //IS_MPI | 
| 758 |  |  |  | 
| 759 |  |  | // first thing first, suspend fatalities. | 
| 760 |  |  | painCave.isEventLoop = 1; | 
| 761 |  |  |  | 
| 762 |  |  | int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone | 
| 763 | chrisfen | 2107 | int haveError; | 
| 764 | chrisfen | 2101 | int which_node; | 
| 765 |  |  |  | 
| 766 |  |  | haveError = 0; | 
| 767 | chrisfen | 2107 |  | 
| 768 | chrisfen | 2903 | for (int i=0 ; i < info_->getNGlobalMolecules(); i++) { | 
| 769 |  |  |  | 
| 770 |  |  | // Get the Node number which has this atom | 
| 771 |  |  | which_node = info_->getMolToProc(i); | 
| 772 | chrisfen | 2107 |  | 
| 773 | chrisfen | 2903 | // each processor zeroes its own integrable objects | 
| 774 |  |  | if (which_node == worldRank) { | 
| 775 |  |  | mol = info_->getMoleculeByGlobalIndex(i); | 
| 776 |  |  |  | 
| 777 |  |  | if(mol == NULL) { | 
| 778 |  |  | sprintf( painCave.errMsg, | 
| 779 |  |  | "Molecule not found on node %i!", | 
| 780 |  |  | which_node ); | 
| 781 |  |  | haveError = 1; | 
| 782 |  |  | simError(); | 
| 783 | chrisfen | 2107 | } | 
| 784 |  |  |  | 
| 785 | chrisfen | 2903 | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 786 |  |  | integrableObject != NULL; | 
| 787 |  |  | integrableObject = mol->nextIntegrableObject(ii)){ | 
| 788 | chrisfen | 2107 |  | 
| 789 | chrisfen | 2903 | integrableObject->setZangle( 0.0 ); | 
| 790 | chrisfen | 2107 |  | 
| 791 | chrisfen | 2903 | } | 
| 792 | chrisfen | 2101 | } | 
| 793 |  |  | } | 
| 794 | chrisfen | 2903 |  | 
| 795 | chrisfen | 2101 | #endif | 
| 796 |  |  | } | 
| 797 |  |  |  | 
| 798 |  |  | } // end namespace oopse |