ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/LJFF.cpp
Revision: 631
Committed: Thu Jul 17 19:25:51 2003 UTC (20 years, 11 months ago) by chuckv
File size: 11882 byte(s)
Log Message:
Added massive changes for eam....

File Contents

# User Rev Content
1 mmeineke 559 #include <cstdlib>
2     #include <cstdio>
3     #include <cstring>
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 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     char errMsg[1000];
127    
128     headAtomType = NULL;
129     currentAtomType = NULL;
130    
131     // do the funtion wrapping
132     wrapMeFF( this );
133    
134     #ifdef IS_MPI
135     int i;
136    
137     // **********************************************************************
138     // Init the atomStruct mpi type
139    
140     atomStruct atomProto; // mpiPrototype
141     int atomBC[3] = {15,3,2}; // block counts
142     MPI_Aint atomDspls[3]; // displacements
143     MPI_Datatype atomMbrTypes[3]; // member mpi types
144    
145     MPI_Address(&atomProto.name, &atomDspls[0]);
146     MPI_Address(&atomProto.mass, &atomDspls[1]);
147     MPI_Address(&atomProto.ident, &atomDspls[2]);
148    
149     atomMbrTypes[0] = MPI_CHAR;
150     atomMbrTypes[1] = MPI_DOUBLE;
151     atomMbrTypes[2] = MPI_INT;
152    
153     for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
154    
155     MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
156     MPI_Type_commit(&mpiAtomStructType);
157    
158     // ***********************************************************************
159    
160     if( worldRank == 0 ){
161     #endif
162    
163     // generate the force file name
164    
165 mmeineke 561 strcpy( fileName, "LJFF.frc" );
166     fprintf( stderr,"Trying to open %s\n", fileName );
167 mmeineke 559
168     // attempt to open the file in the current directory first.
169    
170     frcFile = fopen( fileName, "r" );
171    
172     if( frcFile == NULL ){
173    
174     // next see if the force path enviorment variable is set
175    
176     ffPath = getenv( ffPath_env );
177     if( ffPath == NULL ) {
178     STR_DEFINE(ffPath, FRC_PATH );
179     }
180    
181    
182     strcpy( temp, ffPath );
183     strcat( temp, "/" );
184     strcat( temp, fileName );
185     strcpy( fileName, temp );
186    
187     frcFile = fopen( fileName, "r" );
188    
189     if( frcFile == NULL ){
190    
191     sprintf( painCave.errMsg,
192     "Error opening the force field parameter file: %s\n"
193     "Have you tried setting the FORCE_PARAM_PATH environment "
194     "vairable?\n",
195     fileName );
196     painCave.isFatal = 1;
197     simError();
198     }
199     }
200    
201     #ifdef IS_MPI
202     }
203    
204 mmeineke 561 sprintf( checkPointMsg, "LJFF file opened sucessfully." );
205 mmeineke 559 MPIcheckPoint();
206    
207     #endif // is_mpi
208     }
209    
210    
211 mmeineke 561 LJFF::~LJFF(){
212 mmeineke 559
213     if( headAtomType != NULL ) delete headAtomType;
214    
215     #ifdef IS_MPI
216     if( worldRank == 0 ){
217     #endif // is_mpi
218    
219     fclose( frcFile );
220    
221     #ifdef IS_MPI
222     }
223     #endif // is_mpi
224     }
225    
226 mmeineke 561 void LJFF::initForceField( int ljMixRule ){
227 mmeineke 559
228     initFortran( ljMixRule, 0 );
229     }
230    
231 mmeineke 561 void LJFF::cleanMe( void ){
232 mmeineke 559
233     #ifdef IS_MPI
234    
235     // keep the linked list in the mpi version
236    
237     #else // is_mpi
238    
239     // delete the linked list in the single processor version
240    
241     if( headAtomType != NULL ) delete headAtomType;
242    
243     #endif // is_mpi
244     }
245    
246 mmeineke 561 void LJFF::readParams( void ){
247 mmeineke 559
248     atomStruct info;
249     info.last = 1; // initialize last to have the last set.
250     // if things go well, last will be set to 0
251    
252     int i;
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 mmeineke 561 "LJFF atom structures read successfully." );
307 mmeineke 559 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     recieveFrcStruct( &info, mpiAtomStructType );
337    
338     while( !info.last ){
339    
340    
341    
342     headAtomType->add( info );
343    
344     MPIcheckPoint();
345    
346     recieveFrcStruct( &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 chuckv 631 int isEAM = 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 mmeineke 559 &(currentAtomType->epslon),
375     &(currentAtomType->sigma),
376     &dipole,
377     &isError );
378     if( isError ){
379     sprintf( painCave.errMsg,
380     "Error initializing the \"%s\" atom type in fortran\n",
381     currentAtomType->name );
382     painCave.isFatal = 1;
383     simError();
384     }
385     }
386     currentAtomType = currentAtomType->next;
387     }
388    
389     entry_plug->useLJ = 1;
390    
391     #ifdef IS_MPI
392     sprintf( checkPointMsg,
393 mmeineke 561 "LJFF atom structures successfully sent to fortran\n" );
394 mmeineke 559 MPIcheckPoint();
395     #endif // is_mpi
396    
397     }
398    
399    
400 mmeineke 561 void LJFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
401 mmeineke 559
402     int i;
403    
404     // initialize the atoms
405    
406    
407     Atom* thisAtom;
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]->setEpslon( currentAtomType->epslon );
422     the_atoms[i]->setSigma( currentAtomType->sigma );
423     the_atoms[i]->setIdent( currentAtomType->ident );
424     the_atoms[i]->setLJ();
425    
426     if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
427     }
428     }
429    
430 mmeineke 561 void LJFF::initializeBonds( int nBonds, Bond** BondArray,
431 mmeineke 559 bond_pair* the_bonds ){
432    
433     if( nBonds ){
434     sprintf( painCave.errMsg,
435 mmeineke 561 "LJFF does not support bonds.\n" );
436 mmeineke 559 painCave.isFatal = 1;
437     simError();
438     }
439     }
440    
441 mmeineke 561 void LJFF::initializeBends( int nBends, Bend** bendArray,
442 mmeineke 559 bend_set* the_bends ){
443    
444     if( nBends ){
445     sprintf( painCave.errMsg,
446 mmeineke 561 "LJFF does not support bends.\n" );
447 mmeineke 559 painCave.isFatal = 1;
448     simError();
449     }
450     }
451    
452 mmeineke 561 void LJFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
453 mmeineke 559 torsion_set* the_torsions ){
454    
455     if( nTorsions ){
456     sprintf( painCave.errMsg,
457 mmeineke 561 "LJFF does not support torsions.\n" );
458 mmeineke 559 painCave.isFatal = 1;
459     simError();
460     }
461     }
462    
463 mmeineke 561 void LJFF::fastForward( char* stopText, char* searchOwner ){
464 mmeineke 559
465     int foundText = 0;
466     char* the_token;
467    
468     rewind( frcFile );
469     lineNum = 0;
470    
471     eof_test = fgets( readLine, sizeof(readLine), frcFile );
472     lineNum++;
473     if( eof_test == NULL ){
474     sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
475     " file is empty.\n",
476     searchOwner );
477     painCave.isFatal = 1;
478     simError();
479     }
480    
481    
482     while( !foundText ){
483     while( eof_test != NULL && readLine[0] != '#' ){
484     eof_test = fgets( readLine, sizeof(readLine), frcFile );
485     lineNum++;
486     }
487     if( eof_test == NULL ){
488     sprintf( painCave.errMsg,
489     "Error fast forwarding force file for %s at "
490     "line %d: file ended unexpectedly.\n",
491     searchOwner,
492     lineNum );
493     painCave.isFatal = 1;
494     simError();
495     }
496    
497     the_token = strtok( readLine, " ,;\t#\n" );
498     foundText = !strcmp( stopText, the_token );
499    
500     if( !foundText ){
501     eof_test = fgets( readLine, sizeof(readLine), frcFile );
502     lineNum++;
503    
504     if( eof_test == NULL ){
505     sprintf( painCave.errMsg,
506     "Error fast forwarding force file for %s at "
507     "line %d: file ended unexpectedly.\n",
508     searchOwner,
509     lineNum );
510     painCave.isFatal = 1;
511     simError();
512     }
513     }
514     }
515     }
516    
517    
518    
519     int LJ_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
520    
521     char* the_token;
522    
523     the_token = strtok( lineBuffer, " \n\t,;" );
524     if( the_token != NULL ){
525    
526     strcpy( info.name, the_token );
527    
528     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
529     sprintf( painCave.errMsg,
530     "Error parseing AtomTypes: line %d\n", lineNum );
531     painCave.isFatal = 1;
532     simError();
533     }
534    
535     info.mass = atof( the_token );
536    
537     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
538     sprintf( painCave.errMsg,
539     "Error parseing AtomTypes: line %d\n", lineNum );
540     painCave.isFatal = 1;
541     simError();
542     }
543    
544     info.epslon = atof( the_token );
545    
546     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
547     sprintf( painCave.errMsg,
548     "Error parseing AtomTypes: line %d\n", lineNum );
549     painCave.isFatal = 1;
550     simError();
551     }
552    
553     info.sigma = atof( the_token );
554    
555     return 1;
556     }
557     else return 0;
558     }