1 |
chrisfen |
417 |
/* |
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 |
431 |
|
72 |
chrisfen |
417 |
idealName = "idealCrystal.in"; |
73 |
chrisfen |
996 |
|
74 |
chrisfen |
417 |
#ifdef IS_MPI |
75 |
|
|
if (worldRank == 0) { |
76 |
|
|
#endif |
77 |
chrisfen |
996 |
|
78 |
|
|
inIdealFile = new std::ifstream(idealName.c_str()); |
79 |
|
|
|
80 |
|
|
if(inIdealFile->fail()){ |
81 |
chrisfen |
417 |
sprintf(painCave.errMsg, |
82 |
chrisfen |
996 |
"RestReader: Cannot open file: %s\n", |
83 |
|
|
idealName.c_str()); |
84 |
chrisfen |
417 |
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 |
431 |
|
95 |
chrisfen |
417 |
return; |
96 |
|
|
} |
97 |
|
|
|
98 |
|
|
RestReader :: ~RestReader( ){ |
99 |
|
|
#ifdef IS_MPI |
100 |
|
|
if (worldRank == 0) { |
101 |
|
|
#endif |
102 |
chrisfen |
996 |
|
103 |
|
|
delete inIdealFile; |
104 |
|
|
delete inAngFile; |
105 |
|
|
|
106 |
chrisfen |
417 |
#ifdef IS_MPI |
107 |
|
|
} |
108 |
chrisfen |
996 |
strcpy( checkPointMsg, |
109 |
|
|
"File idealCrystal.in (and .zang0 if present) closed successfully." ); |
110 |
chrisfen |
417 |
MPIcheckPoint(); |
111 |
|
|
#endif |
112 |
|
|
|
113 |
|
|
return; |
114 |
|
|
} |
115 |
|
|
|
116 |
|
|
|
117 |
|
|
void RestReader :: readIdealCrystal(){ |
118 |
chrisfen |
990 |
|
119 |
chrisfen |
417 |
#ifdef IS_MPI |
120 |
chrisfen |
996 |
int which_node; |
121 |
|
|
int i, j; |
122 |
chrisfen |
417 |
#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 |
996 |
inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
140 |
|
|
|
141 |
|
|
if( inIdealFile->eof() ){ |
142 |
chrisfen |
417 |
sprintf( painCave.errMsg, |
143 |
|
|
"RestraintReader error: error reading 1st line of \"%s\"\n", |
144 |
chrisfen |
996 |
idealName.c_str() ); |
145 |
chrisfen |
417 |
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 |
996 |
idealName.c_str(), |
156 |
|
|
nTotObjs, |
157 |
chrisfen |
417 |
info_->getNGlobalIntegrableObjects()); |
158 |
|
|
painCave.isFatal = 1; |
159 |
|
|
simError(); |
160 |
|
|
} |
161 |
|
|
|
162 |
|
|
// skip over the comment line |
163 |
chrisfen |
996 |
inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
164 |
|
|
|
165 |
|
|
if( inIdealFile->eof() ){ |
166 |
chrisfen |
417 |
sprintf( painCave.errMsg, |
167 |
chrisfen |
996 |
"error in reading commment in %s\n", |
168 |
|
|
idealName.c_str()); |
169 |
chrisfen |
417 |
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 |
996 |
|
186 |
|
|
inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
187 |
|
|
|
188 |
|
|
if( inIdealFile->eof() ){ |
189 |
chrisfen |
417 |
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 |
996 |
idealName.c_str(), nTotObjs, |
194 |
chrisfen |
417 |
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 |
423 |
int haveError; |
220 |
|
|
|
221 |
|
|
nTotObjs = info_->getNGlobalIntegrableObjects(); |
222 |
chrisfen |
417 |
haveError = 0; |
223 |
chrisfen |
431 |
|
224 |
chrisfen |
417 |
if (worldRank == masterNode) { |
225 |
chrisfen |
996 |
inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
226 |
|
|
|
227 |
|
|
if( inIdealFile->eof() ){ |
228 |
chrisfen |
417 |
sprintf( painCave.errMsg, |
229 |
chrisfen |
996 |
"Error reading 1st line of %s \n ",idealName.c_str()); |
230 |
chrisfen |
417 |
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 |
996 |
idealName.c_str(), nTotObjs, |
244 |
chrisfen |
417 |
info_->getNGlobalIntegrableObjects()); |
245 |
|
|
painCave.isFatal = 1; |
246 |
|
|
simError(); |
247 |
|
|
} |
248 |
|
|
|
249 |
|
|
// skip over the comment line |
250 |
chrisfen |
996 |
inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
251 |
|
|
|
252 |
|
|
if( inIdealFile->eof() ){ |
253 |
chrisfen |
417 |
sprintf( painCave.errMsg, |
254 |
chrisfen |
996 |
"error in reading commment in %s\n", idealName.c_str()); |
255 |
chrisfen |
417 |
painCave.isFatal = 1; |
256 |
|
|
simError(); |
257 |
|
|
} |
258 |
chrisfen |
431 |
|
259 |
chrisfen |
990 |
MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
260 |
|
|
|
261 |
chrisfen |
417 |
for (i=0 ; i < info_->getNGlobalMolecules(); i++) { |
262 |
chrisfen |
996 |
which_node = info_->getMolToProc(i); |
263 |
chrisfen |
417 |
|
264 |
|
|
if(which_node == masterNode){ |
265 |
|
|
//molecules belong to master node |
266 |
|
|
|
267 |
chrisfen |
423 |
mol = info_->getMoleculeByGlobalIndex(i); |
268 |
chrisfen |
417 |
|
269 |
chrisfen |
423 |
if(mol == NULL) { |
270 |
|
|
sprintf(painCave.errMsg, |
271 |
gezelter |
507 |
"RestReader Error: Molecule not found on node %d!\n", |
272 |
|
|
worldRank); |
273 |
chrisfen |
417 |
painCave.isFatal = 1; |
274 |
|
|
simError(); |
275 |
|
|
} |
276 |
|
|
|
277 |
|
|
for (integrableObject = mol->beginIntegrableObject(ii); |
278 |
|
|
integrableObject != NULL; |
279 |
|
|
integrableObject = mol->nextIntegrableObject(ii)){ |
280 |
|
|
|
281 |
chrisfen |
996 |
inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
282 |
chrisfen |
417 |
|
283 |
chrisfen |
996 |
if( inIdealFile->eof() ){ |
284 |
chrisfen |
417 |
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 |
996 |
idealName.c_str(), nTotObjs, i ); |
289 |
chrisfen |
417 |
painCave.isFatal = 1; |
290 |
|
|
simError(); |
291 |
|
|
} |
292 |
chrisfen |
431 |
|
293 |
|
|
parseIdealLine(read_buffer, integrableObject); |
294 |
|
|
|
295 |
chrisfen |
417 |
} |
296 |
chrisfen |
990 |
|
297 |
chrisfen |
417 |
} 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 |
990 |
for(j = 0; j < nCurObj; j++){ |
304 |
chrisfen |
996 |
inIdealFile->getline(read_buffer, sizeof(read_buffer)); |
305 |
|
|
|
306 |
|
|
if( inIdealFile->eof() ){ |
307 |
chrisfen |
417 |
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 |
996 |
idealName.c_str(), nTotObjs, i ); |
312 |
chrisfen |
417 |
painCave.isFatal = 1; |
313 |
|
|
simError(); |
314 |
|
|
} |
315 |
|
|
|
316 |
chrisfen |
423 |
MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node, |
317 |
chrisfen |
417 |
TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD); |
318 |
|
|
} |
319 |
|
|
} |
320 |
|
|
} |
321 |
|
|
} else { |
322 |
chrisfen |
990 |
//actions taken at slave nodes |
323 |
|
|
MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
324 |
|
|
|
325 |
chrisfen |
417 |
for (i=0 ; i < info_->getNGlobalMolecules(); i++) { |
326 |
chrisfen |
990 |
int which_node = info_->getMolToProc(i); |
327 |
|
|
|
328 |
chrisfen |
417 |
if(which_node == worldRank){ |
329 |
|
|
//molecule with global index i belongs to this processor |
330 |
|
|
|
331 |
chrisfen |
423 |
mol = info_->getMoleculeByGlobalIndex(i); |
332 |
chrisfen |
417 |
|
333 |
chrisfen |
423 |
if(mol == NULL) { |
334 |
chrisfen |
417 |
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 |
423 |
nCurObj = mol->getNIntegrableObjects(); |
342 |
chrisfen |
417 |
|
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 |
963 |
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 |
417 |
StringTokenizer tokenizer(readLine); |
375 |
|
|
int nTokens; |
376 |
|
|
|
377 |
|
|
nTokens = tokenizer.countTokens(); |
378 |
chrisfen |
431 |
|
379 |
chrisfen |
417 |
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 |
431 |
|
388 |
chrisfen |
417 |
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 |
996 |
name.c_str(), idealName.c_str(), |
394 |
chrisfen |
417 |
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 |
431 |
|
403 |
chrisfen |
417 |
// 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 |
431 |
|
409 |
chrisfen |
417 |
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 |
431 |
|
419 |
chrisfen |
417 |
// 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 |
996 |
int which_node; |
481 |
chrisfen |
417 |
MPI_Status istatus; |
482 |
|
|
#endif //is_mpi |
483 |
|
|
|
484 |
|
|
const int BUFFERSIZE = 2000; // size of the read buffer |
485 |
chrisfen |
996 |
unsigned int nTotObjs; // the number of atoms |
486 |
chrisfen |
417 |
char read_buffer[BUFFERSIZE]; //the line buffer for reading |
487 |
|
|
|
488 |
|
|
std::vector<StuntDouble*> vecParticles; |
489 |
tim |
963 |
std::vector<RealType> tempZangs; |
490 |
chrisfen |
417 |
|
491 |
chrisfen |
996 |
angFile = info_->getRestFileName(); |
492 |
chrisfen |
417 |
|
493 |
chrisfen |
996 |
angFile += "0"; |
494 |
chrisfen |
417 |
|
495 |
|
|
// open the omega value file for reading |
496 |
|
|
#ifdef IS_MPI |
497 |
|
|
if (worldRank == 0) { |
498 |
|
|
#endif |
499 |
|
|
isPresent = 1; |
500 |
chrisfen |
996 |
|
501 |
|
|
inAngFile = new std::ifstream(angFile.c_str()); |
502 |
|
|
|
503 |
|
|
if(inAngFile->fail()){ |
504 |
chrisfen |
417 |
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 |
996 |
angFile.c_str()); |
511 |
chrisfen |
417 |
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 |
996 |
if (!isPresent) { |
529 |
|
|
// master node zeroes out its zAngles if .zang0 isn't present |
530 |
|
|
zeroZangle(); |
531 |
|
|
return; |
532 |
|
|
} |
533 |
|
|
|
534 |
chrisfen |
417 |
} |
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 |
996 |
// 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 |
417 |
sprintf( painCave.errMsg, |
556 |
|
|
"RestraintReader error: error reading 1st line of \"%s\"\n", |
557 |
chrisfen |
996 |
angFile.c_str() ); |
558 |
chrisfen |
417 |
painCave.isFatal = 1; |
559 |
|
|
simError(); |
560 |
|
|
} |
561 |
chrisfen |
996 |
|
562 |
|
|
// read the file and load the values into a vector |
563 |
|
|
inAngFile->getline(read_buffer, sizeof(read_buffer)); |
564 |
chrisfen |
417 |
|
565 |
chrisfen |
996 |
while ( !inAngFile->eof() ) { |
566 |
chrisfen |
417 |
// check for and ignore blank lines |
567 |
|
|
if ( read_buffer != NULL ) |
568 |
|
|
tempZangs.push_back( atof(read_buffer) ); |
569 |
chrisfen |
996 |
|
570 |
|
|
inAngFile->getline(read_buffer, sizeof(read_buffer)); |
571 |
chrisfen |
417 |
} |
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 |
996 |
"does not match the meta-data file's nIntegrable, %i.\n", |
579 |
|
|
angFile.c_str(), |
580 |
|
|
tempZangs.size(), |
581 |
|
|
nTotObjs ); |
582 |
chrisfen |
417 |
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 |
423 |
|
606 |
|
|
int masterNode = 0; |
607 |
chrisfen |
996 |
|
608 |
chrisfen |
423 |
int haveError; |
609 |
chrisfen |
990 |
int intObjIndex; |
610 |
|
|
int intObjIndexTransfer; |
611 |
chrisfen |
423 |
|
612 |
chrisfen |
996 |
int j; |
613 |
chrisfen |
417 |
int nCurObj; |
614 |
chrisfen |
996 |
RealType angleTransfer; |
615 |
chrisfen |
417 |
|
616 |
chrisfen |
423 |
nTotObjs = info_->getNGlobalIntegrableObjects(); |
617 |
chrisfen |
417 |
haveError = 0; |
618 |
chrisfen |
423 |
|
619 |
|
|
if (worldRank == masterNode) { |
620 |
chrisfen |
417 |
|
621 |
chrisfen |
996 |
inAngFile->getline(read_buffer, sizeof(read_buffer)); |
622 |
|
|
|
623 |
|
|
if( inAngFile->eof() ){ |
624 |
chrisfen |
417 |
sprintf( painCave.errMsg, |
625 |
chrisfen |
996 |
"Error reading 1st line of %s \n ",angFile.c_str()); |
626 |
chrisfen |
417 |
haveError = 1; |
627 |
|
|
simError(); |
628 |
|
|
} |
629 |
|
|
|
630 |
chrisfen |
996 |
// 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 |
417 |
// check for and ignore blank lines |
635 |
|
|
if ( read_buffer != NULL ) |
636 |
|
|
tempZangs.push_back( atof(read_buffer) ); |
637 |
chrisfen |
996 |
|
638 |
|
|
inAngFile->getline(read_buffer, sizeof(read_buffer)); |
639 |
chrisfen |
417 |
} |
640 |
chrisfen |
990 |
|
641 |
chrisfen |
417 |
// 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 |
996 |
angFile.c_str(), |
649 |
|
|
tempZangs.size(), |
650 |
|
|
nTotObjs); |
651 |
chrisfen |
417 |
haveError= 1; |
652 |
|
|
simError(); |
653 |
|
|
} |
654 |
|
|
|
655 |
chrisfen |
423 |
// At this point, node 0 has a tempZangs vector completed, and |
656 |
|
|
// everyone else has nada |
657 |
chrisfen |
417 |
|
658 |
chrisfen |
423 |
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 |
985 |
if (which_node == masterNode) { |
663 |
chrisfen |
423 |
mol = info_->getMoleculeByGlobalIndex(i); |
664 |
chrisfen |
985 |
|
665 |
chrisfen |
423 |
if(mol == NULL) { |
666 |
|
|
strcpy(painCave.errMsg, "Molecule not found on node 0!"); |
667 |
|
|
haveError = 1; |
668 |
|
|
simError(); |
669 |
|
|
} |
670 |
chrisfen |
985 |
|
671 |
chrisfen |
423 |
for (integrableObject = mol->beginIntegrableObject(ii); |
672 |
|
|
integrableObject != NULL; |
673 |
|
|
integrableObject = mol->nextIntegrableObject(ii)){ |
674 |
chrisfen |
1030 |
intObjIndex = integrableObject->getGlobalIntegrableObjectIndex(); |
675 |
chrisfen |
990 |
integrableObject->setZangle(tempZangs[intObjIndex]); |
676 |
chrisfen |
423 |
} |
677 |
|
|
|
678 |
|
|
} else { |
679 |
|
|
// I am MASTER OF THE UNIVERSE, but I don't own this molecule |
680 |
chrisfen |
990 |
// listen for the number of integrableObjects in the molecule |
681 |
chrisfen |
423 |
MPI_Recv(&nCurObj, 1, MPI_INT, which_node, |
682 |
|
|
TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus); |
683 |
|
|
|
684 |
chrisfen |
990 |
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 |
963 |
MPI_Send(&angleTransfer, 1, MPI_REALTYPE, which_node, |
691 |
chrisfen |
423 |
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 |
996 |
which_node = info_->getMolToProc(i); |
699 |
chrisfen |
423 |
|
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 |
417 |
simError(); |
712 |
|
|
} |
713 |
chrisfen |
423 |
|
714 |
|
|
nCurObj = mol->getNIntegrableObjects(); |
715 |
chrisfen |
990 |
// send the number of integrableObjects in the molecule |
716 |
chrisfen |
423 |
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 |
1030 |
intObjIndexTransfer = integrableObject->getGlobalIntegrableObjectIndex(); |
723 |
chrisfen |
990 |
// 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 |
963 |
MPI_Recv(&angleTransfer, 1, MPI_REALTYPE, 0, |
728 |
chrisfen |
423 |
TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus); |
729 |
|
|
|
730 |
|
|
integrableObject->setZangle(angleTransfer); |
731 |
|
|
} |
732 |
|
|
} |
733 |
chrisfen |
417 |
} |
734 |
chrisfen |
423 |
} |
735 |
chrisfen |
417 |
#endif |
736 |
|
|
} |
737 |
|
|
|
738 |
|
|
void RestReader :: zeroZangle(){ |
739 |
|
|
|
740 |
|
|
Molecule* mol; |
741 |
|
|
StuntDouble* integrableObject; |
742 |
|
|
SimInfo::MoleculeIterator mi; |
743 |
|
|
Molecule::IntegrableObjectIterator ii; |
744 |
chrisfen |
1030 |
|
745 |
chrisfen |
417 |
#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 |
423 |
int haveError; |
764 |
chrisfen |
417 |
int which_node; |
765 |
|
|
|
766 |
|
|
haveError = 0; |
767 |
chrisfen |
423 |
|
768 |
chrisfen |
996 |
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 |
423 |
|
773 |
chrisfen |
996 |
// 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 |
423 |
} |
784 |
|
|
|
785 |
chrisfen |
996 |
for (integrableObject = mol->beginIntegrableObject(ii); |
786 |
|
|
integrableObject != NULL; |
787 |
|
|
integrableObject = mol->nextIntegrableObject(ii)){ |
788 |
chrisfen |
423 |
|
789 |
chrisfen |
996 |
integrableObject->setZangle( 0.0 ); |
790 |
chrisfen |
423 |
|
791 |
chrisfen |
996 |
} |
792 |
chrisfen |
417 |
} |
793 |
|
|
} |
794 |
chrisfen |
996 |
|
795 |
chrisfen |
417 |
#endif |
796 |
|
|
} |
797 |
|
|
|
798 |
|
|
} // end namespace oopse |