ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/LJ_FF.cpp
Revision: 447
Committed: Thu Apr 3 20:21:54 2003 UTC (21 years, 3 months ago) by mmeineke
File size: 12046 byte(s)
Log Message:
fixed some small things with simError.h

File Contents

# User Rev Content
1 mmeineke 377 #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 mmeineke 420
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 LJ_FF param file./n",
69     name );
70     painCave.isFatal = 1;
71     simError();
72     }
73    
74     if( next != NULL ) next->add(info);
75     else{
76 mmeineke 424 next = new LinkedAtomType();
77 mmeineke 420 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 mmeineke 377 }
112    
113     using namespace LJ_NS;
114    
115     //****************************************************************
116     // begins the actual forcefield stuff.
117     //****************************************************************
118    
119    
120     LJ_FF::LJ_FF(){
121    
122     char fileName[200];
123     char* ffPath_env = "FORCE_PARAM_PATH";
124     char* ffPath;
125     char temp[200];
126     char errMsg[1000];
127    
128 mmeineke 420 headAtomType = NULL;
129     currentAtomType = NULL;
130    
131 mmeineke 377 // 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     strcpy( fileName, "LJ_FF.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: %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     sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
205     MPIcheckPoint();
206    
207     #endif // is_mpi
208     }
209    
210    
211     LJ_FF::~LJ_FF(){
212    
213 mmeineke 420 if( headAtomType != NULL ) delete headAtomType;
214    
215 mmeineke 377 #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     void LJ_FF::initForceField( int ljMixRule ){
227    
228     initFortran( ljMixRule, 0 );
229     }
230    
231 mmeineke 420 void LJ_FF::cleanMe( void ){
232 mmeineke 377
233 mmeineke 420 #ifdef IS_MPI
234 mmeineke 377
235 mmeineke 420 // keep the linked list in the mpi version
236 mmeineke 377
237 mmeineke 420 #else // is_mpi
238 mmeineke 377
239 mmeineke 420 // delete the linked list in the single processor version
240 mmeineke 377
241 mmeineke 420 if( headAtomType != NULL ) delete headAtomType;
242 mmeineke 377
243 mmeineke 420 #endif // is_mpi
244     }
245 mmeineke 377
246 mmeineke 420 void LJ_FF::readParams( void ){
247 mmeineke 377
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 mmeineke 420
256     bigSigma = 0.0;
257 mmeineke 377 #ifdef IS_MPI
258     if( worldRank == 0 ){
259     #endif
260    
261     // read in the atom types.
262    
263 mmeineke 420 headAtomType = new LinkedAtomType;
264 mmeineke 377
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     "LJ_FF 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 mmeineke 428 headAtomType = new LinkedAtomType;
336 mmeineke 377 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     double w0 = 0.0;
361     double v0 = 0.0;
362     double dipole = 0.0;
363     double GB_dummy = 0.0;
364    
365    
366     currentAtomType = headAtomType;
367     while( currentAtomType != NULL ){
368    
369     if( currentAtomType->name[0] != '\0' ){
370     isError = 0;
371     makeAtype( &(currentAtomType->ident),
372     &isLJ,
373     &isSSD,
374     &isDipole,
375     &isGB,
376     &(currentAtomType->epslon),
377     &(currentAtomType->sigma),
378     &dipole,
379     &w0,
380     &v0,
381     &GB_dummy,
382     &GB_dummy,
383     &GB_dummy,
384     &GB_dummy,
385     &GB_dummy,
386     &GB_dummy,
387     &isError );
388     if( isError ){
389     sprintf( painCave.errMsg,
390     "Error initializing the \"%s\" atom type in fortran\n",
391     currentAtomType->name );
392     painCave.isFatal = 1;
393     simError();
394     }
395     }
396     currentAtomType = currentAtomType->next;
397     }
398    
399 mmeineke 433 entry_plug->useLJ = 1;
400    
401 mmeineke 377 #ifdef IS_MPI
402     sprintf( checkPointMsg,
403     "LJ_FF atom structures successfully sent to fortran\n" );
404     MPIcheckPoint();
405     #endif // is_mpi
406    
407 mmeineke 420 }
408    
409    
410     void LJ_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
411 mmeineke 377
412 mmeineke 420 int i;
413 mmeineke 377
414     // initialize the atoms
415    
416 mmeineke 420
417 mmeineke 377 Atom* thisAtom;
418    
419     for( i=0; i<nAtoms; i++ ){
420    
421     currentAtomType = headAtomType->find( the_atoms[i]->getType() );
422     if( currentAtomType == NULL ){
423     sprintf( painCave.errMsg,
424     "AtomType error, %s not found in force file.\n",
425     the_atoms[i]->getType() );
426     painCave.isFatal = 1;
427     simError();
428     }
429    
430     the_atoms[i]->setMass( currentAtomType->mass );
431     the_atoms[i]->setEpslon( currentAtomType->epslon );
432     the_atoms[i]->setSigma( currentAtomType->sigma );
433     the_atoms[i]->setIdent( currentAtomType->ident );
434     the_atoms[i]->setLJ();
435    
436     if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
437     }
438     }
439    
440 mmeineke 420 void LJ_FF::initializeBonds( int nBonds, Bond** BondArray,
441     bond_pair* the_bonds ){
442 mmeineke 377
443 mmeineke 420 if( nBonds ){
444 mmeineke 377 sprintf( painCave.errMsg,
445     "LJ_FF does not support bonds.\n" );
446     painCave.isFatal = 1;
447     simError();
448     }
449     }
450    
451 mmeineke 420 void LJ_FF::initializeBends( int nBends, Bend** bendArray,
452     bend_set* the_bends ){
453 mmeineke 377
454 mmeineke 420 if( nBends ){
455 mmeineke 377 sprintf( painCave.errMsg,
456     "LJ_FF does not support bends.\n" );
457     painCave.isFatal = 1;
458     simError();
459     }
460     }
461    
462 mmeineke 420 void LJ_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
463     torsion_set* the_torsions ){
464 mmeineke 377
465 mmeineke 420 if( nTorsions ){
466 mmeineke 377 sprintf( painCave.errMsg,
467     "LJ_FF does not support torsions.\n" );
468     painCave.isFatal = 1;
469     simError();
470     }
471     }
472    
473     void LJ_FF::fastForward( char* stopText, char* searchOwner ){
474    
475     int foundText = 0;
476     char* the_token;
477    
478     rewind( frcFile );
479     lineNum = 0;
480    
481     eof_test = fgets( readLine, sizeof(readLine), frcFile );
482     lineNum++;
483     if( eof_test == NULL ){
484     sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
485     " file is empty.\n",
486     searchOwner );
487     painCave.isFatal = 1;
488     simError();
489     }
490    
491    
492     while( !foundText ){
493     while( eof_test != NULL && readLine[0] != '#' ){
494     eof_test = fgets( readLine, sizeof(readLine), frcFile );
495     lineNum++;
496     }
497     if( eof_test == NULL ){
498     sprintf( painCave.errMsg,
499     "Error fast forwarding force file for %s at "
500     "line %d: file ended unexpectedly.\n",
501     searchOwner,
502     lineNum );
503     painCave.isFatal = 1;
504     simError();
505     }
506    
507     the_token = strtok( readLine, " ,;\t#\n" );
508     foundText = !strcmp( stopText, the_token );
509    
510     if( !foundText ){
511     eof_test = fgets( readLine, sizeof(readLine), frcFile );
512     lineNum++;
513    
514     if( eof_test == NULL ){
515     sprintf( painCave.errMsg,
516     "Error fast forwarding force file for %s at "
517     "line %d: file ended unexpectedly.\n",
518     searchOwner,
519     lineNum );
520     painCave.isFatal = 1;
521     simError();
522     }
523     }
524     }
525     }
526    
527    
528    
529     int LJ_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
530    
531     char* the_token;
532    
533     the_token = strtok( lineBuffer, " \n\t,;" );
534     if( the_token != NULL ){
535    
536     strcpy( info.name, the_token );
537    
538     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
539     sprintf( painCave.errMsg,
540     "Error parseing AtomTypes: line %d\n", lineNum );
541     painCave.isFatal = 1;
542     simError();
543     }
544    
545     info.mass = atof( the_token );
546    
547     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
548     sprintf( painCave.errMsg,
549     "Error parseing AtomTypes: line %d\n", lineNum );
550     painCave.isFatal = 1;
551     simError();
552     }
553    
554     info.epslon = atof( the_token );
555    
556     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
557     sprintf( painCave.errMsg,
558     "Error parseing AtomTypes: line %d\n", lineNum );
559     painCave.isFatal = 1;
560     simError();
561     }
562    
563     info.sigma = atof( the_token );
564    
565     return 1;
566     }
567     else return 0;
568     }