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

Comparing trunk/OOPSE/libmdtools/EAM_FF.cpp (file contents):
Revision 627 by chuckv, Wed Jul 16 21:49:59 2003 UTC vs.
Revision 669 by chuckv, Thu Aug 7 00:47:33 2003 UTC

# Line 33 | Line 33 | namespace EAM_NS{
33      double eam_dr;    // The distance between each of the rho points.    
34      int eam_nrho;  // Number of points indexed by rho
35      int eam_nr;    // The number of points based on r (Both Phi(r) and Rho(r)).
36 <    int eam_rcut;  // The cutoff radius for eam.
36 >    double eam_rcut;  // The cutoff radius for eam.
37      int eam_ident; // Atomic number
38      int ident;
39      int last;      //  0  -> default
# Line 41 | Line 41 | namespace EAM_NS{
41    } atomStruct;
42  
43    int parseAtom( char *lineBuffer, int lineNum, atomStruct &info, char *eamPotFile );
44 <  int parseEAM( atomStruct &info, char *eamPotFile );
44 >  int parseEAM( atomStruct &info, char *eamPotFile, double **eam_rvals,
45 >                double **eam_rhovals, double **eam_Frhovals);
46   #ifdef IS_MPI
47    
48    MPI_Datatype mpiAtomStructType;
# Line 135 | Line 136 | namespace EAM_NS{
136      double eam_drho; // The distance between each of the points indexed by rho.
137      int eam_nr;   // The number of points based on r (Both Phi(r) and Rho(r)).
138      double eam_dr;   // The distance between each of the rho points.
139 <    int eam_rcut; // The cutoff radius for eam.
139 >    double eam_rcut; // The cutoff radius for eam.
140  
141      double *eam_rvals;    // Z of r values
142      double *eam_rhovals;  // rho of r values
# Line 150 | Line 151 | using namespace LJ_NS;
151  
152   }
153  
154 < using namespace LJ_NS;
154 > using namespace EAM_NS;
155  
156   //****************************************************************
157   // begins the actual forcefield stuff.  
# Line 168 | Line 169 | EAM_FF::EAM_FF(){
169    headAtomType = NULL;
170    currentAtomType = NULL;
171  
172 +  // Set eamRcut to 0.0
173 +  eamRcut = 0.0;
174 +
175    // do the funtion wrapping
176    wrapMeFF( this );
177  
# Line 263 | Line 267 | void EAM_FF::initForceField( int ljMixRule ){
267   #endif // is_mpi
268   }
269  
270 < void EAM_FF::initForceField( int ljMixRule ){
270 >
271 > void EAM_FF::calcRcut( void ){
272    
273 +  #ifdef IS_MPI
274 +  double tempEamRcut = eamRcut;
275 +  MPI_Allreduce( &tempEamRcut, &eamRcut, 1, MPI_DOUBLE, MPI_MAX,
276 +                 MPI_COMM_WORLD);
277 + #endif  //is_mpi
278 +  entry_plug->setRcut(eamRcut);
279 + }
280 +
281 +
282 + void EAM_FF::initForceField( int ljMixRule ){
283    initFortran( ljMixRule, 0 );
284   }
285  
# Line 283 | Line 298 | void EAM_FF::readParams( void ){
298   #endif // is_mpi
299   }
300  
301 +
302   void EAM_FF::readParams( void ){
303  
304    atomStruct info;
# Line 335 | Line 351 | void EAM_FF::readParams( void ){
351          // the parser returns 0 if the line was blank
352          if( parseAtom( readLine, lineNum, info, eamPotFile ) ){
353            parseEAM(info,eamPotFile, &eam_rvals,
354 <                   &eam_rhovals, &eam_Frhovals)){
354 >                   &eam_rhovals, &eam_Frhovals);
355            info.ident = identNum;
356            headAtomType->add( info, eam_rvals,
357                               eam_rhovals,eam_Frhovals );
# Line 437 | Line 453 | void EAM_FF::readParams( void ){
453    int isDipole = 0;
454    int isSSD = 0;
455    int isGB = 0;
456 <  int isEam = 1;
456 >  int isEAM= 1;
457    double dipole = 0.0;
458 +  double eamSigma = 0.0;
459 +  double eamEpslon = 0.0;
460    
461    currentAtomType = headAtomType->next;
462    while( currentAtomType != NULL ){
# Line 451 | Line 469 | void EAM_FF::readParams( void ){
469                   &isDipole,
470                   &isGB,
471                   &isEAM,
472 <                 &(currentAtomType->epslon),
473 <                 &(currentAtomType->sigma),
472 >                 &eamEpslon,
473 >                 &eamSigma,
474                   &dipole,
475                   &isError );
476        if( isError ){
# Line 467 | Line 485 | void EAM_FF::readParams( void ){
485    }
486        
487    entry_plug->useLJ = 0;
488 <
488 >  entry_plug->useEAM = 1;
489    // Walk down again and send out EAM type
490    currentAtomType = headAtomType->next;
491    while( currentAtomType != NULL ){
492      
493      if( currentAtomType->name[0] != '\0' ){
494        isError = 0;
495 +
496        newEAMtype( &(currentAtomType->lattice_constant),
497                    &(currentAtomType->eam_nrho),
498                    &(currentAtomType->eam_drho),
499                    &(currentAtomType->eam_nr),
500                    &(currentAtomType->eam_dr),
501 +                  &(currentAtomType->eam_rcut),
502                    currentAtomType->eam_rvals,
503                    currentAtomType->eam_rhovals,
504                    currentAtomType->eam_Frhovals,
505                    &(currentAtomType->eam_ident),
506                    &isError);
507 +
508        if( isError ){
509          sprintf( painCave.errMsg,
510                   "Error initializing the \"%s\" atom type in fortran EAM\n",
# Line 503 | Line 524 | void EAM_FF::readParams( void ){
524    MPIcheckPoint();
525   #endif // is_mpi
526  
527 +
528 +
529   }
530  
531  
532   void EAM_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
533    
534    int i;
535 <
535 >  
536    // initialize the atoms
537    
538  
# Line 529 | Line 552 | void EAM_FF::initializeAtoms( int nAtoms, Atom** the_a
552      the_atoms[i]->setMass( currentAtomType->mass );
553      the_atoms[i]->setIdent( currentAtomType->ident );
554      the_atoms[i]->setEAM();
555 +    the_atoms[i]->setEamRcut( currentAtomType->eam_rcut);
556  
557 +    if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut;
558 +
559    }
560   }
561  
562 < void LJ_FF::initializeBonds( int nBonds, Bond** BondArray,
562 > void EAM_FF::initializeBonds( int nBonds, Bond** BondArray,
563                               bond_pair* the_bonds ){
564    
565      if( nBonds ){
# Line 657 | Line 683 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
683                       double **eam_rvals,
684                       double **eam_rhovals,
685                       double **eam_Frhovals){
686 <  
686 >  double* myEam_rvals;
687 >  double* myEam_rhovals;
688 >  double* myEam_Frhovals;
689 >
690    char* ffPath_env = "FORCE_PARAM_PATH";
691    char* ffPath;
692    char* the_token;
693 +  char* eam_eof_test;
694    FILE *eamFile;
695 +  const int BUFFERSIZE = 3000;
696 +
697    char temp[200];
698    int linenumber;
699    int nReadLines;
700 +  char eam_read_buffer[BUFFERSIZE];
701  
702 +
703    int i,j;
704  
705    linenumber = 0;
# Line 674 | Line 708 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
708    eamFile = fopen( eamPotFile, "r" );
709    
710    
711 <  if( frcFile == NULL ){
711 >  if( eamFile == NULL ){
712      
713        // next see if the force path enviorment variable is set
714      
# Line 693 | Line 727 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
727  
728      
729      
730 <    if( frcFile == NULL ){
730 >    if( eamFile == NULL ){
731        
732        sprintf( painCave.errMsg,
733                 "Error opening the EAM force parameter file: %s\n"
# Line 706 | Line 740 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
740    }
741  
742    // First line is a comment line, read and toss it....
743 <  eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
743 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
744    linenumber++;
745 <  if(eof_test == NULL){
745 >  if(eam_eof_test == NULL){
746      sprintf( painCave.errMsg,
747               "error in reading commment in %s\n", eamPotFile);
748      painCave.isFatal = 1;
# Line 718 | Line 752 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
752  
753  
754    // The Second line contains atomic number, atomic mass and a lattice constant
755 <  eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
755 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
756    linenumber++;
757 <  if(eof_test == NULL){
757 >  if(eam_eof_test == NULL){
758      sprintf( painCave.errMsg,
759               "error in reading Identifier line in %s\n", eamPotFile);
760      painCave.isFatal = 1;
# Line 730 | Line 764 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
764  
765  
766      
767 <  if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
767 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
768      sprintf( painCave.errMsg,
769               "Error parseing EAM ident  line in %s\n", eamPotFile );
770      painCave.isFatal = 1;
# Line 756 | Line 790 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
790    info.lattice_constant = atof( the_token);
791  
792    // Next line is nrho, drho, nr, dr and rcut
793 <  eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile)
794 <  if(eof_test == NULL){
793 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
794 >  if(eam_eof_test == NULL){
795      sprintf( painCave.errMsg,
796               "error in reading number of points line in %s\n", eamPotFile);
797      painCave.isFatal = 1;
798      simError();
799    }
800  
801 <  if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
801 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
802      sprintf( painCave.errMsg,
803               "Error parseing EAM nrho: line in %s\n", eamPotFile );
804      painCave.isFatal = 1;
# Line 806 | Line 840 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
840    info.eam_rcut = atof( the_token);
841  
842  
843 +
844 +
845 +
846    // Ok now we have to allocate point arrays and read in number of points
847    // Index the arrays for fortran, starting at 1
848 <  *eam_Frhovals = new double[info.nrho];
849 <  *eam_rvals    = new double[info.nr];
850 <  *eam_rhovals  = new double[info.nr];
848 >  myEam_Frhovals = new double[info.eam_nrho];
849 >  myEam_rvals    = new double[info.eam_nr];
850 >  myEam_rhovals  = new double[info.eam_nr];
851  
852    // Parse F of rho vals.
853  
854    // Assume for now that we have a complete number of lines
855 <  nReadLines = int(info.nrho/5);
855 >  nReadLines = int(info.eam_nrho/5);
856 >  
857  
858 +
859    for (i=0;i<nReadLines;i++){
860      j = i*5;
861 <    
861 >
862      // Read next line
863 <    eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
863 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
864      linenumber++;
865 <    if(eof_test == NULL){
865 >    if(eam_eof_test == NULL){
866        sprintf( painCave.errMsg,
867                 "error in reading EAM file %s at line %d\n",
868                 eamPotFile,linenumber);
# Line 833 | Line 872 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
872      
873      // Parse 5 values on each line into array
874      // Value 1
875 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
875 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
876        sprintf( painCave.errMsg,
877                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
878        painCave.isFatal = 1;
879        simError();
880      }
881 +
882 +    myEam_Frhovals[j+0] = atof( the_token );
883      
843 –    *eam_Frhovals[j+0] = atof( the_token );
844 –    
884      // Value 2
885      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
886        sprintf( painCave.errMsg,
# Line 849 | Line 888 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
888        painCave.isFatal = 1;
889        simError();
890      }
891 <  
892 <    *eam_Frhovals[j+1] = atof( the_token );
891 >
892 >    myEam_Frhovals[j+1] = atof( the_token );
893      
894      // Value 3
895      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 859 | Line 898 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
898        painCave.isFatal = 1;
899        simError();
900      }
901 +
902 +    myEam_Frhovals[j+2] = atof( the_token );
903      
863 –    *eam_Frhovals[j+2] = atof( the_token );
864 –    
904      // Value 4
905      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
906        sprintf( painCave.errMsg,
# Line 869 | Line 908 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
908        painCave.isFatal = 1;
909        simError();
910      }
872 –    
873 –    *eam_Frhovals[j+3] = atof( the_token );
911  
912 +    myEam_Frhovals[j+3] = atof( the_token );
913 +
914      // Value 5
915      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
916        sprintf( painCave.errMsg,
# Line 879 | Line 918 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
918        painCave.isFatal = 1;
919        simError();
920      }
921 +
922 +    myEam_Frhovals[j+4] = atof( the_token );
923      
883 –    *eam_Frhovals[j+4] = atof( the_token );
884 –    
924    }
925    // Parse Z of r vals
926    
927    // Assume for now that we have a complete number of lines
928 <  nReadLines = int(info.nr/5);
928 >  nReadLines = int(info.eam_nr/5);
929  
930    for (i=0;i<nReadLines;i++){
931      j = i*5;
932  
933      // Read next line
934 <    eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
934 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
935      linenumber++;
936 <    if(eof_test == NULL){
936 >    if(eam_eof_test == NULL){
937        sprintf( painCave.errMsg,
938                 "error in reading EAM file %s at line %d\n",
939                 eamPotFile,linenumber);
# Line 904 | Line 943 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
943      
944      // Parse 5 values on each line into array
945      // Value 1
946 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
946 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
947        sprintf( painCave.errMsg,
948                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
949        painCave.isFatal = 1;
950        simError();
951      }
952      
953 <    *eam_rvals[j+0] = atof( the_token );
953 >    myEam_rvals[j+0] = atof( the_token );
954  
955      // Value 2
956      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 921 | Line 960 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
960        simError();
961      }
962    
963 <    *eam_rvals[j+1] = atof( the_token );
963 >    myEam_rvals[j+1] = atof( the_token );
964  
965      // Value 3
966      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 931 | Line 970 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
970        simError();
971      }
972    
973 <    *eam_rvals[j+2] = atof( the_token );
973 >    myEam_rvals[j+2] = atof( the_token );
974  
975      // Value 4
976      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 941 | Line 980 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
980        simError();
981      }
982    
983 <    *eam_rvals[j+3] = atof( the_token );
983 >    myEam_rvals[j+3] = atof( the_token );
984  
985      // Value 5
986      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 951 | Line 990 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
990        simError();
991      }
992    
993 <    *eam_rvals[j+4] = atof( the_token );
993 >    myEam_rvals[j+4] = atof( the_token );
994  
995    }
996    // Parse rho of r vals
# Line 962 | Line 1001 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
1001      j = i*5;
1002  
1003      // Read next line
1004 <    eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
1004 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
1005      linenumber++;
1006 <    if(eof_test == NULL){
1006 >    if(eam_eof_test == NULL){
1007        sprintf( painCave.errMsg,
1008                 "error in reading EAM file %s at line %d\n",
1009                 eamPotFile,linenumber);
# Line 974 | Line 1013 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
1013    
1014      // Parse 5 values on each line into array
1015      // Value 1
1016 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1016 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
1017        sprintf( painCave.errMsg,
1018                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1019        painCave.isFatal = 1;
1020        simError();
1021      }
1022    
1023 <    *eam_rhovals[j+0] = atof( the_token );
1023 >    myEam_rhovals[j+0] = atof( the_token );
1024  
1025      // Value 2
1026 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1026 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
1027        sprintf( painCave.errMsg,
1028                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1029        painCave.isFatal = 1;
1030        simError();
1031      }
1032    
1033 <    *eam_rhovals[j+1] = atof( the_token );
1033 >    myEam_rhovals[j+1] = atof( the_token );
1034  
1035      // Value 3
1036 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1036 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
1037        sprintf( painCave.errMsg,
1038                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1039        painCave.isFatal = 1;
1040        simError();
1041      }
1042    
1043 <    *eam_rhovals[j+2] = atof( the_token );
1043 >    myEam_rhovals[j+2] = atof( the_token );
1044  
1045      // Value 4
1046 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1046 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
1047        sprintf( painCave.errMsg,
1048                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1049        painCave.isFatal = 1;
1050        simError();
1051      }
1052    
1053 <    *eam_rhovals[j+3] = atof( the_token );
1053 >    myEam_rhovals[j+3] = atof( the_token );
1054  
1055      // Value 5
1056 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1056 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
1057        sprintf( painCave.errMsg,
1058                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1059        painCave.isFatal = 1;
1060        simError();
1061      }
1062    
1063 <    *eam_rhovals[j+4] = atof( the_token );
1063 >    myEam_rhovals[j+4] = atof( the_token );
1064  
1065    }
1066 +  *eam_rvals = myEam_rvals;
1067 +  *eam_rhovals = myEam_rhovals;
1068 +  *eam_Frhovals = myEam_Frhovals;
1069  
1070    fclose(eamFile);
1071    return 0;

Diff Legend

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