Unverified Commit fac7b864 authored by Richard Berger's avatar Richard Berger
Browse files

Update pair eam alloy

parent e08404f5

src/MANYBODY/pair_eam_alloy.cpp

+93 −80
Original line number Diff line number Diff line
@@ -24,11 +24,11 @@ 
#include "memory.h"

#include "error.h"

#include "utils.h"

#include "tokenizer.h"

#include "potential_file_reader.h"



using namespace LAMMPS_NS;



#define MAXLINE 1024



/* ---------------------------------------------------------------------- */



PairEAMAlloy::PairEAMAlloy(LAMMPS *lmp) : PairEAM(lmp)
@@ -61,7 +61,7 @@ void PairEAMAlloy::coeff(int narg, char **arg) 
  if (setfl) {

    for (i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i];

    delete [] setfl->elements;

    delete [] setfl->mass;

    memory->destroy(setfl->mass);

    memory->destroy(setfl->frho);

    memory->destroy(setfl->rhor);

    memory->destroy(setfl->z2r);
@@ -117,98 +117,111 @@ void PairEAMAlloy::read_file(char *filename) 
{

  Setfl *file = setfl;



  // open potential file

  // read potential file

  if(comm->me == 0) {

    PotentialFileReader reader(lmp, filename, "EAM");



  int me = comm->me;

  FILE *fptr;

  char line[MAXLINE];

    try {

      char * line = nullptr;



  if (me == 0) {

    fptr = force->open_potential(filename);

    if (fptr == NULL) {

      char str[128];

      snprintf(str,128,"Cannot open EAM potential file %s",filename);

      error->one(FLERR,str);

    }

  }

      reader.skip_line();

      reader.skip_line();

      reader.skip_line();



  // read and broadcast header

      // extract element names from nelements line

      line = reader.next_line(1);

      ValueTokenizer values(line);

      file->nelements = values.next_int();



  int n;

  if (me == 0) {

    utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);

    utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);

    utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);

    utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);

    n = strlen(line) + 1;

      if (values.count() != file->nelements + 1)

        error->one(FLERR,"Incorrect element names in EAM potential file");



      file->elements = new char*[file->nelements];

      for (int i = 0; i < file->nelements; i++) {

        const std::string word = values.next_string();

        const int n = word.length() + 1;

        file->elements[i] = new char[n];

        strcpy(file->elements[i], word.c_str());

      }

  MPI_Bcast(&n,1,MPI_INT,0,world);

  MPI_Bcast(line,n,MPI_CHAR,0,world);



  sscanf(line,"%d",&file->nelements);

  int nwords = utils::count_words(line);

  if (nwords != file->nelements + 1)

    error->all(FLERR,"Incorrect element names in EAM potential file");

      // 



  char **words = new char*[file->nelements+1];

  nwords = 0;

  strtok(line," \t\n\r\f");

  while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;

      line = reader.next_line(5);

      values = ValueTokenizer(line);

      file->nrho = values.next_int();

      file->drho = values.next_double();

      file->nr   = values.next_int();

      file->dr   = values.next_double();

      file->cut  = values.next_double();

      

      if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))

        error->one(FLERR,"Invalid EAM potential file");

      

      memory->create(file->mass, file->nelements, "pair:mass");

      memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");

      memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");

      memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");



  file->elements = new char*[file->nelements];

      for (int i = 0; i < file->nelements; i++) {

    n = strlen(words[i]) + 1;

    file->elements[i] = new char[n];

    strcpy(file->elements[i],words[i]);

        line = reader.next_line(2);

        values = ValueTokenizer(line);

        values.next_int(); // ignore

        file->mass[i] = values.next_double();



        reader.next_dvector(file->nrho, &file->frho[i][1]);

        reader.next_dvector(file->nr, &file->rhor[i][1]);

      }

  delete [] words;



  if (me == 0) {

    utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);

    nwords = sscanf(line,"%d %lg %d %lg %lg",

           &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);

      for (int i = 0; i < file->nelements; i++) {

        for (int j = 0; j <= i; j++) {

          reader.next_dvector(file->nr, &file->z2r[i][j][1]);

        }

      }

    } catch (TokenizerException & e) {

      error->one(FLERR, e.what());

    }

  }



  // broadcast potential information

  MPI_Bcast(&file->nelements, 1, MPI_INT, 0, world);



  MPI_Bcast(&nwords,1,MPI_INT,0,world);

  MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world);

  MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world);

  MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);

  MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);

  MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);



  if ((nwords != 5) || (file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))

    error->all(FLERR,"Invalid EAM potential file");



  file->mass = new double[file->nelements];

  // allocate memory on other procs

  if (comm->me != 0) {

    file->elements = new char*[file->nelements];

    for (int i = 0; i < file->nelements; i++) file->elements[i] = nullptr;

    memory->create(file->mass, file->nelements, "pair:mass");

    memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");

    memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");

  memory->create(file->z2r,file->nelements,file->nelements,file->nr+1,

                 "pair:z2r");



  int i,j,tmp;

  for (i = 0; i < file->nelements; i++) {

    if (me == 0) {

      utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);

      sscanf(line,"%d %lg",&tmp,&file->mass[i]);

    memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");

  }

    MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);



    if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]);

  // broadcast file->elements string array

  for (int i = 0; i < file->nelements; i++) {

    int n = strlen(file->elements[i]) + 1;

    MPI_Bcast(&n, 1, MPI_INT, 0, world);

    if (comm->me != 0) file->elements[i] = new char[n];

    MPI_Bcast(file->elements[i], n, MPI_CHAR, 0, world);

  }



  // broadcast file->mass, frho, rhor

  for (int i = 0; i < file->nelements; i++) {

    MPI_Bcast(&file->mass[i], 1, MPI_DOUBLE, 0, world);

    MPI_Bcast(&file->frho[i][1], file->nrho, MPI_DOUBLE, 0, world);

    if (me == 0) grab(fptr,file->nr,&file->rhor[i][1]);

    MPI_Bcast(&file->rhor[i][1], file->nr, MPI_DOUBLE, 0, world);

  }



  for (i = 0; i < file->nelements; i++)

    for (j = 0; j <= i; j++) {

      if (me == 0) grab(fptr,file->nr,&file->z2r[i][j][1]);

  // broadcast file->z2r

  for (int i = 0; i < file->nelements; i++) {

    for (int j = 0; j <= i; j++) {

      MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world);

    }



  // close the potential file



  if (me == 0) fclose(fptr);

  }

}



/* ----------------------------------------------------------------------

CRA Git | Maintained and supported by SUSTech CRA and CCSE