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 1097 by gezelter, Mon Apr 12 20:32:20 2004 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 75 | Line 76 | namespace EAM_NS{
76      void add( atomStruct &info, double *the_eam_rvals,
77                double *the_eam_rhovals,double *the_eam_Frhovals ){
78  
78 –      int i;
79 –
79        // check for duplicates
80        
81        if( !strcmp( info.name, name ) ){
# 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->setDefaultRcut(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 395 | Line 407 | void EAM_FF::readParams( void ){
407      MPIcheckPoint();
408  
409      headAtomType = new LinkedAtomType;
410 <    recieveFrcStruct( &info, mpiAtomStructType );
410 >    receiveFrcStruct( &info, mpiAtomStructType );
411  
412      while( !info.last ){
413        
# Line 421 | Line 433 | void EAM_FF::readParams( void ){
433        
434        MPIcheckPoint();
435  
436 <      recieveFrcStruct( &info, mpiAtomStructType );
436 >      receiveFrcStruct( &info, mpiAtomStructType );
437  
438  
439      }
# 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 >  int isCharge = 0;
454    double dipole = 0.0;
455 +  double charge = 0.0;
456 +  double eamSigma = 0.0;
457 +  double eamEpslon = 0.0;
458    
459    currentAtomType = headAtomType->next;
460    while( currentAtomType != NULL ){
# Line 451 | Line 467 | void EAM_FF::readParams( void ){
467                   &isDipole,
468                   &isGB,
469                   &isEAM,
470 <                 &(currentAtomType->epslon),
471 <                 &(currentAtomType->sigma),
470 >                 &isCharge,
471 >                 &eamEpslon,
472 >                 &eamSigma,
473 >                 &charge,
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    
515 –
516 –  Atom* thisAtom;
517 –
538    for( i=0; i<nAtoms; i++ ){
539      
540      currentAtomType = headAtomType->find( the_atoms[i]->getType() );
# Line 528 | Line 548 | void EAM_FF::initializeAtoms( int nAtoms, Atom** the_a
548      
549      the_atoms[i]->setMass( currentAtomType->mass );
550      the_atoms[i]->setIdent( currentAtomType->ident );
531 –    the_atoms[i]->setEAM();
551  
552 +    if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut;
553 +    
554    }
555   }
556  
557 < void LJ_FF::initializeBonds( int nBonds, Bond** BondArray,
557 > void EAM_FF::initializeBonds( int nBonds, Bond** BondArray,
558                               bond_pair* the_bonds ){
559    
560      if( nBonds ){
# Line 657 | Line 678 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
678                       double **eam_rvals,
679                       double **eam_rhovals,
680                       double **eam_Frhovals){
681 <  
681 >  double* myEam_rvals;
682 >  double* myEam_rhovals;
683 >  double* myEam_Frhovals;
684 >
685    char* ffPath_env = "FORCE_PARAM_PATH";
686    char* ffPath;
687    char* the_token;
688 +  char* eam_eof_test;
689    FILE *eamFile;
690 +  const int BUFFERSIZE = 3000;
691 +
692    char temp[200];
693    int linenumber;
694    int nReadLines;
695 +  char eam_read_buffer[BUFFERSIZE];
696  
697 +
698    int i,j;
699  
700    linenumber = 0;
# Line 674 | Line 703 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
703    eamFile = fopen( eamPotFile, "r" );
704    
705    
706 <  if( frcFile == NULL ){
706 >  if( eamFile == NULL ){
707      
708        // next see if the force path enviorment variable is set
709      
# Line 693 | Line 722 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
722  
723      
724      
725 <    if( frcFile == NULL ){
725 >    if( eamFile == NULL ){
726        
727        sprintf( painCave.errMsg,
728                 "Error opening the EAM force parameter file: %s\n"
# Line 706 | Line 735 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
735    }
736  
737    // First line is a comment line, read and toss it....
738 <  eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
738 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
739    linenumber++;
740 <  if(eof_test == NULL){
740 >  if(eam_eof_test == NULL){
741      sprintf( painCave.errMsg,
742               "error in reading commment in %s\n", eamPotFile);
743      painCave.isFatal = 1;
# Line 718 | Line 747 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
747  
748  
749    // The Second line contains atomic number, atomic mass and a lattice constant
750 <  eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
750 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
751    linenumber++;
752 <  if(eof_test == NULL){
752 >  if(eam_eof_test == NULL){
753      sprintf( painCave.errMsg,
754               "error in reading Identifier line in %s\n", eamPotFile);
755      painCave.isFatal = 1;
# Line 730 | Line 759 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
759  
760  
761      
762 <  if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
762 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
763      sprintf( painCave.errMsg,
764               "Error parseing EAM ident  line in %s\n", eamPotFile );
765      painCave.isFatal = 1;
# Line 756 | Line 785 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
785    info.lattice_constant = atof( the_token);
786  
787    // Next line is nrho, drho, nr, dr and rcut
788 <  eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile)
789 <  if(eof_test == NULL){
788 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
789 >  if(eam_eof_test == NULL){
790      sprintf( painCave.errMsg,
791               "error in reading number of points line in %s\n", eamPotFile);
792      painCave.isFatal = 1;
793      simError();
794    }
795  
796 <  if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
796 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
797      sprintf( painCave.errMsg,
798               "Error parseing EAM nrho: line in %s\n", eamPotFile );
799      painCave.isFatal = 1;
# Line 806 | Line 835 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
835    info.eam_rcut = atof( the_token);
836  
837  
838 +
839 +
840 +
841    // Ok now we have to allocate point arrays and read in number of points
842    // Index the arrays for fortran, starting at 1
843 <  *eam_Frhovals = new double[info.nrho];
844 <  *eam_rvals    = new double[info.nr];
845 <  *eam_rhovals  = new double[info.nr];
843 >  myEam_Frhovals = new double[info.eam_nrho];
844 >  myEam_rvals    = new double[info.eam_nr];
845 >  myEam_rhovals  = new double[info.eam_nr];
846  
847    // Parse F of rho vals.
848  
849    // Assume for now that we have a complete number of lines
850 <  nReadLines = int(info.nrho/5);
850 >  nReadLines = int(info.eam_nrho/5);
851 >  
852  
853 +
854    for (i=0;i<nReadLines;i++){
855      j = i*5;
856 <    
856 >
857      // Read next line
858 <    eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
858 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
859      linenumber++;
860 <    if(eof_test == NULL){
860 >    if(eam_eof_test == NULL){
861        sprintf( painCave.errMsg,
862                 "error in reading EAM file %s at line %d\n",
863                 eamPotFile,linenumber);
# Line 833 | Line 867 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
867      
868      // Parse 5 values on each line into array
869      // Value 1
870 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
870 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
871        sprintf( painCave.errMsg,
872                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
873        painCave.isFatal = 1;
874        simError();
875      }
876 +
877 +    myEam_Frhovals[j+0] = atof( the_token );
878      
843 –    *eam_Frhovals[j+0] = atof( the_token );
844 –    
879      // Value 2
880      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
881        sprintf( painCave.errMsg,
# Line 849 | Line 883 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
883        painCave.isFatal = 1;
884        simError();
885      }
886 <  
887 <    *eam_Frhovals[j+1] = atof( the_token );
886 >
887 >    myEam_Frhovals[j+1] = atof( the_token );
888      
889      // Value 3
890      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 859 | Line 893 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
893        painCave.isFatal = 1;
894        simError();
895      }
896 +
897 +    myEam_Frhovals[j+2] = atof( the_token );
898      
863 –    *eam_Frhovals[j+2] = atof( the_token );
864 –    
899      // Value 4
900      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
901        sprintf( painCave.errMsg,
# Line 869 | Line 903 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
903        painCave.isFatal = 1;
904        simError();
905      }
872 –    
873 –    *eam_Frhovals[j+3] = atof( the_token );
906  
907 +    myEam_Frhovals[j+3] = atof( the_token );
908 +
909      // Value 5
910      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
911        sprintf( painCave.errMsg,
# Line 879 | Line 913 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
913        painCave.isFatal = 1;
914        simError();
915      }
916 +
917 +    myEam_Frhovals[j+4] = atof( the_token );
918      
883 –    *eam_Frhovals[j+4] = atof( the_token );
884 –    
919    }
920    // Parse Z of r vals
921    
922    // Assume for now that we have a complete number of lines
923 <  nReadLines = int(info.nr/5);
923 >  nReadLines = int(info.eam_nr/5);
924  
925    for (i=0;i<nReadLines;i++){
926      j = i*5;
927  
928      // Read next line
929 <    eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
929 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
930      linenumber++;
931 <    if(eof_test == NULL){
931 >    if(eam_eof_test == NULL){
932        sprintf( painCave.errMsg,
933                 "error in reading EAM file %s at line %d\n",
934                 eamPotFile,linenumber);
# Line 904 | Line 938 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
938      
939      // Parse 5 values on each line into array
940      // Value 1
941 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
941 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
942        sprintf( painCave.errMsg,
943                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
944        painCave.isFatal = 1;
945        simError();
946      }
947      
948 <    *eam_rvals[j+0] = atof( the_token );
948 >    myEam_rvals[j+0] = atof( the_token );
949  
950      // Value 2
951      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 921 | Line 955 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
955        simError();
956      }
957    
958 <    *eam_rvals[j+1] = atof( the_token );
958 >    myEam_rvals[j+1] = atof( the_token );
959  
960      // Value 3
961      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 931 | Line 965 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
965        simError();
966      }
967    
968 <    *eam_rvals[j+2] = atof( the_token );
968 >    myEam_rvals[j+2] = atof( the_token );
969  
970      // Value 4
971      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 941 | Line 975 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
975        simError();
976      }
977    
978 <    *eam_rvals[j+3] = atof( the_token );
978 >    myEam_rvals[j+3] = atof( the_token );
979  
980      // Value 5
981      if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
# Line 951 | Line 985 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
985        simError();
986      }
987    
988 <    *eam_rvals[j+4] = atof( the_token );
988 >    myEam_rvals[j+4] = atof( the_token );
989  
990    }
991    // Parse rho of r vals
# Line 962 | Line 996 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
996      j = i*5;
997  
998      // Read next line
999 <    eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile);
999 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
1000      linenumber++;
1001 <    if(eof_test == NULL){
1001 >    if(eam_eof_test == NULL){
1002        sprintf( painCave.errMsg,
1003                 "error in reading EAM file %s at line %d\n",
1004                 eamPotFile,linenumber);
# Line 974 | Line 1008 | int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil
1008    
1009      // Parse 5 values on each line into array
1010      // Value 1
1011 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1011 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
1012        sprintf( painCave.errMsg,
1013                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1014        painCave.isFatal = 1;
1015        simError();
1016      }
1017    
1018 <    *eam_rhovals[j+0] = atof( the_token );
1018 >    myEam_rhovals[j+0] = atof( the_token );
1019  
1020      // Value 2
1021 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1021 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1022        sprintf( painCave.errMsg,
1023                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1024        painCave.isFatal = 1;
1025        simError();
1026      }
1027    
1028 <    *eam_rhovals[j+1] = atof( the_token );
1028 >    myEam_rhovals[j+1] = atof( the_token );
1029  
1030      // Value 3
1031 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1031 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1032        sprintf( painCave.errMsg,
1033                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1034        painCave.isFatal = 1;
1035        simError();
1036      }
1037    
1038 <    *eam_rhovals[j+2] = atof( the_token );
1038 >    myEam_rhovals[j+2] = atof( the_token );
1039  
1040      // Value 4
1041 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1041 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1042        sprintf( painCave.errMsg,
1043                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1044        painCave.isFatal = 1;
1045        simError();
1046      }
1047    
1048 <    *eam_rhovals[j+3] = atof( the_token );
1048 >    myEam_rhovals[j+3] = atof( the_token );
1049  
1050      // Value 5
1051 <    if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){
1051 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1052        sprintf( painCave.errMsg,
1053                 "Error parseing EAM nrho: line in %s\n", eamPotFile );
1054        painCave.isFatal = 1;
1055        simError();
1056      }
1057    
1058 <    *eam_rhovals[j+4] = atof( the_token );
1059 <
1058 >    myEam_rhovals[j+4] = atof( the_token );
1059 >
1060    }
1061 +  *eam_rvals = myEam_rvals;
1062 +  *eam_rhovals = myEam_rhovals;
1063 +  *eam_Frhovals = myEam_Frhovals;
1064  
1065    fclose(eamFile);
1066    return 0;

Diff Legend

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