ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/UseTheForce/LJFF.cpp
Revision: 1638
Committed: Fri Oct 22 23:00:06 2004 UTC (19 years, 8 months ago) by chrisfen
Original Path: trunk/OOPSE-4/src/UseTheForce/LJFF.cpp
File size: 12512 byte(s)
Log Message:
Dear god!  It runs and conserves energy!

File Contents

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