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 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC

# Line 1 | Line 1
1 < #include <cstdlib>
2 < #include <cstdio>
3 < #include <cstring>
1 > #include <stdlib.h>
2 > #include <stdio.h>
3 > #include <string.h>
4  
5   #include <iostream>
6   using namespace std;
# Line 30 | Line 30 | namespace EAM_NS{
30      double mass;
31      double lattice_constant;
32      double eam_drho;  // The distance between each of the points indexed by rho.
33 +    double eam_rcut;  // The cutoff radius for eam.
34      double eam_dr;    // The distance between each of the rho points.    
35      int eam_nrho;  // Number of points indexed by rho
36      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.
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 74 | Line 75 | namespace EAM_NS{
75  
76      void add( atomStruct &info, double *the_eam_rvals,
77                double *the_eam_rhovals,double *the_eam_Frhovals ){
77 –
78 –      int i;
78  
79        // check for duplicates
80        
# Line 135 | Line 134 | namespace EAM_NS{
134      double eam_drho; // The distance between each of the points indexed by rho.
135      int eam_nr;   // The number of points based on r (Both Phi(r) and Rho(r)).
136      double eam_dr;   // The distance between each of the rho points.
137 <    int eam_rcut; // The cutoff radius for eam.
137 >    double eam_rcut; // The cutoff radius for eam.
138  
139      double *eam_rvals;    // Z of r values
140      double *eam_rhovals;  // rho of r values
# Line 150 | Line 149 | namespace EAM_NS{
149  
150   }
151  
152 < using namespace LJ_NS;
152 > using namespace EAM_NS;
153  
154   //****************************************************************
155   // begins the actual forcefield stuff.  
# Line 163 | Line 162 | EAM_FF::EAM_FF(){
162    char* ffPath_env = "FORCE_PARAM_PATH";
163    char* ffPath;
164    char temp[200];
166 –  char errMsg[1000];
165  
166    headAtomType = NULL;
167    currentAtomType = NULL;
168  
169 +  // Set eamRcut to 0.0
170 +  eamRcut = 0.0;
171 +
172    // do the funtion wrapping
173    wrapMeFF( this );
174  
# Line 178 | Line 179 | EAM_FF::EAM_FF(){
179    // Init the atomStruct mpi type
180  
181    atomStruct atomProto; // mpiPrototype
182 <  int atomBC[3] = {15,4,6};  // block counts
182 >  int atomBC[3] = {15,5,5};  // block counts
183    MPI_Aint atomDspls[3];           // displacements
184    MPI_Datatype atomMbrTypes[3];    // member mpi types
185  
# Line 263 | Line 264 | EAM_FF::~EAM_FF(){
264   #endif // is_mpi
265   }
266  
267 < void EAM_FF::initForceField( int ljMixRule ){
267 >
268 > void EAM_FF::calcRcut( void ){
269    
270 +  #ifdef IS_MPI
271 +  double tempEamRcut = eamRcut;
272 +  MPI_Allreduce( &tempEamRcut, &eamRcut, 1, MPI_DOUBLE, MPI_MAX,
273 +                 MPI_COMM_WORLD);
274 + #endif  //is_mpi
275 +  entry_plug->setRcut(eamRcut);
276 + }
277 +
278 +
279 + void EAM_FF::initForceField( int ljMixRule ){
280    initFortran( ljMixRule, 0 );
281   }
282  
# Line 283 | Line 295 | void EAM_FF::cleanMe( void ){
295   #endif // is_mpi
296   }
297  
298 +
299   void EAM_FF::readParams( void ){
300  
301    atomStruct info;
302    info.last = 1; // initialize last to have the last set.
303                   // if things go well, last will be set to 0
304  
292 –  int i;
305    int identNum;
306    double *eam_rvals;    // Z of r values
307    double *eam_rhovals;  // rho of r values
# Line 335 | Line 347 | void EAM_FF::readParams( void ){
347          // the parser returns 0 if the line was blank
348          if( parseAtom( readLine, lineNum, info, eamPotFile ) ){
349            parseEAM(info,eamPotFile, &eam_rvals,
350 <                   &eam_rhovals, &eam_Frhovals)){
350 >                   &eam_rhovals, &eam_Frhovals);
351            info.ident = identNum;
352            headAtomType->add( info, eam_rvals,
353                               eam_rhovals,eam_Frhovals );
# Line 437 | Line 449 | void EAM_FF::readParams( void ){
449    int isDipole = 0;
450    int isSSD = 0;
451    int isGB = 0;
452 <  int isEam = 1;
452 >  int isEAM= 1;
453    double dipole = 0.0;
454 +  double eamSigma = 0.0;
455 +  double eamEpslon = 0.0;
456    
457    currentAtomType = headAtomType->next;
458    while( currentAtomType != NULL ){
# Line 451 | Line 465 | void EAM_FF::readParams( void ){
465                   &isDipole,
466                   &isGB,
467                   &isEAM,
468 <                 &(currentAtomType->epslon),
469 <                 &(currentAtomType->sigma),
468 >                 &eamEpslon,
469 >                 &eamSigma,
470                   &dipole,
471                   &isError );
472        if( isError ){
# Line 467 | Line 481 | void EAM_FF::readParams( void ){
481    }
482        
483    entry_plug->useLJ = 0;
484 <
484 >  entry_plug->useEAM = 1;
485    // Walk down again and send out EAM type
486    currentAtomType = headAtomType->next;
487    while( currentAtomType != NULL ){
488      
489      if( currentAtomType->name[0] != '\0' ){
490        isError = 0;
491 +
492        newEAMtype( &(currentAtomType->lattice_constant),
493                    &(currentAtomType->eam_nrho),
494                    &(currentAtomType->eam_drho),
495                    &(currentAtomType->eam_nr),
496                    &(currentAtomType->eam_dr),
497 +                  &(currentAtomType->eam_rcut),
498                    currentAtomType->eam_rvals,
499                    currentAtomType->eam_rhovals,
500                    currentAtomType->eam_Frhovals,
501                    &(currentAtomType->eam_ident),
502                    &isError);
503 +
504        if( isError ){
505          sprintf( painCave.errMsg,
506                   "Error initializing the \"%s\" atom type in fortran EAM\n",
# Line 503 | Line 520 | void EAM_FF::readParams( void ){
520    MPIcheckPoint();
521   #endif // is_mpi
522  
523 +
524 +
525   }
526  
527  
528   void EAM_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
529    
530    int i;
531 <
531 >  
532    // initialize the atoms
533    
515 –
516 –  Atom* thisAtom;
517 –
534    for( i=0; i<nAtoms; i++ ){
535      
536      currentAtomType = headAtomType->find( the_atoms[i]->getType() );
# Line 529 | Line 545 | void EAM_FF::initializeAtoms( int nAtoms, Atom** the_a
545      the_atoms[i]->setMass( currentAtomType->mass );
546      the_atoms[i]->setIdent( currentAtomType->ident );
547      the_atoms[i]->setEAM();
548 +    the_atoms[i]->setEamRcut( currentAtomType->eam_rcut);
549  
550 +    if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut;
551 +
552    }
553   }
554  
555 < void LJ_FF::initializeBonds( int nBonds, Bond** BondArray,
555 > void EAM_FF::initializeBonds( int nBonds, Bond** BondArray,
556                               bond_pair* the_bonds ){
557    
558      if( nBonds ){
# Line 657 | Line 676 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
676                       double **eam_rvals,
677                       double **eam_rhovals,
678                       double **eam_Frhovals){
679 <  
679 >  double* myEam_rvals;
680 >  double* myEam_rhovals;
681 >  double* myEam_Frhovals;
682 >
683    char* ffPath_env = "FORCE_PARAM_PATH";
684    char* ffPath;
685    char* the_token;
686 +  char* eam_eof_test;
687    FILE *eamFile;
688 +  const int BUFFERSIZE = 3000;
689 +
690    char temp[200];
691    int linenumber;
692    int nReadLines;
693 +  char eam_read_buffer[BUFFERSIZE];
694 +
695  
696    int i,j;
697  
# Line 674 | Line 701 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
701    eamFile = fopen( eamPotFile, "r" );
702    
703    
704 <  if( frcFile == NULL ){
704 >  if( eamFile == NULL ){
705      
706        // next see if the force path enviorment variable is set
707      
# Line 693 | Line 720 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
720  
721      
722      
723 <    if( frcFile == NULL ){
723 >    if( eamFile == NULL ){
724        
725        sprintf( painCave.errMsg,
726                 "Error opening the EAM force parameter file: %s\n"
# Line 706 | Line 733 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
733    }
734  
735    // First line is a comment line, read and toss it....
736 <  eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
736 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
737    linenumber++;
738 <  if(eof_test == NULL){
738 >  if(eam_eof_test == NULL){
739      sprintf( painCave.errMsg,
740               "error in reading commment in %s\n", eamPotFile);
741      painCave.isFatal = 1;
# Line 718 | Line 745 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
745  
746  
747    // The Second line contains atomic number, atomic mass and a lattice constant
748 <  eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
748 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
749    linenumber++;
750 <  if(eof_test == NULL){
750 >  if(eam_eof_test == NULL){
751      sprintf( painCave.errMsg,
752               "error in reading Identifier line in %s\n", eamPotFile);
753      painCave.isFatal = 1;
# Line 730 | Line 757 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
757  
758  
759      
760 <  if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
760 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
761      sprintf( painCave.errMsg,
762               "Error parseing EAM ident  line in %s\n", eamPotFile );
763      painCave.isFatal = 1;
# Line 756 | Line 783 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
783    info.lattice_constant = atof( the_token);
784  
785    // Next line is nrho, drho, nr, dr and rcut
786 <  eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile)
787 <  if(eof_test == NULL){
786 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
787 >  if(eam_eof_test == NULL){
788      sprintf( painCave.errMsg,
789               "error in reading number of points line in %s\n", eamPotFile);
790      painCave.isFatal = 1;
791      simError();
792    }
793  
794 <  if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
794 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
795      sprintf( painCave.errMsg,
796               "Error parseing EAM nrho: line in %s\n", eamPotFile );
797      painCave.isFatal = 1;
# Line 806 | Line 833 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
833    info.eam_rcut = atof( the_token);
834  
835  
836 +
837 +
838 +
839    // Ok now we have to allocate point arrays and read in number of points
840    // Index the arrays for fortran, starting at 1
841 <  *eam_Frhovals = new double[info.nrho];
842 <  *eam_rvals    = new double[info.nr];
843 <  *eam_rhovals  = new double[info.nr];
841 >  myEam_Frhovals = new double[info.eam_nrho];
842 >  myEam_rvals    = new double[info.eam_nr];
843 >  myEam_rhovals  = new double[info.eam_nr];
844  
845    // Parse F of rho vals.
846  
847    // Assume for now that we have a complete number of lines
848 <  nReadLines = int(info.nrho/5);
848 >  nReadLines = int(info.eam_nrho/5);
849 >  
850  
851 +
852    for (i=0;i<nReadLines;i++){
853      j = i*5;
854 <    
854 >
855      // Read next line
856 <    eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
856 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
857      linenumber++;
858 <    if(eof_test == NULL){
858 >    if(eam_eof_test == NULL){
859        sprintf( painCave.errMsg,
860                 "error in reading EAM file %s at line %d\n",
861                 eamPotFile,linenumber);
# Line 833 | Line 865 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
865      
866      // Parse 5 values on each line into array
867      // Value 1
868 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
868 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
869        sprintf( painCave.errMsg,
870                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
871        painCave.isFatal = 1;
872        simError();
873      }
874 +
875 +    myEam_Frhovals[j+0] = atof( the_token );
876      
843 –    *eam_Frhovals[j+0] = atof( the_token );
844 –    
877      // Value 2
878      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
879        sprintf( painCave.errMsg,
# Line 849 | Line 881 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
881        painCave.isFatal = 1;
882        simError();
883      }
884 <  
885 <    *eam_Frhovals[j+1] = atof( the_token );
884 >
885 >    myEam_Frhovals[j+1] = atof( the_token );
886      
887      // Value 3
888      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 859 | Line 891 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
891        painCave.isFatal = 1;
892        simError();
893      }
894 +
895 +    myEam_Frhovals[j+2] = atof( the_token );
896      
863 –    *eam_Frhovals[j+2] = atof( the_token );
864 –    
897      // Value 4
898      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
899        sprintf( painCave.errMsg,
# Line 869 | Line 901 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
901        painCave.isFatal = 1;
902        simError();
903      }
872 –    
873 –    *eam_Frhovals[j+3] = atof( the_token );
904  
905 +    myEam_Frhovals[j+3] = atof( the_token );
906 +
907      // Value 5
908      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
909        sprintf( painCave.errMsg,
# Line 879 | Line 911 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
911        painCave.isFatal = 1;
912        simError();
913      }
914 +
915 +    myEam_Frhovals[j+4] = atof( the_token );
916      
883 –    *eam_Frhovals[j+4] = atof( the_token );
884 –    
917    }
918    // Parse Z of r vals
919    
920    // Assume for now that we have a complete number of lines
921 <  nReadLines = int(info.nr/5);
921 >  nReadLines = int(info.eam_nr/5);
922  
923    for (i=0;i<nReadLines;i++){
924      j = i*5;
925  
926      // Read next line
927 <    eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
927 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
928      linenumber++;
929 <    if(eof_test == NULL){
929 >    if(eam_eof_test == NULL){
930        sprintf( painCave.errMsg,
931                 "error in reading EAM file %s at line %d\n",
932                 eamPotFile,linenumber);
# Line 904 | Line 936 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
936      
937      // Parse 5 values on each line into array
938      // Value 1
939 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
939 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
940        sprintf( painCave.errMsg,
941                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
942        painCave.isFatal = 1;
943        simError();
944      }
945      
946 <    *eam_rvals[j+0] = atof( the_token );
946 >    myEam_rvals[j+0] = atof( the_token );
947  
948      // Value 2
949      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 921 | Line 953 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
953        simError();
954      }
955    
956 <    *eam_rvals[j+1] = atof( the_token );
956 >    myEam_rvals[j+1] = atof( the_token );
957  
958      // Value 3
959      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 931 | Line 963 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
963        simError();
964      }
965    
966 <    *eam_rvals[j+2] = atof( the_token );
966 >    myEam_rvals[j+2] = atof( the_token );
967  
968      // Value 4
969      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 941 | Line 973 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
973        simError();
974      }
975    
976 <    *eam_rvals[j+3] = atof( the_token );
976 >    myEam_rvals[j+3] = atof( the_token );
977  
978      // Value 5
979      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 951 | Line 983 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
983        simError();
984      }
985    
986 <    *eam_rvals[j+4] = atof( the_token );
986 >    myEam_rvals[j+4] = atof( the_token );
987  
988    }
989    // Parse rho of r vals
# Line 962 | Line 994 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
994      j = i*5;
995  
996      // Read next line
997 <    eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
997 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
998      linenumber++;
999 <    if(eof_test == NULL){
999 >    if(eam_eof_test == NULL){
1000        sprintf( painCave.errMsg,
1001                 "error in reading EAM file %s at line %d\n",
1002                 eamPotFile,linenumber);
# Line 974 | Line 1006 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
1006    
1007      // Parse 5 values on each line into array
1008      // Value 1
1009 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1009 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
1010        sprintf( painCave.errMsg,
1011                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1012        painCave.isFatal = 1;
1013        simError();
1014      }
1015    
1016 <    *eam_rhovals[j+0] = atof( the_token );
1016 >    myEam_rhovals[j+0] = atof( the_token );
1017  
1018      // Value 2
1019 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1019 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1020        sprintf( painCave.errMsg,
1021                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1022        painCave.isFatal = 1;
1023        simError();
1024      }
1025    
1026 <    *eam_rhovals[j+1] = atof( the_token );
1026 >    myEam_rhovals[j+1] = atof( the_token );
1027  
1028      // Value 3
1029 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1029 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1030        sprintf( painCave.errMsg,
1031                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1032        painCave.isFatal = 1;
1033        simError();
1034      }
1035    
1036 <    *eam_rhovals[j+2] = atof( the_token );
1036 >    myEam_rhovals[j+2] = atof( the_token );
1037  
1038      // Value 4
1039 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1039 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1040        sprintf( painCave.errMsg,
1041                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1042        painCave.isFatal = 1;
1043        simError();
1044      }
1045    
1046 <    *eam_rhovals[j+3] = atof( the_token );
1046 >    myEam_rhovals[j+3] = atof( the_token );
1047  
1048      // Value 5
1049 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1049 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1050        sprintf( painCave.errMsg,
1051                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1052        painCave.isFatal = 1;
1053        simError();
1054      }
1055    
1056 <    *eam_rhovals[j+4] = atof( the_token );
1057 <
1056 >    myEam_rhovals[j+4] = atof( the_token );
1057 >
1058    }
1059 +  *eam_rvals = myEam_rvals;
1060 +  *eam_rhovals = myEam_rhovals;
1061 +  *eam_Frhovals = myEam_Frhovals;
1062  
1063    fclose(eamFile);
1064    return 0;

Diff Legend

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