From fe3c6670f6b42bcc8eeeaa53ebdd4a9a45e33ae4 Mon Sep 17 00:00:00 2001 From: athomps Date: Wed, 23 Sep 2015 00:30:30 +0000 Subject: [PATCH] Modifed ZBL to accomodate mixing git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14047 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/pair_zbl.cpp | 138 ++++++++++++++++++++++++++--------------------- src/pair_zbl.h | 1 + 2 files changed, 79 insertions(+), 60 deletions(-) diff --git a/src/pair_zbl.cpp b/src/pair_zbl.cpp index ba9208ddf2..57278c7a3a 100644 --- a/src/pair_zbl.cpp +++ b/src/pair_zbl.cpp @@ -199,11 +199,14 @@ void PairZBL::settings(int narg, char **arg) /* ---------------------------------------------------------------------- set coeffs for one or more type pairs + ------------------------------------------------------------------------- */ void PairZBL::coeff(int narg, char **arg) { - if (narg != 3) + double z_one, z_two; + + if (narg != 4) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); @@ -213,7 +216,8 @@ void PairZBL::coeff(int narg, char **arg) int jlo,jhi; force->bounds(arg[1],atom->ntypes,jlo,jhi); - double z_one = force->numeric(FLERR,arg[2]); + z_one = force->numeric(FLERR,arg[2]); + z_two = force->numeric(FLERR,arg[3]); // set flag for each i-j pair // set z-parameter only for i-i pairs @@ -221,8 +225,13 @@ void PairZBL::coeff(int narg, char **arg) int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - if (i == j) z[i] = z_one; + if (i == j) { + if (z_one != z_two) + error->all(FLERR,"Incorrect args for pair coefficients"); + z[i] = z_one; + } setflag[i][j] = 1; + set_coeff(i, j, z_one, z_two); count++; } } @@ -248,63 +257,8 @@ void PairZBL::init_style() double PairZBL::init_one(int i, int j) { - - double ainv = (pow(z[i],pzbl) + pow(z[j],pzbl))/(a0*force->angstrom); - d1a[i][j] = d1*ainv; - d2a[i][j] = d2*ainv; - d3a[i][j] = d3*ainv; - d4a[i][j] = d4*ainv; - zze[i][j] = z[i]*z[j]*force->qqr2e*force->qelectron*force->qelectron; - - d1a[j][i] = d1a[i][j]; - d2a[j][i] = d2a[i][j]; - d3a[j][i] = d3a[i][j]; - d4a[j][i] = d4a[i][j]; - zze[j][i] = zze[i][j]; - - // e = t^3 (sw3 + sw4*t) + sw5 - // = A/3*t^3 + B/4*t^4 + C - // sw3 = A/3 - // sw4 = B/4 - // sw5 = C - - // dedr = t^2 (sw1 + sw2*t) - // = A*t^2 + B*t^3 - // sw1 = A - // sw2 = B - - // de2dr2 = 2*A*t + 3*B*t^2 - - // Require that at t = tc: - // e = -Fc - // dedr = -Fc' - // d2edr2 = -Fc'' - - // Hence: - // A = (-3Fc' + tc*Fc'')/tc^2 - // B = ( 2Fc' - tc*Fc'')/tc^3 - // C = -Fc + tc/2*Fc' - tc^2/12*Fc'' - - double tc = cut_global - cut_inner; - double fc = e_zbl(cut_global, i, j); - double fcp = dzbldr(cut_global, i, j); - double fcpp = d2zbldr2(cut_global, i, j); - - double swa = (-3.0*fcp + tc*fcpp)/(tc*tc); - double swb = ( 2.0*fcp - tc*fcpp)/(tc*tc*tc); - double swc = -fc + (tc/2.0)*fcp - (tc*tc/12.0)*fcpp; - - sw1[i][j] = swa; - sw2[i][j] = swb; - sw3[i][j] = swa/3.0; - sw4[i][j] = swb/4.0; - sw5[i][j] = swc; - - sw1[j][i] = sw1[i][j]; - sw2[j][i] = sw2[i][j]; - sw3[j][i] = sw3[i][j]; - sw4[j][i] = sw4[i][j]; - sw5[j][i] = sw5[i][j]; + if (setflag[i][j] == 0) + set_coeff(i, j, z[i], z[j]); return cut_global; } @@ -434,3 +388,67 @@ double PairZBL::d2zbldr2(double r, int i, int j) { return result; } +/* ---------------------------------------------------------------------- + calculate the i,j entries in the various coeff arrays +------------------------------------------------------------------------- */ + +void PairZBL::set_coeff(int i, int j, double zi, double zj) +{ + double ainv = (pow(zi,pzbl) + pow(zj,pzbl))/(a0*force->angstrom); + d1a[i][j] = d1*ainv; + d2a[i][j] = d2*ainv; + d3a[i][j] = d3*ainv; + d4a[i][j] = d4*ainv; + zze[i][j] = zi*zj*force->qqr2e*force->qelectron*force->qelectron; + + d1a[j][i] = d1a[i][j]; + d2a[j][i] = d2a[i][j]; + d3a[j][i] = d3a[i][j]; + d4a[j][i] = d4a[i][j]; + zze[j][i] = zze[i][j]; + + // e = t^3 (sw3 + sw4*t) + sw5 + // = A/3*t^3 + B/4*t^4 + C + // sw3 = A/3 + // sw4 = B/4 + // sw5 = C + + // dedr = t^2 (sw1 + sw2*t) + // = A*t^2 + B*t^3 + // sw1 = A + // sw2 = B + + // de2dr2 = 2*A*t + 3*B*t^2 + + // Require that at t = tc: + // e = -Fc + // dedr = -Fc' + // d2edr2 = -Fc'' + + // Hence: + // A = (-3Fc' + tc*Fc'')/tc^2 + // B = ( 2Fc' - tc*Fc'')/tc^3 + // C = -Fc + tc/2*Fc' - tc^2/12*Fc'' + + double tc = cut_global - cut_inner; + double fc = e_zbl(cut_global, i, j); + double fcp = dzbldr(cut_global, i, j); + double fcpp = d2zbldr2(cut_global, i, j); + + double swa = (-3.0*fcp + tc*fcpp)/(tc*tc); + double swb = ( 2.0*fcp - tc*fcpp)/(tc*tc*tc); + double swc = -fc + (tc/2.0)*fcp - (tc*tc/12.0)*fcpp; + + sw1[i][j] = swa; + sw2[i][j] = swb; + sw3[i][j] = swa/3.0; + sw4[i][j] = swb/4.0; + sw5[i][j] = swc; + + sw1[j][i] = sw1[i][j]; + sw2[j][i] = sw2[i][j]; + sw3[j][i] = sw3[i][j]; + sw4[j][i] = sw4[i][j]; + sw5[j][i] = sw5[i][j]; +} + diff --git a/src/pair_zbl.h b/src/pair_zbl.h index 5aea8ebefe..0dc01731d3 100644 --- a/src/pair_zbl.h +++ b/src/pair_zbl.h @@ -46,6 +46,7 @@ class PairZBL : public Pair { double e_zbl(double, int, int); double dzbldr(double, int, int); double d2zbldr2(double, int, int); + void set_coeff(int, int, double, double); }; namespace PairZBLConstants {