ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/LJ_FF.cpp
(Generate patch)

Comparing trunk/mdtools/interface_implementation/LJ_FF.cpp (file contents):
Revision 215 by chuckv, Thu Dec 19 21:59:51 2002 UTC vs.
Revision 238 by mmeineke, Fri Jan 17 21:53:36 2003 UTC

# Line 10 | Line 10 | using namespace std;
10   #include "simError.h"
11  
12  
13 < #ifdef IS_MPI
13 > // Declare the structures that will be passed by the parser and  MPI
14  
15 – #include "mpiForceField.h"
16 –
17 –
18 – // Declare the structures that will be passed by MPI
19 –
15   typedef struct{
16    char name[15];
17    double mass;
18    double epslon;
19    double sigma;
20 +  int ident;
21    int last;      //  0  -> default
22 <                 //  1  -> tells nodes to stop listening
22 >                 //  1  -> in MPI: tells nodes to stop listening
23   } atomStruct;
24 +
25 + int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info );
26 +
27 + #ifdef IS_MPI
28 +
29 + #include "mpiForceField.h"
30 +
31   MPI_Datatype mpiAtomStructType;
32  
33   #endif
34  
35  
36 + // declaration of functions needed to wrap the fortran module
37  
38 + extern "C" {
39 +  void forcefactory_( char* forceName,
40 +                      int* status,
41 +                      void (*wrapFunction)( void (*p1)( int* ident,
42 +                                                        double* mass,
43 +                                                        double* epslon,
44 +                                                        double* sigma,
45 +                                                        int* status ),
46 +                                            void (*p2)( void ),
47 +                                            void (*p3)( double* positionArray,
48 +                                                        double* forceArray,
49 +                                                        double* potentialEnergy,
50 +                                                        short int* doPotentialCalc )),
51 +                      int forceNameLength );
52 + }
53 +
54 +
55 + void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon,
56 +                                   double* sigma, int* status ),
57 +                        void (*p2)( void ),
58 +                        void (*p3)( double* positionArray,double* forceArray,
59 +                                    double* potentialEnergy,
60 +                                    short int* doPotentialCalc ) );
61 +
62 + void (*newLJtype)( int* ident, double* mass, double* epslon, double* sigma,
63 +                   int* status );
64 +
65 + void (*initLJfortran)( void );
66 +
67 + LJ_FF* currentLJwrap;
68 +
69 +
70 + //****************************************************************
71 + // begins the actual forcefield stuff.  
72 + //****************************************************************
73 +
74   LJ_FF::LJ_FF(){
75  
76    char fileName[200];
# Line 39 | Line 79 | LJ_FF::LJ_FF(){
79    char temp[200];
80    char errMsg[1000];
81  
82 +  // do the funtion wrapping
83 +  currentLJwrap = this;
84 +  wrapMe();
85 +
86   #ifdef IS_MPI
87    int i;
88    
# Line 46 | Line 90 | LJ_FF::LJ_FF(){
90    // Init the atomStruct mpi type
91  
92    atomStruct atomProto; // mpiPrototype
93 <  int atomBC[3] = {15,3,1};  // block counts
93 >  int atomBC[3] = {15,3,2};  // block counts
94    MPI_Aint atomDspls[3];           // displacements
95    MPI_Datatype atomMbrTypes[3];    // member mpi types
96  
97    MPI_Address(&atomProto.name, &atomDspls[0]);
98    MPI_Address(&atomProto.mass, &atomDspls[1]);
99 <  MPI_Address(&atomProto.last, &atomDspls[2]);
99 >  MPI_Address(&atomProto.ident, &atomDspls[2]);
100    
101    atomMbrTypes[0] = MPI_CHAR;
102    atomMbrTypes[1] = MPI_DOUBLE;
# Line 135 | Line 179 | void LJ_FF::initializeAtoms( void ){
179   #endif // is_mpi
180   }
181  
182 +
183 + void LJ_FF::wrapMe( void ){
184 +  
185 +  char* currentFF = "LJ";
186 +  int isError = 0;
187 +  
188 +  forcefactory_( currentFF, &isError, LJfunctionWrapper, strlen(currentFF) );
189 +
190 +  if( isError ){
191 +    
192 +    sprintf( painCave.errMsg,
193 +             "LJ_FF error: an error was returned from fortran when the "
194 +             "the functions were being wrapped.\n" );
195 +    painCave.isFatal = 1;
196 +    simError();
197 +  }
198 +
199 + #ifdef IS_MPI
200 +  sprintf( checkPointMsg, "LJ_FF functions succesfully wrapped." );
201 +  MPIcheckPoint();
202 + #endif // is_mpi
203 + }
204 +  
205 +
206 + void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon,
207 +                                   double* sigma, int* status ),
208 +                        void (*p2)( void ),
209 +                        void (*p3)( double* positionArray,double* forceArray,
210 +                                    double* potentialEnergy,
211 +                                    short int* doPotentialCalc ) ){
212 +  
213 +  
214 +  p1 = newLJtype;
215 +  p2 = initLJfortran;
216 +  this->setLJfortran( p3 );
217 + }
218 +
219 +
220 +
221   void LJ_FF::initializeAtoms( void ){
222    
223    class LinkedType {
# Line 151 | Line 234 | void LJ_FF::initializeAtoms( void ){
234        return NULL;
235      }
236      
237 < #ifdef IS_MPI
237 >
238      void add( atomStruct &info ){
239 +    
240 +      // check for duplicates
241 +      
242 +      if( !strcmp( info.name, name ) ){
243 +        sprintf( simError.painCave,
244 +                 "Duplicate LJ atom type \"%s\" found in "
245 +                 "the LJ_FF param file./n",
246 +                 name );
247 +        painCave.isFatal = 1;
248 +        simError();
249 +      }
250 +      
251        if( next != NULL ) next->add(info);
252        else{
253          next = new LinkedType();
# Line 160 | Line 255 | void LJ_FF::initializeAtoms( void ){
255          next->mass     = info.mass;
256          next->epslon   = info.epslon;
257          next->sigma    = info.sigma;
258 +        next->ident    = info.ident;
259        }
260      }
261      
262 +
263 + #ifdef IS_MPI
264 +    
265      void duplicate( atomStruct &info ){
266        strcpy(info.name, name);
267        info.mass     = mass;
268        info.epslon   = epslon;
269        info.sigma    = sigma;
270 +      info.ident    = ident;
271        info.last     = 0;
272      }
273  
# Line 178 | Line 278 | void LJ_FF::initializeAtoms( void ){
278      double mass;
279      double epslon;
280      double sigma;
281 +    int ident;
282      LinkedType* next;
283    };
284    
285    LinkedType* headAtomType;
286    LinkedType* currentAtomType;
186 –  LinkedType* tempAtomType;
187 –
188 – #ifdef IS_MPI
287    atomStruct info;
288    info.last = 1; // initialize last to have the last set.
289                   // if things go well, last will be set to 0
192 – #endif
193 –  
290  
195 –  char readLine[500];
196 –  char* the_token;
197 –  char* eof_test;
198 –  int foundAtom = 0;
199 –  int lineNum = 0;
291    int i;
292 <
292 >  int identNum;
293    
203 –  //////////////////////////////////////////////////
204 –  // a quick water fix
205 –
206 –  double waterI[3][3];
207 –  waterI[0][0] = 1.76958347772500;
208 –  waterI[0][1] = 0.0;
209 –  waterI[0][2] = 0.0;
210 –
211 –  waterI[1][0] = 0.0;
212 –  waterI[1][1] = 0.614537057924513;
213 –  waterI[1][2] = 0.0;
214 –
215 –  waterI[2][0] = 0.0;
216 –  waterI[2][1] = 0.0;
217 –  waterI[2][2] = 1.15504641980049;
218 –
219 –
220 –  double headI[3][3];
221 –  headI[0][0] = 1125;
222 –  headI[0][1] = 0.0;
223 –  headI[0][2] = 0.0;
224 –
225 –  headI[1][0] = 0.0;
226 –  headI[1][1] = 1125;
227 –  headI[1][2] = 0.0;
228 –
229 –  headI[2][0] = 0.0;
230 –  headI[2][1] = 0.0;
231 –  headI[2][2] = 250;
232 –
233 –  
234 –
235 –  //////////////////////////////////////////////////
236 –
294    Atom** the_atoms;
295    int nAtoms;
296    the_atoms = entry_plug->atoms;
# Line 245 | Line 302 | void LJ_FF::initializeAtoms( void ){
302   #endif
303      
304      // read in the atom types.
305 <    
249 <    rewind( frcFile );
305 >
306      headAtomType = new LinkedType;
251 –    currentAtomType = headAtomType;
307      
308 <    eof_test = fgets( readLine, sizeof(readLine), frcFile );
309 <    lineNum++;
255 <    if( eof_test == NULL ){
256 <      sprintf( painCave.errMsg, "Error in reading Atoms from force file.\n" );
257 <      painCave.isFatal = 1;
258 <      simError();
259 <    }
260 <    
261 <    
262 <    while( !foundAtom ){
263 <      while( eof_test != NULL && readLine[0] != '#' ){
264 <        eof_test = fgets( readLine, sizeof(readLine), frcFile );
265 <        lineNum++;
266 <      }
267 <      if( eof_test == NULL ){
268 <        sprintf( painCave.errMsg,
269 <                 "Error in reading Atoms from force file at line %d.\n",
270 <                 lineNum );
271 <        painCave.isFatal = 1;
272 <        simError();
273 <      }
274 <      
275 <      the_token = strtok( readLine, " ,;\t#\n" );
276 <      foundAtom = !strcmp( "AtomTypes", the_token );
277 <      
278 <      if( !foundAtom ){
279 <        eof_test = fgets( readLine, sizeof(readLine), frcFile );
280 <        lineNum++;
281 <        
282 <        if( eof_test == NULL ){
283 <          sprintf( painCave.errMsg,
284 <                   "Error in reading Atoms from force file at line %d.\n",
285 <                   lineNum );
286 <          painCave.isFatal = 1;
287 <          simError();
288 <        }
289 <      }
290 <    }
291 <    
308 >    fastFoward( "AtomTypes", "initializeAtoms" );
309 >
310      // we are now at the AtomTypes section.
311      
312      eof_test = fgets( readLine, sizeof(readLine), frcFile );
313      lineNum++;
314      
315 +    
316 +    // read a line, and start parseing out the atom types
317 +
318      if( eof_test == NULL ){
319        sprintf( painCave.errMsg,
320                 "Error in reading Atoms from force file at line %d.\n",
# Line 302 | Line 323 | void LJ_FF::initializeAtoms( void ){
323        simError();
324      }
325      
326 +    identNum = 1;
327 +    // stop reading at end of file, or at next section
328      while( readLine[0] != '#' && eof_test != NULL ){
329 <      
329 >
330 >      // toss comment lines
331        if( readLine[0] != '!' ){
332          
333 <        the_token = strtok( readLine, " \n\t,;" );
334 <        if( the_token != NULL ){
335 <          
336 <          strcpy( currentAtomType->name, the_token );
337 <          
314 <          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
315 <            sprintf( painCave.errMsg,
316 <                     "Error parseing AtomTypes: line %d\n", lineNum );
317 <            painCave.isFatal = 1;
318 <            simError();
319 <          }
320 <          
321 <          sscanf( the_token, "%lf", &currentAtomType->mass );
322 <          
323 <          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
324 <            sprintf( painCave.errMsg,
325 <                     "Error parseing AtomTypes: line %d\n", lineNum );
326 <            painCave.isFatal = 1;
327 <            simError();
328 <          }
329 <          
330 <          sscanf( the_token, "%lf", &currentAtomType->epslon );
331 <          
332 <          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
333 <            sprintf( painCave.errMsg,
334 <                     "Error parseing AtomTypes: line %d\n", lineNum );
335 <            painCave.isFatal = 1;
336 <            simError();
337 <          }
338 <          
339 <          sscanf( the_token, "%lf", &currentAtomType->sigma );
333 >        // the parser returns 0 if the line was blank
334 >        if( parseAtomLJ( readLine, lineNum, info ) ){
335 >          info.ident = identNum;
336 >          headAtomType->add( info );;
337 >          identNum++;
338          }
339        }
342 –      
343 –      tempAtomType = new LinkedType;
344 –      currentAtomType->next = tempAtomType;
345 –      currentAtomType = tempAtomType;
346 –      
340        eof_test = fgets( readLine, sizeof(readLine), frcFile );
341        lineNum++;
342      }
# Line 383 | Line 376 | void LJ_FF::initializeAtoms( void ){
376    }
377   #endif // is_mpi
378  
379 +  // call new A_types in fortran
380    
381 +  int isError;
382 +  currentAtomType = headAtomType;
383 +  while( currentAtomType != NULL ){
384 +    
385 +    if( currentAtomType->name[0] != NULL ){
386 +      isError = 0;
387 +      newLJtype( &(currentAtomType->ident),
388 +                 &(currentAtomType->mass),
389 +                 &(currentAtomType->epslon),
390 +                 &(currentAtomType->sigma),
391 +                 isError );
392 +      if( isError ){
393 +        sprintf( painCave.errMsg,
394 +                 "Error initializing the \"%s\" atom type in fortran\n",
395 +                 currentAtomType->name );
396 +        painCave.isFatal = 1;
397 +        simError();
398 +      }
399 +    }
400 +    currentAtomType = currentAtomType->next;
401 +  }
402 +      
403 + #ifdef IS_MPI
404 +  sprintf( checkPointMsg,
405 +           "LJ_FF atom structures successfully sent to fortran\n" );
406 +  MPIcheckPoint();
407 + #endif // is_mpi
408 +
409 +  
410 +
411    // initialize the atoms
412    
413    Atom* thisAtom;
# Line 402 | Line 426 | void LJ_FF::initializeAtoms( void ){
426      the_atoms[i]->setMass( currentAtomType->mass );
427      the_atoms[i]->setEpslon( currentAtomType->epslon );
428      the_atoms[i]->setSigma( currentAtomType->sigma );
429 +    the_atoms[i]->setIdent( currentAtomType->ident );
430      the_atoms[i]->setLJ();
431    }
432  
# Line 458 | Line 483 | void LJ_FF::initializeTorsions( torsion_set* the_torsi
483   #endif // is_mpi
484  
485   }
486 +
487 +
488 + void LJ_FF::fastForward( char* stopText, char* searchOwner ){
489 +
490 +  int foundText = 0;
491 +  char* the_token;
492 +
493 +  rewind( frcFile );
494 +  lineNum = 0;
495 +
496 +  eof_test = fgets( readLine, sizeof(readLine), frcFile );
497 +  lineNum++;
498 +  if( eof_test == NULL ){
499 +    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
500 +             " file is empty.\n",
501 +             searchOwner );
502 +    painCave.isFatal = 1;
503 +    simError();
504 +  }
505 +  
506 +  
507 +  while( !foundText ){
508 +    while( eof_test != NULL && readLine[0] != '#' ){
509 +      eof_test = fgets( readLine, sizeof(readLine), frcFile );
510 +      lineNum++;
511 +    }
512 +    if( eof_test == NULL ){
513 +      sprintf( painCave.errMsg,
514 +               "Error fast forwarding force file for %s at "
515 +               "line %d: file ended unexpectedly.\n",
516 +               searchOwner,
517 +               lineNum );
518 +      painCave.isFatal = 1;
519 +      simError();
520 +    }
521 +    
522 +    the_token = strtok( readLine, " ,;\t#\n" );
523 +    foundText = !strcmp( stopText, the_token );
524 +    
525 +    if( !foundText ){
526 +      eof_test = fgets( readLine, sizeof(readLine), frcFile );
527 +      lineNum++;
528 +      
529 +      if( eof_test == NULL ){
530 +        sprintf( painCave.errMsg,
531 +                 "Error fast forwarding force file for %s at "
532 +                 "line %d: file ended unexpectedly.\n",
533 +                 searchOwner,
534 +                 lineNum );
535 +        painCave.isFatal = 1;
536 +        simError();
537 +      }
538 +    }
539 +  }  
540 + }
541 +
542 +
543 +
544 + int parseAtomLJ( char *lineBuffer, int lineNum,  atomStruct &info ){
545 +
546 +  char* the_token;
547 +  
548 +  the_token = strtok( lineBuffer, " \n\t,;" );
549 +  if( the_token != NULL ){
550 +    
551 +    strcpy( info.name, the_token );
552 +          
553 +    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
554 +      sprintf( painCave.errMsg,
555 +               "Error parseing AtomTypes: line %d\n", lineNum );
556 +      painCave.isFatal = 1;
557 +      simError();
558 +    }
559 +    
560 +    info.mass = atof( the_token );
561 +    
562 +    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
563 +      sprintf( painCave.errMsg,
564 +               "Error parseing AtomTypes: line %d\n", lineNum );
565 +      painCave.isFatal = 1;
566 +      simError();
567 +    }
568 +        
569 +    info.epslon = atof( the_token );
570 +          
571 +    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
572 +      sprintf( painCave.errMsg,
573 +               "Error parseing AtomTypes: line %d\n", lineNum );
574 +      painCave.isFatal = 1;
575 +      simError();
576 +    }
577 +        
578 +    info.sigma = atof( the_token );
579 +    
580 +    return 1;
581 +  }
582 +  else return 0;
583 + }
584 +
585 +
586 + void LJ_FF::doForces( void ){
587 +
588 +  int i;
589 +  double* frc;
590 +  double* pos;
591 +  double potE;
592 +  short int calcPot = 0;
593 +
594 +  // forces are zeroed here, before any are acumulated.
595 +  // NOTE: do not rezero the forces in Fortran.
596 +
597 +  for(i=0; i<entry_plug->n_atoms; i++){
598 +    entry_plug->atoms[i]->zeroForces();
599 +  }
600 +
601 +  frc = Atom::getFrcArray();
602 +  pos = Atom::getPosArray();
603 +
604 +  doLJfortran( pos, frc, potE, calcPot );
605 + }
606 +  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines