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 263 by chuckv, Tue Feb 4 20:15:48 2003 UTC

# Line 9 | Line 9 | using namespace std;
9   #include "SRI.hpp"
10   #include "simError.h"
11  
12 + // Declare the structures that will be passed by the parser and  MPI
13  
13 – #ifdef IS_MPI
14 –
15 – #include "mpiForceField.h"
16 –
17 –
18 – // Declare the structures that will be passed by MPI
19 –
14   typedef struct{
15    char name[15];
16    double mass;
17    double epslon;
18    double sigma;
19 +  int ident;
20    int last;      //  0  -> default
21 <                 //  1  -> tells nodes to stop listening
21 >                 //  1  -> in MPI: tells nodes to stop listening
22   } atomStruct;
23 +
24 + int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info );
25 +
26 + #ifdef IS_MPI
27 + #include "mpiForceField.h"
28 +
29   MPI_Datatype mpiAtomStructType;
30  
31   #endif
32  
33  
34 + // declaration of functions needed to wrap the fortran module
35  
36 + extern "C" {
37 +  void forcefactory_( char* forceName,
38 +                      int* status,
39 +                      void (*wrapFunction)( void (*p1)( int* ident,
40 +                                                        double* mass,
41 +                                                        double* epslon,
42 +                                                        double* sigma,
43 +                                                        int* status ),
44 +                                            void (*p2)( int *nLocal,
45 +                                                        int *identArray,
46 +                                                        int *isError ),
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)( int *nLocal, int *identArray, int *isError ),
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) ( int *nLocal, int *identArray, int *isError );
66 +
67 + LJ_FF* currentLJwrap;
68 +
69 +
70 + //****************************************************************
71 + // begins the actual forcefield stuff.  
72 + //****************************************************************
73 +
74 +
75   LJ_FF::LJ_FF(){
76  
77    char fileName[200];
# Line 39 | Line 80 | LJ_FF::LJ_FF(){
80    char temp[200];
81    char errMsg[1000];
82  
83 +  // do the funtion wrapping
84 +  currentLJwrap = this;
85 +  wrapMe();
86 +
87   #ifdef IS_MPI
88    int i;
89    
# Line 46 | Line 91 | LJ_FF::LJ_FF(){
91    // Init the atomStruct mpi type
92  
93    atomStruct atomProto; // mpiPrototype
94 <  int atomBC[3] = {15,3,1};  // block counts
94 >  int atomBC[3] = {15,3,2};  // block counts
95    MPI_Aint atomDspls[3];           // displacements
96    MPI_Datatype atomMbrTypes[3];    // member mpi types
97  
98    MPI_Address(&atomProto.name, &atomDspls[0]);
99    MPI_Address(&atomProto.mass, &atomDspls[1]);
100 <  MPI_Address(&atomProto.last, &atomDspls[2]);
100 >  MPI_Address(&atomProto.ident, &atomDspls[2]);
101    
102    atomMbrTypes[0] = MPI_CHAR;
103    atomMbrTypes[1] = MPI_DOUBLE;
# Line 135 | Line 180 | void LJ_FF::initializeAtoms( void ){
180   #endif // is_mpi
181   }
182  
183 +
184 + void LJ_FF::wrapMe( void ){
185 +  
186 +  char* currentFF = "LJ";
187 +  int isError = 0;
188 +  
189 +  forcefactory_( currentFF, &isError, LJfunctionWrapper, strlen(currentFF) );
190 +
191 +  if( isError ){
192 +    
193 +    sprintf( painCave.errMsg,
194 +             "LJ_FF error: an error was returned from fortran when the "
195 +             "the functions were being wrapped.\n" );
196 +    painCave.isFatal = 1;
197 +    simError();
198 +  }
199 +
200 + #ifdef IS_MPI
201 +  sprintf( checkPointMsg, "LJ_FF functions succesfully wrapped." );
202 +  MPIcheckPoint();
203 + #endif // is_mpi
204 + }
205 +  
206 +
207 + void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon,
208 +                                   double* sigma, int* status ),
209 +                        void (*p2)( int*, int*, int* ),
210 +                        void (*p3)( double* positionArray,double* forceArray,
211 +                                    double* potentialEnergy,
212 +                                    short int* doPotentialCalc ) ){
213 +  
214 +  
215 +  newLJtype = p1;
216 +  initLJfortran = p2;
217 +  currentLJwrap->setLJfortran( p3 );
218 + }
219 +
220 +
221 +
222   void LJ_FF::initializeAtoms( void ){
223    
224    class LinkedType {
# Line 151 | Line 235 | void LJ_FF::initializeAtoms( void ){
235        return NULL;
236      }
237      
238 < #ifdef IS_MPI
238 >
239      void add( atomStruct &info ){
240 +    
241 +      // check for duplicates
242 +      
243 +      if( !strcmp( info.name, name ) ){
244 +        sprintf( painCave.errMsg,
245 +                 "Duplicate LJ atom type \"%s\" found in "
246 +                 "the LJ_FF param file./n",
247 +                 name );
248 +        painCave.isFatal = 1;
249 +        simError();
250 +      }
251 +      
252        if( next != NULL ) next->add(info);
253        else{
254          next = new LinkedType();
# Line 160 | Line 256 | void LJ_FF::initializeAtoms( void ){
256          next->mass     = info.mass;
257          next->epslon   = info.epslon;
258          next->sigma    = info.sigma;
259 +        next->ident    = info.ident;
260        }
261      }
262      
263 +
264 + #ifdef IS_MPI
265 +    
266      void duplicate( atomStruct &info ){
267        strcpy(info.name, name);
268        info.mass     = mass;
269        info.epslon   = epslon;
270        info.sigma    = sigma;
271 +      info.ident    = ident;
272        info.last     = 0;
273      }
274  
# Line 178 | Line 279 | void LJ_FF::initializeAtoms( void ){
279      double mass;
280      double epslon;
281      double sigma;
282 +    int ident;
283      LinkedType* next;
284    };
285    
286    LinkedType* headAtomType;
287    LinkedType* currentAtomType;
186 –  LinkedType* tempAtomType;
187 –
188 – #ifdef IS_MPI
288    atomStruct info;
289    info.last = 1; // initialize last to have the last set.
290                   // if things go well, last will be set to 0
192 – #endif
193 –  
291  
195 –  char readLine[500];
196 –  char* the_token;
197 –  char* eof_test;
198 –  int foundAtom = 0;
199 –  int lineNum = 0;
292    int i;
293 <
293 >  int identNum;
294    
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 –
295    Atom** the_atoms;
296    int nAtoms;
297    the_atoms = entry_plug->atoms;
# Line 245 | Line 303 | void LJ_FF::initializeAtoms( void ){
303   #endif
304      
305      // read in the atom types.
306 <    
249 <    rewind( frcFile );
306 >
307      headAtomType = new LinkedType;
251 –    currentAtomType = headAtomType;
308      
309 <    eof_test = fgets( readLine, sizeof(readLine), frcFile );
310 <    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 <    
309 >    fastForward( "AtomTypes", "initializeAtoms" );
310 >
311      // we are now at the AtomTypes section.
312      
313      eof_test = fgets( readLine, sizeof(readLine), frcFile );
314      lineNum++;
315      
316 +    
317 +    // read a line, and start parseing out the atom types
318 +
319      if( eof_test == NULL ){
320        sprintf( painCave.errMsg,
321                 "Error in reading Atoms from force file at line %d.\n",
# Line 302 | Line 324 | void LJ_FF::initializeAtoms( void ){
324        simError();
325      }
326      
327 +    identNum = 1;
328 +    // stop reading at end of file, or at next section
329      while( readLine[0] != '#' && eof_test != NULL ){
330 <      
330 >
331 >      // toss comment lines
332        if( readLine[0] != '!' ){
333          
334 <        the_token = strtok( readLine, " \n\t,;" );
335 <        if( the_token != NULL ){
336 <          
337 <          strcpy( currentAtomType->name, the_token );
338 <          
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 );
334 >        // the parser returns 0 if the line was blank
335 >        if( parseAtomLJ( readLine, lineNum, info ) ){
336 >          info.ident = identNum;
337 >          headAtomType->add( info );;
338 >          identNum++;
339          }
340        }
342 –      
343 –      tempAtomType = new LinkedType;
344 –      currentAtomType->next = tempAtomType;
345 –      currentAtomType = tempAtomType;
346 –      
341        eof_test = fgets( readLine, sizeof(readLine), frcFile );
342        lineNum++;
343      }
# Line 356 | Line 350 | void LJ_FF::initializeAtoms( void ){
350               "LJ_FF atom structures read successfully." );
351      MPIcheckPoint();
352  
353 <    currentAtomType = headAtomType;
353 >    currentAtomType = headAtomType->next; //skip the first element who is a place holder.
354      while( currentAtomType != NULL ){
355        currentAtomType->duplicate( info );
356 +
357 +
358 +
359        sendFrcStruct( &info, mpiAtomStructType );
360 +
361 +      sprintf( checkPointMsg,
362 +               "successfully sent lJ force type: \"%s\"\n",
363 +               info.name );
364 +      MPIcheckPoint();
365 +
366        currentAtomType = currentAtomType->next;
367      }
368      info.last = 1;
# Line 375 | Line 378 | void LJ_FF::initializeAtoms( void ){
378  
379      headAtomType = new LinkedType;
380      recieveFrcStruct( &info, mpiAtomStructType );
381 +    
382      while( !info.last ){
383  
384 +
385 +
386        headAtomType->add( info );
387 +      
388 +      MPIcheckPoint();
389 +
390        recieveFrcStruct( &info, mpiAtomStructType );
391      }
392    }
393   #endif // is_mpi
394  
395 +  // call new A_types in fortran
396    
397 +  int isError;
398 +  currentAtomType = headAtomType;
399 +  while( currentAtomType != NULL ){
400 +    
401 +    if( currentAtomType->name[0] != '\0' ){
402 +      isError = 0;
403 +          newLJtype( &(currentAtomType->ident),
404 +                 &(currentAtomType->mass),
405 +                 &(currentAtomType->epslon),
406 +                 &(currentAtomType->sigma),
407 +                 &isError );
408 +      if( isError ){
409 +        sprintf( painCave.errMsg,
410 +                 "Error initializing the \"%s\" atom type in fortran\n",
411 +                 currentAtomType->name );
412 +        painCave.isFatal = 1;
413 +        simError();
414 +      }
415 +    }
416 +    currentAtomType = currentAtomType->next;
417 +  }
418 +      
419 + #ifdef IS_MPI
420 +  sprintf( checkPointMsg,
421 +           "LJ_FF atom structures successfully sent to fortran\n" );
422 +  MPIcheckPoint();
423 + #endif // is_mpi
424 +
425 +  
426 +
427    // initialize the atoms
428    
429 +  double bigSigma = 0.0;
430    Atom* thisAtom;
431  
432    for( i=0; i<nAtoms; i++ ){
# Line 402 | Line 443 | void LJ_FF::initializeAtoms( void ){
443      the_atoms[i]->setMass( currentAtomType->mass );
444      the_atoms[i]->setEpslon( currentAtomType->epslon );
445      the_atoms[i]->setSigma( currentAtomType->sigma );
446 +    the_atoms[i]->setIdent( currentAtomType->ident );
447      the_atoms[i]->setLJ();
448 +
449 +    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
450    }
451  
452 +  
453 + #ifdef IS_MPI
454 +  double tempBig = bigSigma;
455 +  MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
456 + #endif  //is_mpi
457  
458 +  //calc rCut and rList
459 +
460 +  entry_plug->rCut = 2.5 * bigSigma;
461 +  if(entry_plug->rCut > (entry_plug->box_x / 2.0)) entry_plug->rCut = entry_plug->box_x / 2.0;
462 +  if(entry_plug->rCut > (entry_plug->box_y / 2.0)) entry_plug->rCut = entry_plug->box_y / 2.0;
463 +  if(entry_plug->rCut > (entry_plug->box_z / 2.0)) entry_plug->rCut = entry_plug->box_z / 2.0;
464 +
465 +  entry_plug->rList = entry_plug->rCut + 1.0;
466 +
467    // clean up the memory
468    
469    delete headAtomType;
# Line 415 | Line 473 | void LJ_FF::initializeAtoms( void ){
473    MPIcheckPoint();
474   #endif // is_mpi
475  
476 +  initFortran();
477 +  entry_plug->refreshSim();
478 +
479   }
480  
481   void LJ_FF::initializeBonds( bond_pair* the_bonds ){
# Line 458 | Line 519 | void LJ_FF::initializeTorsions( torsion_set* the_torsi
519   #endif // is_mpi
520  
521   }
522 +
523 +
524 + void LJ_FF::fastForward( char* stopText, char* searchOwner ){
525 +
526 +  int foundText = 0;
527 +  char* the_token;
528 +
529 +  rewind( frcFile );
530 +  lineNum = 0;
531 +
532 +  eof_test = fgets( readLine, sizeof(readLine), frcFile );
533 +  lineNum++;
534 +  if( eof_test == NULL ){
535 +    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
536 +             " file is empty.\n",
537 +             searchOwner );
538 +    painCave.isFatal = 1;
539 +    simError();
540 +  }
541 +  
542 +  
543 +  while( !foundText ){
544 +    while( eof_test != NULL && readLine[0] != '#' ){
545 +      eof_test = fgets( readLine, sizeof(readLine), frcFile );
546 +      lineNum++;
547 +    }
548 +    if( eof_test == NULL ){
549 +      sprintf( painCave.errMsg,
550 +               "Error fast forwarding force file for %s at "
551 +               "line %d: file ended unexpectedly.\n",
552 +               searchOwner,
553 +               lineNum );
554 +      painCave.isFatal = 1;
555 +      simError();
556 +    }
557 +    
558 +    the_token = strtok( readLine, " ,;\t#\n" );
559 +    foundText = !strcmp( stopText, the_token );
560 +    
561 +    if( !foundText ){
562 +      eof_test = fgets( readLine, sizeof(readLine), frcFile );
563 +      lineNum++;
564 +      
565 +      if( eof_test == NULL ){
566 +        sprintf( painCave.errMsg,
567 +                 "Error fast forwarding force file for %s at "
568 +                 "line %d: file ended unexpectedly.\n",
569 +                 searchOwner,
570 +                 lineNum );
571 +        painCave.isFatal = 1;
572 +        simError();
573 +      }
574 +    }
575 +  }  
576 + }
577 +
578 +
579 +
580 + int parseAtomLJ( char *lineBuffer, int lineNum,  atomStruct &info ){
581 +
582 +  char* the_token;
583 +  
584 +  the_token = strtok( lineBuffer, " \n\t,;" );
585 +  if( the_token != NULL ){
586 +    
587 +    strcpy( info.name, the_token );
588 +          
589 +    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
590 +      sprintf( painCave.errMsg,
591 +               "Error parseing AtomTypes: line %d\n", lineNum );
592 +      painCave.isFatal = 1;
593 +      simError();
594 +    }
595 +    
596 +    info.mass = atof( the_token );
597 +    
598 +    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
599 +      sprintf( painCave.errMsg,
600 +               "Error parseing AtomTypes: line %d\n", lineNum );
601 +      painCave.isFatal = 1;
602 +      simError();
603 +    }
604 +        
605 +    info.epslon = atof( the_token );
606 +          
607 +    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
608 +      sprintf( painCave.errMsg,
609 +               "Error parseing AtomTypes: line %d\n", lineNum );
610 +      painCave.isFatal = 1;
611 +      simError();
612 +    }
613 +        
614 +    info.sigma = atof( the_token );
615 +    
616 +    return 1;
617 +  }
618 +  else return 0;
619 + }
620 +
621 +
622 + void LJ_FF::doForces( int calcPot ){
623 +
624 +  int i;
625 +  double* frc;
626 +  double* pos;
627 +  short int passedCalcPot = (short int)calcPot;
628 +
629 +  // forces are zeroed here, before any are acumulated.
630 +  // NOTE: do not rezero the forces in Fortran.
631 +
632 +  for(i=0; i<entry_plug->n_atoms; i++){
633 +    entry_plug->atoms[i]->zeroForces();
634 +  }
635 +
636 +  frc = Atom::getFrcArray();
637 +  pos = Atom::getPosArray();
638 +
639 + //   entry_plug->lrPot = -1;
640 +  doLJfortran( pos, frc, &(entry_plug->lrPot), &passedCalcPot );
641 +
642 +
643 +  //   fprintf( stderr,
644 +  //   "lrPot =  %lf\n", entry_plug->lrPot );
645 +  
646 + }
647 +  
648 + void LJ_FF::initFortran( void ){
649 +  
650 +  int nLocal = entry_plug->n_atoms;
651 +  int *ident;
652 +  int isError;
653 +  int i;
654 +
655 +  ident = new int[nLocal];
656 +
657 +  for(i=0; i<nLocal; i++){
658 +    ident[i] = entry_plug->atoms[i]->getIdent();
659 +  }
660 +
661 +  isError = 0;
662 +  initLJfortran( &nLocal, ident, &isError );
663 +  
664 +  if(isError){
665 +    sprintf( painCave.errMsg,
666 +             "LJ_FF error: There was an error initializing the component list in fortran.\n" );
667 +    painCave.isFatal = 1;
668 +    simError();
669 +  }
670 +
671 +  
672 + #ifdef IS_MPI
673 +  sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" );
674 +  MPIcheckPoint();
675 + #endif // is_mpi
676 +  
677 +  delete[] ident;
678 +
679 + }
680 +  

Diff Legend

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