ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/LJ_FF.cpp
Revision: 433
Committed: Fri Mar 28 15:28:53 2003 UTC (21 years, 3 months ago) by mmeineke
File size: 12341 byte(s)
Log Message:
fixed long range interactions in Trappe

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     #include <mpi++.h>
11     #endif //is_mpi
12    
13     #include "ForceFields.hpp"
14     #include "SRI.hpp"
15     #include "simError.h"
16    
17     #include "fortranWrappers.hpp"
18    
19     #ifdef IS_MPI
20     #include "mpiForceField.h"
21     #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 mmeineke 420
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 LJ_FF param file./n",
70     name );
71     painCave.isFatal = 1;
72     simError();
73     }
74    
75     if( next != NULL ) next->add(info);
76     else{
77 mmeineke 424 next = new LinkedAtomType();
78 mmeineke 420 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 mmeineke 377 }
113    
114     using namespace LJ_NS;
115    
116     //****************************************************************
117     // begins the actual forcefield stuff.
118     //****************************************************************
119    
120    
121     LJ_FF::LJ_FF(){
122    
123     char fileName[200];
124     char* ffPath_env = "FORCE_PARAM_PATH";
125     char* ffPath;
126     char temp[200];
127     char errMsg[1000];
128    
129 mmeineke 420 headAtomType = NULL;
130     currentAtomType = NULL;
131    
132 mmeineke 377 // do the funtion wrapping
133     wrapMeFF( this );
134    
135     #ifdef IS_MPI
136     int i;
137    
138     // **********************************************************************
139     // Init the atomStruct mpi type
140    
141     atomStruct atomProto; // mpiPrototype
142     int atomBC[3] = {15,3,2}; // block counts
143     MPI_Aint atomDspls[3]; // displacements
144     MPI_Datatype atomMbrTypes[3]; // member mpi types
145    
146     MPI_Address(&atomProto.name, &atomDspls[0]);
147     MPI_Address(&atomProto.mass, &atomDspls[1]);
148     MPI_Address(&atomProto.ident, &atomDspls[2]);
149    
150     atomMbrTypes[0] = MPI_CHAR;
151     atomMbrTypes[1] = MPI_DOUBLE;
152     atomMbrTypes[2] = MPI_INT;
153    
154     for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
155    
156     MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
157     MPI_Type_commit(&mpiAtomStructType);
158    
159     // ***********************************************************************
160    
161     if( worldRank == 0 ){
162     #endif
163    
164     // generate the force file name
165    
166     strcpy( fileName, "LJ_FF.frc" );
167     // fprintf( stderr,"Trying to open %s\n", fileName );
168    
169     // attempt to open the file in the current directory first.
170    
171     frcFile = fopen( fileName, "r" );
172    
173     if( frcFile == NULL ){
174    
175     // next see if the force path enviorment variable is set
176    
177     ffPath = getenv( ffPath_env );
178     if( ffPath == NULL ) {
179     STR_DEFINE(ffPath, FRC_PATH );
180     }
181    
182    
183     strcpy( temp, ffPath );
184     strcat( temp, "/" );
185     strcat( temp, fileName );
186     strcpy( fileName, temp );
187    
188     frcFile = fopen( fileName, "r" );
189    
190     if( frcFile == NULL ){
191    
192     sprintf( painCave.errMsg,
193     "Error opening the force field parameter file: %s\n"
194     "Have you tried setting the FORCE_PARAM_PATH environment "
195     "vairable?\n",
196     fileName );
197     painCave.isFatal = 1;
198     simError();
199     }
200     }
201    
202     #ifdef IS_MPI
203     }
204    
205     sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
206     MPIcheckPoint();
207    
208     #endif // is_mpi
209     }
210    
211    
212     LJ_FF::~LJ_FF(){
213    
214 mmeineke 420 if( headAtomType != NULL ) delete headAtomType;
215    
216 mmeineke 377 #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 LJ_FF::initForceField( int ljMixRule ){
228    
229     initFortran( ljMixRule, 0 );
230     }
231    
232 mmeineke 420 void LJ_FF::cleanMe( void ){
233 mmeineke 377
234 mmeineke 420 #ifdef IS_MPI
235 mmeineke 377
236 mmeineke 420 // keep the linked list in the mpi version
237 mmeineke 377
238 mmeineke 420 #else // is_mpi
239 mmeineke 377
240 mmeineke 420 // delete the linked list in the single processor version
241 mmeineke 377
242 mmeineke 420 if( headAtomType != NULL ) delete headAtomType;
243 mmeineke 377
244 mmeineke 420 #endif // is_mpi
245     }
246 mmeineke 377
247 mmeineke 420 void LJ_FF::readParams( void ){
248 mmeineke 377
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 i;
254     int identNum;
255    
256 mmeineke 420
257     bigSigma = 0.0;
258 mmeineke 377 #ifdef IS_MPI
259     if( worldRank == 0 ){
260     #endif
261    
262     // read in the atom types.
263    
264 mmeineke 420 headAtomType = new LinkedAtomType;
265 mmeineke 377
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     "LJ_FF 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 mmeineke 428 headAtomType = new LinkedAtomType;
337 mmeineke 377 recieveFrcStruct( &info, mpiAtomStructType );
338    
339     while( !info.last ){
340    
341    
342    
343     headAtomType->add( info );
344    
345     MPIcheckPoint();
346    
347     recieveFrcStruct( &info, mpiAtomStructType );
348     }
349     }
350     #endif // is_mpi
351    
352     // call new A_types in fortran
353    
354     int isError;
355    
356     // dummy variables
357     int isLJ = 1;
358     int isDipole = 0;
359     int isSSD = 0;
360     int isGB = 0;
361     double w0 = 0.0;
362     double v0 = 0.0;
363     double dipole = 0.0;
364     double GB_dummy = 0.0;
365    
366    
367     currentAtomType = headAtomType;
368     while( currentAtomType != NULL ){
369    
370     if( currentAtomType->name[0] != '\0' ){
371     isError = 0;
372     makeAtype( &(currentAtomType->ident),
373     &isLJ,
374     &isSSD,
375     &isDipole,
376     &isGB,
377     &(currentAtomType->epslon),
378     &(currentAtomType->sigma),
379     &dipole,
380     &w0,
381     &v0,
382     &GB_dummy,
383     &GB_dummy,
384     &GB_dummy,
385     &GB_dummy,
386     &GB_dummy,
387     &GB_dummy,
388     &isError );
389     if( isError ){
390     sprintf( painCave.errMsg,
391     "Error initializing the \"%s\" atom type in fortran\n",
392     currentAtomType->name );
393     painCave.isFatal = 1;
394     simError();
395     }
396     }
397     currentAtomType = currentAtomType->next;
398     }
399    
400 mmeineke 433 entry_plug->useLJ = 1;
401    
402 mmeineke 377 #ifdef IS_MPI
403     sprintf( checkPointMsg,
404     "LJ_FF atom structures successfully sent to fortran\n" );
405     MPIcheckPoint();
406     #endif // is_mpi
407    
408 mmeineke 420 }
409    
410    
411     void LJ_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
412 mmeineke 377
413 mmeineke 420 int i;
414 mmeineke 377
415     // initialize the atoms
416    
417 mmeineke 420
418 mmeineke 377 Atom* thisAtom;
419    
420     for( i=0; i<nAtoms; i++ ){
421    
422     currentAtomType = headAtomType->find( the_atoms[i]->getType() );
423     if( currentAtomType == NULL ){
424     sprintf( painCave.errMsg,
425     "AtomType error, %s not found in force file.\n",
426     the_atoms[i]->getType() );
427     painCave.isFatal = 1;
428     simError();
429     }
430    
431     the_atoms[i]->setMass( currentAtomType->mass );
432     the_atoms[i]->setEpslon( currentAtomType->epslon );
433     the_atoms[i]->setSigma( currentAtomType->sigma );
434     the_atoms[i]->setIdent( currentAtomType->ident );
435     the_atoms[i]->setLJ();
436    
437     if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
438     }
439    
440    
441 mmeineke 433
442 mmeineke 377 #ifdef IS_MPI
443     sprintf( checkPointMsg, "LJ_FF atoms initialized succesfully" );
444     MPIcheckPoint();
445     #endif // is_mpi
446    
447     }
448    
449 mmeineke 420 void LJ_FF::initializeBonds( int nBonds, Bond** BondArray,
450     bond_pair* the_bonds ){
451 mmeineke 377
452 mmeineke 420 if( nBonds ){
453 mmeineke 377 sprintf( painCave.errMsg,
454     "LJ_FF does not support bonds.\n" );
455     painCave.isFatal = 1;
456     simError();
457     }
458     #ifdef IS_MPI
459     MPIcheckPoint();
460     #endif // is_mpi
461    
462     }
463    
464 mmeineke 420 void LJ_FF::initializeBends( int nBends, Bend** bendArray,
465     bend_set* the_bends ){
466 mmeineke 377
467 mmeineke 420 if( nBends ){
468 mmeineke 377 sprintf( painCave.errMsg,
469     "LJ_FF does not support bends.\n" );
470     painCave.isFatal = 1;
471     simError();
472     }
473     #ifdef IS_MPI
474     MPIcheckPoint();
475     #endif // is_mpi
476    
477     }
478    
479 mmeineke 420 void LJ_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
480     torsion_set* the_torsions ){
481 mmeineke 377
482 mmeineke 420 if( nTorsions ){
483 mmeineke 377 sprintf( painCave.errMsg,
484     "LJ_FF does not support torsions.\n" );
485     painCave.isFatal = 1;
486     simError();
487     }
488     #ifdef IS_MPI
489     MPIcheckPoint();
490     #endif // is_mpi
491    
492     }
493    
494     void LJ_FF::fastForward( char* stopText, char* searchOwner ){
495    
496     int foundText = 0;
497     char* the_token;
498    
499     rewind( frcFile );
500     lineNum = 0;
501    
502     eof_test = fgets( readLine, sizeof(readLine), frcFile );
503     lineNum++;
504     if( eof_test == NULL ){
505     sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
506     " file is empty.\n",
507     searchOwner );
508     painCave.isFatal = 1;
509     simError();
510     }
511    
512    
513     while( !foundText ){
514     while( eof_test != NULL && readLine[0] != '#' ){
515     eof_test = fgets( readLine, sizeof(readLine), frcFile );
516     lineNum++;
517     }
518     if( eof_test == NULL ){
519     sprintf( painCave.errMsg,
520     "Error fast forwarding force file for %s at "
521     "line %d: file ended unexpectedly.\n",
522     searchOwner,
523     lineNum );
524     painCave.isFatal = 1;
525     simError();
526     }
527    
528     the_token = strtok( readLine, " ,;\t#\n" );
529     foundText = !strcmp( stopText, the_token );
530    
531     if( !foundText ){
532     eof_test = fgets( readLine, sizeof(readLine), frcFile );
533     lineNum++;
534    
535     if( eof_test == NULL ){
536     sprintf( painCave.errMsg,
537     "Error fast forwarding force file for %s at "
538     "line %d: file ended unexpectedly.\n",
539     searchOwner,
540     lineNum );
541     painCave.isFatal = 1;
542     simError();
543     }
544     }
545     }
546     }
547    
548    
549    
550     int LJ_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
551    
552     char* the_token;
553    
554     the_token = strtok( lineBuffer, " \n\t,;" );
555     if( the_token != NULL ){
556    
557     strcpy( info.name, the_token );
558    
559     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
560     sprintf( painCave.errMsg,
561     "Error parseing AtomTypes: line %d\n", lineNum );
562     painCave.isFatal = 1;
563     simError();
564     }
565    
566     info.mass = atof( the_token );
567    
568     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
569     sprintf( painCave.errMsg,
570     "Error parseing AtomTypes: line %d\n", lineNum );
571     painCave.isFatal = 1;
572     simError();
573     }
574    
575     info.epslon = atof( the_token );
576    
577     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
578     sprintf( painCave.errMsg,
579     "Error parseing AtomTypes: line %d\n", lineNum );
580     painCave.isFatal = 1;
581     simError();
582     }
583    
584     info.sigma = atof( the_token );
585    
586     return 1;
587     }
588     else return 0;
589     }