ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-1.0/libmdtools/LJFF.cpp
Revision: 1426
Committed: Wed Jul 28 16:26:33 2004 UTC (19 years, 11 months ago) by gezelter
File size: 12197 byte(s)
Log Message:
Added utility functions to grab the mass from the force field.

File Contents

# User Rev Content
1 gezelter 1334 #include <stdlib.h>
2     #include <stdio.h>
3     #include <string.h>
4    
5     #include <iostream>
6     using namespace std;
7    
8     #ifdef IS_MPI
9     #include <mpi.h>
10     #endif //is_mpi
11    
12     #include "ForceFields.hpp"
13     #include "SRI.hpp"
14     #include "simError.h"
15    
16     #include "fortranWrappers.hpp"
17    
18     #ifdef IS_MPI
19     #include "mpiForceField.h"
20     #endif // is_mpi
21    
22    
23    
24     namespace LJ_NS{
25    
26     // Declare the structures that will be passed by the parser and MPI
27    
28     typedef struct{
29     char name[15];
30     double mass;
31     double epslon;
32     double sigma;
33     int ident;
34     int last; // 0 -> default
35     // 1 -> in MPI: tells nodes to stop listening
36     } atomStruct;
37    
38     int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
39    
40     #ifdef IS_MPI
41    
42     MPI_Datatype mpiAtomStructType;
43    
44     #endif
45    
46     class LinkedAtomType {
47     public:
48     LinkedAtomType(){
49     next = NULL;
50     name[0] = '\0';
51     }
52     ~LinkedAtomType(){ if( next != NULL ) delete next; }
53    
54     LinkedAtomType* find(char* key){
55     if( !strcmp(name, key) ) return this;
56     if( next != NULL ) return next->find(key);
57     return NULL;
58     }
59    
60    
61     void add( atomStruct &info ){
62    
63     // check for duplicates
64    
65     if( !strcmp( info.name, name ) ){
66     sprintf( painCave.errMsg,
67     "Duplicate LJ atom type \"%s\" found in "
68     "the LJFF param file./n",
69     name );
70     painCave.isFatal = 1;
71     simError();
72     }
73    
74     if( next != NULL ) next->add(info);
75     else{
76     next = new LinkedAtomType();
77     strcpy(next->name, info.name);
78     next->mass = info.mass;
79     next->epslon = info.epslon;
80     next->sigma = info.sigma;
81     next->ident = info.ident;
82     }
83     }
84    
85    
86     #ifdef IS_MPI
87    
88     void duplicate( atomStruct &info ){
89     strcpy(info.name, name);
90     info.mass = mass;
91     info.epslon = epslon;
92     info.sigma = sigma;
93     info.ident = ident;
94     info.last = 0;
95     }
96    
97    
98     #endif
99    
100     char name[15];
101     double mass;
102     double epslon;
103     double sigma;
104     int ident;
105     LinkedAtomType* next;
106     };
107    
108     LinkedAtomType* headAtomType;
109     LinkedAtomType* currentAtomType;
110    
111     }
112    
113     using namespace LJ_NS;
114    
115     //****************************************************************
116     // begins the actual forcefield stuff.
117     //****************************************************************
118    
119    
120     LJFF::LJFF(){
121    
122     char fileName[200];
123     char* ffPath_env = "FORCE_PARAM_PATH";
124     char* ffPath;
125     char temp[200];
126    
127     headAtomType = NULL;
128     currentAtomType = NULL;
129    
130     // do the funtion wrapping
131     wrapMeFF( this );
132    
133     #ifdef IS_MPI
134     int i;
135    
136     // **********************************************************************
137     // Init the atomStruct mpi type
138    
139     atomStruct atomProto; // mpiPrototype
140     int atomBC[3] = {15,3,2}; // block counts
141     MPI_Aint atomDspls[3]; // displacements
142     MPI_Datatype atomMbrTypes[3]; // member mpi types
143    
144     MPI_Address(&atomProto.name, &atomDspls[0]);
145     MPI_Address(&atomProto.mass, &atomDspls[1]);
146     MPI_Address(&atomProto.ident, &atomDspls[2]);
147    
148     atomMbrTypes[0] = MPI_CHAR;
149     atomMbrTypes[1] = MPI_DOUBLE;
150     atomMbrTypes[2] = MPI_INT;
151    
152     for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
153    
154     MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
155     MPI_Type_commit(&mpiAtomStructType);
156    
157     // ***********************************************************************
158    
159     if( worldRank == 0 ){
160     #endif
161    
162     // generate the force file name
163    
164     strcpy( fileName, "LJFF.frc" );
165     // fprintf( stderr,"Trying to open %s\n", fileName );
166    
167     // attempt to open the file in the current directory first.
168    
169     frcFile = fopen( fileName, "r" );
170    
171     if( frcFile == NULL ){
172    
173     // next see if the force path enviorment variable is set
174    
175     ffPath = getenv( ffPath_env );
176     if( ffPath == NULL ) {
177     STR_DEFINE(ffPath, FRC_PATH );
178     }
179    
180    
181     strcpy( temp, ffPath );
182     strcat( temp, "/" );
183     strcat( temp, fileName );
184     strcpy( fileName, temp );
185    
186     frcFile = fopen( fileName, "r" );
187    
188     if( frcFile == NULL ){
189    
190     sprintf( painCave.errMsg,
191     "Error opening the force field parameter file:\n"
192     "\t%s\n"
193     "\tHave you tried setting the FORCE_PARAM_PATH environment "
194     "variable?\n",
195     fileName );
196     painCave.severity = OOPSE_ERROR;
197     painCave.isFatal = 1;
198     simError();
199     }
200     }
201    
202     #ifdef IS_MPI
203     }
204    
205     sprintf( checkPointMsg, "LJFF file opened sucessfully." );
206     MPIcheckPoint();
207    
208     #endif // is_mpi
209     }
210    
211    
212     LJFF::~LJFF(){
213    
214     if( headAtomType != NULL ) delete headAtomType;
215    
216     #ifdef IS_MPI
217     if( worldRank == 0 ){
218     #endif // is_mpi
219    
220     fclose( frcFile );
221    
222     #ifdef IS_MPI
223     }
224     #endif // is_mpi
225     }
226    
227     void LJFF::initForceField( int ljMixRule ){
228    
229     initFortran( ljMixRule, 0 );
230     }
231    
232     void LJFF::cleanMe( void ){
233    
234     #ifdef IS_MPI
235    
236     // keep the linked list in the mpi version
237    
238     #else // is_mpi
239    
240     // delete the linked list in the single processor version
241    
242     if( headAtomType != NULL ) delete headAtomType;
243    
244     #endif // is_mpi
245     }
246    
247     void LJFF::readParams( void ){
248    
249     atomStruct info;
250     info.last = 1; // initialize last to have the last set.
251     // if things go well, last will be set to 0
252    
253     int identNum;
254    
255    
256     bigSigma = 0.0;
257     #ifdef IS_MPI
258     if( worldRank == 0 ){
259     #endif
260    
261     // read in the atom types.
262    
263     headAtomType = new LinkedAtomType;
264    
265     fastForward( "AtomTypes", "initializeAtoms" );
266    
267     // we are now at the AtomTypes section.
268    
269     eof_test = fgets( readLine, sizeof(readLine), frcFile );
270     lineNum++;
271    
272    
273     // read a line, and start parseing out the atom types
274    
275     if( eof_test == NULL ){
276     sprintf( painCave.errMsg,
277     "Error in reading Atoms from force file at line %d.\n",
278     lineNum );
279     painCave.isFatal = 1;
280     simError();
281     }
282    
283     identNum = 1;
284     // stop reading at end of file, or at next section
285     while( readLine[0] != '#' && eof_test != NULL ){
286    
287     // toss comment lines
288     if( readLine[0] != '!' ){
289    
290     // the parser returns 0 if the line was blank
291     if( parseAtom( readLine, lineNum, info ) ){
292     info.ident = identNum;
293     headAtomType->add( info );;
294     identNum++;
295     }
296     }
297     eof_test = fgets( readLine, sizeof(readLine), frcFile );
298     lineNum++;
299     }
300    
301     #ifdef IS_MPI
302    
303     // send out the linked list to all the other processes
304    
305     sprintf( checkPointMsg,
306     "LJFF atom structures read successfully." );
307     MPIcheckPoint();
308    
309     currentAtomType = headAtomType->next; //skip the first element who is a place holder.
310     while( currentAtomType != NULL ){
311     currentAtomType->duplicate( info );
312    
313    
314    
315     sendFrcStruct( &info, mpiAtomStructType );
316    
317     sprintf( checkPointMsg,
318     "successfully sent lJ force type: \"%s\"\n",
319     info.name );
320     MPIcheckPoint();
321    
322     currentAtomType = currentAtomType->next;
323     }
324     info.last = 1;
325     sendFrcStruct( &info, mpiAtomStructType );
326    
327     }
328    
329     else{
330    
331     // listen for node 0 to send out the force params
332    
333     MPIcheckPoint();
334    
335     headAtomType = new LinkedAtomType;
336     receiveFrcStruct( &info, mpiAtomStructType );
337    
338     while( !info.last ){
339    
340    
341    
342     headAtomType->add( info );
343    
344     MPIcheckPoint();
345    
346     receiveFrcStruct( &info, mpiAtomStructType );
347     }
348     }
349     #endif // is_mpi
350    
351     // call new A_types in fortran
352    
353     int isError;
354    
355     // dummy variables
356     int isLJ = 1;
357     int isDipole = 0;
358     int isSSD = 0;
359     int isGB = 0;
360     int isEAM = 0;
361     int isCharge = 0;
362     double charge = 0.0;
363     double dipole = 0.0;
364    
365     currentAtomType = headAtomType;
366     while( currentAtomType != NULL ){
367    
368     if( currentAtomType->name[0] != '\0' ){
369     isError = 0;
370     makeAtype( &(currentAtomType->ident),
371     &isLJ,
372     &isSSD,
373     &isDipole,
374     &isGB,
375     &isEAM,
376     &isCharge,
377     &(currentAtomType->epslon),
378     &(currentAtomType->sigma),
379     &charge,
380     &dipole,
381     &isError );
382     if( isError ){
383     sprintf( painCave.errMsg,
384     "Error initializing the \"%s\" atom type in fortran\n",
385     currentAtomType->name );
386     painCave.isFatal = 1;
387     simError();
388     }
389     }
390     currentAtomType = currentAtomType->next;
391     }
392    
393     entry_plug->useLJ = 1;
394    
395     #ifdef IS_MPI
396     sprintf( checkPointMsg,
397     "LJFF atom structures successfully sent to fortran\n" );
398     MPIcheckPoint();
399     #endif // is_mpi
400    
401     }
402    
403 gezelter 1426 double LJFF::getAtomTypeMass (char* atomType) {
404 gezelter 1334
405 gezelter 1426 currentAtomType = headAtomType->find( atomType );
406     if( currentAtomType == NULL ){
407     sprintf( painCave.errMsg,
408     "AtomType error, %s not found in force file.\n",
409     atomType );
410     painCave.isFatal = 1;
411     simError();
412     }
413    
414     return currentAtomType->mass;
415     }
416    
417 gezelter 1334 void LJFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
418    
419     int i;
420    
421     // initialize the atoms
422    
423    
424     for( i=0; i<nAtoms; i++ ){
425    
426     currentAtomType = headAtomType->find( the_atoms[i]->getType() );
427     if( currentAtomType == NULL ){
428     sprintf( painCave.errMsg,
429     "AtomType error, %s not found in force file.\n",
430     the_atoms[i]->getType() );
431     painCave.isFatal = 1;
432     simError();
433     }
434    
435     the_atoms[i]->setMass( currentAtomType->mass );
436     the_atoms[i]->setIdent( currentAtomType->ident );
437    
438     if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
439     }
440     }
441    
442     void LJFF::initializeBonds( int nBonds, Bond** BondArray,
443     bond_pair* the_bonds ){
444    
445     if( nBonds ){
446     sprintf( painCave.errMsg,
447     "LJFF does not support bonds.\n" );
448     painCave.isFatal = 1;
449     simError();
450     }
451     }
452    
453     void LJFF::initializeBends( int nBends, Bend** bendArray,
454     bend_set* the_bends ){
455    
456     if( nBends ){
457     sprintf( painCave.errMsg,
458     "LJFF does not support bends.\n" );
459     painCave.isFatal = 1;
460     simError();
461     }
462     }
463    
464     void LJFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
465     torsion_set* the_torsions ){
466    
467     if( nTorsions ){
468     sprintf( painCave.errMsg,
469     "LJFF does not support torsions.\n" );
470     painCave.isFatal = 1;
471     simError();
472     }
473     }
474    
475     void LJFF::fastForward( char* stopText, char* searchOwner ){
476    
477     int foundText = 0;
478     char* the_token;
479    
480     rewind( frcFile );
481     lineNum = 0;
482    
483     eof_test = fgets( readLine, sizeof(readLine), frcFile );
484     lineNum++;
485     if( eof_test == NULL ){
486     sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
487     " file is empty.\n",
488     searchOwner );
489     painCave.isFatal = 1;
490     simError();
491     }
492    
493    
494     while( !foundText ){
495     while( eof_test != NULL && readLine[0] != '#' ){
496     eof_test = fgets( readLine, sizeof(readLine), frcFile );
497     lineNum++;
498     }
499     if( eof_test == NULL ){
500     sprintf( painCave.errMsg,
501     "Error fast forwarding force file for %s at "
502     "line %d: file ended unexpectedly.\n",
503     searchOwner,
504     lineNum );
505     painCave.isFatal = 1;
506     simError();
507     }
508    
509     the_token = strtok( readLine, " ,;\t#\n" );
510     foundText = !strcmp( stopText, the_token );
511    
512     if( !foundText ){
513     eof_test = fgets( readLine, sizeof(readLine), frcFile );
514     lineNum++;
515    
516     if( eof_test == NULL ){
517     sprintf( painCave.errMsg,
518     "Error fast forwarding force file for %s at "
519     "line %d: file ended unexpectedly.\n",
520     searchOwner,
521     lineNum );
522     painCave.isFatal = 1;
523     simError();
524     }
525     }
526     }
527     }
528    
529    
530    
531     int LJ_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
532    
533     char* the_token;
534    
535     the_token = strtok( lineBuffer, " \n\t,;" );
536     if( the_token != NULL ){
537    
538     strcpy( info.name, the_token );
539    
540     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
541     sprintf( painCave.errMsg,
542     "Error parseing AtomTypes: line %d\n", lineNum );
543     painCave.isFatal = 1;
544     simError();
545     }
546    
547     info.mass = atof( the_token );
548    
549     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
550     sprintf( painCave.errMsg,
551     "Error parseing AtomTypes: line %d\n", lineNum );
552     painCave.isFatal = 1;
553     simError();
554     }
555    
556     info.epslon = atof( the_token );
557    
558     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
559     sprintf( painCave.errMsg,
560     "Error parseing AtomTypes: line %d\n", lineNum );
561     painCave.isFatal = 1;
562     simError();
563     }
564    
565     info.sigma = atof( the_token );
566    
567     return 1;
568     }
569     else return 0;
570     }
571 gezelter 1426
572