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 253 by chuckv, Thu Jan 30 15:20:21 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)( int *nLocal,
47 +                                                        int *identArray,
48 +                                                        int *isError ),
49 +                                            void (*p3)( double* positionArray,
50 +                                                        double* forceArray,
51 +                                                        double* potentialEnergy,
52 +                                                        short int* doPotentialCalc )),
53 +                      int forceNameLength );
54 + }
55 +
56 +
57 + void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon,
58 +                                   double* sigma, int* status ),
59 +                        void (*p2)( int *nLocal, int *identArray, int *isError ),
60 +                        void (*p3)( double* positionArray,double* forceArray,
61 +                                    double* potentialEnergy,
62 +                                    short int* doPotentialCalc ) );
63 +
64 + void (*newLJtype)( int* ident, double* mass, double* epslon, double* sigma,
65 +                   int* status );
66 +
67 + void (*initLJfortran) ( int *nLocal, int *identArray, int *isError );
68 +
69 + LJ_FF* currentLJwrap;
70 +
71 +
72 + //****************************************************************
73 + // begins the actual forcefield stuff.  
74 + //****************************************************************
75 +
76   LJ_FF::LJ_FF(){
77  
78    char fileName[200];
# Line 39 | Line 81 | LJ_FF::LJ_FF(){
81    char temp[200];
82    char errMsg[1000];
83  
84 +  // do the funtion wrapping
85 +  currentLJwrap = this;
86 +  wrapMe();
87 +
88   #ifdef IS_MPI
89    int i;
90    
# Line 46 | Line 92 | LJ_FF::LJ_FF(){
92    // Init the atomStruct mpi type
93  
94    atomStruct atomProto; // mpiPrototype
95 <  int atomBC[3] = {15,3,1};  // block counts
95 >  int atomBC[3] = {15,3,2};  // block counts
96    MPI_Aint atomDspls[3];           // displacements
97    MPI_Datatype atomMbrTypes[3];    // member mpi types
98  
99    MPI_Address(&atomProto.name, &atomDspls[0]);
100    MPI_Address(&atomProto.mass, &atomDspls[1]);
101 <  MPI_Address(&atomProto.last, &atomDspls[2]);
101 >  MPI_Address(&atomProto.ident, &atomDspls[2]);
102    
103    atomMbrTypes[0] = MPI_CHAR;
104    atomMbrTypes[1] = MPI_DOUBLE;
# Line 135 | Line 181 | void LJ_FF::initializeAtoms( void ){
181   #endif // is_mpi
182   }
183  
184 +
185 + void LJ_FF::wrapMe( void ){
186 +  
187 +  char* currentFF = "LJ";
188 +  int isError = 0;
189 +  
190 +  forcefactory_( currentFF, &isError, LJfunctionWrapper, strlen(currentFF) );
191 +
192 +  if( isError ){
193 +    
194 +    sprintf( painCave.errMsg,
195 +             "LJ_FF error: an error was returned from fortran when the "
196 +             "the functions were being wrapped.\n" );
197 +    painCave.isFatal = 1;
198 +    simError();
199 +  }
200 +
201 + #ifdef IS_MPI
202 +  sprintf( checkPointMsg, "LJ_FF functions succesfully wrapped." );
203 +  MPIcheckPoint();
204 + #endif // is_mpi
205 + }
206 +  
207 +
208 + void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon,
209 +                                   double* sigma, int* status ),
210 +                        void (*p2)( int*, int*, int* ),
211 +                        void (*p3)( double* positionArray,double* forceArray,
212 +                                    double* potentialEnergy,
213 +                                    short int* doPotentialCalc ) ){
214 +  
215 +  
216 +  newLJtype = p1;
217 +  initLJfortran = p2;
218 +  currentLJwrap->setLJfortran( p3 );
219 + }
220 +
221 +
222 +
223   void LJ_FF::initializeAtoms( void ){
224    
225    class LinkedType {
# Line 151 | Line 236 | void LJ_FF::initializeAtoms( void ){
236        return NULL;
237      }
238      
239 < #ifdef IS_MPI
239 >
240      void add( atomStruct &info ){
241 +    
242 +      // check for duplicates
243 +      
244 +      if( !strcmp( info.name, name ) ){
245 +        sprintf( painCave.errMsg,
246 +                 "Duplicate LJ atom type \"%s\" found in "
247 +                 "the LJ_FF param file./n",
248 +                 name );
249 +        painCave.isFatal = 1;
250 +        simError();
251 +      }
252 +      
253        if( next != NULL ) next->add(info);
254        else{
255          next = new LinkedType();
# Line 160 | Line 257 | void LJ_FF::initializeAtoms( void ){
257          next->mass     = info.mass;
258          next->epslon   = info.epslon;
259          next->sigma    = info.sigma;
260 +        next->ident    = info.ident;
261        }
262      }
263      
264 +
265 + #ifdef IS_MPI
266 +    
267      void duplicate( atomStruct &info ){
268        strcpy(info.name, name);
269        info.mass     = mass;
270        info.epslon   = epslon;
271        info.sigma    = sigma;
272 +      info.ident    = ident;
273        info.last     = 0;
274      }
275  
# Line 178 | Line 280 | void LJ_FF::initializeAtoms( void ){
280      double mass;
281      double epslon;
282      double sigma;
283 +    int ident;
284      LinkedType* next;
285    };
286    
287    LinkedType* headAtomType;
288    LinkedType* currentAtomType;
186 –  LinkedType* tempAtomType;
187 –
188 – #ifdef IS_MPI
289    atomStruct info;
290    info.last = 1; // initialize last to have the last set.
291                   // if things go well, last will be set to 0
192 – #endif
193 –  
292  
195 –  char readLine[500];
196 –  char* the_token;
197 –  char* eof_test;
198 –  int foundAtom = 0;
199 –  int lineNum = 0;
293    int i;
294 <
294 >  int identNum;
295    
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 –
296    Atom** the_atoms;
297    int nAtoms;
298    the_atoms = entry_plug->atoms;
# Line 245 | Line 304 | void LJ_FF::initializeAtoms( void ){
304   #endif
305      
306      // read in the atom types.
307 <    
249 <    rewind( frcFile );
307 >
308      headAtomType = new LinkedType;
251 –    currentAtomType = headAtomType;
309      
310 <    eof_test = fgets( readLine, sizeof(readLine), frcFile );
311 <    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 <    
310 >    fastForward( "AtomTypes", "initializeAtoms" );
311 >
312      // we are now at the AtomTypes section.
313      
314      eof_test = fgets( readLine, sizeof(readLine), frcFile );
315      lineNum++;
316      
317 +    
318 +    // read a line, and start parseing out the atom types
319 +
320      if( eof_test == NULL ){
321        sprintf( painCave.errMsg,
322                 "Error in reading Atoms from force file at line %d.\n",
# Line 302 | Line 325 | void LJ_FF::initializeAtoms( void ){
325        simError();
326      }
327      
328 +    identNum = 1;
329 +    // stop reading at end of file, or at next section
330      while( readLine[0] != '#' && eof_test != NULL ){
331 <      
331 >
332 >      // toss comment lines
333        if( readLine[0] != '!' ){
334          
335 <        the_token = strtok( readLine, " \n\t,;" );
336 <        if( the_token != NULL ){
337 <          
338 <          strcpy( currentAtomType->name, the_token );
339 <          
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 );
335 >        // the parser returns 0 if the line was blank
336 >        if( parseAtomLJ( readLine, lineNum, info ) ){
337 >          info.ident = identNum;
338 >          headAtomType->add( info );;
339 >          identNum++;
340          }
341        }
342 –      
343 –      tempAtomType = new LinkedType;
344 –      currentAtomType->next = tempAtomType;
345 –      currentAtomType = tempAtomType;
346 –      
342        eof_test = fgets( readLine, sizeof(readLine), frcFile );
343        lineNum++;
344      }
# Line 383 | Line 378 | void LJ_FF::initializeAtoms( void ){
378    }
379   #endif // is_mpi
380  
381 +  // call new A_types in fortran
382    
383 +  int isError;
384 +  currentAtomType = headAtomType;
385 +  while( currentAtomType != NULL ){
386 +    
387 +    if( currentAtomType->name[0] != '\0' ){
388 +      isError = 0;
389 +          newLJtype( &(currentAtomType->ident),
390 +                 &(currentAtomType->mass),
391 +                 &(currentAtomType->epslon),
392 +                 &(currentAtomType->sigma),
393 +                 &isError );
394 +      if( isError ){
395 +        sprintf( painCave.errMsg,
396 +                 "Error initializing the \"%s\" atom type in fortran\n",
397 +                 currentAtomType->name );
398 +        painCave.isFatal = 1;
399 +        simError();
400 +      }
401 +    }
402 +    currentAtomType = currentAtomType->next;
403 +  }
404 +      
405 + #ifdef IS_MPI
406 +  sprintf( checkPointMsg,
407 +           "LJ_FF atom structures successfully sent to fortran\n" );
408 +  MPIcheckPoint();
409 + #endif // is_mpi
410 +
411 +  
412 +
413    // initialize the atoms
414    
415 +  double bigSigma = 0.0;
416    Atom* thisAtom;
417  
418    for( i=0; i<nAtoms; i++ ){
# Line 402 | Line 429 | void LJ_FF::initializeAtoms( void ){
429      the_atoms[i]->setMass( currentAtomType->mass );
430      the_atoms[i]->setEpslon( currentAtomType->epslon );
431      the_atoms[i]->setSigma( currentAtomType->sigma );
432 +    the_atoms[i]->setIdent( currentAtomType->ident );
433      the_atoms[i]->setLJ();
434 +
435 +    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
436    }
437  
438 +  
439 + #ifdef IS_MPI
440 +  double tempBig = bigSigma;
441 +  MPI::COMM_WORLD::Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
442 + #endif  //is_mpi
443  
444 +  //calc rCut and rList
445 +
446 +  entry_plug->rCut = 2.5 * bigSigma;
447 +  if(entry_plug->rCut > (entry_plug->box_x / 2.0)) entry_plug->rCut = entry_plug->box_x / 2.0;
448 +  if(entry_plug->rCut > (entry_plug->box_y / 2.0)) entry_plug->rCut = entry_plug->box_y / 2.0;
449 +  if(entry_plug->rCut > (entry_plug->box_z / 2.0)) entry_plug->rCut = entry_plug->box_z / 2.0;
450 +
451 +  entry_plug->rList = entry_plug->rCut + 1.0;
452 +
453    // clean up the memory
454    
455    delete headAtomType;
# Line 415 | Line 459 | void LJ_FF::initializeAtoms( void ){
459    MPIcheckPoint();
460   #endif // is_mpi
461  
462 +  initFortran();
463 +  entry_plug->refreshSim();
464 +
465   }
466  
467   void LJ_FF::initializeBonds( bond_pair* the_bonds ){
# Line 458 | Line 505 | void LJ_FF::initializeTorsions( torsion_set* the_torsi
505   #endif // is_mpi
506  
507   }
508 +
509 +
510 + void LJ_FF::fastForward( char* stopText, char* searchOwner ){
511 +
512 +  int foundText = 0;
513 +  char* the_token;
514 +
515 +  rewind( frcFile );
516 +  lineNum = 0;
517 +
518 +  eof_test = fgets( readLine, sizeof(readLine), frcFile );
519 +  lineNum++;
520 +  if( eof_test == NULL ){
521 +    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
522 +             " file is empty.\n",
523 +             searchOwner );
524 +    painCave.isFatal = 1;
525 +    simError();
526 +  }
527 +  
528 +  
529 +  while( !foundText ){
530 +    while( eof_test != NULL && readLine[0] != '#' ){
531 +      eof_test = fgets( readLine, sizeof(readLine), frcFile );
532 +      lineNum++;
533 +    }
534 +    if( eof_test == NULL ){
535 +      sprintf( painCave.errMsg,
536 +               "Error fast forwarding force file for %s at "
537 +               "line %d: file ended unexpectedly.\n",
538 +               searchOwner,
539 +               lineNum );
540 +      painCave.isFatal = 1;
541 +      simError();
542 +    }
543 +    
544 +    the_token = strtok( readLine, " ,;\t#\n" );
545 +    foundText = !strcmp( stopText, the_token );
546 +    
547 +    if( !foundText ){
548 +      eof_test = fgets( readLine, sizeof(readLine), frcFile );
549 +      lineNum++;
550 +      
551 +      if( eof_test == NULL ){
552 +        sprintf( painCave.errMsg,
553 +                 "Error fast forwarding force file for %s at "
554 +                 "line %d: file ended unexpectedly.\n",
555 +                 searchOwner,
556 +                 lineNum );
557 +        painCave.isFatal = 1;
558 +        simError();
559 +      }
560 +    }
561 +  }  
562 + }
563 +
564 +
565 +
566 + int parseAtomLJ( char *lineBuffer, int lineNum,  atomStruct &info ){
567 +
568 +  char* the_token;
569 +  
570 +  the_token = strtok( lineBuffer, " \n\t,;" );
571 +  if( the_token != NULL ){
572 +    
573 +    strcpy( info.name, the_token );
574 +          
575 +    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
576 +      sprintf( painCave.errMsg,
577 +               "Error parseing AtomTypes: line %d\n", lineNum );
578 +      painCave.isFatal = 1;
579 +      simError();
580 +    }
581 +    
582 +    info.mass = atof( the_token );
583 +    
584 +    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
585 +      sprintf( painCave.errMsg,
586 +               "Error parseing AtomTypes: line %d\n", lineNum );
587 +      painCave.isFatal = 1;
588 +      simError();
589 +    }
590 +        
591 +    info.epslon = atof( the_token );
592 +          
593 +    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
594 +      sprintf( painCave.errMsg,
595 +               "Error parseing AtomTypes: line %d\n", lineNum );
596 +      painCave.isFatal = 1;
597 +      simError();
598 +    }
599 +        
600 +    info.sigma = atof( the_token );
601 +    
602 +    return 1;
603 +  }
604 +  else return 0;
605 + }
606 +
607 +
608 + void LJ_FF::doForces( int calcPot ){
609 +
610 +  int i;
611 +  double* frc;
612 +  double* pos;
613 +  short int passedCalcPot = (short int)calcPot;
614 +
615 +  // forces are zeroed here, before any are acumulated.
616 +  // NOTE: do not rezero the forces in Fortran.
617 +
618 +  for(i=0; i<entry_plug->n_atoms; i++){
619 +    entry_plug->atoms[i]->zeroForces();
620 +  }
621 +
622 +  frc = Atom::getFrcArray();
623 +  pos = Atom::getPosArray();
624 +
625 + //   entry_plug->lrPot = -1;
626 +  doLJfortran( pos, frc, &(entry_plug->lrPot), &passedCalcPot );
627 +
628 +
629 + //  fprintf( stderr,
630 + //         "lrPot =  %lf\n", entry_plug->lrPot );
631 +  
632 + }
633 +  
634 + void LJ_FF::initFortran( void ){
635 +  
636 +  int nLocal = entry_plug->n_atoms;
637 +  int *ident;
638 +  int isError;
639 +  int i;
640 +
641 +  ident = new int[nLocal];
642 +
643 +  for(i=0; i<nLocal; i++){
644 +    ident[i] = entry_plug->atoms[i]->getIdent();
645 +  }
646 +
647 +  isError = 0;
648 +  initLJfortran( &nLocal, ident, &isError );
649 +  
650 +  if(isError){
651 +    sprintf( painCave.errMsg,
652 +             "LJ_FF error: There was an error initializing the component list in fortran.\n" );
653 +    painCave.isFatal = 1;
654 +    simError();
655 +  }
656 +
657 +  
658 + #ifdef IS_MPI
659 +  sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" );
660 +  MPIcheckPoint();
661 + #endif // is_mpi
662 +  
663 +  delete[] ident;
664 +
665 + }
666 +  

Diff Legend

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