ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/Shapes_FF.cpp
Revision: 1646
Committed: Tue Oct 26 18:02:46 2004 UTC (19 years, 8 months ago) by chrisfen
File size: 13714 byte(s)
Log Message:
Changes to Shapes force field reader

File Contents

# Content
1 /*
2 * Shapes_FF.cpp
3 * oopse
4 *
5 * Created by Chris Fennell on 10/20/04.
6 * Copyright 2004 __MyCompanyName__. All rights reserved.
7 *
8 */
9
10 #include <stdlib.h>
11 #include <stdio.h>
12 #include <string.h>
13 #include <map>
14 #include <iostream>
15
16 using namespace std;
17 using namespace oopse;
18
19 #ifdef IS_MPI
20 #include <mpi.h>
21 #endif //is_mpi
22
23 #include "UseTheForce/ForceFields.hpp"
24 #include "primitives/SRI.hpp"
25 #include "utils/simError.h"
26 #include "io/basic_ifstrstream.hpp"
27 #include "math/RealSphericalHarmonic.hpp"
28 #include "math/SquareMatrix3.hpp"
29
30 #include "UseTheForce/DarkSide/atype_interface.h"
31 #include "UseTheForce/DarkSide/shapes_interface.h"
32
33 #ifdef IS_MPI
34 #include "UseTheForce/mpiForceField.h"
35 #endif // is_mpi
36
37 Shapes_FF::Shapes_FF() {
38 Shapes_FF("");
39 }
40
41 Shapes_FF::Shapes_FF(char* the_variant){
42 ffPath_env = "FORCE_PARAM_PATH";
43
44 }
45
46 Shapes_FF::~Shapes_FF(){
47
48 #ifdef IS_MPI
49 if( worldRank == 0 ){
50 #endif // is_mpi
51
52 fclose( frcFile );
53
54 #ifdef IS_MPI
55 }
56 #endif // is_mpi
57 }
58
59
60 void Shapes_FF::calcRcut( void ){
61
62 #ifdef IS_MPI
63 double tempShapesRcut = shapesRcut;
64 MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
65 MPI_COMM_WORLD);
66 #endif //is_mpi
67 entry_plug->setDefaultRcut(shapesRcut);
68 }
69
70
71 void Shapes_FF::initForceField(){
72 initFortran(0);
73 }
74
75
76 void Shapes_FF::readForceFile( void ){
77
78 char readLine[1024];
79 char fileName[200];
80 char temp[200];
81 char* ffPath;
82 char *shapeFileName;
83 char *nameToken;
84 char *delim = " ,;\t\n";
85 int nTokens, i;
86 int nContact = 0;
87 int nRange = 0;
88 int nStrength = 0;
89 int myATID;
90 string nameString;
91 ShapeType* st;
92 map<string, AtomType*> atomTypeMap;
93 map<string, AtomType*>::iterator iter;
94
95 // vectors for shape transfer to fortran
96 vector<RealSphericalHarmonic> tempSHVector;
97 vector<int> contactL;
98 vector<int> contactM;
99 vector<int> contactFunc;
100 vector<double> contactCoeff;
101 vector<int> rangeL;
102 vector<int> rangeM;
103 vector<int> rangeFunc;
104 vector<double> rangeCoeff;
105 vector<int> strengthL;
106 vector<int> strengthM;
107 vector<int> strengthFunc;
108 vector<double> strengthCoeff;
109
110 // build a basic file reader
111 ifstrsteam frcReader;
112
113 // generate the force file name
114 strcpy( fileName, "Shapes.frc" );
115
116 // attempt to open the file in the current directory first.
117 frcReader.open( fileName );
118
119 if( frcReader == NULL ){
120
121 // next see if the force path enviorment variable is set
122
123 ffPath = getenv( ffPath_env );
124 if( ffPath == NULL ) {
125 STR_DEFINE(ffPath, FRC_PATH );
126 }
127
128 strcpy( temp, ffPath );
129 strcat( temp, "/" );
130 strcat( temp, fileName );
131 strcpy( fileName, temp );
132
133 frcReader.open( fileName );
134
135 if( frcFile == NULL ){
136
137 sprintf( painCave.errMsg,
138 "Error opening the force field parameter file:\n"
139 "\t%s\n"
140 "\tHave you tried setting the FORCE_PARAM_PATH environment "
141 "variable?\n",
142 fileName );
143 painCave.severity = OOPSE_ERROR;
144 painCave.isFatal = 1;
145 simError();
146 }
147 }
148
149 // read in the shape types.
150
151 findBegin( "ShapeTypes" );
152
153 while( !frcReader.eof() ){
154 frcReader.getline( readLine, sizeof(readLine) );
155
156 // toss comment lines
157 if( readLine[0] != '!' && readLine[0] != '#' ){
158
159 if (isEndLine(readLine)) break;
160
161 nTokens = count_tokens(readLine, " ,;\t");
162 if (nTokens != 0) {
163 if (nTokens < 2) {
164 sprintf( painCave.errMsg,
165 "Not enough data on a ShapeTypes line in file: %s\n"
166 fileName );
167 painCave.severity = OOPSE_ERROR;
168 painCave.isFatal = 1;
169 simError();
170 }
171
172 nameToken = strtok( readLine, delim );
173 shapeFileName = strtok( NULL, delim );
174
175 // strings are not char arrays!
176 nameString = nameToken;
177
178 // does this AtomType name already exist in the map?
179 iter = atomTypeMap.find(nameString);
180 if (iter == atomTypeMap.end()) {
181 // no, it doesn't, so we may proceed:
182
183 st = new ShapeType();
184
185 st->setName(nameString);
186 myATID = atomTypeMap.size();
187 st->setIdent(myATID);
188 parseShapeFile(shapeFileName, st);
189 st->complete();
190 atomTypeMap.insert(make_pair(nameString, st));
191
192 } else {
193 // atomType map already contained this string (i.e. it was
194 // declared in a previous block, and we just need to add
195 // the shape-specific information for this AtomType:
196
197 st = (ShapeType*)(iter->second);
198 parseShapeFile(shapeFileName, st);
199 }
200 }
201 }
202 }
203
204 #ifdef IS_MPI
205
206 // looks like all the processors have their ShapeType vectors ready...
207 sprintf( checkPointMsg,
208 "Shapes_FF shape objects read successfully." );
209 MPIcheckPoint();
210
211 #endif // is_mpi
212
213 // pack up and send the necessary info to fortran
214 for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
215
216 at = (AtomType*)(iter.second);
217
218 if (at->isDirectional()) {
219
220 dat = (DirectionalAtomType*)at;
221
222 if (dat->isShape()) {
223
224 st = (ShapeAtomType*)at;
225
226 contactL.clear();
227 contactM.clear();
228 contactFunc.clear();
229 contactCoeff.clear();
230
231 tempSHVector = st->getContactFuncs();
232
233 nContact = tempSHVector.size();
234 for (i=0; i<nContact; i++){
235 contactL.push_back(tempSHVector[i].getL());
236 contactM.push_back(tempSHVector[i].getM());
237 contactFunc.push_back(tempSHVector[i].getFunctionType());
238 contactCoeff.push_back(tempSHVector[i].getCoefficient());
239 }
240
241 rangeL.clear();
242 rangeM.clear();
243 rangeFunc.clear();
244 rangeCoeff.clear();
245
246 tempSHVector = st->getRangeFuncs();
247
248 nRange = tempSHVector.size();
249 for (i=0; i<nRange; i++){
250 rangeL.push_back(tempSHVector[i].getL());
251 rangeM.push_back(tempSHVector[i].getM());
252 rangeFunc.push_back(tempSHVector[i].getFunctionType());
253 rangeCoeff.push_back(tempSHVector[i].getCoefficient());
254 }
255
256 strengthL.clear();
257 strengthM.clear();
258 strengthFunc.clear();
259 strengthCoeff.clear();
260
261 tempSHVector = st->getStrengthFuncs();
262
263 nStrength = tempSHVector.size();
264 for (i=0; i<nStrength; i++){
265 strengthL.push_back(tempSHVector[i].getL());
266 strengthM.push_back(tempSHVector[i].getM());
267 strengthFunc.push_back(tempSHVector[i].getFunctionType());
268 strengthCoeff.push_back(tempSHVector[i].getCoefficient());
269 }
270
271 isError = 0;
272 myATID = at->getIdent();
273 makeShape( &nContact, &contactL, &contactM, &contactFunc,
274 &contactCoeff,
275 &nRange, &rangeL, &rangeM, &rangeFunc, &rangeCoeff,
276 &nStrength, &strengthL, &strengthM,
277 &strengthFunc, &strengthCoeff,
278 &myATID,
279 &isError);
280 if( isError ){
281 sprintf( painCave.errMsg,
282 "Error initializing the \"%s\" shape in fortran\n",
283 marker->first );
284 painCave.isFatal = 1;
285 simError();
286 }
287 }
288 }
289 }
290
291 #ifdef IS_MPI
292 sprintf( checkPointMsg,
293 "Shapes_FF atom structures successfully sent to fortran\n" );
294 MPIcheckPoint();
295 #endif // is_mpi
296
297 }
298
299 void SHAPES_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
300
301 int i,j,k;
302
303 // initialize the atoms
304 DirectionalAtom* dAtom;
305 AtomType* at;
306 DirectionalAtomType* dat;
307 double mySigma;
308 double ji[3];
309 double inertialMat[3][3];
310 Mat3x3d momInt;
311 string myTypeString;
312
313 for( i=0; i<nAtoms; i++ ){
314
315 myTypeString = the_atoms[i]->getType();
316 iter = atomTypeMap.find(myTypeString);
317
318 if (iter == atomTypeMap.end()) {
319 sprintf( painCave.errMsg,
320 "AtomType error, %s not found in force file.\n",
321 the_atoms[i]->getType() );
322 painCave.isFatal = 1;
323 simError();
324 } else {
325
326 at = (AtomType*)(iter->second);
327
328 the_atoms[i]->setMass( at->getMass() );
329 the_atoms[i]->setIdent( at->getIdent() );
330
331 if( at->isLennardJones() ) {
332 mySigma = at->properties->getPropertyByName("sigma");
333 if( bigSigma < mySigma ) bigSigma = mySigma;
334 }
335
336 the_atoms[i]->setHasCharge(at->isCharge());
337
338 if( at->isDirectional() ){
339
340 dat = (DirectionalAtomType*)at;
341 dAtom = (DirectionalAtom *) the_atoms[i];
342
343 momInt = dat->getI();
344
345 // zero out the moments of inertia matrix
346 for( j=0; j<3; j++ )
347 for( k=0; k<3; k++ )
348 inertialMat[j][k] = momInt(j,k);
349 dAtom->setI( inertialMat );
350
351 ji[0] = 0.0;
352 ji[1] = 0.0;
353 ji[2] = 0.0;
354 dAtom->setJ( ji );
355
356 if (dat->isDipole()) {
357 dAtom->setHasDipole( dat->isDipole() );
358 entry_plug->n_dipoles++;
359 }
360 }
361 }
362 }
363 }
364
365 void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
366 bond_pair* the_bonds ){
367
368 if( nBonds ){
369 sprintf( painCave.errMsg,
370 "Shapes_FF does not support bonds.\n" );
371 painCave.isFatal = 1;
372 simError();
373 }
374 }
375
376 void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
377 bend_set* the_bends ){
378
379 if( nBends ){
380 sprintf( painCave.errMsg,
381 "Shapes_FF does not support bends.\n" );
382 painCave.isFatal = 1;
383 simError();
384 }
385 }
386
387 void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
388 torsion_set* the_torsions ){
389
390 if( nTorsions ){
391 sprintf( painCave.errMsg,
392 "Shapes_FF does not support torsions.\n" );
393 painCave.isFatal = 1;
394 simError();
395 }
396 }
397
398 int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAtomType* st){
399 const int MAXLEN = 1024;
400 char inLine[MAXLEN];
401 char *token;
402 char *delim = " ,;\t\n";
403 Mat3x3d momInert;
404 RealSphericalHarmonic tempHarmonic;
405 vector<RealSphericalHarmonic> functionVector;
406
407 ifstrstream shapeFile;
408
409 // first grab the values in the ShapeInfo section
410 findBegin(shapeFile, "ShapeInfo");
411
412 shapeFile.getline(inLine, MAXLEN);
413 while( !shapeFile.eof() ) {
414 // toss comment lines
415 if( inLine[0] != '!' && inLine[0] != '#' ){
416 // end marks section completion
417 if (isEndLine(inLine)) break;
418
419 nTokens = count_tokens(inLine, delim);
420 if (nTokens != 0) {
421 if (nTokens < 5) {
422 sprintf( painCave.errMsg,
423 "Not enough data on a ShapeInfo line in file: %s\n"
424 fileName );
425 painCave.severity = OOPSE_ERROR;
426 painCave.isFatal = 1;
427 simError();
428
429 token = strtok(inLine, delim);
430 token = strtok(NULL, delim);
431 st->setMass(atof(token));
432 token = strtok(NULL, delim);
433 momInert(0,0) = atof(token);
434 token = strtok(NULL, delim);
435 momInert(1,1) = atof(token);
436 token = strtok(NULL, delim);
437 momInert(2,2) = atof(token);
438 st->setI(momInert);
439 }
440 }
441 }
442 }
443
444 // now grab the contact functions
445 findBegin(shapeFile, "ContactFunctions");
446 functionVector.clear();
447
448 shapeFile.getline(inLine, MAXLEN);
449 while( !shapeFile.eof() ) {
450 // toss comment lines
451 if( inLine[0] != '!' && inLine[0] != '#' ){
452 // end marks section completion
453 if (isEndLine(inLine)) break;
454
455 nTokens = count_tokens(inLine, delim);
456 if (nTokens != 0) {
457 if (nTokens < 4) {
458 sprintf( painCave.errMsg,
459 "Not enough data on a ContactFunctions line in file: %s\n"
460 fileName );
461 painCave.severity = OOPSE_ERROR;
462 painCave.isFatal = 1;
463 simError();
464
465 // read in a spherical harmonic function
466 token = strtok(inLine, delim);
467 tempHarmonic.setL(atoi(token));
468 token = strtok(NULL, delim);
469 tempHarmonic.setM(atoi(token));
470 token = strtok(NULL, delim);
471 if (!strcasecmp("sin",token))
472 tempHarmonic.makeSinFunction();
473 else
474 tempHarmonic.makeCosFunction();
475 token = strtok(NULL, delim);
476 tempHarmonic.setCoefficient(atof(token));
477
478 functionVector.push_back(tempHarmonic);
479 }
480 }
481 }
482 }
483
484 // pass contact functions to ShapeType
485 st->setContactFuncs(functionVector);
486
487 // now grab the range functions
488 findBegin(shapeFile, "RangeFunctions");
489 functionVector.clear();
490
491 shapeFile.getline(inLine, MAXLEN);
492 while( !shapeFile.eof() ) {
493 // toss comment lines
494 if( inLine[0] != '!' && inLine[0] != '#' ){
495 // end marks section completion
496 if (isEndLine(inLine)) break;
497
498 nTokens = count_tokens(inLine, delim);
499 if (nTokens != 0) {
500 if (nTokens < 4) {
501 sprintf( painCave.errMsg,
502 "Not enough data on a RangeFunctions line in file: %s\n"
503 fileName );
504 painCave.severity = OOPSE_ERROR;
505 painCave.isFatal = 1;
506 simError();
507
508 // read in a spherical harmonic function
509 token = strtok(inLine, delim);
510 tempHarmonic.setL(atoi(token));
511 token = strtok(NULL, delim);
512 tempHarmonic.setM(atoi(token));
513 token = strtok(NULL, delim);
514 if (!strcasecmp("sin",token))
515 tempHarmonic.makeSinFunction();
516 else
517 tempHarmonic.makeCosFunction();
518 token = strtok(NULL, delim);
519 tempHarmonic.setCoefficient(atof(token));
520
521 functionVector.push_back(tempHarmonic);
522 }
523 }
524 }
525 }
526
527 // pass range functions to ShapeType
528 st->setRangeFuncs(functionVector);
529
530 // finally grab the strength functions
531 findBegin(shapeFile, "StrengthFunctions");
532 functionVector.clear();
533
534 shapeFile.getline(inLine, MAXLEN);
535 while( !shapeFile.eof() ) {
536 // toss comment lines
537 if( inLine[0] != '!' && inLine[0] != '#' ){
538 // end marks section completion
539 if (isEndLine(inLine)) break;
540
541 nTokens = count_tokens(inLine, delim);
542 if (nTokens != 0) {
543 if (nTokens < 4) {
544 sprintf( painCave.errMsg,
545 "Not enough data on a StrengthFunctions line in file: %s\n"
546 fileName );
547 painCave.severity = OOPSE_ERROR;
548 painCave.isFatal = 1;
549 simError();
550
551 // read in a spherical harmonic function
552 token = strtok(inLine, delim);
553 tempHarmonic.setL(atoi(token));
554 token = strtok(NULL, delim);
555 tempHarmonic.setM(atoi(token));
556 token = strtok(NULL, delim);
557 if (!strcasecmp("sin",token))
558 tempHarmonic.makeSinFunction();
559 else
560 tempHarmonic.makeCosFunction();
561 token = strtok(NULL, delim);
562 tempHarmonic.setCoefficient(atof(token));
563
564 functionVector.push_back(tempHarmonic);
565 }
566 }
567 }
568 }
569
570 // pass strength functions to ShapeType
571 st->setStrengthFuncs(functionVector);
572
573 // we're done reading from this file
574 shapeFile.close();
575 return 0;
576 }
577