From 609842916cc08041127da136b0a767e25f42d232 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Jul 2013 23:52:14 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10392 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-OMP/pair_eam_alloy_omp.cpp.orig | 323 -------------------- src/USER-OMP/pair_eam_fs_omp.cpp.orig | 332 --------------------- src/USER-OMP/pair_tersoff_zbl_omp.cpp.orig | 297 ------------------ 3 files changed, 952 deletions(-) delete mode 100644 src/USER-OMP/pair_eam_alloy_omp.cpp.orig delete mode 100644 src/USER-OMP/pair_eam_fs_omp.cpp.orig delete mode 100644 src/USER-OMP/pair_tersoff_zbl_omp.cpp.orig diff --git a/src/USER-OMP/pair_eam_alloy_omp.cpp.orig b/src/USER-OMP/pair_eam_alloy_omp.cpp.orig deleted file mode 100644 index 3dff868b6a..0000000000 --- a/src/USER-OMP/pair_eam_alloy_omp.cpp.orig +++ /dev/null @@ -1,323 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Stephen Foiles (SNL), Murray Daw (SNL) -------------------------------------------------------------------------- */ - -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "pair_eam_alloy_omp.h" -#include "atom.h" -#include "comm.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define MAXLINE 1024 - -/* ---------------------------------------------------------------------- */ - -PairEAMAlloyOMP::PairEAMAlloyOMP(LAMMPS *lmp) : PairEAMOMP(lmp) -{ - one_coeff = 1; -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs - read DYNAMO setfl file -------------------------------------------------------------------------- */ - -void PairEAMAlloyOMP::coeff(int narg, char **arg) -{ - int i,j; - - if (!allocated) allocate(); - - if (narg != 3 + atom->ntypes) - error->all(FLERR,"Incorrect args for pair coefficients"); - - // insure I,J args are * * - - if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) - error->all(FLERR,"Incorrect args for pair coefficients"); - - // read EAM setfl file - - if (setfl) { - for (i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i]; - delete [] setfl->elements; - delete [] setfl->mass; - memory->destroy(setfl->frho); - memory->destroy(setfl->rhor); - memory->destroy(setfl->z2r); - delete setfl; - } - setfl = new Setfl(); - read_file(arg[2]); - - // read args that map atom types to elements in potential file - // map[i] = which element the Ith atom type is, -1 if NULL - - for (i = 3; i < narg; i++) { - if (strcmp(arg[i],"NULL") == 0) { - map[i-2] = -1; - continue; - } - for (j = 0; j < setfl->nelements; j++) - if (strcmp(arg[i],setfl->elements[j]) == 0) break; - if (j < setfl->nelements) map[i-2] = j; - else error->all(FLERR,"No matching element in EAM potential file"); - } - - // clear setflag since coeff() called once with I,J = * * - - int n = atom->ntypes; - for (i = 1; i <= n; i++) - for (j = i; j <= n; j++) - setflag[i][j] = 0; - - // set setflag i,j for type pairs where both are mapped to elements - // set mass of atom type if i = j - - int count = 0; - for (i = 1; i <= n; i++) { - for (j = i; j <= n; j++) { - if (map[i] >= 0 && map[j] >= 0) { - setflag[i][j] = 1; - if (i == j) atom->set_mass(i,setfl->mass[map[i]]); - count++; - } - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - read a multi-element DYNAMO setfl file -------------------------------------------------------------------------- */ - -void PairEAMAlloyOMP::read_file(char *filename) -{ - Setfl *file = setfl; - - // open potential file - - int me = comm->me; - FILE *fptr; - char line[MAXLINE]; - - if (me == 0) { - fptr = fopen(filename,"r"); - if (fptr == NULL) { - char str[128]; - sprintf(str,"Cannot open EAM potential file %s",filename); - error->one(FLERR,str); - } - } - - // read and broadcast header - // extract element names from nelements line - - int n; - if (me == 0) { - fgets(line,MAXLINE,fptr); - fgets(line,MAXLINE,fptr); - fgets(line,MAXLINE,fptr); - fgets(line,MAXLINE,fptr); - n = strlen(line) + 1; - } - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - - sscanf(line,"%d",&file->nelements); - int nwords = atom->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; - - 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]); - } - delete [] words; - - if (me == 0) { - fgets(line,MAXLINE,fptr); - sscanf(line,"%d %lg %d %lg %lg", - &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut); - } - - 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); - - file->mass = new double[file->nelements]; - 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) { - fgets(line,MAXLINE,fptr); - sscanf(line,"%d %lg",&tmp,&file->mass[i]); - } - MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world); - - if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]); - 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]); - MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world); - } - - // close the potential file - - if (me == 0) fclose(fptr); -} - -/* ---------------------------------------------------------------------- - copy read-in setfl potential to standard array format -------------------------------------------------------------------------- */ - -void PairEAMAlloyOMP::file2array() -{ - int i,j,m,n; - int ntypes = atom->ntypes; - - // set function params directly from setfl file - - nrho = setfl->nrho; - nr = setfl->nr; - drho = setfl->drho; - dr = setfl->dr; - - // ------------------------------------------------------------------ - // setup frho arrays - // ------------------------------------------------------------------ - - // allocate frho arrays - // nfrho = # of setfl elements + 1 for zero array - - nfrho = setfl->nelements + 1; - memory->destroy(frho); - memory->create(frho,nfrho,nrho+1,"pair:frho"); - - // copy each element's frho to global frho - - for (i = 0; i < setfl->nelements; i++) - for (m = 1; m <= nrho; m++) frho[i][m] = setfl->frho[i][m]; - - // add extra frho of zeroes for non-EAM types to point to (pair hybrid) - // this is necessary b/c fp is still computed for non-EAM atoms - - for (m = 1; m <= nrho; m++) frho[nfrho-1][m] = 0.0; - - // type2frho[i] = which frho array (0 to nfrho-1) each atom type maps to - // if atom type doesn't point to element (non-EAM atom in pair hybrid) - // then map it to last frho array of zeroes - - for (i = 1; i <= ntypes; i++) - if (map[i] >= 0) type2frho[i] = map[i]; - else type2frho[i] = nfrho-1; - - // ------------------------------------------------------------------ - // setup rhor arrays - // ------------------------------------------------------------------ - - // allocate rhor arrays - // nrhor = # of setfl elements - - nrhor = setfl->nelements; - memory->destroy(rhor); - memory->create(rhor,nrhor,nr+1,"pair:rhor"); - - // copy each element's rhor to global rhor - - for (i = 0; i < setfl->nelements; i++) - for (m = 1; m <= nr; m++) rhor[i][m] = setfl->rhor[i][m]; - - // type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to - // for setfl files, I,J mapping only depends on I - // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used - - for (i = 1; i <= ntypes; i++) - for (j = 1; j <= ntypes; j++) - type2rhor[i][j] = map[i]; - - // ------------------------------------------------------------------ - // setup z2r arrays - // ------------------------------------------------------------------ - - // allocate z2r arrays - // nz2r = N*(N+1)/2 where N = # of setfl elements - - nz2r = setfl->nelements * (setfl->nelements+1) / 2; - memory->destroy(z2r); - memory->create(z2r,nz2r,nr+1,"pair:z2r"); - - // copy each element pair z2r to global z2r, only for I >= J - - n = 0; - for (i = 0; i < setfl->nelements; i++) - for (j = 0; j <= i; j++) { - for (m = 1; m <= nr; m++) z2r[n][m] = setfl->z2r[i][j][m]; - n++; - } - - // type2z2r[i][j] = which z2r array (0 to nz2r-1) each type pair maps to - // set of z2r arrays only fill lower triangular Nelement matrix - // value = n = sum over rows of lower-triangular matrix until reach irow,icol - // swap indices when irow < icol to stay lower triangular - // if map = -1 (non-EAM atom in pair hybrid): - // type2z2r is not used by non-opt - // but set type2z2r to 0 since accessed by opt - - int irow,icol; - for (i = 1; i <= ntypes; i++) { - for (j = 1; j <= ntypes; j++) { - irow = map[i]; - icol = map[j]; - if (irow == -1 || icol == -1) { - type2z2r[i][j] = 0; - continue; - } - if (irow < icol) { - irow = map[j]; - icol = map[i]; - } - n = 0; - for (m = 0; m < irow; m++) n += m + 1; - n += icol; - type2z2r[i][j] = n; - } - } -} diff --git a/src/USER-OMP/pair_eam_fs_omp.cpp.orig b/src/USER-OMP/pair_eam_fs_omp.cpp.orig deleted file mode 100644 index 84101913c5..0000000000 --- a/src/USER-OMP/pair_eam_fs_omp.cpp.orig +++ /dev/null @@ -1,332 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Tim Lau (MIT) -------------------------------------------------------------------------- */ - -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "pair_eam_fs_omp.h" -#include "atom.h" -#include "comm.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define MAXLINE 1024 - -/* ---------------------------------------------------------------------- */ - -PairEAMFSOMP::PairEAMFSOMP(LAMMPS *lmp) : PairEAMOMP(lmp) -{ - one_coeff = 1; -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs - read EAM Finnis-Sinclair file -------------------------------------------------------------------------- */ - -void PairEAMFSOMP::coeff(int narg, char **arg) -{ - int i,j; - - if (!allocated) allocate(); - - if (narg != 3 + atom->ntypes) - error->all(FLERR,"Incorrect args for pair coefficients"); - - // insure I,J args are * * - - if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) - error->all(FLERR,"Incorrect args for pair coefficients"); - - // read EAM Finnis-Sinclair file - - if (fs) { - for (i = 0; i < fs->nelements; i++) delete [] fs->elements[i]; - delete [] fs->elements; - delete [] fs->mass; - memory->destroy(fs->frho); - memory->destroy(fs->rhor); - memory->destroy(fs->z2r); - delete fs; - } - fs = new Fs(); - read_file(arg[2]); - - // read args that map atom types to elements in potential file - // map[i] = which element the Ith atom type is, -1 if NULL - - for (i = 3; i < narg; i++) { - if (strcmp(arg[i],"NULL") == 0) { - map[i-2] = -1; - continue; - } - for (j = 0; j < fs->nelements; j++) - if (strcmp(arg[i],fs->elements[j]) == 0) break; - if (j < fs->nelements) map[i-2] = j; - else error->all(FLERR,"No matching element in EAM potential file"); - } - - // clear setflag since coeff() called once with I,J = * * - - int n = atom->ntypes; - for (i = 1; i <= n; i++) - for (j = i; j <= n; j++) - setflag[i][j] = 0; - - // set setflag i,j for type pairs where both are mapped to elements - // set mass of atom type if i = j - - int count = 0; - for (i = 1; i <= n; i++) { - for (j = i; j <= n; j++) { - if (map[i] >= 0 && map[j] >= 0) { - setflag[i][j] = 1; - if (i == j) atom->set_mass(i,fs->mass[map[i]]); - count++; - } - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - read a multi-element DYNAMO setfl file -------------------------------------------------------------------------- */ - -void PairEAMFSOMP::read_file(char *filename) -{ - Fs *file = fs; - - // open potential file - - int me = comm->me; - FILE *fptr; - char line[MAXLINE]; - - if (me == 0) { - fptr = fopen(filename,"r"); - if (fptr == NULL) { - char str[128]; - sprintf(str,"Cannot open EAM potential file %s",filename); - error->one(FLERR,str); - } - } - - // read and broadcast header - // extract element names from nelements line - - int n; - if (me == 0) { - fgets(line,MAXLINE,fptr); - fgets(line,MAXLINE,fptr); - fgets(line,MAXLINE,fptr); - fgets(line,MAXLINE,fptr); - n = strlen(line) + 1; - } - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - - sscanf(line,"%d",&file->nelements); - int nwords = atom->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; - - 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]); - } - delete [] words; - - if (me == 0) { - fgets(line,MAXLINE,fptr); - sscanf(line,"%d %lg %d %lg %lg", - &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut); - } - - 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); - - file->mass = new double[file->nelements]; - memory->create(file->frho,file->nelements,file->nrho+1, - "pair:frho"); - memory->create(file->rhor,file->nelements,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) { - fgets(line,MAXLINE,fptr); - sscanf(line,"%d %lg",&tmp,&file->mass[i]); - } - MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world); - - if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]); - MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world); - - for (j = 0; j < file->nelements; j++) { - if (me == 0) grab(fptr,file->nr,&file->rhor[i][j][1]); - MPI_Bcast(&file->rhor[i][j][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]); - MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world); - } - - // close the potential file - - if (me == 0) fclose(fptr); -} - -/* ---------------------------------------------------------------------- - copy read-in setfl potential to standard array format -------------------------------------------------------------------------- */ - -void PairEAMFSOMP::file2array() -{ - int i,j,m,n; - int ntypes = atom->ntypes; - - // set function params directly from fs file - - nrho = fs->nrho; - nr = fs->nr; - drho = fs->drho; - dr = fs->dr; - - // ------------------------------------------------------------------ - // setup frho arrays - // ------------------------------------------------------------------ - - // allocate frho arrays - // nfrho = # of fs elements + 1 for zero array - - nfrho = fs->nelements + 1; - memory->destroy(frho); - memory->create(frho,nfrho,nrho+1,"pair:frho"); - - // copy each element's frho to global frho - - for (i = 0; i < fs->nelements; i++) - for (m = 1; m <= nrho; m++) frho[i][m] = fs->frho[i][m]; - - // add extra frho of zeroes for non-EAM types to point to (pair hybrid) - // this is necessary b/c fp is still computed for non-EAM atoms - - for (m = 1; m <= nrho; m++) frho[nfrho-1][m] = 0.0; - - // type2frho[i] = which frho array (0 to nfrho-1) each atom type maps to - // if atom type doesn't point to element (non-EAM atom in pair hybrid) - // then map it to last frho array of zeroes - - for (i = 1; i <= ntypes; i++) - if (map[i] >= 0) type2frho[i] = map[i]; - else type2frho[i] = nfrho-1; - - // ------------------------------------------------------------------ - // setup rhor arrays - // ------------------------------------------------------------------ - - // allocate rhor arrays - // nrhor = square of # of fs elements - - nrhor = fs->nelements * fs->nelements; - memory->destroy(rhor); - memory->create(rhor,nrhor,nr+1,"pair:rhor"); - - // copy each element pair rhor to global rhor - - n = 0; - for (i = 0; i < fs->nelements; i++) - for (j = 0; j < fs->nelements; j++) { - for (m = 1; m <= nr; m++) rhor[n][m] = fs->rhor[i][j][m]; - n++; - } - - // type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to - // for fs files, there is a full NxN set of rhor arrays - // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used - - for (i = 1; i <= ntypes; i++) - for (j = 1; j <= ntypes; j++) - type2rhor[i][j] = map[i] * fs->nelements + map[j]; - - // ------------------------------------------------------------------ - // setup z2r arrays - // ------------------------------------------------------------------ - - // allocate z2r arrays - // nz2r = N*(N+1)/2 where N = # of fs elements - - nz2r = fs->nelements * (fs->nelements+1) / 2; - memory->destroy(z2r); - memory->create(z2r,nz2r,nr+1,"pair:z2r"); - - // copy each element pair z2r to global z2r, only for I >= J - - n = 0; - for (i = 0; i < fs->nelements; i++) - for (j = 0; j <= i; j++) { - for (m = 1; m <= nr; m++) z2r[n][m] = fs->z2r[i][j][m]; - n++; - } - - // type2z2r[i][j] = which z2r array (0 to nz2r-1) each type pair maps to - // set of z2r arrays only fill lower triangular Nelement matrix - // value = n = sum over rows of lower-triangular matrix until reach irow,icol - // swap indices when irow < icol to stay lower triangular - // if map = -1 (non-EAM atom in pair hybrid): - // type2z2r is not used by non-opt - // but set type2z2r to 0 since accessed by opt - - int irow,icol; - for (i = 1; i <= ntypes; i++) { - for (j = 1; j <= ntypes; j++) { - irow = map[i]; - icol = map[j]; - if (irow == -1 || icol == -1) { - type2z2r[i][j] = 0; - continue; - } - if (irow < icol) { - irow = map[j]; - icol = map[i]; - } - n = 0; - for (m = 0; m < irow; m++) n += m + 1; - n += icol; - type2z2r[i][j] = n; - } - } -} diff --git a/src/USER-OMP/pair_tersoff_zbl_omp.cpp.orig b/src/USER-OMP/pair_tersoff_zbl_omp.cpp.orig deleted file mode 100644 index 9c0f397118..0000000000 --- a/src/USER-OMP/pair_tersoff_zbl_omp.cpp.orig +++ /dev/null @@ -1,297 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Aidan Thompson (SNL) - original Tersoff implementation - David Farrell (NWU) - ZBL addition -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "pair_tersoff_zbl_omp.h" -#include "atom.h" -#include "update.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "force.h" -#include "comm.h" -#include "memory.h" -#include "error.h" - -#include "math_const.h" -#include "math_special.h" -using namespace LAMMPS_NS; -using namespace MathConst; -using namespace MathSpecial; - -#define MAXLINE 1024 -#define DELTA 4 - -/* ---------------------------------------------------------------------- - Fermi-like smoothing function -------------------------------------------------------------------------- */ - -static double F_fermi(const double r, const double expsc, const double cut) -{ - return 1.0 / (1.0 + exp(-expsc*(r-cut))); -} - -/* ---------------------------------------------------------------------- - Fermi-like smoothing function derivative with respect to r -------------------------------------------------------------------------- */ - -static double F_fermi_d(const double r, const double expsc, const double cut) -{ - return expsc*exp(-expsc*(r-cut)) / square(1.0 + exp(-expsc*(r-cut))); -} - -/* ---------------------------------------------------------------------- */ - -PairTersoffZBLOMP::PairTersoffZBLOMP(LAMMPS *lmp) : PairTersoffOMP(lmp) -{ - // hard-wired constants in metal or real units - // a0 = Bohr radius - // epsilon0 = permittivity of vacuum = q / energy-distance units - // e = unit charge - // 1 Kcal/mole = 0.043365121 eV - - if (strcmp(update->unit_style,"metal") == 0) { - global_a_0 = 0.529; - global_epsilon_0 = 0.00552635; - global_e = 1.0; - } else if (strcmp(update->unit_style,"real") == 0) { - global_a_0 = 0.529; - global_epsilon_0 = 0.00552635 * 0.043365121; - global_e = 1.0; - } else error->all(FLERR,"Pair tersoff/zbl requires metal or real units"); -} - -/* ---------------------------------------------------------------------- */ - -void PairTersoffZBLOMP::read_file(char *file) -{ - int params_per_line = 21; - char **words = new char*[params_per_line+1]; - - memory->sfree(params); - params = NULL; - nparams = 0; - - // open file on proc 0 - - FILE *fp; - if (comm->me == 0) { - fp = fopen(file,"r"); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open Tersoff potential file %s",file); - error->one(FLERR,str); - } - } - - // read each line out of file, skipping blank lines or leading '#' - // store line of params if all 3 element tags are in element list - - int n,nwords,ielement,jelement,kelement; - char line[MAXLINE],*ptr; - int eof = 0; - - while (1) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fp); - if (ptr == NULL) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; - } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - - // strip comment, skip line if blank - - if (ptr = strchr(line,'#')) *ptr = '\0'; - nwords = atom->count_words(line); - if (nwords == 0) continue; - - // concatenate additional lines until have params_per_line words - - while (nwords < params_per_line) { - n = strlen(line); - if (comm->me == 0) { - ptr = fgets(&line[n],MAXLINE-n,fp); - if (ptr == NULL) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; - } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - if (ptr = strchr(line,'#')) *ptr = '\0'; - nwords = atom->count_words(line); - } - - if (nwords != params_per_line) - error->all(FLERR,"Incorrect format in Tersoff potential file"); - - // words = ptrs to all words in line - - nwords = 0; - words[nwords++] = strtok(line," \t\n\r\f"); - while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; - - // ielement,jelement,kelement = 1st args - // if all 3 args are in element list, then parse this line - // else skip to next line - - for (ielement = 0; ielement < nelements; ielement++) - if (strcmp(words[0],elements[ielement]) == 0) break; - if (ielement == nelements) continue; - for (jelement = 0; jelement < nelements; jelement++) - if (strcmp(words[1],elements[jelement]) == 0) break; - if (jelement == nelements) continue; - for (kelement = 0; kelement < nelements; kelement++) - if (strcmp(words[2],elements[kelement]) == 0) break; - if (kelement == nelements) continue; - - // load up parameter settings and error check their values - - if (nparams == maxparam) { - maxparam += DELTA; - params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), - "pair:params"); - } - - params[nparams].ielement = ielement; - params[nparams].jelement = jelement; - params[nparams].kelement = kelement; - params[nparams].powerm = atof(words[3]); - params[nparams].gamma = atof(words[4]); - params[nparams].lam3 = atof(words[5]); - params[nparams].c = atof(words[6]); - params[nparams].d = atof(words[7]); - params[nparams].h = atof(words[8]); - params[nparams].powern = atof(words[9]); - params[nparams].beta = atof(words[10]); - params[nparams].lam2 = atof(words[11]); - params[nparams].bigb = atof(words[12]); - params[nparams].bigr = atof(words[13]); - params[nparams].bigd = atof(words[14]); - params[nparams].lam1 = atof(words[15]); - params[nparams].biga = atof(words[16]); - params[nparams].Z_i = atof(words[17]); - params[nparams].Z_j = atof(words[18]); - params[nparams].ZBLcut = atof(words[19]); - params[nparams].ZBLexpscale = atof(words[20]); - - // currently only allow m exponent of 1 or 3 - - params[nparams].powermint = int(params[nparams].powerm); - - if ( - params[nparams].lam3 < 0.0 || params[nparams].c < 0.0 || - params[nparams].d < 0.0 || params[nparams].powern < 0.0 || - params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 || - params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 || - params[nparams].bigd < 0.0 || - params[nparams].bigd > params[nparams].bigr || - params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 || - params[nparams].powerm - params[nparams].powermint != 0.0 || - (params[nparams].powermint != 3 && params[nparams].powermint != 1) || - params[nparams].gamma < 0.0 || - params[nparams].Z_i < 1.0 || params[nparams].Z_j < 1.0 || - params[nparams].ZBLcut < 0.0 || params[nparams].ZBLexpscale < 0.0) - error->all(FLERR,"Illegal Tersoff parameter"); - - nparams++; - } - - delete [] words; -} - -/* ---------------------------------------------------------------------- */ - -void PairTersoffZBLOMP::force_zeta(Param *param, double rsq, double zeta_ij, - double &fforce, double &prefactor, - int eflag, double &eng) -{ - double r,fa,fa_d,bij; - - r = sqrt(rsq); - - fa = (r > param->bigr + param->bigd) ? 0.0 : - -param->bigb * exp(-param->lam2 * r) * ters_fc(r,param) * - F_fermi(r,param->ZBLexpscale,param->ZBLcut); - - fa_d = (r > param->bigr + param->bigd) ? 0.0 : - param->bigb * exp(-param->lam2 * r) * - (param->lam2 * ters_fc(r,param) * - F_fermi(r,param->ZBLexpscale,param->ZBLcut) - - ters_fc_d(r,param) * F_fermi(r,param->ZBLexpscale,param->ZBLcut) - - ters_fc(r,param) * F_fermi_d(r,param->ZBLexpscale,param->ZBLcut)); - - bij = ters_bij(zeta_ij,param); - fforce = 0.5*bij*fa_d / r; - prefactor = -0.5*fa * ters_bij_d(zeta_ij,param); - if (eflag) eng = 0.5*bij*fa; -} - -/* ---------------------------------------------------------------------- */ - -void PairTersoffZBLOMP::repulsive(Param *param, double rsq, double &fforce, - int eflag, double &eng) -{ - double r,tmp_fc,tmp_fc_d,tmp_exp; - - // Tersoff repulsive portion - - r = sqrt(rsq); - tmp_fc = ters_fc(r,param); - tmp_fc_d = ters_fc_d(r,param); - tmp_exp = exp(-param->lam1 * r); - double fforce_ters = param->biga * tmp_exp * (tmp_fc_d - tmp_fc*param->lam1); - double eng_ters = tmp_fc * param->biga * tmp_exp; - - // ZBL repulsive portion - - double esq = square(global_e); - double a_ij = (0.8854*global_a_0) / - (pow(param->Z_i,0.23) + pow(param->Z_j,0.23)); - double premult = (param->Z_i * param->Z_j * esq)/(4.0*MY_PI*global_epsilon_0); - double r_ov_a = r/a_ij; - double phi = 0.1818*exp(-3.2*r_ov_a) + 0.5099*exp(-0.9423*r_ov_a) + - 0.2802*exp(-0.4029*r_ov_a) + 0.02817*exp(-0.2016*r_ov_a); - double dphi = (1.0/a_ij) * (-3.2*0.1818*exp(-3.2*r_ov_a) - - 0.9423*0.5099*exp(-0.9423*r_ov_a) - - 0.4029*0.2802*exp(-0.4029*r_ov_a) - - 0.2016*0.02817*exp(-0.2016*r_ov_a)); - double fforce_ZBL = premult*-rsq* phi + premult/r*dphi; - double eng_ZBL = premult/r*phi; - - // combine two parts with smoothing by Fermi-like function - - fforce = -(-F_fermi_d(r,param->ZBLexpscale,param->ZBLcut) * eng_ZBL + - (1.0 - F_fermi(r,param->ZBLexpscale,param->ZBLcut))*fforce_ZBL + - F_fermi_d(r,param->ZBLexpscale,param->ZBLcut)*eng_ters + - F_fermi(r,param->ZBLexpscale,param->ZBLcut)*fforce_ters) / r; - - if (eflag) - eng = (1.0 - F_fermi(r,param->ZBLexpscale,param->ZBLcut))*eng_ZBL + - F_fermi(r,param->ZBLexpscale,param->ZBLcut)*eng_ters; -}