# | Line 1 | Line 1 | |
---|---|---|
1 | < | /* |
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 | < | */ |
1 | > | /* |
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 | > | #include "io/DumpReader.hpp" |
43 | > | #include "io/RestReader.hpp" |
44 | > | #include "primitives/Molecule.hpp" |
45 | > | #include "restraints/ObjectRestraint.hpp" |
46 | > | #include "restraints/MolecularRestraint.hpp" |
47 | ||
48 | < | #define _LARGEFILE_SOURCE64 |
49 | < | #define _FILE_OFFSET_BITS 64 |
48 | > | namespace oopse { |
49 | > | |
50 | > | void RestReader::readReferenceStructure() { |
51 | ||
52 | < | #include <sys/types.h> |
46 | < | #include <sys/stat.h> |
52 | > | // some of this comes directly from DumpReader. |
53 | ||
54 | < | #include <iostream> |
55 | < | #include <math.h> |
54 | > | if (!isScanned_) |
55 | > | scanFile(); |
56 | > | |
57 | > | int storageLayout = info_->getSnapshotManager()->getStorageLayout(); |
58 | > | |
59 | > | if (storageLayout & DataStorage::dslPosition) { |
60 | > | needPos_ = true; |
61 | > | } else { |
62 | > | needPos_ = false; |
63 | > | } |
64 | > | |
65 | > | needVel_ = false; |
66 | > | |
67 | > | if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) { |
68 | > | needQuaternion_ = true; |
69 | > | } else { |
70 | > | needQuaternion_ = false; |
71 | > | } |
72 | ||
73 | < | #include <stdio.h> |
52 | < | #include <stdlib.h> |
53 | < | #include <string.h> |
73 | > | needAngMom_ = false; |
74 | ||
75 | < | #include "primitives/Molecule.hpp" |
76 | < | #include "utils/MemoryUtils.hpp" |
77 | < | #include "utils/StringTokenizer.hpp" |
58 | < | #include "io/RestReader.hpp" |
59 | < | #include "utils/simError.h" |
75 | > | // We need temporary storage to keep track of all StuntDouble positions |
76 | > | // in case some of the restraints are molecular (i.e. if they use |
77 | > | // multiple SD positions to determine restrained orientations or positions: |
78 | ||
79 | < | #ifdef IS_MPI |
80 | < | #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 |
79 | > | all_pos_.clear(); |
80 | > | all_pos_.resize(info_->getNGlobalIntegrableObjects() ); |
81 | ||
82 | < | namespace oopse { |
83 | < | |
84 | < | RestReader::RestReader( SimInfo* info ) : info_(info){ |
82 | > | |
83 | > | // Restraint files are just standard dump files, but with the reference |
84 | > | // structure stored in the first frame (frame 0). |
85 | > | // RestReader overloads readSet and explicitly handles all of the |
86 | > | // ObjectRestraints in that method: |
87 | > | |
88 | > | readSet(0); |
89 | ||
90 | < | idealName = "idealCrystal.in"; |
91 | < | |
74 | < | #ifdef IS_MPI |
75 | < | if (worldRank == 0) { |
76 | < | #endif |
90 | > | // all ObjectRestraints have been handled, now we have to worry about |
91 | > | // molecular restraints: |
92 | ||
93 | < | inIdealFile = new std::ifstream(idealName.c_str()); |
93 | > | SimInfo::MoleculeIterator i; |
94 | > | Molecule::IntegrableObjectIterator j; |
95 | > | Molecule* mol; |
96 | > | StuntDouble* sd; |
97 | ||
98 | < | if(inIdealFile->fail()){ |
99 | < | sprintf(painCave.errMsg, |
100 | < | "RestReader: Cannot open file: %s\n", |
101 | < | idealName.c_str()); |
102 | < | painCave.isFatal = 1; |
103 | < | simError(); |
104 | < | } |
98 | > | // no need to worry about parallel molecules, as molecules are not |
99 | > | // split across processor boundaries. Just loop over all molecules |
100 | > | // we know about: |
101 | > | |
102 | > | for (mol = info_->beginMolecule(i); mol != NULL; |
103 | > | mol = info_->nextMolecule(i)) { |
104 | > | |
105 | > | // is this molecule restrained? |
106 | > | GenericData* data = mol->getPropertyByName("Restraint"); |
107 | ||
108 | < | #ifdef IS_MPI |
89 | < | } |
90 | < | strcpy( checkPointMsg, |
91 | < | "File \"idealCrystal.in\" opened successfully for reading." ); |
92 | < | MPIcheckPoint(); |
93 | < | #endif |
108 | > | if (data != NULL) { |
109 | ||
110 | < | return; |
96 | < | } |
97 | < | |
98 | < | RestReader :: ~RestReader( ){ |
99 | < | #ifdef IS_MPI |
100 | < | if (worldRank == 0) { |
101 | < | #endif |
110 | > | // make sure we can reinterpret the generic data as restraint data: |
111 | ||
112 | < | delete inIdealFile; |
104 | < | delete inAngFile; |
112 | > | RestraintData* restData= dynamic_cast<RestraintData*>(data); |
113 | ||
114 | < | #ifdef IS_MPI |
107 | < | } |
108 | < | strcpy( checkPointMsg, |
109 | < | "File idealCrystal.in (and .zang0 if present) closed successfully." ); |
110 | < | MPIcheckPoint(); |
111 | < | #endif |
112 | < | |
113 | < | return; |
114 | < | } |
115 | < | |
116 | < | |
117 | < | void RestReader :: readIdealCrystal(){ |
118 | < | |
119 | < | #ifdef IS_MPI |
120 | < | int which_node; |
121 | < | int i, j; |
122 | < | #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 | < | inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
114 | > | if (restData != NULL) { |
115 | ||
116 | < | if( inIdealFile->eof() ){ |
117 | < | sprintf( painCave.errMsg, |
143 | < | "RestraintReader error: error reading 1st line of \"%s\"\n", |
144 | < | idealName.c_str() ); |
145 | < | 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 | < | idealName.c_str(), |
156 | < | nTotObjs, |
157 | < | info_->getNGlobalIntegrableObjects()); |
158 | < | painCave.isFatal = 1; |
159 | < | simError(); |
160 | < | } |
161 | < | |
162 | < | // skip over the comment line |
163 | < | inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
116 | > | // make sure we can reinterpet the restraint data as a |
117 | > | // pointer to a MolecularRestraint: |
118 | ||
119 | < | if( inIdealFile->eof() ){ |
166 | < | sprintf( painCave.errMsg, |
167 | < | "error in reading commment in %s\n", |
168 | < | idealName.c_str()); |
169 | < | 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)) { |
119 | > | MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData()); |
120 | ||
121 | < | inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
121 | > | if (mRest != NULL) { |
122 | ||
123 | < | if( inIdealFile->eof() ){ |
124 | < | 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 | < | idealName.c_str(), nTotObjs, |
194 | < | 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 | < | int haveError; |
123 | > | // now we need to pack the stunt doubles for the reference |
124 | > | // structure: |
125 | ||
126 | < | nTotObjs = info_->getNGlobalIntegrableObjects(); |
127 | < | haveError = 0; |
128 | < | |
129 | < | if (worldRank == masterNode) { |
225 | < | inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
226 | < | |
227 | < | if( inIdealFile->eof() ){ |
228 | < | sprintf( painCave.errMsg, |
229 | < | "Error reading 1st line of %s \n ",idealName.c_str()); |
230 | < | 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 | < | idealName.c_str(), nTotObjs, |
244 | < | info_->getNGlobalIntegrableObjects()); |
245 | < | painCave.isFatal = 1; |
246 | < | simError(); |
247 | < | } |
248 | < | |
249 | < | // skip over the comment line |
250 | < | inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
251 | < | |
252 | < | if( inIdealFile->eof() ){ |
253 | < | sprintf( painCave.errMsg, |
254 | < | "error in reading commment in %s\n", idealName.c_str()); |
255 | < | painCave.isFatal = 1; |
256 | < | simError(); |
257 | < | } |
258 | < | |
259 | < | MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
260 | < | |
261 | < | for (i=0 ; i < info_->getNGlobalMolecules(); i++) { |
262 | < | which_node = info_->getMolToProc(i); |
263 | < | |
264 | < | if(which_node == masterNode){ |
265 | < | //molecules belong to master node |
266 | < | |
267 | < | mol = info_->getMoleculeByGlobalIndex(i); |
268 | < | |
269 | < | if(mol == NULL) { |
270 | < | sprintf(painCave.errMsg, |
271 | < | "RestReader Error: Molecule not found on node %d!\n", |
272 | < | worldRank); |
273 | < | painCave.isFatal = 1; |
274 | < | simError(); |
275 | < | } |
276 | < | |
277 | < | for (integrableObject = mol->beginIntegrableObject(ii); |
278 | < | integrableObject != NULL; |
279 | < | integrableObject = mol->nextIntegrableObject(ii)){ |
126 | > | std::vector<Vector3d> ref; |
127 | > | int count = 0; |
128 | > | RealType mass, mTot; |
129 | > | Vector3d COM(0.0); |
130 | ||
131 | < | inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
131 | > | mTot = 0.0; |
132 | ||
133 | < | if( inIdealFile->eof() ){ |
134 | < | sprintf(painCave.errMsg, |
135 | < | "RestReader Error: error in reading file %s\n" |
136 | < | "natoms = %d; index = %d\n" |
137 | < | "error reading the line from the file.\n", |
138 | < | idealName.c_str(), nTotObjs, i ); |
139 | < | painCave.isFatal = 1; |
140 | < | simError(); |
133 | > | // loop over the stunt doubles in this molecule in the order we |
134 | > | // will be looping them in the restraint code: |
135 | > | |
136 | > | for (sd = mol->beginIntegrableObject(j); sd != NULL; |
137 | > | sd = mol->nextIntegrableObject(j)) { |
138 | > | |
139 | > | // push back the reference positions of the stunt |
140 | > | // doubles from the *globally* sorted array of |
141 | > | // positions: |
142 | > | |
143 | > | ref.push_back( all_pos_[sd->getGlobalIndex()] ); |
144 | > | |
145 | > | mass = sd->getMass(); |
146 | > | |
147 | > | COM = COM + mass * ref[count]; |
148 | > | mTot = mTot + mass; |
149 | > | count = count + 1; |
150 | } | |
292 | – | |
293 | – | parseIdealLine(read_buffer, integrableObject); |
294 | – | |
295 | – | } |
296 | – | |
297 | – | } 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 | – | for(j = 0; j < nCurObj; j++){ |
304 | – | inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
305 | – | |
306 | – | if( inIdealFile->eof() ){ |
307 | – | 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 | – | idealName.c_str(), nTotObjs, i ); |
312 | – | painCave.isFatal = 1; |
313 | – | simError(); |
314 | – | } |
151 | ||
152 | < | MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node, |
317 | < | TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD); |
318 | < | } |
319 | < | } |
320 | < | } |
321 | < | } else { |
322 | < | //actions taken at slave nodes |
323 | < | MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
152 | > | COM /= mTot; |
153 | ||
154 | < | for (i=0 ; i < info_->getNGlobalMolecules(); i++) { |
326 | < | int which_node = info_->getMolToProc(i); |
154 | > | mRest->setReferenceStructure(ref, COM); |
155 | ||
328 | – | if(which_node == worldRank){ |
329 | – | //molecule with global index i belongs to this processor |
330 | – | |
331 | – | mol = info_->getMoleculeByGlobalIndex(i); |
332 | – | |
333 | – | if(mol == NULL) { |
334 | – | sprintf(painCave.errMsg, |
335 | – | "RestReader Error: molecule not found on node %d\n", |
336 | – | worldRank); |
337 | – | painCave.isFatal = 1; |
338 | – | simError(); |
156 | } | |
340 | – | |
341 | – | nCurObj = mol->getNIntegrableObjects(); |
342 | – | |
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 | – | } |
157 | } | |
158 | } | |
159 | } | |
363 | – | #endif |
160 | } | |
365 | – | |
366 | – | char* RestReader::parseIdealLine(char* readLine, StuntDouble* sd){ |
367 | – | |
368 | – | 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 | – | StringTokenizer tokenizer(readLine); |
375 | – | int nTokens; |
376 | – | |
377 | – | nTokens = tokenizer.countTokens(); |
161 | ||
162 | < | if (nTokens < 14) { |
163 | < | sprintf(painCave.errMsg, |
164 | < | "RestReader Error: Not enough Tokens.\n"); |
165 | < | painCave.isFatal = 1; |
166 | < | simError(); |
167 | < | } |
168 | < | |
169 | < | std::string name = tokenizer.nextToken(); |
162 | > | |
163 | > | |
164 | > | void RestReader::parseDumpLine(const std::string& line) { |
165 | > | StringTokenizer tokenizer(line); |
166 | > | int nTokens; |
167 | > | |
168 | > | nTokens = tokenizer.countTokens(); |
169 | > | |
170 | > | if (nTokens < 2) { |
171 | > | sprintf(painCave.errMsg, |
172 | > | "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); |
173 | > | painCave.isFatal = 1; |
174 | > | simError(); |
175 | > | } |
176 | ||
177 | < | if (name != sd->getType()) { |
178 | < | |
179 | < | sprintf(painCave.errMsg, |
391 | < | "RestReader Error: Atom type [%s] in %s does not " |
392 | < | "match Atom Type [%s] in .md file.\n", |
393 | < | name.c_str(), idealName.c_str(), |
394 | < | 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(); |
177 | > | int index = tokenizer.nextTokenAsInt(); |
178 | > | |
179 | > | StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index); |
180 | ||
181 | < | // store the positions in the stuntdouble as generic data doubles |
182 | < | DoubleGenericData* refPosX = new DoubleGenericData(); |
183 | < | refPosX->setID("refPosX"); |
406 | < | refPosX->setData(pos[0]); |
407 | < | sd->addProperty(refPosX); |
181 | > | if (integrableObject == NULL) { |
182 | > | return; |
183 | > | } |
184 | ||
185 | < | DoubleGenericData* refPosY = new DoubleGenericData(); |
186 | < | refPosY->setID("refPosY"); |
187 | < | refPosY->setData(pos[1]); |
188 | < | sd->addProperty(refPosY); |
189 | < | |
190 | < | DoubleGenericData* refPosZ = new DoubleGenericData(); |
191 | < | refPosZ->setID("refPosZ"); |
192 | < | refPosZ->setData(pos[2]); |
193 | < | sd->addProperty(refPosZ); |
185 | > | std::string type = tokenizer.nextToken(); |
186 | > | int size = type.size(); |
187 | > | Vector3d pos; |
188 | > | Vector3d vel; |
189 | > | Quat4d q; |
190 | > | Vector3d ji; |
191 | > | Vector3d force; |
192 | > | Vector3d torque; |
193 | > | |
194 | > | for(int i = 0; i < size; ++i) { |
195 | > | switch(type[i]) { |
196 | > | |
197 | > | case 'p': { |
198 | > | pos[0] = tokenizer.nextTokenAsDouble(); |
199 | > | pos[1] = tokenizer.nextTokenAsDouble(); |
200 | > | pos[2] = tokenizer.nextTokenAsDouble(); |
201 | > | break; |
202 | > | } |
203 | > | case 'v' : { |
204 | > | vel[0] = tokenizer.nextTokenAsDouble(); |
205 | > | vel[1] = tokenizer.nextTokenAsDouble(); |
206 | > | vel[2] = tokenizer.nextTokenAsDouble(); |
207 | > | break; |
208 | > | } |
209 | ||
210 | < | // we don't need the velocities |
211 | < | uselessToken = tokenizer.nextTokenAsDouble(); |
212 | < | uselessToken = tokenizer.nextTokenAsDouble(); |
213 | < | uselessToken = tokenizer.nextTokenAsDouble(); |
214 | < | |
215 | < | if (sd->isDirectional()) { |
216 | < | |
217 | < | q[0] = tokenizer.nextTokenAsDouble(); |
218 | < | q[1] = tokenizer.nextTokenAsDouble(); |
219 | < | q[2] = tokenizer.nextTokenAsDouble(); |
220 | < | q[3] = tokenizer.nextTokenAsDouble(); |
221 | < | |
222 | < | // now build the rotation matrix and find the unit vectors |
223 | < | RfromQ[0][0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3]; |
224 | < | RfromQ[0][1] = 2*(q[1]*q[2] + q[0]*q[3]); |
225 | < | RfromQ[0][2] = 2*(q[1]*q[3] - q[0]*q[2]); |
226 | < | RfromQ[1][0] = 2*(q[1]*q[2] - q[0]*q[3]); |
227 | < | RfromQ[1][1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3]; |
228 | < | RfromQ[1][2] = 2*(q[2]*q[3] + q[0]*q[1]); |
229 | < | RfromQ[2][0] = 2*(q[1]*q[3] + q[0]*q[2]); |
230 | < | RfromQ[2][1] = 2*(q[2]*q[3] - q[0]*q[1]); |
231 | < | RfromQ[2][2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3]; |
232 | < | |
233 | < | normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1] |
234 | < | + RfromQ[2][2]*RfromQ[2][2]); |
235 | < | uX = RfromQ[2][0]/normalize; |
236 | < | uY = RfromQ[2][1]/normalize; |
237 | < | uZ = RfromQ[2][2]/normalize; |
238 | < | |
239 | < | // store reference unit vectors as generic data in the stuntdouble |
240 | < | DoubleGenericData* refVectorX = new DoubleGenericData(); |
241 | < | refVectorX->setID("refVectorX"); |
242 | < | refVectorX->setData(uX); |
243 | < | 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 | < | int which_node; |
481 | < | MPI_Status istatus; |
482 | < | #endif //is_mpi |
483 | < | |
484 | < | const int BUFFERSIZE = 2000; // size of the read buffer |
485 | < | unsigned int nTotObjs; // the number of atoms |
486 | < | char read_buffer[BUFFERSIZE]; //the line buffer for reading |
487 | < | |
488 | < | std::vector<StuntDouble*> vecParticles; |
489 | < | std::vector<RealType> tempZangs; |
490 | < | |
491 | < | angFile = info_->getRestFileName(); |
492 | < | |
493 | < | angFile += "0"; |
494 | < | |
495 | < | // open the omega value file for reading |
496 | < | #ifdef IS_MPI |
497 | < | if (worldRank == 0) { |
498 | < | #endif |
499 | < | isPresent = 1; |
500 | < | |
501 | < | inAngFile = new std::ifstream(angFile.c_str()); |
502 | < | |
503 | < | if(inAngFile->fail()){ |
504 | < | 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 | < | angFile.c_str()); |
511 | < | 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 | < | if (!isPresent) { |
529 | < | // master node zeroes out its zAngles if .zang0 isn't present |
530 | < | zeroZangle(); |
531 | < | return; |
532 | < | } |
210 | > | case 'q' : { |
211 | > | if (integrableObject->isDirectional()) { |
212 | > | |
213 | > | q[0] = tokenizer.nextTokenAsDouble(); |
214 | > | q[1] = tokenizer.nextTokenAsDouble(); |
215 | > | q[2] = tokenizer.nextTokenAsDouble(); |
216 | > | q[3] = tokenizer.nextTokenAsDouble(); |
217 | > | |
218 | > | RealType qlen = q.length(); |
219 | > | if (qlen < oopse::epsilon) { //check quaternion is not equal to 0 |
220 | > | |
221 | > | sprintf(painCave.errMsg, |
222 | > | "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n"); |
223 | > | painCave.isFatal = 1; |
224 | > | simError(); |
225 | > | |
226 | > | } |
227 | > | |
228 | > | q.normalize(); |
229 | > | } |
230 | > | break; |
231 | > | } |
232 | > | case 'j' : { |
233 | > | if (integrableObject->isDirectional()) { |
234 | > | ji[0] = tokenizer.nextTokenAsDouble(); |
235 | > | ji[1] = tokenizer.nextTokenAsDouble(); |
236 | > | ji[2] = tokenizer.nextTokenAsDouble(); |
237 | > | } |
238 | > | break; |
239 | > | } |
240 | > | case 'f': { |
241 | > | force[0] = tokenizer.nextTokenAsDouble(); |
242 | > | force[1] = tokenizer.nextTokenAsDouble(); |
243 | > | force[2] = tokenizer.nextTokenAsDouble(); |
244 | ||
245 | < | } |
246 | < | |
247 | < | // listen to node 0 to see if we should exit this function |
248 | < | if (worldRank != 0) { |
249 | < | MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD); |
250 | < | if (!isPresent) { |
251 | < | zeroZangle(); |
252 | < | return; |
245 | > | break; |
246 | > | } |
247 | > | case 't' : { |
248 | > | torque[0] = tokenizer.nextTokenAsDouble(); |
249 | > | torque[1] = tokenizer.nextTokenAsDouble(); |
250 | > | torque[2] = tokenizer.nextTokenAsDouble(); |
251 | > | |
252 | > | break; |
253 | > | } |
254 | > | default: { |
255 | > | sprintf(painCave.errMsg, |
256 | > | "DumpReader Error: %s is an unrecognized type\n", type.c_str()); |
257 | > | painCave.isFatal = 1; |
258 | > | simError(); |
259 | > | break; |
260 | > | } |
261 | > | |
262 | } | |
263 | } | |
264 | + | |
265 | + | // keep the position in case we need it for a molecular restraint: |
266 | ||
267 | < | strcpy( checkPointMsg, "zAngle file opened successfully for reading." ); |
546 | < | MPIcheckPoint(); |
547 | < | #endif |
548 | < | |
549 | < | #ifndef IS_MPI |
550 | < | |
551 | < | // read the first line and die if there is a failure |
552 | < | inAngFile->getline(read_buffer, sizeof(read_buffer)); |
267 | > | all_pos_[index] = pos; |
268 | ||
269 | < | if( inAngFile->eof() ){ |
270 | < | sprintf( painCave.errMsg, |
271 | < | "RestraintReader error: error reading 1st line of \"%s\"\n", |
557 | < | angFile.c_str() ); |
558 | < | painCave.isFatal = 1; |
559 | < | simError(); |
560 | < | } |
269 | > | // is this io restrained? |
270 | > | GenericData* data = integrableObject->getPropertyByName("Restraint"); |
271 | > | ObjectRestraint* oRest; |
272 | ||
273 | < | // read the file and load the values into a vector |
274 | < | inAngFile->getline(read_buffer, sizeof(read_buffer)); |
275 | < | |
276 | < | while ( !inAngFile->eof() ) { |
277 | < | // check for and ignore blank lines |
278 | < | if ( read_buffer != NULL ) |
279 | < | tempZangs.push_back( atof(read_buffer) ); |
280 | < | |
281 | < | inAngFile->getline(read_buffer, sizeof(read_buffer)); |
282 | < | } |
283 | < | |
284 | < | nTotObjs = info_->getNGlobalIntegrableObjects(); |
285 | < | |
286 | < | if( nTotObjs != tempZangs.size() ){ |
576 | < | sprintf( painCave.errMsg, |
577 | < | "RestraintReader zAngle reading error. %s nIntegrable, %d, " |
578 | < | "does not match the meta-data file's nIntegrable, %i.\n", |
579 | < | angFile.c_str(), |
580 | < | tempZangs.size(), |
581 | < | nTotObjs ); |
582 | < | 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++; |
273 | > | if (data != NULL) { |
274 | > | // make sure we can reinterpret the generic data as restraint data: |
275 | > | RestraintData* restData= dynamic_cast<RestraintData*>(data); |
276 | > | if (restData != NULL) { |
277 | > | // make sure we can reinterpet the restraint data as a pointer to |
278 | > | // an ObjectRestraint: |
279 | > | oRest = dynamic_cast<ObjectRestraint*>(restData->getData()); |
280 | > | if (oRest != NULL) { |
281 | > | if (integrableObject->isDirectional()) { |
282 | > | oRest->setReferenceStructure(pos, q.toRotationMatrix3()); |
283 | > | } else { |
284 | > | oRest->setReferenceStructure(pos); |
285 | > | } |
286 | > | } |
287 | } | |
288 | } | |
599 | – | |
600 | – | // MPI Section of code.......... |
601 | – | #else //IS_MPI |
602 | – | |
603 | – | // first thing first, suspend fatalities. |
604 | – | painCave.isEventLoop = 1; |
289 | ||
290 | < | int masterNode = 0; |
290 | > | } |
291 | > | |
292 | ||
608 | – | int haveError; |
609 | – | int intObjIndex; |
610 | – | int intObjIndexTransfer; |
293 | ||
294 | < | int j; |
295 | < | int nCurObj; |
296 | < | RealType angleTransfer; |
615 | < | |
616 | < | nTotObjs = info_->getNGlobalIntegrableObjects(); |
617 | < | haveError = 0; |
294 | > | void RestReader::readFrameProperties(std::istream& inputStream) { |
295 | > | inputStream.getline(buffer, bufferSize); |
296 | > | std::string line(buffer); |
297 | ||
298 | < | if (worldRank == masterNode) { |
299 | < | |
300 | < | inAngFile->getline(read_buffer, sizeof(read_buffer)); |
298 | > | if (line.find("<FrameData>") == std::string::npos) { |
299 | > | sprintf(painCave.errMsg, |
300 | > | "DumpReader Error: Missing <FrameData>\n"); |
301 | > | painCave.isFatal = 1; |
302 | > | simError(); |
303 | > | } |
304 | ||
305 | < | if( inAngFile->eof() ){ |
306 | < | sprintf( painCave.errMsg, |
307 | < | "Error reading 1st line of %s \n ",angFile.c_str()); |
626 | < | haveError = 1; |
627 | < | simError(); |
628 | < | } |
629 | < | |
630 | < | // let the master node read the file and load the temporary angle vector |
631 | < | inAngFile->getline(read_buffer, sizeof(read_buffer)); |
305 | > | // restraints don't care about frame data (unless we need to wrap |
306 | > | // coordinates, but we'll worry about that later), so |
307 | > | // we'll just scan ahead until the end of the frame data: |
308 | ||
309 | < | while ( !inAngFile->eof() ) { |
310 | < | // check for and ignore blank lines |
635 | < | if ( read_buffer != NULL ) |
636 | < | tempZangs.push_back( atof(read_buffer) ); |
637 | < | |
638 | < | inAngFile->getline(read_buffer, sizeof(read_buffer)); |
639 | < | } |
640 | < | |
641 | < | // 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 | < | angFile.c_str(), |
649 | < | tempZangs.size(), |
650 | < | nTotObjs); |
651 | < | haveError= 1; |
652 | < | simError(); |
653 | < | } |
309 | > | while(inputStream.getline(buffer, bufferSize)) { |
310 | > | line = buffer; |
311 | ||
312 | < | // At this point, node 0 has a tempZangs vector completed, and |
313 | < | // everyone else has nada |
657 | < | |
658 | < | for (i=0 ; i < info_->getNGlobalMolecules(); i++) { |
659 | < | // Get the Node number which has this atom |
660 | < | which_node = info_->getMolToProc(i); |
661 | < | |
662 | < | if (which_node == masterNode) { |
663 | < | mol = info_->getMoleculeByGlobalIndex(i); |
664 | < | |
665 | < | if(mol == NULL) { |
666 | < | strcpy(painCave.errMsg, "Molecule not found on node 0!"); |
667 | < | haveError = 1; |
668 | < | simError(); |
669 | < | } |
670 | < | |
671 | < | for (integrableObject = mol->beginIntegrableObject(ii); |
672 | < | integrableObject != NULL; |
673 | < | integrableObject = mol->nextIntegrableObject(ii)){ |
674 | < | intObjIndex = integrableObject->getGlobalIndex(); |
675 | < | integrableObject->setZangle(tempZangs[intObjIndex]); |
676 | < | } |
677 | < | |
678 | < | } else { |
679 | < | // I am MASTER OF THE UNIVERSE, but I don't own this molecule |
680 | < | // listen for the number of integrableObjects in the molecule |
681 | < | MPI_Recv(&nCurObj, 1, MPI_INT, which_node, |
682 | < | TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus); |
683 | < | |
684 | < | 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 | < | MPI_Send(&angleTransfer, 1, MPI_REALTYPE, which_node, |
691 | < | TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD); |
692 | < | } |
693 | < | } |
312 | > | if(line.find("</FrameData>") != std::string::npos) { |
313 | > | break; |
314 | } | |
695 | – | } else { |
696 | – | // I am SLAVE TO THE MASTER |
697 | – | for (i=0 ; i < info_->getNGlobalMolecules(); i++) { |
698 | – | which_node = info_->getMolToProc(i); |
699 | – | |
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 | – | simError(); |
712 | – | } |
713 | – | |
714 | – | nCurObj = mol->getNIntegrableObjects(); |
715 | – | // send the number of integrableObjects in the molecule |
716 | – | 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 | – | intObjIndexTransfer = integrableObject->getGlobalIndex(); |
723 | – | // 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 | – | MPI_Recv(&angleTransfer, 1, MPI_REALTYPE, 0, |
728 | – | TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus); |
729 | – | |
730 | – | integrableObject->setZangle(angleTransfer); |
731 | – | } |
732 | – | } |
733 | – | } |
734 | – | } |
735 | – | #endif |
736 | – | } |
737 | – | |
738 | – | void RestReader :: zeroZangle(){ |
739 | – | |
740 | – | Molecule* mol; |
741 | – | StuntDouble* integrableObject; |
742 | – | SimInfo::MoleculeIterator mi; |
743 | – | Molecule::IntegrableObjectIterator ii; |
744 | – | |
745 | – | #ifndef IS_MPI |
746 | – | // set all zAngles to 0.0 |
747 | – | for (mol = info_->beginMolecule(mi); mol != NULL; |
748 | – | mol = info_->nextMolecule(mi)) |
315 | ||
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 | – | int haveError; |
764 | – | int which_node; |
765 | – | |
766 | – | haveError = 0; |
767 | – | |
768 | – | 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 | – | |
773 | – | // 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 | – | } |
784 | – | |
785 | – | for (integrableObject = mol->beginIntegrableObject(ii); |
786 | – | integrableObject != NULL; |
787 | – | integrableObject = mol->nextIntegrableObject(ii)){ |
788 | – | |
789 | – | integrableObject->setZangle( 0.0 ); |
790 | – | |
791 | – | } |
792 | – | } |
316 | } | |
317 | < | |
795 | < | #endif |
317 | > | |
318 | } | |
319 | < | |
320 | < | } // end namespace oopse |
319 | > | |
320 | > | |
321 | > | }//end namespace oopse |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |