ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

enable and apply clang-format, remove redundant code, improve error message for EAM potential pair styles

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-03-25 06:48:59 -04:00
parent 7edb42f190
commit 22a1cf935e
16 changed files with 385 additions and 424 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
src/GPU/pair_eam_alloy_gpu.cpp
View File

@ -138,8 +138,7 @@ void PairEAMAlloyGPU::compute(int eflag, int vflag)
eam_alloy_gpu_compute_force(nullptr, eflag, vflag, eflag_atom, vflag_atom);
else
eam_alloy_gpu_compute_force(ilist, eflag, vflag, eflag_atom, vflag_atom);
if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0)
neighbor->build_topology();
if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0) neighbor->build_topology();
}
/* ----------------------------------------------------------------------
@ -282,12 +281,8 @@ void PairEAMAlloyGPU::coeff(int narg, char **arg)
if (!allocated) allocate();
if (narg != 3 + atom->ntypes) error->all(FLERR, "Incorrect args for pair coefficients");
// ensure I,J args are * *
if (strcmp(arg[0], "*") != 0 || strcmp(arg[1], "*") != 0)
error->all(FLERR, "Incorrect args for pair coefficients");
if (narg != 3 + atom->ntypes)
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM setfl file

11
src/GPU/pair_eam_fs_gpu.cpp
View File

@ -138,8 +138,7 @@ void PairEAMFSGPU::compute(int eflag, int vflag)
eam_fs_gpu_compute_force(nullptr, eflag, vflag, eflag_atom, vflag_atom);
else
eam_fs_gpu_compute_force(ilist, eflag, vflag, eflag_atom, vflag_atom);
if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0)
neighbor->build_topology();
if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0) neighbor->build_topology();
}
/* ----------------------------------------------------------------------
@ -282,12 +281,8 @@ void PairEAMFSGPU::coeff(int narg, char **arg)
if (!allocated) allocate();
if (narg != 3 + atom->ntypes) error->all(FLERR, "Incorrect args for pair coefficients");
// ensure I,J args are * *
if (strcmp(arg[0], "*") != 0 || strcmp(arg[1], "*") != 0)
error->all(FLERR, "Incorrect args for pair coefficients");
if (narg != 3 + atom->ntypes)
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM Finnis-Sinclair file

3
src/GPU/pair_eam_gpu.cpp
View File

@ -136,8 +136,7 @@ void PairEAMGPU::compute(int eflag, int vflag)
eam_gpu_compute_force(nullptr, eflag, vflag, eflag_atom, vflag_atom);
else
eam_gpu_compute_force(ilist, eflag, vflag, eflag_atom, vflag_atom);
if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0)
neighbor->build_topology();
if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0) neighbor->build_topology();
}
/* ----------------------------------------------------------------------

89
src/INTEL/pair_eam_alloy_intel.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -43,24 +42,19 @@ PairEAMAlloyIntel::PairEAMAlloyIntel(LAMMPS *lmp) : PairEAMIntel(lmp)
void PairEAMAlloyIntel::coeff(int narg, char **arg)
{
int i,j;
int i, j;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// ensure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM setfl file
if (setfl) {
for (i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i];
delete [] setfl->elements;
delete [] setfl->mass;
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);
@ -73,22 +67,23 @@ void PairEAMAlloyIntel::coeff(int narg, char **arg)
// 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;
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");
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;
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
@ -98,14 +93,14 @@ void PairEAMAlloyIntel::coeff(int narg, char **arg)
for (j = i; j <= n; j++) {
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(FLERR,i,setfl->mass[map[i]]);
if (i == j) atom->set_mass(FLERR, i, setfl->mass[map[i]]);
count++;
}
scale[i][j] = 1.0;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -123,8 +118,7 @@ void PairEAMAlloyIntel::read_file(char *filename)
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
double conversion_factor = utils::get_conversion_factor(utils::ENERGY, unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -135,21 +129,21 @@ void PairEAMAlloyIntel::read_file(char *filename)
file->nelements = values.next_int();
if (values.count() != file->nelements + 1)
error->one(FLERR,"Incorrect element names in EAM potential file");
error->one(FLERR, "Incorrect element names in EAM potential file");
file->elements = new char*[file->nelements];
file->elements = new char *[file->nelements];
for (int i = 0; i < file->nelements; i++)
file->elements[i] = utils::strdup(values.next_string());
values = reader.next_values(5);
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();
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");
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");
@ -158,14 +152,13 @@ void PairEAMAlloyIntel::read_file(char *filename)
for (int i = 0; i < file->nelements; i++) {
values = reader.next_values(2);
values.next_int(); // ignore
values.next_int(); // ignore
file->mass[i] = values.next_double();
reader.next_dvector(&file->frho[i][1], file->nrho);
reader.next_dvector(&file->rhor[i][1], file->nr);
if (unit_convert) {
for (int j = 1; j < file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
for (int j = 1; j < file->nrho; ++j) file->frho[i][j] *= conversion_factor;
}
}
@ -173,8 +166,7 @@ void PairEAMAlloyIntel::read_file(char *filename)
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k < file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
for (int k = 1; k < file->nr; ++k) file->z2r[i][j][k] *= conversion_factor;
}
}
}
@ -194,7 +186,7 @@ void PairEAMAlloyIntel::read_file(char *filename)
// allocate memory on other procs
if (comm->me != 0) {
file->elements = new char*[file->nelements];
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");
@ -220,9 +212,7 @@ void PairEAMAlloyIntel::read_file(char *filename)
// 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);
}
for (int j = 0; j <= i; j++) { MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world); }
}
}
@ -232,7 +222,7 @@ void PairEAMAlloyIntel::read_file(char *filename)
void PairEAMAlloyIntel::file2array()
{
int i,j,m,n;
int i, j, m, n;
int ntypes = atom->ntypes;
// set function params directly from setfl file
@ -241,7 +231,7 @@ void PairEAMAlloyIntel::file2array()
nr = setfl->nr;
drho = setfl->drho;
dr = setfl->dr;
rhomax = (nrho-1) * drho;
rhomax = (nrho - 1) * drho;
// ------------------------------------------------------------------
// setup frho arrays
@ -252,7 +242,7 @@ void PairEAMAlloyIntel::file2array()
nfrho = setfl->nelements + 1;
memory->destroy(frho);
memory->create(frho,nfrho,nrho+1,"pair:frho");
memory->create(frho, nfrho, nrho + 1, "pair:frho");
// copy each element's frho to global frho
@ -262,15 +252,17 @@ void PairEAMAlloyIntel::file2array()
// 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;
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;
if (map[i] >= 0)
type2frho[i] = map[i];
else
type2frho[i] = nfrho - 1;
// ------------------------------------------------------------------
// setup rhor arrays
@ -281,7 +273,7 @@ void PairEAMAlloyIntel::file2array()
nrhor = setfl->nelements;
memory->destroy(rhor);
memory->create(rhor,nrhor,nr+1,"pair:rhor");
memory->create(rhor, nrhor, nr + 1, "pair:rhor");
// copy each element's rhor to global rhor
@ -293,8 +285,7 @@ void PairEAMAlloyIntel::file2array()
// 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];
for (j = 1; j <= ntypes; j++) type2rhor[i][j] = map[i];
// ------------------------------------------------------------------
// setup z2r arrays
@ -303,9 +294,9 @@ void PairEAMAlloyIntel::file2array()
// allocate z2r arrays
// nz2r = N*(N+1)/2 where N = # of setfl elements
nz2r = setfl->nelements * (setfl->nelements+1) / 2;
nz2r = setfl->nelements * (setfl->nelements + 1) / 2;
memory->destroy(z2r);
memory->create(z2r,nz2r,nr+1,"pair:z2r");
memory->create(z2r, nz2r, nr + 1, "pair:z2r");
// copy each element pair z2r to global z2r, only for I >= J
@ -324,7 +315,7 @@ void PairEAMAlloyIntel::file2array()
// type2z2r is not used by non-opt
// but set type2z2r to 0 since accessed by opt
int irow,icol;
int irow, icol;
for (i = 1; i <= ntypes; i++) {
for (j = 1; j <= ntypes; j++) {
irow = map[i];

89
src/INTEL/pair_eam_fs_intel.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -43,24 +42,19 @@ PairEAMFSIntel::PairEAMFSIntel(LAMMPS *lmp) : PairEAMIntel(lmp)
void PairEAMFSIntel::coeff(int narg, char **arg)
{
int i,j;
int i, j;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// ensure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM Finnis-Sinclair file
if (fs) {
for (i = 0; i < fs->nelements; i++) delete [] fs->elements[i];
delete [] fs->elements;
delete [] fs->mass;
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);
@ -73,22 +67,23 @@ void PairEAMFSIntel::coeff(int narg, char **arg)
// 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;
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");
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;
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
@ -98,14 +93,14 @@ void PairEAMFSIntel::coeff(int narg, char **arg)
for (j = i; j <= n; j++) {
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(FLERR,i,fs->mass[map[i]]);
if (i == j) atom->set_mass(FLERR, i, fs->mass[map[i]]);
count++;
}
scale[i][j] = 1.0;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -123,8 +118,7 @@ void PairEAMFSIntel::read_file(char *filename)
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
double conversion_factor = utils::get_conversion_factor(utils::ENERGY, unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -135,21 +129,21 @@ void PairEAMFSIntel::read_file(char *filename)
file->nelements = values.next_int();
if (values.count() != file->nelements + 1)
error->one(FLERR,"Incorrect element names in EAM potential file");
error->one(FLERR, "Incorrect element names in EAM potential file");
file->elements = new char*[file->nelements];
file->elements = new char *[file->nelements];
for (int i = 0; i < file->nelements; i++)
file->elements[i] = utils::strdup(values.next_string());
values = reader.next_values(5);
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();
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");
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");
@ -158,13 +152,12 @@ void PairEAMFSIntel::read_file(char *filename)
for (int i = 0; i < file->nelements; i++) {
values = reader.next_values(2);
values.next_int(); // ignore
values.next_int(); // ignore
file->mass[i] = values.next_double();
reader.next_dvector(&file->frho[i][1], file->nrho);
if (unit_convert) {
for (int j = 1; j <= file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
for (int j = 1; j <= file->nrho; ++j) file->frho[i][j] *= conversion_factor;
}
for (int j = 0; j < file->nelements; j++) {
@ -176,8 +169,7 @@ void PairEAMFSIntel::read_file(char *filename)
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k <= file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
for (int k = 1; k <= file->nr; ++k) file->z2r[i][j][k] *= conversion_factor;
}
}
}
@ -197,7 +189,7 @@ void PairEAMFSIntel::read_file(char *filename)
// allocate memory on other procs
if (comm->me != 0) {
file->elements = new char*[file->nelements];
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");
@ -226,9 +218,7 @@ void PairEAMFSIntel::read_file(char *filename)
// 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);
}
for (int j = 0; j <= i; j++) { MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world); }
}
}
@ -238,7 +228,7 @@ void PairEAMFSIntel::read_file(char *filename)
void PairEAMFSIntel::file2array()
{
int i,j,m,n;
int i, j, m, n;
int ntypes = atom->ntypes;
// set function params directly from fs file
@ -247,7 +237,7 @@ void PairEAMFSIntel::file2array()
nr = fs->nr;
drho = fs->drho;
dr = fs->dr;
rhomax = (nrho-1) * drho;
rhomax = (nrho - 1) * drho;
// ------------------------------------------------------------------
// setup frho arrays
@ -258,7 +248,7 @@ void PairEAMFSIntel::file2array()
nfrho = fs->nelements + 1;
memory->destroy(frho);
memory->create(frho,nfrho,nrho+1,"pair:frho");
memory->create(frho, nfrho, nrho + 1, "pair:frho");
// copy each element's frho to global frho
@ -268,15 +258,17 @@ void PairEAMFSIntel::file2array()
// 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;
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;
if (map[i] >= 0)
type2frho[i] = map[i];
else
type2frho[i] = nfrho - 1;
// ------------------------------------------------------------------
// setup rhor arrays
@ -287,7 +279,7 @@ void PairEAMFSIntel::file2array()
nrhor = fs->nelements * fs->nelements;
memory->destroy(rhor);
memory->create(rhor,nrhor,nr+1,"pair:rhor");
memory->create(rhor, nrhor, nr + 1, "pair:rhor");
// copy each element pair rhor to global rhor
@ -303,8 +295,7 @@ void PairEAMFSIntel::file2array()
// 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];
for (j = 1; j <= ntypes; j++) type2rhor[i][j] = map[i] * fs->nelements + map[j];
// ------------------------------------------------------------------
// setup z2r arrays
@ -313,9 +304,9 @@ void PairEAMFSIntel::file2array()
// allocate z2r arrays
// nz2r = N*(N+1)/2 where N = # of fs elements
nz2r = fs->nelements * (fs->nelements+1) / 2;
nz2r = fs->nelements * (fs->nelements + 1) / 2;
memory->destroy(z2r);
memory->create(z2r,nz2r,nr+1,"pair:z2r");
memory->create(z2r, nz2r, nr + 1, "pair:z2r");
// copy each element pair z2r to global z2r, only for I >= J
@ -334,7 +325,7 @@ void PairEAMFSIntel::file2array()
// type2z2r is not used by non-opt
// but set type2z2r to 0 since accessed by opt
int irow,icol;
int irow, icol;
for (i = 1; i <= ntypes; i++) {
for (j = 1; j <= ntypes; j++) {
irow = map[i];

13
src/KOKKOS/pair_eam_alloy_kokkos.cpp
View File

@ -919,19 +919,14 @@ void PairEAMAlloyKokkos<DeviceType>::coeff(int narg, char **arg)
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// ensure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM setfl file
if (setfl) {
for (i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i];
delete [] setfl->elements;
delete [] setfl->mass;
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);

13
src/KOKKOS/pair_eam_fs_kokkos.cpp
View File

@ -920,19 +920,14 @@ void PairEAMFSKokkos<DeviceType>::coeff(int narg, char **arg)
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// ensure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM Finnis-Sinclair file
if (fs) {
for (i = 0; i < fs->nelements; i++) delete [] fs->elements[i];
delete [] fs->elements;
delete [] fs->mass;
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);

84
src/MANYBODY/pair_eam_alloy.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -43,18 +42,18 @@ PairEAMAlloy::PairEAMAlloy(LAMMPS *lmp) : PairEAM(lmp)
void PairEAMAlloy::coeff(int narg, char **arg)
{
int i,j;
int i, j;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM setfl file
if (setfl) {
for (i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i];
delete [] setfl->elements;
for (i = 0; i < setfl->nelements; i++) delete[] setfl->elements[i];
delete[] setfl->elements;
memory->destroy(setfl->mass);
memory->destroy(setfl->frho);
memory->destroy(setfl->rhor);
@ -68,22 +67,23 @@ void PairEAMAlloy::coeff(int narg, char **arg)
// 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;
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");
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;
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
@ -93,14 +93,14 @@ void PairEAMAlloy::coeff(int narg, char **arg)
for (j = i; j <= n; j++) {
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(FLERR,i,setfl->mass[map[i]]);
if (i == j) atom->set_mass(FLERR, i, setfl->mass[map[i]]);
count++;
}
scale[i][j] = 1.0;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -118,8 +118,7 @@ void PairEAMAlloy::read_file(char *filename)
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
double conversion_factor = utils::get_conversion_factor(utils::ENERGY, unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -129,22 +128,22 @@ void PairEAMAlloy::read_file(char *filename)
ValueTokenizer values = reader.next_values(1);
file->nelements = values.next_int();
if ((int)values.count() != file->nelements + 1)
error->one(FLERR,"Incorrect element names in EAM potential file");
if ((int) values.count() != file->nelements + 1)
error->one(FLERR, "Incorrect element names in EAM potential file");
file->elements = new char*[file->nelements];
file->elements = new char *[file->nelements];
for (int i = 0; i < file->nelements; i++)
file->elements[i] = utils::strdup(values.next_string());
values = reader.next_values(5);
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();
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");
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");
@ -153,14 +152,13 @@ void PairEAMAlloy::read_file(char *filename)
for (int i = 0; i < file->nelements; i++) {
values = reader.next_values(2);
values.next_int(); // ignore
values.next_int(); // ignore
file->mass[i] = values.next_double();
reader.next_dvector(&file->frho[i][1], file->nrho);
reader.next_dvector(&file->rhor[i][1], file->nr);
if (unit_convert) {
for (int j = 1; j < file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
for (int j = 1; j < file->nrho; ++j) file->frho[i][j] *= conversion_factor;
}
}
@ -168,8 +166,7 @@ void PairEAMAlloy::read_file(char *filename)
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k < file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
for (int k = 1; k < file->nr; ++k) file->z2r[i][j][k] *= conversion_factor;
}
}
}
@ -189,7 +186,7 @@ void PairEAMAlloy::read_file(char *filename)
// allocate memory on other procs
if (comm->me != 0) {
file->elements = new char*[file->nelements];
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");
@ -215,9 +212,7 @@ void PairEAMAlloy::read_file(char *filename)
// 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);
}
for (int j = 0; j <= i; j++) { MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world); }
}
}
@ -227,7 +222,7 @@ void PairEAMAlloy::read_file(char *filename)
void PairEAMAlloy::file2array()
{
int i,j,m,n;
int i, j, m, n;
int ntypes = atom->ntypes;
// set function params directly from setfl file
@ -236,7 +231,7 @@ void PairEAMAlloy::file2array()
nr = setfl->nr;
drho = setfl->drho;
dr = setfl->dr;
rhomax = (nrho-1) * drho;
rhomax = (nrho - 1) * drho;
// ------------------------------------------------------------------
// setup frho arrays
@ -247,7 +242,7 @@ void PairEAMAlloy::file2array()
nfrho = setfl->nelements + 1;
memory->destroy(frho);
memory->create(frho,nfrho,nrho+1,"pair:frho");
memory->create(frho, nfrho, nrho + 1, "pair:frho");
// copy each element's frho to global frho
@ -257,15 +252,17 @@ void PairEAMAlloy::file2array()
// 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;
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;
if (map[i] >= 0)
type2frho[i] = map[i];
else
type2frho[i] = nfrho - 1;
// ------------------------------------------------------------------
// setup rhor arrays
@ -276,7 +273,7 @@ void PairEAMAlloy::file2array()
nrhor = setfl->nelements;
memory->destroy(rhor);
memory->create(rhor,nrhor,nr+1,"pair:rhor");
memory->create(rhor, nrhor, nr + 1, "pair:rhor");
// copy each element's rhor to global rhor
@ -288,8 +285,7 @@ void PairEAMAlloy::file2array()
// 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];
for (j = 1; j <= ntypes; j++) type2rhor[i][j] = map[i];
// ------------------------------------------------------------------
// setup z2r arrays
@ -298,9 +294,9 @@ void PairEAMAlloy::file2array()
// allocate z2r arrays
// nz2r = N*(N+1)/2 where N = # of setfl elements
nz2r = setfl->nelements * (setfl->nelements+1) / 2;
nz2r = setfl->nelements * (setfl->nelements + 1) / 2;
memory->destroy(z2r);
memory->create(z2r,nz2r,nr+1,"pair:z2r");
memory->create(z2r, nz2r, nr + 1, "pair:z2r");
// copy each element pair z2r to global z2r, only for I >= J
@ -319,7 +315,7 @@ void PairEAMAlloy::file2array()
// type2z2r is not used by non-opt
// but set type2z2r to 0 since accessed by opt
int irow,icol;
int irow, icol;
for (i = 1; i <= ntypes; i++) {
for (j = 1; j <= ntypes; j++) {
irow = map[i];

97
src/MANYBODY/pair_eam_fs.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -44,18 +43,18 @@ PairEAMFS::PairEAMFS(LAMMPS *lmp) : PairEAM(lmp)
void PairEAMFS::coeff(int narg, char **arg)
{
int i,j;
int i, j;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM Finnis-Sinclair file
if (fs) {
for (i = 0; i < fs->nelements; i++) delete [] fs->elements[i];
delete [] fs->elements;
for (i = 0; i < fs->nelements; i++) delete[] fs->elements[i];
delete[] fs->elements;
memory->destroy(fs->mass);
memory->destroy(fs->frho);
memory->destroy(fs->rhor);
@ -69,22 +68,23 @@ void PairEAMFS::coeff(int narg, char **arg)
// 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;
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");
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;
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
@ -94,14 +94,14 @@ void PairEAMFS::coeff(int narg, char **arg)
for (j = i; j <= n; j++) {
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(FLERR,i,fs->mass[map[i]]);
if (i == j) atom->set_mass(FLERR, i, fs->mass[map[i]]);
count++;
}
scale[i][j] = 1.0;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -114,14 +114,12 @@ void PairEAMFS::read_file(char *filename)
// read potential file
if (comm->me == 0) {
PotentialFileReader reader(lmp, filename, he_flag ? "eam/he" : "eam/fs",
unit_convert_flag);
PotentialFileReader reader(lmp, filename, he_flag ? "eam/he" : "eam/fs", unit_convert_flag);
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
double conversion_factor = utils::get_conversion_factor(utils::ENERGY, unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -131,24 +129,26 @@ void PairEAMFS::read_file(char *filename)
ValueTokenizer values = reader.next_values(1);
file->nelements = values.next_int();
if ((int)values.count() != file->nelements + 1)
error->one(FLERR,"Incorrect element names in EAM potential file");
if ((int) values.count() != file->nelements + 1)
error->one(FLERR, "Incorrect element names in EAM potential file");
file->elements = new char*[file->nelements];
file->elements = new char *[file->nelements];
for (int i = 0; i < file->nelements; i++)
file->elements[i] = utils::strdup(values.next_string());
if (he_flag) values = reader.next_values(6);
else values = reader.next_values(5);
if (he_flag)
values = reader.next_values(6);
else
values = reader.next_values(5);
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();
file->nr = values.next_int();
file->dr = values.next_double();
file->cut = values.next_double();
if (he_flag) rhomax = values.next_double();
if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
error->one(FLERR,"Invalid EAM potential file");
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");
@ -157,13 +157,12 @@ void PairEAMFS::read_file(char *filename)
for (int i = 0; i < file->nelements; i++) {
values = reader.next_values(2);
values.next_int(); // ignore
values.next_int(); // ignore
file->mass[i] = values.next_double();
reader.next_dvector(&file->frho[i][1], file->nrho);
if (unit_convert) {
for (int j = 1; j <= file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
for (int j = 1; j <= file->nrho; ++j) file->frho[i][j] *= conversion_factor;
}
for (int j = 0; j < file->nelements; j++) {
@ -175,8 +174,7 @@ void PairEAMFS::read_file(char *filename)
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k <= file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
for (int k = 1; k <= file->nr; ++k) file->z2r[i][j][k] *= conversion_factor;
}
}
}
@ -197,7 +195,7 @@ void PairEAMFS::read_file(char *filename)
// allocate memory on other procs
if (comm->me != 0) {
file->elements = new char*[file->nelements];
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");
@ -226,9 +224,7 @@ void PairEAMFS::read_file(char *filename)
// 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);
}
for (int j = 0; j <= i; j++) { MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world); }
}
}
@ -238,7 +234,7 @@ void PairEAMFS::read_file(char *filename)
void PairEAMFS::file2array()
{
int i,j,m,n;
int i, j, m, n;
int ntypes = atom->ntypes;
// set function params directly from fs file
@ -247,8 +243,10 @@ void PairEAMFS::file2array()
nr = fs->nr;
drho = fs->drho;
dr = fs->dr;
if (he_flag) rhomin = rhomax - (nrho-1) * drho;
else rhomax = (nrho-1) * drho;
if (he_flag)
rhomin = rhomax - (nrho - 1) * drho;
else
rhomax = (nrho - 1) * drho;
// ------------------------------------------------------------------
// setup frho arrays
@ -259,7 +257,7 @@ void PairEAMFS::file2array()
nfrho = fs->nelements + 1;
memory->destroy(frho);
memory->create(frho,nfrho,nrho+1,"pair:frho");
memory->create(frho, nfrho, nrho + 1, "pair:frho");
// copy each element's frho to global frho
@ -269,15 +267,17 @@ void PairEAMFS::file2array()
// 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;
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;
if (map[i] >= 0)
type2frho[i] = map[i];
else
type2frho[i] = nfrho - 1;
// ------------------------------------------------------------------
// setup rhor arrays
@ -288,7 +288,7 @@ void PairEAMFS::file2array()
nrhor = fs->nelements * fs->nelements;
memory->destroy(rhor);
memory->create(rhor,nrhor,nr+1,"pair:rhor");
memory->create(rhor, nrhor, nr + 1, "pair:rhor");
// copy each element pair rhor to global rhor
@ -304,8 +304,7 @@ void PairEAMFS::file2array()
// 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];
for (j = 1; j <= ntypes; j++) type2rhor[i][j] = map[i] * fs->nelements + map[j];
// ------------------------------------------------------------------
// setup z2r arrays
@ -314,9 +313,9 @@ void PairEAMFS::file2array()
// allocate z2r arrays
// nz2r = N*(N+1)/2 where N = # of fs elements
nz2r = fs->nelements * (fs->nelements+1) / 2;
nz2r = fs->nelements * (fs->nelements + 1) / 2;
memory->destroy(z2r);
memory->create(z2r,nz2r,nr+1,"pair:z2r");
memory->create(z2r, nz2r, nr + 1, "pair:z2r");
// copy each element pair z2r to global z2r, only for I >= J
@ -335,7 +334,7 @@ void PairEAMFS::file2array()
// type2z2r is not used by non-opt
// but set type2z2r to 0 since accessed by opt
int irow,icol;
int irow, icol;
for (i = 1; i <= ntypes; i++) {
for (j = 1; j <= ntypes; j++) {
irow = map[i];

96
src/MANYBODY/pair_eam_he.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -38,14 +37,14 @@ PairEAMHE::PairEAMHE(LAMMPS *lmp) : PairEAM(lmp), PairEAMFS(lmp)
void PairEAMHE::compute(int eflag, int vflag)
{
int i,j,ii,jj,m,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r,p,rhoip,rhojp,z2,z2p,recip,phip,psip,phi;
int i, j, ii, jj, m, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double rsq, r, p, rhoip, rhojp, z2, z2p, recip, phip, psip, phi;
double *coeff;
int *ilist,*jlist,*numneigh,**firstneigh;
int *ilist, *jlist, *numneigh, **firstneigh;
evdwl = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
// grow energy and fp arrays if necessary
// need to be atom->nmax in length
@ -55,9 +54,9 @@ void PairEAMHE::compute(int eflag, int vflag)
memory->destroy(fp);
memory->destroy(numforce);
nmax = atom->nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(fp,nmax,"pair:fp");
memory->create(numforce,nmax,"pair:numforce");
memory->create(rho, nmax, "pair:rho");
memory->create(fp, nmax, "pair:fp");
memory->create(numforce, nmax, "pair:numforce");
}
double **x = atom->x;
@ -76,7 +75,8 @@ void PairEAMHE::compute(int eflag, int vflag)
if (newton_pair) {
for (i = 0; i < nall; i++) rho[i] = 0.0;
} else for (i = 0; i < nlocal; i++) rho[i] = 0.0;
} else
for (i = 0; i < nlocal; i++) rho[i] = 0.0;
// rho = density at each atom
// loop over neighbors of my atoms
@ -97,20 +97,20 @@ void PairEAMHE::compute(int eflag, int vflag)
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cutforcesq) {
jtype = type[j];
p = sqrt(rsq)*rdr + 1.0;
m = static_cast<int> (p);
m = MIN(m,nr-1);
p = sqrt(rsq) * rdr + 1.0;
m = static_cast<int>(p);
m = MIN(m, nr - 1);
p -= m;
p = MIN(p,1.0);
p = MIN(p, 1.0);
coeff = rhor_spline[type2rhor[jtype][itype]][m];
rho[i] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
rho[i] += ((coeff[3] * p + coeff[4]) * p + coeff[5]) * p + coeff[6];
if (newton_pair || j < nlocal) {
coeff = rhor_spline[type2rhor[itype][jtype]][m];
rho[j] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
rho[j] += ((coeff[3] * p + coeff[4]) * p + coeff[5]) * p + coeff[6];
}
}
}
@ -127,12 +127,12 @@ void PairEAMHE::compute(int eflag, int vflag)
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
p = (rho[i]-rhomin)*rdrho + 1.0;
m = static_cast<int> (p);
p = (rho[i] - rhomin) * rdrho + 1.0;
m = static_cast<int>(p);
if (m < 2) {
m = 2;
} else if (m > nrho-1) {
m = nrho-1;
} else if (m > nrho - 1) {
m = nrho - 1;
}
p -= m;
if (p < -1.0) {
@ -141,13 +141,13 @@ void PairEAMHE::compute(int eflag, int vflag)
p = 1.0;
}
coeff = frho_spline[type2frho[type[i]]][m];
fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2];
fp[i] = (coeff[0] * p + coeff[1]) * p + coeff[2];
if (eflag) {
phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
phi = ((coeff[3] * p + coeff[4]) * p + coeff[5]) * p + coeff[6];
if (rho[i] < rhomin) {
phi += fp[i] * (rho[i]-rhomin);
phi += fp[i] * (rho[i] - rhomin);
} else if (rho[i] > rhomax) {
phi += fp[i] * (rho[i]-rhomax);
phi += fp[i] * (rho[i] - rhomax);
}
phi *= scale[type[i]][type[i]];
if (eflag_global) eng_vdwl += phi;
@ -181,17 +181,17 @@ void PairEAMHE::compute(int eflag, int vflag)
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cutforcesq) {
++numforce[i];
jtype = type[j];
r = sqrt(rsq);
p = r*rdr + 1.0;
m = static_cast<int> (p);
m = MIN(m,nr-1);
p = r * rdr + 1.0;
m = static_cast<int>(p);
m = MIN(m, nr - 1);
p -= m;
p = MIN(p,1.0);
p = MIN(p, 1.0);
// rhoip = derivative of (density at atom j due to atom i)
// rhojp = derivative of (density at atom i due to atom j)
@ -205,35 +205,33 @@ void PairEAMHE::compute(int eflag, int vflag)
// scale factor can be applied by thermodynamic integration
coeff = rhor_spline[type2rhor[itype][jtype]][m];
rhoip = (coeff[0]*p + coeff[1])*p + coeff[2];
rhoip = (coeff[0] * p + coeff[1]) * p + coeff[2];
coeff = rhor_spline[type2rhor[jtype][itype]][m];
rhojp = (coeff[0]*p + coeff[1])*p + coeff[2];
rhojp = (coeff[0] * p + coeff[1]) * p + coeff[2];
coeff = z2r_spline[type2z2r[itype][jtype]][m];
z2p = (coeff[0]*p + coeff[1])*p + coeff[2];
z2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
z2p = (coeff[0] * p + coeff[1]) * p + coeff[2];
z2 = ((coeff[3] * p + coeff[4]) * p + coeff[5]) * p + coeff[6];
recip = 1.0/r;
phi = z2*recip;
phip = z2p*recip - phi*recip;
psip = fp[i]*rhojp + fp[j]*rhoip + phip;
fpair = -scale[itype][jtype]*psip*recip;
recip = 1.0 / r;
phi = z2 * recip;
phip = z2p * recip - phi * recip;
psip = fp[i] * rhojp + fp[j] * rhoip + phip;
fpair = -scale[itype][jtype] * psip * recip;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (eflag) evdwl = scale[itype][jtype]*phi;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
if (eflag) evdwl = scale[itype][jtype] * phi;
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}

3
src/OPENMP/pair_eam_alloy_omp.cpp
View File

@ -45,7 +45,8 @@ void PairEAMAlloyOMP::coeff(int narg, char **arg)
if (!allocated) allocate();
if (narg != 3 + atom->ntypes) error->all(FLERR, "Incorrect args for pair coefficients");
if (narg != 3 + atom->ntypes)
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM setfl file

3
src/OPENMP/pair_eam_fs_omp.cpp
View File

@ -45,7 +45,8 @@ void PairEAMFSOMP::coeff(int narg, char **arg)
if (!allocated) allocate();
if (narg != 3 + atom->ntypes) error->all(FLERR, "Incorrect args for pair coefficients");
if (narg != 3 + atom->ntypes)
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
// read EAM Finnis-Sinclair file

4
src/OPT/pair_eam_alloy_opt.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -31,5 +30,4 @@ using namespace LAMMPS_NS;
inherit everything else from PairEAMAlloy
------------------------------------------------------------------------- */
PairEAMAlloyOpt::PairEAMAlloyOpt(LAMMPS *lmp) :
PairEAM(lmp), PairEAMAlloy(lmp), PairEAMOpt(lmp) {}
PairEAMAlloyOpt::PairEAMAlloyOpt(LAMMPS *lmp) : PairEAM(lmp), PairEAMAlloy(lmp), PairEAMOpt(lmp) {}

4
src/OPT/pair_eam_fs_opt.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -31,5 +30,4 @@ using namespace LAMMPS_NS;
inherit everything else from PairEAMFS
------------------------------------------------------------------------- */
PairEAMFSOpt::PairEAMFSOpt(LAMMPS *lmp) :
PairEAM(lmp), PairEAMFS(lmp), PairEAMOpt(lmp) {}
PairEAMFSOpt::PairEAMFSOpt(LAMMPS *lmp) : PairEAM(lmp), PairEAMFS(lmp), PairEAMOpt(lmp) {}

287
src/OPT/pair_eam_opt.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -21,15 +20,16 @@
------------------------------------------------------------------------- */
#include "pair_eam_opt.h"
#include <cmath>
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neigh_list.h"
#include "memory.h"
#include "neigh_list.h"
#include "update.h"
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
@ -40,44 +40,51 @@ PairEAMOpt::PairEAMOpt(LAMMPS *lmp) : PairEAM(lmp) {}
void PairEAMOpt::compute(int eflag, int vflag)
{
ev_init(eflag,vflag);
ev_init(eflag, vflag);
if (evflag) {
if (eflag) {
if (force->newton_pair) return eval<1,1,1>();
else return eval<1,1,0>();
if (force->newton_pair)
return eval<1, 1, 1>();
else
return eval<1, 1, 0>();
} else {
if (force->newton_pair) return eval<1,0,1>();
else return eval<1,0,0>();
if (force->newton_pair)
return eval<1, 0, 1>();
else
return eval<1, 0, 0>();
}
} else {
if (force->newton_pair) return eval<0,0,1>();
else return eval<0,0,0>();
if (force->newton_pair)
return eval<0, 0, 1>();
else
return eval<0, 0, 0>();
}
}
/* ---------------------------------------------------------------------- */
template < int EVFLAG, int EFLAG, int NEWTON_PAIR >
void PairEAMOpt::eval()
template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairEAMOpt::eval()
{
typedef struct { double x,y,z; } vec3_t;
typedef struct {
double x, y, z;
} vec3_t;
typedef struct {
double rhor0i,rhor1i,rhor2i,rhor3i;
double rhor0j,rhor1j,rhor2j,rhor3j;
double rhor0i, rhor1i, rhor2i, rhor3i;
double rhor0j, rhor1j, rhor2j, rhor3j;
} fast_alpha_t;
typedef struct {
double rhor4i,rhor5i,rhor6i;
double rhor4j,rhor5j,rhor6j;
double z2r0,z2r1,z2r2,z2r3,z2r4,z2r5,z2r6;
double rhor4i, rhor5i, rhor6i;
double rhor4j, rhor5j, rhor6j;
double z2r0, z2r1, z2r2, z2r3, z2r4, z2r5, z2r6;
double _pad[3];
} fast_gamma_t;
int i,j,ii,jj,inum,jnum,itype,jtype;
int i, j, ii, jj, inum, jnum, itype, jtype;
double evdwl = 0.0;
double* _noalias coeff;
double *_noalias coeff;
// grow energy array if necessary
@ -86,94 +93,97 @@ void PairEAMOpt::eval()
memory->destroy(fp);
memory->destroy(numforce);
nmax = atom->nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(fp,nmax,"pair:fp");
memory->create(numforce,nmax,"pair:numforce");
memory->create(rho, nmax, "pair:rho");
memory->create(fp, nmax, "pair:fp");
memory->create(numforce, nmax, "pair:numforce");
}
double** _noalias x = atom->x;
double** _noalias f = atom->f;
int* _noalias type = atom->type;
double **_noalias x = atom->x;
double **_noalias f = atom->f;
int *_noalias type = atom->type;
int nlocal = atom->nlocal;
auto * _noalias xx = (vec3_t*)x[0];
auto * _noalias ff = (vec3_t*)f[0];
auto *_noalias xx = (vec3_t *) x[0];
auto *_noalias ff = (vec3_t *) f[0];
double tmp_cutforcesq = cutforcesq;
double tmp_rdr = rdr;
int nr2 = nr-2;
int nr1 = nr-1;
int nr2 = nr - 2;
int nr1 = nr - 1;
inum = list->inum;
int* _noalias ilist = list->ilist;
int** _noalias firstneigh = list->firstneigh;
int* _noalias numneigh = list->numneigh;
int *_noalias ilist = list->ilist;
int **_noalias firstneigh = list->firstneigh;
int *_noalias numneigh = list->numneigh;
int ntypes = atom->ntypes;
int ntypes2 = ntypes*ntypes;
int ntypes2 = ntypes * ntypes;
auto * _noalias fast_alpha =
(fast_alpha_t*) malloc((size_t)ntypes2*(nr+1)*sizeof(fast_alpha_t));
for (i = 0; i < ntypes; i++) for (j = 0; j < ntypes; j++) {
auto * _noalias tab = &fast_alpha[i*ntypes*nr+j*nr];
if (type2rhor[i+1][j+1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].rhor0i = rhor_spline[type2rhor[i+1][j+1]][m][6];
tab[m].rhor1i = rhor_spline[type2rhor[i+1][j+1]][m][5];
tab[m].rhor2i = rhor_spline[type2rhor[i+1][j+1]][m][4];
tab[m].rhor3i = rhor_spline[type2rhor[i+1][j+1]][m][3];
auto *_noalias fast_alpha =
(fast_alpha_t *) malloc((size_t) ntypes2 * (nr + 1) * sizeof(fast_alpha_t));
for (i = 0; i < ntypes; i++)
for (j = 0; j < ntypes; j++) {
auto *_noalias tab = &fast_alpha[i * ntypes * nr + j * nr];
if (type2rhor[i + 1][j + 1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].rhor0i = rhor_spline[type2rhor[i + 1][j + 1]][m][6];
tab[m].rhor1i = rhor_spline[type2rhor[i + 1][j + 1]][m][5];
tab[m].rhor2i = rhor_spline[type2rhor[i + 1][j + 1]][m][4];
tab[m].rhor3i = rhor_spline[type2rhor[i + 1][j + 1]][m][3];
}
}
if (type2rhor[j + 1][i + 1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].rhor0j = rhor_spline[type2rhor[j + 1][i + 1]][m][6];
tab[m].rhor1j = rhor_spline[type2rhor[j + 1][i + 1]][m][5];
tab[m].rhor2j = rhor_spline[type2rhor[j + 1][i + 1]][m][4];
tab[m].rhor3j = rhor_spline[type2rhor[j + 1][i + 1]][m][3];
}
}
}
if (type2rhor[j+1][i+1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].rhor0j = rhor_spline[type2rhor[j+1][i+1]][m][6];
tab[m].rhor1j = rhor_spline[type2rhor[j+1][i+1]][m][5];
tab[m].rhor2j = rhor_spline[type2rhor[j+1][i+1]][m][4];
tab[m].rhor3j = rhor_spline[type2rhor[j+1][i+1]][m][3];
}
}
}
auto * _noalias tabeight = fast_alpha;
auto *_noalias tabeight = fast_alpha;
auto * _noalias fast_gamma =
(fast_gamma_t*) malloc((size_t)ntypes2*(nr+1)*sizeof(fast_gamma_t));
for (i = 0; i < ntypes; i++) for (j = 0; j < ntypes; j++) {
auto * _noalias tab = &fast_gamma[i*ntypes*nr+j*nr];
if (type2rhor[i+1][j+1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].rhor4i = rhor_spline[type2rhor[i+1][j+1]][m][2];
tab[m].rhor5i = rhor_spline[type2rhor[i+1][j+1]][m][1];
tab[m].rhor6i = rhor_spline[type2rhor[i+1][j+1]][m][0];
auto *_noalias fast_gamma =
(fast_gamma_t *) malloc((size_t) ntypes2 * (nr + 1) * sizeof(fast_gamma_t));
for (i = 0; i < ntypes; i++)
for (j = 0; j < ntypes; j++) {
auto *_noalias tab = &fast_gamma[i * ntypes * nr + j * nr];
if (type2rhor[i + 1][j + 1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].rhor4i = rhor_spline[type2rhor[i + 1][j + 1]][m][2];
tab[m].rhor5i = rhor_spline[type2rhor[i + 1][j + 1]][m][1];
tab[m].rhor6i = rhor_spline[type2rhor[i + 1][j + 1]][m][0];
}
}
if (type2rhor[j + 1][i + 1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].rhor4j = rhor_spline[type2rhor[j + 1][i + 1]][m][2];
tab[m].rhor5j = rhor_spline[type2rhor[j + 1][i + 1]][m][1];
tab[m].rhor6j = rhor_spline[type2rhor[j + 1][i + 1]][m][0];
tab[m].z2r6 = z2r_spline[type2z2r[i + 1][j + 1]][m][0];
}
}
if (type2z2r[i + 1][j + 1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].z2r0 = z2r_spline[type2z2r[i + 1][j + 1]][m][6];
tab[m].z2r1 = z2r_spline[type2z2r[i + 1][j + 1]][m][5];
tab[m].z2r2 = z2r_spline[type2z2r[i + 1][j + 1]][m][4];
tab[m].z2r3 = z2r_spline[type2z2r[i + 1][j + 1]][m][3];
tab[m].z2r4 = z2r_spline[type2z2r[i + 1][j + 1]][m][2];
tab[m].z2r5 = z2r_spline[type2z2r[i + 1][j + 1]][m][1];
tab[m].z2r6 = z2r_spline[type2z2r[i + 1][j + 1]][m][0];
}
}
}
if (type2rhor[j+1][i+1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].rhor4j = rhor_spline[type2rhor[j+1][i+1]][m][2];
tab[m].rhor5j = rhor_spline[type2rhor[j+1][i+1]][m][1];
tab[m].rhor6j = rhor_spline[type2rhor[j+1][i+1]][m][0];
tab[m].z2r6 = z2r_spline[type2z2r[i+1][j+1]][m][0];
}
}
if (type2z2r[i+1][j+1] >= 0) {
for (int m = 1; m <= nr; m++) {
tab[m].z2r0 = z2r_spline[type2z2r[i+1][j+1]][m][6];
tab[m].z2r1 = z2r_spline[type2z2r[i+1][j+1]][m][5];
tab[m].z2r2 = z2r_spline[type2z2r[i+1][j+1]][m][4];
tab[m].z2r3 = z2r_spline[type2z2r[i+1][j+1]][m][3];
tab[m].z2r4 = z2r_spline[type2z2r[i+1][j+1]][m][2];
tab[m].z2r5 = z2r_spline[type2z2r[i+1][j+1]][m][1];
tab[m].z2r6 = z2r_spline[type2z2r[i+1][j+1]][m][0];
}
}
}
auto * _noalias tabss = fast_gamma;
auto *_noalias tabss = fast_gamma;
// zero out density
if (NEWTON_PAIR) {
int m = nlocal + atom->nghost;
for (i = 0; i < m; i++) rho[i] = 0.0;
} else for (i = 0; i < nlocal; i++) rho[i] = 0.0;
} else
for (i = 0; i < nlocal; i++) rho[i] = 0.0;
// rho = density at each atom
// loop over neighbors of my atoms
@ -184,11 +194,11 @@ void PairEAMOpt::eval()
double ytmp = xx[i].y;
double ztmp = xx[i].z;
itype = type[i] - 1;
int* _noalias jlist = firstneigh[i];
int *_noalias jlist = firstneigh[i];
jnum = numneigh[i];
double tmprho = rho[i];
auto * _noalias tabeighti = &tabeight[itype*ntypes*nr];
auto *_noalias tabeighti = &tabeight[itype * ntypes * nr];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
@ -197,26 +207,24 @@ void PairEAMOpt::eval()
double delx = xtmp - xx[j].x;
double dely = ytmp - xx[j].y;
double delz = ztmp - xx[j].z;
double rsq = delx*delx + dely*dely + delz*delz;
double rsq = delx * delx + dely * dely + delz * delz;
if (rsq < tmp_cutforcesq) {
jtype = type[j] - 1;
double p = sqrt(rsq)*tmp_rdr;
if ((int)p <= nr2) {
int m = (int)p + 1;
p -= (double)((int)p);
fast_alpha_t& a = tabeighti[jtype*nr+m];
tmprho += ((a.rhor3j*p+a.rhor2j)*p+a.rhor1j)*p+a.rhor0j;
double p = sqrt(rsq) * tmp_rdr;
if ((int) p <= nr2) {
int m = (int) p + 1;
p -= (double) ((int) p);
fast_alpha_t &a = tabeighti[jtype * nr + m];
tmprho += ((a.rhor3j * p + a.rhor2j) * p + a.rhor1j) * p + a.rhor0j;
if (NEWTON_PAIR || j < nlocal) {
rho[j] += ((a.rhor3i*p+a.rhor2i)*p+a.rhor1i)*p+a.rhor0i;
rho[j] += ((a.rhor3i * p + a.rhor2i) * p + a.rhor1i) * p + a.rhor0i;
}
} else {
fast_alpha_t& a = tabeighti[jtype*nr+nr1];
tmprho += a.rhor3j+a.rhor2j+a.rhor1j+a.rhor0j;
if (NEWTON_PAIR || j < nlocal) {
rho[j] += a.rhor3i+a.rhor2i+a.rhor1i+a.rhor0i;
}
fast_alpha_t &a = tabeighti[jtype * nr + nr1];
tmprho += a.rhor3j + a.rhor2j + a.rhor1j + a.rhor0j;
if (NEWTON_PAIR || j < nlocal) { rho[j] += a.rhor3i + a.rhor2i + a.rhor1i + a.rhor0i; }
}
}
}
@ -234,16 +242,16 @@ void PairEAMOpt::eval()
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
double p = rho[i]*rdrho + 1.0;
int m = static_cast<int> (p);
m = MAX(1,MIN(m,nrho-1));
double p = rho[i] * rdrho + 1.0;
int m = static_cast<int>(p);
m = MAX(1, MIN(m, nrho - 1));
p -= m;
p = MIN(p,1.0);
p = MIN(p, 1.0);
coeff = frho_spline[type2frho[type[i]]][m];
fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2];
fp[i] = (coeff[0] * p + coeff[1]) * p + coeff[2];
if (EFLAG) {
double phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
double phi = ((coeff[3] * p + coeff[4]) * p + coeff[5]) * p + coeff[6];
if (rho[i] > rhomax) phi += fp[i] * (rho[i] - rhomax);
phi *= scale[type[i]][type[i]];
if (eflag_global) eng_vdwl += phi;
if (eflag_atom) eatom[i] += phi;
@ -264,15 +272,15 @@ void PairEAMOpt::eval()
double ytmp = xx[i].y;
double ztmp = xx[i].z;
int itype1 = type[i] - 1;
int* _noalias jlist = firstneigh[i];
int *_noalias jlist = firstneigh[i];
jnum = numneigh[i];
double tmpfx = 0.0;
double tmpfy = 0.0;
double tmpfz = 0.0;
auto * _noalias tabssi = &tabss[itype1*ntypes*nr];
double* _noalias scale_i = scale[itype1+1]+1;
auto *_noalias tabssi = &tabss[itype1 * ntypes * nr];
double *_noalias scale_i = scale[itype1 + 1] + 1;
numforce[i] = 0;
for (jj = 0; jj < jnum; jj++) {
@ -282,29 +290,29 @@ void PairEAMOpt::eval()
double delx = xtmp - xx[j].x;
double dely = ytmp - xx[j].y;
double delz = ztmp - xx[j].z;
double rsq = delx*delx + dely*dely + delz*delz;
double rsq = delx * delx + dely * dely + delz * delz;
if (rsq < tmp_cutforcesq) {
++numforce[i];
jtype = type[j] - 1;
double r = sqrt(rsq);
double rhoip,rhojp,z2,z2p;
double p = r*tmp_rdr;
if ((int)p <= nr2) {
double rhoip, rhojp, z2, z2p;
double p = r * tmp_rdr;
if ((int) p <= nr2) {
int m = (int) p + 1;
m = MIN(m,nr-1);
p -= (double)((int) p);
p = MIN(p,1.0);
m = MIN(m, nr - 1);
p -= (double) ((int) p);
p = MIN(p, 1.0);
fast_gamma_t& a = tabssi[jtype*nr+m];
rhoip = (a.rhor6i*p + a.rhor5i)*p + a.rhor4i;
rhojp = (a.rhor6j*p + a.rhor5j)*p + a.rhor4j;
z2 = ((a.z2r3*p + a.z2r2)*p + a.z2r1)*p + a.z2r0;
z2p = (a.z2r6*p + a.z2r5)*p + a.z2r4;
fast_gamma_t &a = tabssi[jtype * nr + m];
rhoip = (a.rhor6i * p + a.rhor5i) * p + a.rhor4i;
rhojp = (a.rhor6j * p + a.rhor5j) * p + a.rhor4j;
z2 = ((a.z2r3 * p + a.z2r2) * p + a.z2r1) * p + a.z2r0;
z2p = (a.z2r6 * p + a.z2r5) * p + a.z2r4;
} else {
fast_gamma_t& a = tabssi[jtype*nr+nr1];
fast_gamma_t &a = tabssi[jtype * nr + nr1];
rhoip = a.rhor6i + a.rhor5i + a.rhor4i;
rhojp = a.rhor6j + a.rhor5j + a.rhor4j;
z2 = a.z2r3 + a.z2r2 + a.z2r1 + a.z2r0;
@ -322,25 +330,24 @@ void PairEAMOpt::eval()
// hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
// scale factor can be applied by thermodynamic integration
double recip = 1.0/r;
double phi = z2*recip;
double phip = z2p*recip - phi*recip;
double psip = fp[i]*rhojp + fp[j]*rhoip + phip;
double fpair = -scale_i[jtype]*psip*recip;
double recip = 1.0 / r;
double phi = z2 * recip;
double phip = z2p * recip - phi * recip;
double psip = fp[i] * rhojp + fp[j] * rhoip + phip;
double fpair = -scale_i[jtype] * psip * recip;
tmpfx += delx*fpair;
tmpfy += dely*fpair;
tmpfz += delz*fpair;
tmpfx += delx * fpair;
tmpfy += dely * fpair;
tmpfz += delz * fpair;
if (NEWTON_PAIR || j < nlocal) {
ff[j].x -= delx*fpair;
ff[j].y -= dely*fpair;
ff[j].z -= delz*fpair;
ff[j].x -= delx * fpair;
ff[j].y -= dely * fpair;
ff[j].z -= delz * fpair;
}
if (EFLAG) evdwl = scale_i[jtype]*phi;
if (EFLAG) evdwl = scale_i[jtype] * phi;
if (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR,
evdwl,0.0,fpair,delx,dely,delz);
if (EVFLAG) ev_tally(i, j, nlocal, NEWTON_PAIR, evdwl, 0.0, fpair, delx, dely, delz);
}
}
@ -349,8 +356,10 @@ void PairEAMOpt::eval()
ff[i].z += tmpfz;
}
free(fast_alpha); fast_alpha = nullptr;
free(fast_gamma); fast_gamma = nullptr;
free(fast_alpha);
fast_alpha = nullptr;
free(fast_gamma);
fast_gamma = nullptr;
if (vflag_fdotr) virial_fdotr_compute();
}

2
src/pair.cpp
View File

@ -821,7 +821,7 @@ void Pair::map_element2type(int narg, char **arg, bool update_setflag)
// elements = list of element names
if (narg != ntypes)
error->all(FLERR,"Number of element to type mappings does not match number of atom types");
error->all(FLERR, "Number of element to type mappings does not match number of atom types");
if (elements) {
for (i = 0; i < nelements; i++) delete[] elements[i];