ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/LJFF.cpp
Revision: 1172
Committed: Wed May 12 20:14:09 2004 UTC (20 years, 2 months ago) by gezelter
File size: 11799 byte(s)
Log Message:
Removed an extraneous write

File Contents

# User Rev Content
1 gezelter 829 #include <stdlib.h>
2     #include <stdio.h>
3     #include <string.h>
4 mmeineke 559
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 mmeineke 561 "the LJFF param file./n",
69 mmeineke 559 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 mmeineke 561 LJFF::LJFF(){
121 mmeineke 559
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 mmeineke 561 strcpy( fileName, "LJFF.frc" );
165 gezelter 1172 // fprintf( stderr,"Trying to open %s\n", fileName );
166 mmeineke 559
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: %s\n"
192     "Have you tried setting the FORCE_PARAM_PATH environment "
193     "vairable?\n",
194     fileName );
195     painCave.isFatal = 1;
196     simError();
197     }
198     }
199    
200     #ifdef IS_MPI
201     }
202    
203 mmeineke 561 sprintf( checkPointMsg, "LJFF file opened sucessfully." );
204 mmeineke 559 MPIcheckPoint();
205    
206     #endif // is_mpi
207     }
208    
209    
210 mmeineke 561 LJFF::~LJFF(){
211 mmeineke 559
212     if( headAtomType != NULL ) delete headAtomType;
213    
214     #ifdef IS_MPI
215     if( worldRank == 0 ){
216     #endif // is_mpi
217    
218     fclose( frcFile );
219    
220     #ifdef IS_MPI
221     }
222     #endif // is_mpi
223     }
224    
225 mmeineke 561 void LJFF::initForceField( int ljMixRule ){
226 mmeineke 559
227     initFortran( ljMixRule, 0 );
228     }
229    
230 mmeineke 561 void LJFF::cleanMe( void ){
231 mmeineke 559
232     #ifdef IS_MPI
233    
234     // keep the linked list in the mpi version
235    
236     #else // is_mpi
237    
238     // delete the linked list in the single processor version
239    
240     if( headAtomType != NULL ) delete headAtomType;
241    
242     #endif // is_mpi
243     }
244    
245 mmeineke 561 void LJFF::readParams( void ){
246 mmeineke 559
247     atomStruct info;
248     info.last = 1; // initialize last to have the last set.
249     // if things go well, last will be set to 0
250    
251     int identNum;
252    
253    
254     bigSigma = 0.0;
255     #ifdef IS_MPI
256     if( worldRank == 0 ){
257     #endif
258    
259     // read in the atom types.
260    
261     headAtomType = new LinkedAtomType;
262    
263     fastForward( "AtomTypes", "initializeAtoms" );
264    
265     // we are now at the AtomTypes section.
266    
267     eof_test = fgets( readLine, sizeof(readLine), frcFile );
268     lineNum++;
269    
270    
271     // read a line, and start parseing out the atom types
272    
273     if( eof_test == NULL ){
274     sprintf( painCave.errMsg,
275     "Error in reading Atoms from force file at line %d.\n",
276     lineNum );
277     painCave.isFatal = 1;
278     simError();
279     }
280    
281     identNum = 1;
282     // stop reading at end of file, or at next section
283     while( readLine[0] != '#' && eof_test != NULL ){
284    
285     // toss comment lines
286     if( readLine[0] != '!' ){
287    
288     // the parser returns 0 if the line was blank
289     if( parseAtom( readLine, lineNum, info ) ){
290     info.ident = identNum;
291     headAtomType->add( info );;
292     identNum++;
293     }
294     }
295     eof_test = fgets( readLine, sizeof(readLine), frcFile );
296     lineNum++;
297     }
298    
299     #ifdef IS_MPI
300    
301     // send out the linked list to all the other processes
302    
303     sprintf( checkPointMsg,
304 mmeineke 561 "LJFF atom structures read successfully." );
305 mmeineke 559 MPIcheckPoint();
306    
307     currentAtomType = headAtomType->next; //skip the first element who is a place holder.
308     while( currentAtomType != NULL ){
309     currentAtomType->duplicate( info );
310    
311    
312    
313     sendFrcStruct( &info, mpiAtomStructType );
314    
315     sprintf( checkPointMsg,
316     "successfully sent lJ force type: \"%s\"\n",
317     info.name );
318     MPIcheckPoint();
319    
320     currentAtomType = currentAtomType->next;
321     }
322     info.last = 1;
323     sendFrcStruct( &info, mpiAtomStructType );
324    
325     }
326    
327     else{
328    
329     // listen for node 0 to send out the force params
330    
331     MPIcheckPoint();
332    
333     headAtomType = new LinkedAtomType;
334 chrisfen 976 receiveFrcStruct( &info, mpiAtomStructType );
335 mmeineke 559
336     while( !info.last ){
337    
338    
339    
340     headAtomType->add( info );
341    
342     MPIcheckPoint();
343    
344 chrisfen 976 receiveFrcStruct( &info, mpiAtomStructType );
345 mmeineke 559 }
346     }
347     #endif // is_mpi
348    
349     // call new A_types in fortran
350    
351     int isError;
352    
353     // dummy variables
354     int isLJ = 1;
355     int isDipole = 0;
356     int isSSD = 0;
357     int isGB = 0;
358 chuckv 631 int isEAM = 0;
359 gezelter 941 int isCharge = 0;
360     double charge = 0.0;
361 mmeineke 559 double dipole = 0.0;
362    
363     currentAtomType = headAtomType;
364     while( currentAtomType != NULL ){
365    
366     if( currentAtomType->name[0] != '\0' ){
367     isError = 0;
368     makeAtype( &(currentAtomType->ident),
369     &isLJ,
370     &isSSD,
371     &isDipole,
372     &isGB,
373 chuckv 631 &isEAM,
374 gezelter 941 &isCharge,
375 mmeineke 559 &(currentAtomType->epslon),
376     &(currentAtomType->sigma),
377 gezelter 941 &charge,
378 mmeineke 559 &dipole,
379     &isError );
380     if( isError ){
381     sprintf( painCave.errMsg,
382     "Error initializing the \"%s\" atom type in fortran\n",
383     currentAtomType->name );
384     painCave.isFatal = 1;
385     simError();
386     }
387     }
388     currentAtomType = currentAtomType->next;
389     }
390    
391     entry_plug->useLJ = 1;
392    
393     #ifdef IS_MPI
394     sprintf( checkPointMsg,
395 mmeineke 561 "LJFF atom structures successfully sent to fortran\n" );
396 mmeineke 559 MPIcheckPoint();
397     #endif // is_mpi
398    
399     }
400    
401    
402 mmeineke 561 void LJFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
403 mmeineke 559
404     int i;
405    
406     // initialize the atoms
407    
408    
409     for( i=0; i<nAtoms; i++ ){
410    
411     currentAtomType = headAtomType->find( the_atoms[i]->getType() );
412     if( currentAtomType == NULL ){
413     sprintf( painCave.errMsg,
414     "AtomType error, %s not found in force file.\n",
415     the_atoms[i]->getType() );
416     painCave.isFatal = 1;
417     simError();
418     }
419    
420     the_atoms[i]->setMass( currentAtomType->mass );
421     the_atoms[i]->setIdent( currentAtomType->ident );
422    
423     if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
424     }
425     }
426    
427 mmeineke 561 void LJFF::initializeBonds( int nBonds, Bond** BondArray,
428 mmeineke 559 bond_pair* the_bonds ){
429    
430     if( nBonds ){
431     sprintf( painCave.errMsg,
432 mmeineke 561 "LJFF does not support bonds.\n" );
433 mmeineke 559 painCave.isFatal = 1;
434     simError();
435     }
436     }
437    
438 mmeineke 561 void LJFF::initializeBends( int nBends, Bend** bendArray,
439 mmeineke 559 bend_set* the_bends ){
440    
441     if( nBends ){
442     sprintf( painCave.errMsg,
443 mmeineke 561 "LJFF does not support bends.\n" );
444 mmeineke 559 painCave.isFatal = 1;
445     simError();
446     }
447     }
448    
449 mmeineke 561 void LJFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
450 mmeineke 559 torsion_set* the_torsions ){
451    
452     if( nTorsions ){
453     sprintf( painCave.errMsg,
454 mmeineke 561 "LJFF does not support torsions.\n" );
455 mmeineke 559 painCave.isFatal = 1;
456     simError();
457     }
458     }
459    
460 mmeineke 561 void LJFF::fastForward( char* stopText, char* searchOwner ){
461 mmeineke 559
462     int foundText = 0;
463     char* the_token;
464    
465     rewind( frcFile );
466     lineNum = 0;
467    
468     eof_test = fgets( readLine, sizeof(readLine), frcFile );
469     lineNum++;
470     if( eof_test == NULL ){
471     sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
472     " file is empty.\n",
473     searchOwner );
474     painCave.isFatal = 1;
475     simError();
476     }
477    
478    
479     while( !foundText ){
480     while( eof_test != NULL && readLine[0] != '#' ){
481     eof_test = fgets( readLine, sizeof(readLine), frcFile );
482     lineNum++;
483     }
484     if( eof_test == NULL ){
485     sprintf( painCave.errMsg,
486     "Error fast forwarding force file for %s at "
487     "line %d: file ended unexpectedly.\n",
488     searchOwner,
489     lineNum );
490     painCave.isFatal = 1;
491     simError();
492     }
493    
494     the_token = strtok( readLine, " ,;\t#\n" );
495     foundText = !strcmp( stopText, the_token );
496    
497     if( !foundText ){
498     eof_test = fgets( readLine, sizeof(readLine), frcFile );
499     lineNum++;
500    
501     if( eof_test == NULL ){
502     sprintf( painCave.errMsg,
503     "Error fast forwarding force file for %s at "
504     "line %d: file ended unexpectedly.\n",
505     searchOwner,
506     lineNum );
507     painCave.isFatal = 1;
508     simError();
509     }
510     }
511     }
512     }
513    
514    
515    
516     int LJ_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
517    
518     char* the_token;
519    
520     the_token = strtok( lineBuffer, " \n\t,;" );
521     if( the_token != NULL ){
522    
523     strcpy( info.name, the_token );
524    
525     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
526     sprintf( painCave.errMsg,
527     "Error parseing AtomTypes: line %d\n", lineNum );
528     painCave.isFatal = 1;
529     simError();
530     }
531    
532     info.mass = atof( the_token );
533    
534     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
535     sprintf( painCave.errMsg,
536     "Error parseing AtomTypes: line %d\n", lineNum );
537     painCave.isFatal = 1;
538     simError();
539     }
540    
541     info.epslon = atof( the_token );
542    
543     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
544     sprintf( painCave.errMsg,
545     "Error parseing AtomTypes: line %d\n", lineNum );
546     painCave.isFatal = 1;
547     simError();
548     }
549    
550     info.sigma = atof( the_token );
551    
552     return 1;
553     }
554     else return 0;
555     }