Commit 7efb42f0 authored by Sebastian Hütter's avatar Sebastian Hütter
Browse files

MEAM/C: code optimization for faster codegen

parent 0d4bb861
Loading
Loading
Loading
Loading
+79 −79
Original line number Original line Diff line number Diff line
@@ -127,33 +127,39 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
    delyij = yjtmp - yitmp;
    delyij = yjtmp - yitmp;
    delzij = zjtmp - zitmp;
    delzij = zjtmp - zitmp;
    rij2 = delxij * delxij + delyij * delyij + delzij * delzij;
    rij2 = delxij * delxij + delyij * delyij + delzij * delzij;
    rij = sqrt(rij2);

    if (rij2 > this->cutforcesq) {
      dscrfcn[jn] = 0.0;
      scrfcn[jn] = 0.0;
      fcpair[jn] = 0.0;
      continue;
    }


    const double rbound = this->ebound_meam[elti][eltj] * rij2;
    const double rbound = this->ebound_meam[elti][eltj] * rij2;
    if (rij > this->rc_meam) {
    rij = sqrt(rij2);
      fcij = 0.0;
    rnorm = (this->cutforce - rij) * drinv;
      dfcij = 0.0;
      sij = 0.0;
    } else {
      rnorm = (this->rc_meam - rij) * drinv;
    sij = 1.0;
    sij = 1.0;


    //     if rjk2 > ebound*rijsq, atom k is definitely outside the ellipse
    //     if rjk2 > ebound*rijsq, atom k is definitely outside the ellipse
    for (kn = 0; kn < numneigh_full; kn++) {
    for (kn = 0; kn < numneigh_full; kn++) {
      k = firstneigh_full[kn];
      k = firstneigh_full[kn];
      if (k == j) continue;
      eltk = fmap[type[k]];
      eltk = fmap[type[k]];
      if (eltk < 0) continue;
      if (eltk < 0) continue;
        if (k == j) continue;


        delxjk = x[k][0] - xjtmp;
      xktmp = x[k][0];
        delyjk = x[k][1] - yjtmp;
      yktmp = x[k][1];
        delzjk = x[k][2] - zjtmp;
      zktmp = x[k][2];

      delxjk = xktmp - xjtmp;
      delyjk = yktmp - yjtmp;
      delzjk = zktmp - zjtmp;
      rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
      rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
      if (rjk2 > rbound) continue;
      if (rjk2 > rbound) continue;


        delxik = x[k][0] - xitmp;
      delxik = xktmp - xitmp;
        delyik = x[k][1] - yitmp;
      delyik = yktmp - yitmp;
        delzik = x[k][2] - zitmp;
      delzik = zktmp - zitmp;
      rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
      rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
      if (rik2 > rbound) continue;
      if (rik2 > rbound) continue;


@@ -185,32 +191,26 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
    fc = dfcut(rnorm, dfc);
    fc = dfcut(rnorm, dfc);
    fcij = fc;
    fcij = fc;
    dfcij = dfc * drinv;
    dfcij = dfc * drinv;
    }


    //     Now compute derivatives
    //     Now compute derivatives
    dscrfcn[jn] = 0.0;
    dscrfcn[jn] = 0.0;
    sfcij = sij * fcij;
    sfcij = sij * fcij;
    if (iszero(sfcij) || iszero(sfcij - 1.0))
    if (!iszero(sfcij) && !iszero(sfcij - 1.0)) {
      goto LABEL_100;

      for (kn = 0; kn < numneigh_full; kn++) {
      for (kn = 0; kn < numneigh_full; kn++) {
        k = firstneigh_full[kn];
        k = firstneigh_full[kn];
        if (k == j) continue;
        if (k == j) continue;
        eltk = fmap[type[k]];
        eltk = fmap[type[k]];
        if (eltk < 0) continue;
        if (eltk < 0) continue;


      xktmp = x[k][0];
        delxjk = x[k][0] - xjtmp;
      yktmp = x[k][1];
        delyjk = x[k][1] - yjtmp;
      zktmp = x[k][2];
        delzjk = x[k][2] - zjtmp;
      delxjk = xktmp - xjtmp;
      delyjk = yktmp - yjtmp;
      delzjk = zktmp - zjtmp;
        rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
        rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
        if (rjk2 > rbound) continue;
        if (rjk2 > rbound) continue;


      delxik = xktmp - xitmp;
        delxik = x[k][0] - xitmp;
      delyik = yktmp - yitmp;
        delyik = x[k][1] - yitmp;
      delzik = zktmp - zitmp;
        delzik = x[k][2] - zitmp;
        rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
        rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
        if (rik2 > rbound) continue;
        if (rik2 > rbound) continue;


@@ -245,8 +245,8 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
      coef1 = sfcij;
      coef1 = sfcij;
      coef2 = sij * dfcij / rij;
      coef2 = sij * dfcij / rij;
      dscrfcn[jn] = dscrfcn[jn] * coef1 - coef2;
      dscrfcn[jn] = dscrfcn[jn] * coef1 - coef2;
    }


    LABEL_100:
    scrfcn[jn] = sij;
    scrfcn[jn] = sij;
    fcpair[jn] = fcij;
    fcpair[jn] = fcij;
  }
  }
+4 −3
Original line number Original line Diff line number Diff line
@@ -52,7 +52,7 @@ MEAM::G_gam(const double gamma, const int ibar, int& errorflag) const
    case 1:
    case 1:
      return MathSpecial::fm_exp(gamma / 2.0);
      return MathSpecial::fm_exp(gamma / 2.0);
    case 3:
    case 3:
      return 2.0 / (1.0 + exp(-gamma));
      return 2.0 / (1.0 + MathSpecial::fm_exp(-gamma));
    case -5:
    case -5:
      if ((1.0 + gamma) >= 0) {
      if ((1.0 + gamma) >= 0) {
        return sqrt(1.0 + gamma);
        return sqrt(1.0 + gamma);
@@ -152,8 +152,9 @@ MEAM::embedding(const double A, const double Ec, const double rhobar, double& dF
  const double AEc = A * Ec;
  const double AEc = A * Ec;


  if (rhobar > 0.0) {
  if (rhobar > 0.0) {
      dF = AEc * (1.0 + log(rhobar));
      const double lrb = log(rhobar);
      return AEc * rhobar * log(rhobar);
      dF = AEc * (1.0 + lrb);
      return AEc * rhobar * lrb;
  } else {
  } else {
    if (this->emb_lin_neg == 0) {
    if (this->emb_lin_neg == 0) {
      dF = 0.0;
      dF = 0.0;
+1 −1
Original line number Original line Diff line number Diff line
@@ -697,7 +697,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
      get_sijk(C, a, a, b, &s112);
      get_sijk(C, a, a, b, &s112);
      get_sijk(C, b, b, a, &s221);
      get_sijk(C, b, b, a, &s221);
      S11 = s111 * s111 * s112 * s112;
      S11 = s111 * s111 * s112 * s112;
      S22 = pow(s221, 4);
      S22 = s221 * s221 * s221 * s221;
      *rho01 = *rho01 + 6 * S11 * rhoa01nn;
      *rho01 = *rho01 + 6 * S11 * rhoa01nn;
      *rho02 = *rho02 + 6 * S22 * rhoa02nn;
      *rho02 = *rho02 + 6 * S22 * rhoa02nn;


+1 −1

File changed.

Contains only whitespace changes.