ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/WATER.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/WATER.cpp (file contents):
Revision 988 by chrisfen, Thu Jan 22 17:34:20 2004 UTC vs.
Revision 989 by gezelter, Tue Jan 27 19:38:27 2004 UTC

# Line 23 | Line 23 | namespace WATER_NS{
23  
24   namespace WATER_NS{
25  
26 <  // Declare the structures that will be passed by the parser and  MPI
26 >  // Declare the structures that will be passed by the parser and MPI
27    
28    typedef struct{
29      char name[15];
# Line 41 | Line 41 | namespace WATER_NS{
41  
42    typedef struct{
43      char name[15];
44 +    double Ixx;
45 +    double Iyy;
46 +    double Izz;
47      double dipole;
48      double w0;
49      double v0;
# Line 172 | Line 175 | namespace WATER_NS{
175          strcpy(next->name, info.name);
176          next->isDipole = info.isDipole;
177          next->isSticky = info.isSticky;
178 <        next->dipole   = info.dipole;
178 >        next->Ixx      = info.Ixx;
179 >        next->Iyy      = info.Iyy;
180 >        next->Izz      = info.Izz;
181 >        next->dipole   = info.dipole;
182          next->w0       = info.w0;
183          next->v0       = info.v0;
184          next->v0p      = info.v0p;
# Line 189 | Line 195 | namespace WATER_NS{
195        strcpy(info.name, name);
196        info.isDipole = isDipole;
197        info.isSticky = isSticky;
198 +      info.Ixx      = Ixx;
199 +      info.Iyy      = Iyy;
200 +      info.Izz      = Izz;
201        info.dipole   = dipole;
202        info.w0       = w0;
203        info.v0       = v0;
# Line 205 | Line 214 | namespace WATER_NS{
214      char name[15];
215      int isDipole;
216      int isSticky;
217 +    double Ixx;
218 +    double Iyy;
219 +    double Izz;
220      double dipole;
221      double w0;
222      double v0;
# Line 220 | Line 232 | namespace WATER_NS{
232    LinkedAtomType* currentAtomType;
233    LinkedDirectionalType* headDirectionalType;
234    LinkedDirectionalType* currentDirectionalType;
235 < }
235 > } // namespace
236  
237   using namespace WATER_NS;
238  
# Line 257 | Line 269 | WATER::WATER(){
269  
270    MPI_Address(&atomProto.name, &atomDspls[0]);
271    MPI_Address(&atomProto.mass, &atomDspls[1]);
272 <  MPI_Address(&atomProto.ident, &atomDspls[2]);
272 >  MPI_Address(&atomProto.isDirectional, &atomDspls[2]);
273    
274    atomMbrTypes[0] = MPI_CHAR;
275    atomMbrTypes[1] = MPI_DOUBLE;
# Line 274 | Line 286 | WATER::WATER(){
286    // Init the directionalStruct mpi type
287  
288    directionalStruct directionalProto; // mpiPrototype
289 <  int directionalBC[3] = {15,8,3};  // block counts
289 >  int directionalBC[3] = {15,11,3};  // block counts
290    MPI_Aint directionalDspls[3];           // displacements
291    MPI_Datatype directionalMbrTypes[3];    // member mpi types
292  
293    MPI_Address(&directionalProto.name, &directionalDspls[0]);
294 <  MPI_Address(&directionalProto.mass, &directionalDspls[1]);
295 <  MPI_Address(&directionalProto.ident, &directionalDspls[2]);
294 >  MPI_Address(&directionalProto.Ixx, &directionalDspls[1]);
295 >  MPI_Address(&directionalProto.isDipole, &directionalDspls[2]);
296    
297    directionalMbrTypes[0] = MPI_CHAR;
298    directionalMbrTypes[1] = MPI_DOUBLE;
# Line 390 | Line 402 | void WATER::readParams( void ){
402  
403    atomStruct atomInfo;
404    directionalStruct directionalInfo;
405 +  fpos_t *atomPos;
406  
407    atomInfo.last = 1;         // initialize last to have the last set.
408    directionalInfo.last = 1;  // if things go well, last will be set to 0
# Line 404 | Line 417 | void WATER::readParams( void ){
417      // read in the atom types.
418  
419      headAtomType = new LinkedAtomType;
420 +    headDirectionalType = new LinkedDirectionalType;
421      
422      fastForward( "AtomTypes", "initializeAtoms" );
423  
# Line 442 | Line 456 | void WATER::readParams( void ){
456              sectionSearch( "DirectionalTypes", atomInfo.name,
457                             "initializeDirectional" );
458              parseDirectional( readLine, lineNum, directionalInfo );
459 +            headDirectionalType->add( directionalInfo );
460  
461              // return to the AtomTypes section
462              fsetpos( frcFile, atomPos );
# Line 460 | Line 475 | void WATER::readParams( void ){
475      sprintf( checkPointMsg,
476               "WATER atom structures read successfully." );
477      MPIcheckPoint();
463 –    
478      currentAtomType = headAtomType->next; //skip the first element who is a place holder.
479      while( currentAtomType != NULL ){
480        currentAtomType->duplicate( atomInfo );
# Line 474 | Line 488 | void WATER::readParams( void ){
488        
489        currentAtomType = currentAtomType->next;
490      }
491 +
492      atomInfo.last = 1;
493      sendFrcStruct( &atomInfo, mpiAtomStructType );
494  
# Line 494 | Line 509 | void WATER::readParams( void ){
509        sendFrcStruct( &directionalInfo, mpiDirectionalStructType );
510      }
511    }
512 <  
512 >
513    else{
514      
515      // listen for node 0 to send out the force params
516      
517      MPIcheckPoint();
518 <    
518 >
519      headAtomType = new LinkedAtomType;
520      receiveFrcStruct( &atomInfo, mpiAtomStructType );
521      
# Line 533 | Line 548 | void WATER::readParams( void ){
548    }
549  
550   #endif // is_mpi
536 –
537 –
551    
552    // call new A_types in fortran
553    
# Line 543 | Line 556 | void WATER::readParams( void ){
556    // dummy variables
557    int isGB = 0;
558    int isEAM = 0;
559 <  
560 <  currentAtomType = headAtomType;
559 >  int notDipole = 0;
560 >  int notSSD = 0;
561 >  double noDipMoment = 0.0;
562 >
563 >
564 >  currentAtomType = headAtomType->next;
565 >  currentDirectionalType = headDirectionalType->next;
566 >
567    while( currentAtomType != NULL ){
568  
569      if( currentAtomType->isLJ ) entry_plug->useLJ = 1;
570 <    if( currentAtomType->isCharge ) entry_plug->useCharge = 1;
570 >    if( currentAtomType->isCharge ) entry_plug->useCharges = 1;
571      if( currentAtomType->isDirectional ){
572 <      if ( currentAtomType->isDipole ){
573 <        entry_plug->useDipoles = 1;
574 <        &(currentAtomType->dipole);
556 <      }
557 <      if ( currentAtomType->isSticky ) {
572 >      // only directional atoms can have dipoles or be sticky
573 >      if ( currentDirectionalType->isDipole ) entry_plug->useDipoles = 1;
574 >      if ( currentDirectionalType->isSticky ) {
575          entry_plug->useSticky = 1;
576 <        set_sticky_params( &(currentAtomType->w0), &(currentAtomType->v0),
577 <                           &(currentAtomType->v0p),
578 <                           &(currentAtomType->rl), &(currentAtomType->ru),
579 <                           &(currentAtomType->rlp), &(currentAtomType->rup));
576 >        set_sticky_params( &(currentDirectionalType->w0),
577 >                           &(currentDirectionalType->v0),
578 >                           &(currentDirectionalType->v0p),
579 >                           &(currentDirectionalType->rl),
580 >                           &(currentDirectionalType->ru),
581 >                           &(currentDirectionalType->rlp),
582 >                           &(currentDirectionalType->rup));
583        }
584 +      if( currentAtomType->name[0] != '\0' ){
585 +        isError = 0;
586 +        makeAtype( &(currentAtomType->ident),
587 +                   &(currentAtomType->isLJ),
588 +                   &(currentDirectionalType->isSticky),
589 +                   &(currentDirectionalType->isDipole),
590 +                   &isGB,
591 +                   &isEAM,
592 +                   &(currentAtomType->isCharge),
593 +                   &(currentAtomType->epslon),
594 +                   &(currentAtomType->sigma),
595 +                   &(currentAtomType->charge),
596 +                   &(currentDirectionalType->dipole),
597 +                   &isError );
598 +        if( isError ){
599 +          sprintf( painCave.errMsg,
600 +                   "Error initializing the \"%s\" atom type in fortran\n",
601 +                   currentAtomType->name );
602 +          painCave.isFatal = 1;
603 +          simError();
604 +        }
605 +      }
606 +      currentDirectionalType->next;
607      }
608 <    if( currentAtomType->name[0] != '\0' ){
609 <      isError = 0;
610 <      makeAtype( &(currentAtomType->ident),
611 <                 &isGB,
612 <                 &isEAM,
613 <                 &(currentAtomType->isDirectional),
614 <                 &(currentAtomType->isLJ),
615 <                 &(currentAtomType->isCharge),
616 <                 &(currentAtomType->epslon),
617 <                 &(currentAtomType->sigma),
618 <                 &(currentAtomType->charge),
619 <                 &isError );
620 <      if( isError ){
621 <        sprintf( painCave.errMsg,
622 <                 "Error initializing the \"%s\" atom type in fortran\n",
623 <                 currentAtomType->name );
624 <        painCave.isFatal = 1;
625 <        simError();
608 >
609 >    else {
610 >      // use all dummy variables if this is not a directional atom
611 >      if( currentAtomType->name[0] != '\0' ){
612 >        isError = 0;
613 >        makeAtype( &(currentAtomType->ident),
614 >                   &(currentAtomType->isLJ),
615 >                   &notSSD,
616 >                   &notDipole,
617 >                   &isGB,
618 >                   &isEAM,
619 >                   &(currentAtomType->isCharge),
620 >                   &(currentAtomType->epslon),
621 >                   &(currentAtomType->sigma),
622 >                   &(currentAtomType->charge),
623 >                   &noDipMoment,
624 >                   &isError );
625 >        if( isError ){
626 >          sprintf( painCave.errMsg,
627 >                   "Error initializing the \"%s\" atom type in fortran\n",
628 >                   currentAtomType->name );
629 >          painCave.isFatal = 1;
630 >          simError();
631 >        }
632        }
633      }
634      currentAtomType = currentAtomType->next;
# Line 590 | Line 639 | void WATER::readParams( void ){
639             "WATER atom structures successfully sent to fortran\n" );
640    MPIcheckPoint();
641   #endif // is_mpi
593 –
594 –
595 –
596 –  // now read in the directional stuff
597 –
598 – #ifdef IS_MPI
599 –  if( worldRank == 0 ) {
600 – #endif
601 –
602 –    // read in the directional types
603 –
604 –    headDirectionalType = new LinkedDirectionalType;
605 –    
606 –    fastForward( "DirectionalTypes", "initializeDirectionals" );
607 –  
608 –    // we are now at the directionalTypes section
609 –    
610 –    eof_test =  fgets( readLine, sizeof(readLine), frcFile );
611 –    lineNum++;
612 –    
613 –    
614 –    // read a line, and start parsing out the atom types
615 –
616 –    if( eof_test == NULL ){
617 –      sprintf( painCave.errMsg,
618 –               "Error in reading directionals from force file at line %d.\n",
619 –               lineNum );
620 –      painCave.isFatal = 1;
621 –      simError();
622 –    }
623 –    
624 –    // stop reading at end of file, or at next section
625 –    while( readLine[0] != '#' && eof_test != NULL ){
626 –
627 –      // toss comment lines
628 –      if( readLine[0] != '!' ){
629 –        
630 –        // the parser returns 0 if the line was blank
631 –        if( parseDirectional( readLine, lineNum, directionalInfo ) ){
632 –          headDirectionalType->add( directionalInfo );
633 –        }
634 –      }
635 –      eof_test = fgets( readLine, sizeof(readLine), frcFile );
636 –      lineNum++;
637 –    }
638 –        
639 – #ifdef IS_MPI
640 –    
641 –    // send out the linked list to all the other processes
642 –    
643 –    sprintf( checkPointMsg,
644 –             "DUFF directional structures read successfully." );
645 –    MPIcheckPoint();
646 –    
647 –    currentDirectionalType = headDirectionalType->next;
648 –    while( currentDirectionalType != NULL ){
649 –      currentDirectionalType->duplicate( directionalInfo );
650 –      sendFrcStruct( &directionalInfo, mpiDirectionalStructType );
651 –      currentDirectionalType = currentDirectionalType->next;
652 –    }
653 –    directionalInfo.last = 1;
654 –    sendFrcStruct( &directionalInfo, mpiDirectionalStructType );
655 –    
656 –  }
642  
658 –  else{
659 –    
660 –    // listen for node 0 to send out the force params
661 –    
662 –    MPIcheckPoint();
663 –
664 –    headDirectionalType = new LinkedDirectionalType;
665 –    receiveFrcStruct( &directionalInfo, mpiDirectionalStructType );
666 –    while( !directionalInfo.last ){
667 –
668 –      headDirectionalType->add( directionalInfo );
669 –      receiveFrcStruct( &directionalInfo, mpiDirectionalStructType );
670 –    }
671 –  }
672 –
673 –  sprintf( checkPointMsg,
674 –           "WATER directional structures broadcast successfully." );
675 –  MPIcheckPoint();
676 –
677 – #endif // is_mpi
643   }
644  
645  
646   void WATER::initializeAtoms( int nAtoms, Atom** the_atoms ){
647    
648 <  int i;
648 >  int i,j;
649  
650    // initialize the atoms
651 <  
651 >  DirectionalAtom* dAtom;
652 >  double inertialMat[3][3];
653  
654    for( i=0; i<nAtoms; i++ ){
655 <    
655 >    fprintf(stderr, "flag 1\n");
656      currentAtomType = headAtomType->find( the_atoms[i]->getType() );
657 +    fprintf(stderr, "%s is the type\n", the_atoms[i]->getType());
658      if( currentAtomType == NULL ){
659        sprintf( painCave.errMsg,
660                 "AtomType error, %s not found in force file.\n",
# Line 695 | Line 662 | void WATER::initializeAtoms( int nAtoms, Atom** the_at
662        painCave.isFatal = 1;
663        simError();
664      }
665 <
666 <    the_atoms[i]->setisLJ( currentAtomType->isLJ);
667 <    the_atoms[i]->setisCharge( currentAtomType->isCharge);
665 >    fprintf(stderr, "flag 2\n");
666 >    if( currentAtomType->isLJ ) the_atoms[i]->setLJ();
667 >    if( currentAtomType->isCharge ) the_atoms[i]->setCharged();
668      the_atoms[i]->setMass( currentAtomType->mass );
702 –    the_atoms[i]->setEpslon( currentAtomType->epslon );
703 –    the_atoms[i]->setSigma( currentAtomType->sigma );
704 –    the_atoms[i]->setCharge( currentAtomType->charge );
669      the_atoms[i]->setIdent( currentAtomType->ident );
670 +    fprintf(stderr, "flag 3\n");
671 +    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
672  
673 <  }
674 < }
673 >    if( currentAtomType->isDirectional ){
674 >      fprintf(stderr, "flag 4\n");
675 >      currentDirectionalType =
676 >        headDirectionalType->find( the_atoms[i]->getType() );
677 >      fprintf(stderr, "%s is the type\n", the_atoms[i]->getType());
678 >      if( currentDirectionalType == NULL ){
679 >        sprintf( painCave.errMsg,
680 >                 "DirectionalType error, %s not found in force file.\n",
681 >                 the_atoms[i]->getType() );
682 >        painCave.isFatal = 1;
683 >        simError();
684 >      }
685  
686 < void WATER::initializeDirectional( int nAtoms, Atom** the_atoms ){
687 <  
688 <  int i;
686 >      // zero out the moments of inertia matrix
687 >      for( i=0; i<3; i++ )
688 >        for( j=0; j<3; j++ )
689 >          inertialMat[i][j] = 0.0;
690  
691 <  // initialize the atoms
692 <  
693 <
694 <  for( i=0; i<nAtoms; i++ ){
695 <    
696 <    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
697 <    if( currentAtomType == NULL ){
698 <      sprintf( painCave.errMsg,
699 <               "AtomType error, %s not found in force file.\n",
700 <               the_atoms[i]->getType() );
691 >      // load the force file moments of inertia
692 >      inertialMat[0][0] = currentDirectionalType->Ixx;
693 >      inertialMat[1][1] = currentDirectionalType->Iyy;
694 >      inertialMat[2][2] = currentDirectionalType->Izz;
695 >      fprintf(stderr, "Let's try pointing to isDirectional\n");
696 >      fprintf(stderr, "%i what is this\n",the_atoms[i]->isDirectional());
697 >      dAtom = (DirectionalAtom *) the_atoms[i];
698 >      fprintf(stderr, "%i is isDipole\n", currentDirectionalType->isDipole);
699 >      dAtom->setHasDipole( currentDirectionalType->isDipole );
700 >      dAtom->setMu( currentDirectionalType->dipole );
701 >      fprintf(stderr, "flag 5\n");
702 >      dAtom->setMu( currentDirectionalType->dipole );
703 >      fprintf(stderr,"flag 6\n");
704 >      // if it's sticky then it's an SSD type
705 >      dAtom->setSSD( currentDirectionalType->isSticky );
706 >      dAtom->setJx( 0.0 );
707 >      dAtom->setJy( 0.0 );
708 >      dAtom->setJz( 0.0 );
709 >      dAtom->setI( inertialMat );
710 >      fprintf(stderr, "flag 7\n");
711 >      entry_plug->n_dipoles++;
712 >      fprintf(stderr, "flag 8\n");
713 >    }
714 >    else{
715 >      sprintf( painCave.errMsg,
716 >               "WATER error: Atom \"%s\" is directional, yet no standard"
717 >               " orientation was specifed in the BASS file.\n",
718 >               currentAtomType->name );
719        painCave.isFatal = 1;
720        simError();
721      }
727 –
728 –    the_atoms[i]->setisLJ( currentAtomType->isLJ);
729 –    the_atoms[i]->setisCharge( currentAtomType->isCharge);
730 –    the_atoms[i]->setMass( currentAtomType->mass );
731 –    the_atoms[i]->setEpslon( currentAtomType->epslon );
732 –    the_atoms[i]->setSigma( currentAtomType->sigma );
733 –    the_atoms[i]->setCharge( currentAtomType->charge );
734 –    the_atoms[i]->setIdent( currentAtomType->ident );
735 –
722    }
723   }
724  
# Line 828 | Line 814 | void WATER::sectionSearch( char* secHead, char* stopTe
814    int foundSection = 0;
815    int foundText = 0;
816    char* the_token;
817 <  tempPos = new fpos_t;
817 >  fpos_t *tempPos;
818  
819    rewind( frcFile );
820    lineNum = 0;
821 +  tempPos = new fpos_t;
822  
823    eof_test = fgets( readLine, sizeof(readLine), frcFile );
824    lineNum++;
# Line 843 | Line 830 | void WATER::sectionSearch( char* secHead, char* stopTe
830      simError();
831    }
832    
846 –  
833    while( !foundSection ){
834      while( eof_test != NULL && readLine[0] != '#' ){
835        eof_test = fgets( readLine, sizeof(readLine), frcFile );
# Line 858 | Line 844 | void WATER::sectionSearch( char* secHead, char* stopTe
844        painCave.isFatal = 1;
845        simError();
846      }
861 –    
847      the_token = strtok( readLine, " ,;\t#\n" );
848      foundSection = !strcmp( secHead, the_token );
849      
865 –    
850      if( !foundSection ){
851        eof_test = fgets( readLine, sizeof(readLine), frcFile );
852        lineNum++;
# Line 870 | Line 854 | void WATER::sectionSearch( char* secHead, char* stopTe
854        if( eof_test == NULL ){
855          sprintf( painCave.errMsg,
856                   "Error section searching force file for %s at "
873 –                 "line %d: file ended unexpectedly.\n",
874 –                 searchOwner,
875 –                 lineNum );
876 –        painCave.isFatal = 1;
877 –        simError();
878 –      }
879 –    }
880 –  }  
881 –
882 –  while( !foundText ){
883 –    if( !foundText ){
884 –      fgetpos( frcFile, tempPos );
885 –      eof_test = fgets( readLine, sizeof(readLine), frcFile );
886 –      lineNum++;
887 –      
888 –      if( eof_test == NULL ){
889 –        sprintf( painCave.errMsg,
890 –                 "Error fast forwarding force file for %s at "
857                   "line %d: file ended unexpectedly.\n",
858                   searchOwner,
859                   lineNum );
# Line 895 | Line 861 | void WATER::sectionSearch( char* secHead, char* stopTe
861          simError();
862        }
863      }
864 +  }
865 +
866 +  while( !foundText ){
867      
868 +    fgetpos( frcFile, tempPos );
869 +    eof_test = fgets( readLine, sizeof(readLine), frcFile );
870 +    lineNum++;
871 +    
872 +    if( eof_test == NULL ){
873 +      sprintf( painCave.errMsg,
874 +               "Error fast forwarding force file for %s at "
875 +               "line %d: file ended unexpectedly.\n",
876 +               searchOwner,
877 +               lineNum );
878 +      painCave.isFatal = 1;
879 +        simError();
880 +    }
881 +    
882      the_token = strtok( readLine, " ,;\t#\n" );
883 <    foundText = !strcmp( stopText, the_token );
883 >    if( the_token != NULL ){
884 >      foundText = !strcmp( stopText, the_token );
885 >    }
886    }  
887 <  
888 <  fsetPos( frcFile, tempPos );
887 >  fsetpos( frcFile, tempPos );
888 >  eof_test = fgets( readLine, sizeof(readLine), frcFile );
889 >  lineNum++;
890   }
891  
892  
# Line 1015 | Line 1001 | int WATER_NS::parseDirectional( char *lineBuffer, int
1001        simError();
1002      }
1003  
1004 <    info.I_xx = atof( the_token );
1004 >    info.Ixx = atof( the_token );
1005      
1006      if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1007        sprintf( painCave.errMsg,
# Line 1024 | Line 1010 | int WATER_NS::parseDirectional( char *lineBuffer, int
1010        simError();
1011      }
1012  
1013 <    info.I_yy = atof( the_token );
1013 >    info.Iyy = atof( the_token );
1014      
1015      if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1016        sprintf( painCave.errMsg,
# Line 1033 | Line 1019 | int WATER_NS::parseDirectional( char *lineBuffer, int
1019        simError();
1020      }
1021  
1022 <    info.I_zz = atof( the_token );
1022 >    info.Izz = atof( the_token );
1023      
1024      if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1025        sprintf( painCave.errMsg,

Diff Legend

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