whitespace cleanup
This commit is contained in:
@ -102,7 +102,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
const double* const * const x = atom->x;
|
const double* const * const x = atom->x;
|
||||||
double* const * const forces = atom->f;
|
double* const * const forces = atom->f;
|
||||||
const int ntypes = atom->ntypes;
|
const int ntypes = atom->ntypes;
|
||||||
|
|
||||||
if (eflag || vflag) {
|
if (eflag || vflag) {
|
||||||
ev_setup(eflag, vflag);
|
ev_setup(eflag, vflag);
|
||||||
} else {
|
} else {
|
||||||
@ -123,7 +123,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
int newMaxNeighbors = 0;
|
int newMaxNeighbors = 0;
|
||||||
for(int ii = 0; ii < listfull->inum; ii++) {
|
for(int ii = 0; ii < listfull->inum; ii++) {
|
||||||
int jnum = listfull->numneigh[listfull->ilist[ii]];
|
int jnum = listfull->numneigh[listfull->ilist[ii]];
|
||||||
if(jnum > newMaxNeighbors)
|
if(jnum > newMaxNeighbors)
|
||||||
newMaxNeighbors = jnum;
|
newMaxNeighbors = jnum;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -151,24 +151,24 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
for(int jj = 0; jj < listfull->numneigh[i]; jj++) {
|
for(int jj = 0; jj < listfull->numneigh[i]; jj++) {
|
||||||
int j = listfull->firstneigh[i][jj];
|
int j = listfull->firstneigh[i][jj];
|
||||||
j &= NEIGHMASK;
|
j &= NEIGHMASK;
|
||||||
|
|
||||||
double jdelx = x[j][0] - x[i][0];
|
double jdelx = x[j][0] - x[i][0];
|
||||||
double jdely = x[j][1] - x[i][1];
|
double jdely = x[j][1] - x[i][1];
|
||||||
double jdelz = x[j][2] - x[i][2];
|
double jdelz = x[j][2] - x[i][2];
|
||||||
double rij_sq = jdelx*jdelx + jdely*jdely + jdelz*jdelz;
|
double rij_sq = jdelx*jdelx + jdely*jdely + jdelz*jdelz;
|
||||||
|
|
||||||
if(rij_sq < cutoff*cutoff) {
|
if(rij_sq < cutoff*cutoff) {
|
||||||
double rij = sqrt(rij_sq);
|
double rij = sqrt(rij_sq);
|
||||||
double partial_sum = 0;
|
double partial_sum = 0;
|
||||||
const int jtype = atom->type[j];
|
const int jtype = atom->type[j];
|
||||||
|
|
||||||
nextTwoBodyInfo->tag = j;
|
nextTwoBodyInfo->tag = j;
|
||||||
nextTwoBodyInfo->r = rij;
|
nextTwoBodyInfo->r = rij;
|
||||||
nextTwoBodyInfo->f = fs[i_to_potl(jtype)].eval(rij, nextTwoBodyInfo->fprime);
|
nextTwoBodyInfo->f = fs[i_to_potl(jtype)].eval(rij, nextTwoBodyInfo->fprime);
|
||||||
nextTwoBodyInfo->del[0] = jdelx / rij;
|
nextTwoBodyInfo->del[0] = jdelx / rij;
|
||||||
nextTwoBodyInfo->del[1] = jdely / rij;
|
nextTwoBodyInfo->del[1] = jdely / rij;
|
||||||
nextTwoBodyInfo->del[2] = jdelz / rij;
|
nextTwoBodyInfo->del[2] = jdelz / rij;
|
||||||
|
|
||||||
for(int kk = 0; kk < numBonds; kk++) {
|
for(int kk = 0; kk < numBonds; kk++) {
|
||||||
const MEAM2Body& bondk = twoBodyInfo[kk];
|
const MEAM2Body& bondk = twoBodyInfo[kk];
|
||||||
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
|
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
|
||||||
@ -176,10 +176,10 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
nextTwoBodyInfo->del[2]*bondk.del[2]);
|
nextTwoBodyInfo->del[2]*bondk.del[2]);
|
||||||
partial_sum += bondk.f * gs[ij_to_potl(jtype,atom->type[bondk.tag],ntypes)].eval(cos_theta);
|
partial_sum += bondk.f * gs[ij_to_potl(jtype,atom->type[bondk.tag],ntypes)].eval(cos_theta);
|
||||||
}
|
}
|
||||||
|
|
||||||
rho_value += nextTwoBodyInfo->f * partial_sum;
|
rho_value += nextTwoBodyInfo->f * partial_sum;
|
||||||
rho_value += rhos[i_to_potl(jtype)].eval(rij);
|
rho_value += rhos[i_to_potl(jtype)].eval(rij);
|
||||||
|
|
||||||
numBonds++;
|
numBonds++;
|
||||||
nextTwoBodyInfo++;
|
nextTwoBodyInfo++;
|
||||||
}
|
}
|
||||||
@ -189,7 +189,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
double Uprime_i;
|
double Uprime_i;
|
||||||
double embeddingEnergy = Us[i_to_potl(itype)].eval(rho_value, Uprime_i)
|
double embeddingEnergy = Us[i_to_potl(itype)].eval(rho_value, Uprime_i)
|
||||||
- zero_atom_energies[i_to_potl(itype)];
|
- zero_atom_energies[i_to_potl(itype)];
|
||||||
|
|
||||||
Uprime_values[i] = Uprime_i;
|
Uprime_values[i] = Uprime_i;
|
||||||
if(eflag) {
|
if(eflag) {
|
||||||
if(eflag_global)
|
if(eflag_global)
|
||||||
@ -204,17 +204,17 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
const MEAM2Body bondj = twoBodyInfo[jj];
|
const MEAM2Body bondj = twoBodyInfo[jj];
|
||||||
double rij = bondj.r;
|
double rij = bondj.r;
|
||||||
int j = bondj.tag;
|
int j = bondj.tag;
|
||||||
|
|
||||||
double f_rij_prime = bondj.fprime;
|
double f_rij_prime = bondj.fprime;
|
||||||
double f_rij = bondj.f;
|
double f_rij = bondj.f;
|
||||||
|
|
||||||
double forces_j[3] = {0, 0, 0};
|
double forces_j[3] = {0, 0, 0};
|
||||||
const int jtype = atom->type[j];
|
const int jtype = atom->type[j];
|
||||||
|
|
||||||
MEAM2Body const* bondk = twoBodyInfo;
|
MEAM2Body const* bondk = twoBodyInfo;
|
||||||
for(int kk = 0; kk < jj; kk++, ++bondk) {
|
for(int kk = 0; kk < jj; kk++, ++bondk) {
|
||||||
double rik = bondk->r;
|
double rik = bondk->r;
|
||||||
|
|
||||||
double cos_theta = (bondj.del[0]*bondk->del[0] +
|
double cos_theta = (bondj.del[0]*bondk->del[0] +
|
||||||
bondj.del[1]*bondk->del[1] +
|
bondj.del[1]*bondk->del[1] +
|
||||||
bondj.del[2]*bondk->del[2]);
|
bondj.del[2]*bondk->del[2]);
|
||||||
@ -222,37 +222,37 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
double g_value = gs[ij_to_potl(jtype,atom->type[bondk->tag],ntypes)].eval(cos_theta, g_prime);
|
double g_value = gs[ij_to_potl(jtype,atom->type[bondk->tag],ntypes)].eval(cos_theta, g_prime);
|
||||||
double f_rik_prime = bondk->fprime;
|
double f_rik_prime = bondk->fprime;
|
||||||
double f_rik = bondk->f;
|
double f_rik = bondk->f;
|
||||||
|
|
||||||
double fij = -Uprime_i * g_value * f_rik * f_rij_prime;
|
double fij = -Uprime_i * g_value * f_rik * f_rij_prime;
|
||||||
double fik = -Uprime_i * g_value * f_rij * f_rik_prime;
|
double fik = -Uprime_i * g_value * f_rij * f_rik_prime;
|
||||||
|
|
||||||
double prefactor = Uprime_i * f_rij * f_rik * g_prime;
|
double prefactor = Uprime_i * f_rij * f_rik * g_prime;
|
||||||
double prefactor_ij = prefactor / rij;
|
double prefactor_ij = prefactor / rij;
|
||||||
double prefactor_ik = prefactor / rik;
|
double prefactor_ik = prefactor / rik;
|
||||||
fij += prefactor_ij * cos_theta;
|
fij += prefactor_ij * cos_theta;
|
||||||
fik += prefactor_ik * cos_theta;
|
fik += prefactor_ik * cos_theta;
|
||||||
|
|
||||||
double fj[3], fk[3];
|
double fj[3], fk[3];
|
||||||
|
|
||||||
fj[0] = bondj.del[0] * fij - bondk->del[0] * prefactor_ij;
|
fj[0] = bondj.del[0] * fij - bondk->del[0] * prefactor_ij;
|
||||||
fj[1] = bondj.del[1] * fij - bondk->del[1] * prefactor_ij;
|
fj[1] = bondj.del[1] * fij - bondk->del[1] * prefactor_ij;
|
||||||
fj[2] = bondj.del[2] * fij - bondk->del[2] * prefactor_ij;
|
fj[2] = bondj.del[2] * fij - bondk->del[2] * prefactor_ij;
|
||||||
forces_j[0] += fj[0];
|
forces_j[0] += fj[0];
|
||||||
forces_j[1] += fj[1];
|
forces_j[1] += fj[1];
|
||||||
forces_j[2] += fj[2];
|
forces_j[2] += fj[2];
|
||||||
|
|
||||||
fk[0] = bondk->del[0] * fik - bondj.del[0] * prefactor_ik;
|
fk[0] = bondk->del[0] * fik - bondj.del[0] * prefactor_ik;
|
||||||
fk[1] = bondk->del[1] * fik - bondj.del[1] * prefactor_ik;
|
fk[1] = bondk->del[1] * fik - bondj.del[1] * prefactor_ik;
|
||||||
fk[2] = bondk->del[2] * fik - bondj.del[2] * prefactor_ik;
|
fk[2] = bondk->del[2] * fik - bondj.del[2] * prefactor_ik;
|
||||||
forces_i[0] -= fk[0];
|
forces_i[0] -= fk[0];
|
||||||
forces_i[1] -= fk[1];
|
forces_i[1] -= fk[1];
|
||||||
forces_i[2] -= fk[2];
|
forces_i[2] -= fk[2];
|
||||||
|
|
||||||
int k = bondk->tag;
|
int k = bondk->tag;
|
||||||
forces[k][0] += fk[0];
|
forces[k][0] += fk[0];
|
||||||
forces[k][1] += fk[1];
|
forces[k][1] += fk[1];
|
||||||
forces[k][2] += fk[2];
|
forces[k][2] += fk[2];
|
||||||
|
|
||||||
if(evflag) {
|
if(evflag) {
|
||||||
double delta_ij[3];
|
double delta_ij[3];
|
||||||
double delta_ik[3];
|
double delta_ik[3];
|
||||||
@ -265,7 +265,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
ev_tally3(i, j, k, 0.0, 0.0, fj, fk, delta_ij, delta_ik);
|
ev_tally3(i, j, k, 0.0, 0.0, fj, fk, delta_ij, delta_ik);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
forces[i][0] -= forces_j[0];
|
forces[i][0] -= forces_j[0];
|
||||||
forces[i][1] -= forces_j[1];
|
forces[i][1] -= forces_j[1];
|
||||||
forces[i][2] -= forces_j[2];
|
forces[i][2] -= forces_j[2];
|
||||||
@ -273,7 +273,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
forces[j][1] += forces_j[1];
|
forces[j][1] += forces_j[1];
|
||||||
forces[j][2] += forces_j[2];
|
forces[j][2] += forces_j[2];
|
||||||
}
|
}
|
||||||
|
|
||||||
forces[i][0] += forces_i[0];
|
forces[i][0] += forces_i[0];
|
||||||
forces[i][1] += forces_i[1];
|
forces[i][1] += forces_i[1];
|
||||||
forces[i][2] += forces_i[2];
|
forces[i][2] += forces_i[2];
|
||||||
@ -287,17 +287,17 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
for(int ii = 0; ii < listhalf->inum; ii++) {
|
for(int ii = 0; ii < listhalf->inum; ii++) {
|
||||||
int i = listhalf->ilist[ii];
|
int i = listhalf->ilist[ii];
|
||||||
const int itype = atom->type[i];
|
const int itype = atom->type[i];
|
||||||
|
|
||||||
for(int jj = 0; jj < listhalf->numneigh[i]; jj++) {
|
for(int jj = 0; jj < listhalf->numneigh[i]; jj++) {
|
||||||
int j = listhalf->firstneigh[i][jj];
|
int j = listhalf->firstneigh[i][jj];
|
||||||
j &= NEIGHMASK;
|
j &= NEIGHMASK;
|
||||||
|
|
||||||
double jdel[3];
|
double jdel[3];
|
||||||
jdel[0] = x[j][0] - x[i][0];
|
jdel[0] = x[j][0] - x[i][0];
|
||||||
jdel[1] = x[j][1] - x[i][1];
|
jdel[1] = x[j][1] - x[i][1];
|
||||||
jdel[2] = x[j][2] - x[i][2];
|
jdel[2] = x[j][2] - x[i][2];
|
||||||
double rij_sq = jdel[0]*jdel[0] + jdel[1]*jdel[1] + jdel[2]*jdel[2];
|
double rij_sq = jdel[0]*jdel[0] + jdel[1]*jdel[1] + jdel[2]*jdel[2];
|
||||||
|
|
||||||
if(rij_sq < cutoff*cutoff) {
|
if(rij_sq < cutoff*cutoff) {
|
||||||
double rij = sqrt(rij_sq);
|
double rij = sqrt(rij_sq);
|
||||||
const int jtype = atom->type[j];
|
const int jtype = atom->type[j];
|
||||||
@ -327,7 +327,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
if(vflag_fdotr)
|
if(vflag_fdotr)
|
||||||
virial_fdotr_compute();
|
virial_fdotr_compute();
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -410,12 +410,12 @@ void PairMEAMSpline::coeff(int narg, char **arg)
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
// clear setflag since coeff() called once with I,J = * *
|
// clear setflag since coeff() called once with I,J = * *
|
||||||
|
|
||||||
n = atom->ntypes;
|
n = atom->ntypes;
|
||||||
for (int i = 1; i <= n; i++)
|
for (int i = 1; i <= n; i++)
|
||||||
for (int j = i; j <= n; j++)
|
for (int j = i; j <= n; j++)
|
||||||
setflag[i][j] = 0;
|
setflag[i][j] = 0;
|
||||||
|
|
||||||
// set setflag i,j for type pairs where both are mapped to elements
|
// set setflag i,j for type pairs where both are mapped to elements
|
||||||
|
|
||||||
int count = 0;
|
int count = 0;
|
||||||
@ -442,19 +442,19 @@ void PairMEAMSpline::read_file(const char* filename)
|
|||||||
sprintf(str,"Cannot open spline MEAM potential file %s", filename);
|
sprintf(str,"Cannot open spline MEAM potential file %s", filename);
|
||||||
error->one(FLERR,str);
|
error->one(FLERR,str);
|
||||||
}
|
}
|
||||||
|
|
||||||
// Skip first line of file. It's a comment.
|
// Skip first line of file. It's a comment.
|
||||||
char line[MAXLINE];
|
char line[MAXLINE];
|
||||||
char *ptr;
|
char *ptr;
|
||||||
fgets(line, MAXLINE, fp);
|
fgets(line, MAXLINE, fp);
|
||||||
|
|
||||||
// Second line holds potential type ("meam/spline")
|
// Second line holds potential type ("meam/spline")
|
||||||
// in new potential format.
|
// in new potential format.
|
||||||
|
|
||||||
bool isNewFormat = false;
|
bool isNewFormat = false;
|
||||||
fgets(line, MAXLINE, fp);
|
fgets(line, MAXLINE, fp);
|
||||||
ptr = strtok(line, " \t\n\r\f");
|
ptr = strtok(line, " \t\n\r\f");
|
||||||
|
|
||||||
if (strcmp(ptr, "meam/spline") == 0) {
|
if (strcmp(ptr, "meam/spline") == 0) {
|
||||||
isNewFormat = true;
|
isNewFormat = true;
|
||||||
// parse the rest of the line!
|
// parse the rest of the line!
|
||||||
@ -484,13 +484,13 @@ void PairMEAMSpline::read_file(const char* filename)
|
|||||||
rewind(fp);
|
rewind(fp);
|
||||||
fgets(line, MAXLINE, fp);
|
fgets(line, MAXLINE, fp);
|
||||||
}
|
}
|
||||||
|
|
||||||
nmultichoose2 = ((nelements+1)*nelements)/2;
|
nmultichoose2 = ((nelements+1)*nelements)/2;
|
||||||
// allocate!!
|
// allocate!!
|
||||||
allocate();
|
allocate();
|
||||||
|
|
||||||
// Parse spline functions.
|
// Parse spline functions.
|
||||||
|
|
||||||
for (int i = 0; i < nmultichoose2; i++)
|
for (int i = 0; i < nmultichoose2; i++)
|
||||||
phis[i].parse(fp, error, isNewFormat);
|
phis[i].parse(fp, error, isNewFormat);
|
||||||
for (int i = 0; i < nelements; i++)
|
for (int i = 0; i < nelements; i++)
|
||||||
@ -501,7 +501,7 @@ void PairMEAMSpline::read_file(const char* filename)
|
|||||||
fs[i].parse(fp, error, isNewFormat);
|
fs[i].parse(fp, error, isNewFormat);
|
||||||
for (int i = 0; i < nmultichoose2; i++)
|
for (int i = 0; i < nmultichoose2; i++)
|
||||||
gs[i].parse(fp, error, isNewFormat);
|
gs[i].parse(fp, error, isNewFormat);
|
||||||
|
|
||||||
fclose(fp);
|
fclose(fp);
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -532,11 +532,11 @@ void PairMEAMSpline::read_file(const char* filename)
|
|||||||
Us[i].communicate(world, comm->me);
|
Us[i].communicate(world, comm->me);
|
||||||
for (int i = 0; i < nmultichoose2; i++)
|
for (int i = 0; i < nmultichoose2; i++)
|
||||||
gs[i].communicate(world, comm->me);
|
gs[i].communicate(world, comm->me);
|
||||||
|
|
||||||
// Calculate 'zero-point energy' of single atom in vacuum.
|
// Calculate 'zero-point energy' of single atom in vacuum.
|
||||||
for (int i = 0; i < nelements; i++)
|
for (int i = 0; i < nelements; i++)
|
||||||
zero_atom_energies[i] = Us[i].eval(0.0);
|
zero_atom_energies[i] = Us[i].eval(0.0);
|
||||||
|
|
||||||
// Determine maximum cutoff radius of all relevant spline functions.
|
// Determine maximum cutoff radius of all relevant spline functions.
|
||||||
cutoff = 0.0;
|
cutoff = 0.0;
|
||||||
for (int i = 0; i < nmultichoose2; i++)
|
for (int i = 0; i < nmultichoose2; i++)
|
||||||
@ -548,7 +548,7 @@ void PairMEAMSpline::read_file(const char* filename)
|
|||||||
for (int i = 0; i < nelements; i++)
|
for (int i = 0; i < nelements; i++)
|
||||||
if(fs[i].cutoff() > cutoff)
|
if(fs[i].cutoff() > cutoff)
|
||||||
cutoff = fs[i].cutoff();
|
cutoff = fs[i].cutoff();
|
||||||
|
|
||||||
// Set LAMMPS pair interaction flags.
|
// Set LAMMPS pair interaction flags.
|
||||||
for(int i = 1; i <= atom->ntypes; i++) {
|
for(int i = 1; i <= atom->ntypes; i++) {
|
||||||
for(int j = 1; j <= atom->ntypes; j++) {
|
for(int j = 1; j <= atom->ntypes; j++) {
|
||||||
@ -556,7 +556,7 @@ void PairMEAMSpline::read_file(const char* filename)
|
|||||||
cutsq[i][j] = cutoff;
|
cutsq[i][j] = cutoff;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/* ----------------------------------------------------------------------
|
/* ----------------------------------------------------------------------
|
||||||
@ -643,27 +643,27 @@ void PairMEAMSpline::SplineFunction::parse(FILE* fp, Error* error,
|
|||||||
bool isNewFormat)
|
bool isNewFormat)
|
||||||
{
|
{
|
||||||
char line[MAXLINE];
|
char line[MAXLINE];
|
||||||
|
|
||||||
// If new format, read the spline format. Should always be "spline3eq" for now.
|
// If new format, read the spline format. Should always be "spline3eq" for now.
|
||||||
if (isNewFormat)
|
if (isNewFormat)
|
||||||
fgets(line, MAXLINE, fp);
|
fgets(line, MAXLINE, fp);
|
||||||
|
|
||||||
// Parse number of spline knots.
|
// Parse number of spline knots.
|
||||||
fgets(line, MAXLINE, fp);
|
fgets(line, MAXLINE, fp);
|
||||||
int n = atoi(line);
|
int n = atoi(line);
|
||||||
if(n < 2)
|
if(n < 2)
|
||||||
error->one(FLERR,"Invalid number of spline knots in MEAM potential file");
|
error->one(FLERR,"Invalid number of spline knots in MEAM potential file");
|
||||||
|
|
||||||
// Parse first derivatives at beginning and end of spline.
|
// Parse first derivatives at beginning and end of spline.
|
||||||
fgets(line, MAXLINE, fp);
|
fgets(line, MAXLINE, fp);
|
||||||
double d0 = atof(strtok(line, " \t\n\r\f"));
|
double d0 = atof(strtok(line, " \t\n\r\f"));
|
||||||
double dN = atof(strtok(NULL, " \t\n\r\f"));
|
double dN = atof(strtok(NULL, " \t\n\r\f"));
|
||||||
init(n, d0, dN);
|
init(n, d0, dN);
|
||||||
|
|
||||||
// Skip line in old format
|
// Skip line in old format
|
||||||
if (!isNewFormat)
|
if (!isNewFormat)
|
||||||
fgets(line, MAXLINE, fp);
|
fgets(line, MAXLINE, fp);
|
||||||
|
|
||||||
// Parse knot coordinates.
|
// Parse knot coordinates.
|
||||||
for(int i=0; i<n; i++) {
|
for(int i=0; i<n; i++) {
|
||||||
fgets(line, MAXLINE, fp);
|
fgets(line, MAXLINE, fp);
|
||||||
@ -673,7 +673,7 @@ void PairMEAMSpline::SplineFunction::parse(FILE* fp, Error* error,
|
|||||||
}
|
}
|
||||||
setKnot(i, x, y);
|
setKnot(i, x, y);
|
||||||
}
|
}
|
||||||
|
|
||||||
prepareSpline(error);
|
prepareSpline(error);
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -682,11 +682,11 @@ void PairMEAMSpline::SplineFunction::prepareSpline(Error* error)
|
|||||||
{
|
{
|
||||||
xmin = X[0];
|
xmin = X[0];
|
||||||
xmax = X[N-1];
|
xmax = X[N-1];
|
||||||
|
|
||||||
isGridSpline = true;
|
isGridSpline = true;
|
||||||
h = (xmax-xmin)/(N-1);
|
h = (xmax-xmin)/(N-1);
|
||||||
hsq = h*h;
|
hsq = h*h;
|
||||||
|
|
||||||
double* u = new double[N];
|
double* u = new double[N];
|
||||||
Y2[0] = -0.5;
|
Y2[0] = -0.5;
|
||||||
u[0] = (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - deriv0);
|
u[0] = (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - deriv0);
|
||||||
@ -696,20 +696,20 @@ void PairMEAMSpline::SplineFunction::prepareSpline(Error* error)
|
|||||||
Y2[i] = (sig - 1.0) / p;
|
Y2[i] = (sig - 1.0) / p;
|
||||||
u[i] = (Y[i+1]-Y[i]) / (X[i+1]-X[i]) - (Y[i]-Y[i-1])/(X[i]-X[i-1]);
|
u[i] = (Y[i+1]-Y[i]) / (X[i+1]-X[i]) - (Y[i]-Y[i-1])/(X[i]-X[i-1]);
|
||||||
u[i] = (6.0 * u[i]/(X[i+1]-X[i-1]) - sig*u[i-1])/p;
|
u[i] = (6.0 * u[i]/(X[i+1]-X[i-1]) - sig*u[i-1])/p;
|
||||||
|
|
||||||
if(fabs(h*i+xmin - X[i]) > 1e-8)
|
if(fabs(h*i+xmin - X[i]) > 1e-8)
|
||||||
isGridSpline = false;
|
isGridSpline = false;
|
||||||
}
|
}
|
||||||
|
|
||||||
double qn = 0.5;
|
double qn = 0.5;
|
||||||
double un = (3.0/(X[N-1]-X[N-2])) * (derivN - (Y[N-1]-Y[N-2])/(X[N-1]-X[N-2]));
|
double un = (3.0/(X[N-1]-X[N-2])) * (derivN - (Y[N-1]-Y[N-2])/(X[N-1]-X[N-2]));
|
||||||
Y2[N-1] = (un - qn*u[N-2]) / (qn * Y2[N-2] + 1.0);
|
Y2[N-1] = (un - qn*u[N-2]) / (qn * Y2[N-2] + 1.0);
|
||||||
for(int k = N-2; k >= 0; k--) {
|
for(int k = N-2; k >= 0; k--) {
|
||||||
Y2[k] = Y2[k] * Y2[k+1] + u[k];
|
Y2[k] = Y2[k] * Y2[k+1] + u[k];
|
||||||
}
|
}
|
||||||
|
|
||||||
delete[] u;
|
delete[] u;
|
||||||
|
|
||||||
#if !SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES
|
#if !SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES
|
||||||
if(!isGridSpline)
|
if(!isGridSpline)
|
||||||
error->one(FLERR,"Support for MEAM potentials with non-uniform cubic splines has not been enabled in the MEAM potential code. Set SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES in pair_spline_meam.h to 1 to enable it");
|
error->one(FLERR,"Support for MEAM potentials with non-uniform cubic splines has not been enabled in the MEAM potential code. Set SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES in pair_spline_meam.h to 1 to enable it");
|
||||||
|
|||||||
@ -51,21 +51,21 @@ public:
|
|||||||
void init_style();
|
void init_style();
|
||||||
void init_list(int, class NeighList *);
|
void init_list(int, class NeighList *);
|
||||||
double init_one(int, int);
|
double init_one(int, int);
|
||||||
|
|
||||||
// helper functions for compute()
|
// helper functions for compute()
|
||||||
|
|
||||||
int ij_to_potl(const int itype, const int jtype, const int ntypes) const {
|
int ij_to_potl(const int itype, const int jtype, const int ntypes) const {
|
||||||
return jtype - 1 + (itype-1)*ntypes - (itype-1)*itype/2;
|
return jtype - 1 + (itype-1)*ntypes - (itype-1)*itype/2;
|
||||||
}
|
}
|
||||||
int i_to_potl(const int itype) const { return itype-1; }
|
int i_to_potl(const int itype) const { return itype-1; }
|
||||||
|
|
||||||
|
|
||||||
int pack_forward_comm(int, int *, double *, int, int *);
|
int pack_forward_comm(int, int *, double *, int, int *);
|
||||||
void unpack_forward_comm(int, int, double *);
|
void unpack_forward_comm(int, int, double *);
|
||||||
int pack_reverse_comm(int, int, double *);
|
int pack_reverse_comm(int, int, double *);
|
||||||
void unpack_reverse_comm(int, int *, double *);
|
void unpack_reverse_comm(int, int *, double *);
|
||||||
double memory_usage();
|
double memory_usage();
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
char **elements; // names of unique elements
|
char **elements; // names of unique elements
|
||||||
int *map; // mapping from atom types to elements
|
int *map; // mapping from atom types to elements
|
||||||
@ -75,7 +75,7 @@ protected:
|
|||||||
public:
|
public:
|
||||||
/// Default constructor.
|
/// Default constructor.
|
||||||
SplineFunction() : X(NULL), Xs(NULL), Y(NULL), Y2(NULL), Ydelta(NULL), N(0) {}
|
SplineFunction() : X(NULL), Xs(NULL), Y(NULL), Y2(NULL), Ydelta(NULL), N(0) {}
|
||||||
|
|
||||||
/// Destructor.
|
/// Destructor.
|
||||||
~SplineFunction() {
|
~SplineFunction() {
|
||||||
delete[] X;
|
delete[] X;
|
||||||
@ -84,7 +84,7 @@ protected:
|
|||||||
delete[] Y2;
|
delete[] Y2;
|
||||||
delete[] Ydelta;
|
delete[] Ydelta;
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Initialization of spline function.
|
/// Initialization of spline function.
|
||||||
void init(int _N, double _deriv0, double _derivN) {
|
void init(int _N, double _deriv0, double _derivN) {
|
||||||
N = _N;
|
N = _N;
|
||||||
@ -101,19 +101,19 @@ protected:
|
|||||||
Y2 = new double[N];
|
Y2 = new double[N];
|
||||||
Ydelta = new double[N];
|
Ydelta = new double[N];
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Adds a knot to the spline.
|
/// Adds a knot to the spline.
|
||||||
void setKnot(int n, double x, double y) { X[n] = x; Y[n] = y; }
|
void setKnot(int n, double x, double y) { X[n] = x; Y[n] = y; }
|
||||||
|
|
||||||
/// Returns the number of knots.
|
/// Returns the number of knots.
|
||||||
int numKnots() const { return N; }
|
int numKnots() const { return N; }
|
||||||
|
|
||||||
/// Parses the spline knots from a text file.
|
/// Parses the spline knots from a text file.
|
||||||
void parse(FILE* fp, Error* error, bool isNewFormat);
|
void parse(FILE* fp, Error* error, bool isNewFormat);
|
||||||
|
|
||||||
/// Calculates the second derivatives of the cubic spline.
|
/// Calculates the second derivatives of the cubic spline.
|
||||||
void prepareSpline(Error* error);
|
void prepareSpline(Error* error);
|
||||||
|
|
||||||
/// Evaluates the spline function at position x.
|
/// Evaluates the spline function at position x.
|
||||||
inline double eval(double x) const
|
inline double eval(double x) const
|
||||||
{
|
{
|
||||||
@ -151,7 +151,7 @@ protected:
|
|||||||
#endif
|
#endif
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Evaluates the spline function and its first derivative at position x.
|
/// Evaluates the spline function and its first derivative at position x.
|
||||||
inline double eval(double x, double& deriv) const
|
inline double eval(double x, double& deriv) const
|
||||||
{
|
{
|
||||||
@ -197,16 +197,16 @@ protected:
|
|||||||
#endif
|
#endif
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Returns the number of bytes used by this function object.
|
/// Returns the number of bytes used by this function object.
|
||||||
double memory_usage() const { return sizeof(*this) + sizeof(X[0]) * N * 3; }
|
double memory_usage() const { return sizeof(*this) + sizeof(X[0]) * N * 3; }
|
||||||
|
|
||||||
/// Returns the cutoff radius of this function.
|
/// Returns the cutoff radius of this function.
|
||||||
double cutoff() const { return X[N-1]; }
|
double cutoff() const { return X[N-1]; }
|
||||||
|
|
||||||
/// Writes a Gnuplot script that plots the spline function.
|
/// Writes a Gnuplot script that plots the spline function.
|
||||||
void writeGnuplot(const char* filename, const char* title = NULL) const;
|
void writeGnuplot(const char* filename, const char* title = NULL) const;
|
||||||
|
|
||||||
/// Broadcasts the spline function parameters to all processors.
|
/// Broadcasts the spline function parameters to all processors.
|
||||||
void communicate(MPI_Comm& world, int me);
|
void communicate(MPI_Comm& world, int me);
|
||||||
|
|
||||||
@ -226,7 +226,7 @@ protected:
|
|||||||
double hsq; // The squared distance between knots if this is a grid spline with equidistant knots.
|
double hsq; // The squared distance between knots if this is a grid spline with equidistant knots.
|
||||||
double xmax_shifted; // The end of the spline interval after it has been shifted to begin at X=0.
|
double xmax_shifted; // The end of the spline interval after it has been shifted to begin at X=0.
|
||||||
};
|
};
|
||||||
|
|
||||||
/// Helper data structure for potential routine.
|
/// Helper data structure for potential routine.
|
||||||
struct MEAM2Body {
|
struct MEAM2Body {
|
||||||
int tag; // holds the index of the second atom (j)
|
int tag; // holds the index of the second atom (j)
|
||||||
@ -234,26 +234,26 @@ protected:
|
|||||||
double f, fprime;
|
double f, fprime;
|
||||||
double del[3];
|
double del[3];
|
||||||
};
|
};
|
||||||
|
|
||||||
SplineFunction* phis; // Phi_i(r_ij)
|
SplineFunction* phis; // Phi_i(r_ij)
|
||||||
SplineFunction* rhos; // Rho_ij(r_ij)
|
SplineFunction* rhos; // Rho_ij(r_ij)
|
||||||
SplineFunction* fs; // f_i(r_ij)
|
SplineFunction* fs; // f_i(r_ij)
|
||||||
SplineFunction* Us; // U_i(rho)
|
SplineFunction* Us; // U_i(rho)
|
||||||
SplineFunction* gs; // g_ij(cos_theta)
|
SplineFunction* gs; // g_ij(cos_theta)
|
||||||
double* zero_atom_energies; // Shift embedding energy by this value to make it zero for a single atom in vacuum.
|
double* zero_atom_energies; // Shift embedding energy by this value to make it zero for a single atom in vacuum.
|
||||||
|
|
||||||
double cutoff; // The cutoff radius
|
double cutoff; // The cutoff radius
|
||||||
|
|
||||||
double* Uprime_values; // Used for temporary storage of U'(rho) values
|
double* Uprime_values; // Used for temporary storage of U'(rho) values
|
||||||
int nmax; // Size of temporary array.
|
int nmax; // Size of temporary array.
|
||||||
int maxNeighbors; // The last maximum number of neighbors a single atoms has.
|
int maxNeighbors; // The last maximum number of neighbors a single atoms has.
|
||||||
MEAM2Body* twoBodyInfo; // Temporary array.
|
MEAM2Body* twoBodyInfo; // Temporary array.
|
||||||
|
|
||||||
void read_file(const char* filename);
|
void read_file(const char* filename);
|
||||||
void allocate();
|
void allocate();
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
|||||||
@ -305,7 +305,7 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
|
|||||||
rhos[i_to_potl(itype)].eval(rij,rho_prime_i);
|
rhos[i_to_potl(itype)].eval(rij,rho_prime_i);
|
||||||
rhos[i_to_potl(jtype)].eval(rij,rho_prime_j);
|
rhos[i_to_potl(jtype)].eval(rij,rho_prime_j);
|
||||||
double fpair = rho_prime_j * Uprime_values[i] + rho_prime_i*Uprime_values[j];
|
double fpair = rho_prime_j * Uprime_values[i] + rho_prime_i*Uprime_values[j];
|
||||||
|
|
||||||
double pair_pot_deriv;
|
double pair_pot_deriv;
|
||||||
double pair_pot = phis[ij_to_potl(itype,jtype,ntypes)].eval(rij, pair_pot_deriv);
|
double pair_pot = phis[ij_to_potl(itype,jtype,ntypes)].eval(rij, pair_pot_deriv);
|
||||||
|
|
||||||
|
|||||||
Reference in New Issue
Block a user