Commit 64cf52d3 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

address spline out-of-bounds bug reported in issue #59 and refactor high-level... 

address spline out-of-bounds bug reported in issue #59 and refactor high-level spline code for consistency and efficiency
parent 06c15142

src/MANYBODY/pair_airebo.cpp

+117 −174
Original line number Diff line number Diff line
@@ -3094,72 +3094,54 @@ double PairAIREBO::gSpline(double costh, double Nij, int typei, 
double PairAIREBO::PijSpline(double NijC, double NijH, int typei, int typej,

                             double dN2[2])

{

  int x,y,i,done;

  double Pij,coeffs[16];



  for (i = 0; i < 16; i++) coeffs[i]=0.0;

  int x,y;

  double Pij;



  x = 0;

  y = 0;

  dN2[0] = 0.0;

  dN2[1] = 0.0;

  done = 0;

  Pij = 0.0;



  if (typei == 1) return Pij;



  if (typej == 0) {



    // if inputs are out of bounds set them back to a point in bounds



  if (typei == 0 && typej == 0) {

    if (NijC <  pCCdom[0][0]) NijC=pCCdom[0][0];

    if (NijC > pCCdom[0][1]) NijC=pCCdom[0][1];

    if (NijC >= pCCdom[0][1]) NijC=pCCdom[0][1]-TOL;

    if (NijH <  pCCdom[1][0]) NijH=pCCdom[1][0];

    if (NijH > pCCdom[1][1]) NijH=pCCdom[1][1];

    if (NijH >= pCCdom[1][1]) NijH=pCCdom[1][1]-TOL;

    x = (int) floor(NijC);

    y = (int) floor(NijH);



    if (fabs(NijC-floor(NijC)) < TOL && fabs(NijH-floor(NijH)) < TOL) {

      Pij = PCCf[(int) NijC][(int) NijH];

      dN2[0] = PCCdfdx[(int) NijC][(int) NijH];

      dN2[1] = PCCdfdy[(int) NijC][(int) NijH];

      done = 1;

    }

    if (done == 0) {

      x = (int) (floor(NijC));

      y = (int) (floor(NijH));

      for (i = 0; i<16; i++) coeffs[i] = pCC[x][y][i];

      Pij = Spbicubic(NijC,NijH,coeffs,dN2);

    }

      Pij    = PCCf[x][y];

      dN2[0] = PCCdfdx[x][y];

      dN2[1] = PCCdfdy[x][y];

    } else {

      Pij = Spbicubic(NijC,NijH,pCC[x][y],dN2);

    }



  } else if (typej == 1) {



    // if inputs are out of bounds set them back to a point in bounds



   if (typei == 0 && typej == 1){

    if (NijC <  pCHdom[0][0]) NijC=pCHdom[0][0];

     if (NijC > pCHdom[0][1]) NijC=pCHdom[0][1];

    if (NijC >= pCHdom[0][1]) NijC=pCHdom[0][1]-TOL;

    if (NijH <  pCHdom[1][0]) NijH=pCHdom[1][0];

      if (NijH > pCHdom[1][1]) NijH=pCHdom[1][1];

    if (NijH >= pCHdom[1][1]) NijH=pCHdom[1][1]-TOL;

    x = (int) floor(NijC);

    y = (int) floor(NijH);



    if (fabs(NijC-floor(NijC)) < TOL && fabs(NijH-floor(NijH)) < TOL) {

      Pij = PCHf[(int) NijC][(int) NijH];

      dN2[0] = PCHdfdx[(int) NijC][(int) NijH];

      dN2[1] = PCHdfdy[(int) NijC][(int) NijH];

      done = 1;

    }

    if (done == 0) {

      x = (int) (floor(NijC));

      y = (int) (floor(NijH));

      for (i = 0; i<16; i++) coeffs[i] = pCH[x][y][i];

      Pij = Spbicubic(NijC,NijH,coeffs,dN2);

    }

  }



  if (typei == 1 && typej == 0) {

    Pij = 0.0;

    dN2[0] = 0.0;

    dN2[1] = 0.0;

      Pij = PCHf[x][y];

      dN2[0] = PCHdfdx[x][y];

      dN2[1] = PCHdfdy[x][y];

    } else {

      Pij = Spbicubic(NijC,NijH,pCH[x][y],dN2);

    }





  if (typei == 1 && typej == 1) {

    Pij = 0.0;

    dN2[0] = 0.0;

    dN2[1] = 0.0;

  }

  return Pij;

}
@@ -3171,115 +3153,84 @@ double PairAIREBO::PijSpline(double NijC, double NijH, int typei, int typej, 
double PairAIREBO::piRCSpline(double Nij, double Nji, double Nijconj,

                              int typei, int typej, double dN3[3])

{

  int x,y,z,i,done;

  double piRC,coeffs[64];

  int x,y,z;

  double piRC;

  x=0;

  y=0;

  z=0;

  i=0;



  done=0;



  for (i=0; i<64; i++) coeffs[i]=0.0;

  dN3[0]=0.0;

  dN3[1]=0.0;

  dN3[2]=0.0;



  if (typei==0 && typej==0) {



    // CC interaction



    // if the inputs are out of bounds set them back to a point in bounds



    if (Nij <  piCCdom[0][0]) Nij=piCCdom[0][0];

    if (Nij>piCCdom[0][1]) Nij=piCCdom[0][1];

    if (Nij >= piCCdom[0][1]) Nij=piCCdom[0][1]-TOL;

    if (Nji <  piCCdom[1][0]) Nji=piCCdom[1][0];

    if (Nji>piCCdom[1][1]) Nji=piCCdom[1][1];

    if (Nji >= piCCdom[1][1]) Nji=piCCdom[1][1]-TOL;

    if (Nijconj <  piCCdom[2][0]) Nijconj=piCCdom[2][0];

    if (Nijconj>piCCdom[2][1]) Nijconj=piCCdom[2][1];



    if (fabs(Nij-floor(Nij))<TOL && fabs(Nji-floor(Nji))<TOL &&

        fabs(Nijconj-floor(Nijconj))<TOL) {

      piRC=piCCf[(int) Nij][(int) Nji][(int) Nijconj];

      dN3[0]=piCCdfdx[(int) Nij][(int) Nji][(int) Nijconj];

      dN3[1]=piCCdfdy[(int) Nij][(int) Nji][(int) Nijconj];

      dN3[2]=piCCdfdz[(int) Nij][(int) Nji][(int) Nijconj];

      done=1;

    }



    if (done==0) {

      for (i=0; i<piCCdom[0][1]; i++)

        if (Nij>=(double) i && Nij<=(double) i+1) x=i;

      for (i=0; i<piCCdom[1][1]; i++)

        if (Nji>=(double) i && Nji<=(double) i+1) y=i;

      for (i=0; i<piCCdom[2][1]; i++)

        if (Nijconj>=(double) i && Nijconj<=(double) i+1) z=i;



      for (i=0; i<64; i++) coeffs[i]=piCC[x][y][z][i];

      piRC=Sptricubic(Nij,Nji,Nijconj,coeffs,dN3);

    }

    if (Nijconj >= piCCdom[2][1]) Nijconj=piCCdom[2][1]-TOL;

    x = (int) floor(Nij);

    y = (int) floor(Nji);

    z = (int) floor(Nijconj);



    if (fabs(Nij-floor(Nij)) < TOL && fabs(Nji-floor(Nji)) < TOL

        && fabs(Nijconj-floor(Nijconj)) < TOL) {

      piRC=piCCf[x][y][z];

      dN3[0]=piCCdfdx[x][y][z];

      dN3[1]=piCCdfdy[x][y][z];

      dN3[2]=piCCdfdz[x][y][z];

    } else {

      piRC=Sptricubic(Nij,Nji,Nijconj,piCC[x][y][z],dN3);

    }

  } else if ((typei==0 && typej==1) || (typei==1 && typej==0)) {



    // CH interaction



  if ((typei==0 && typej==1) || (typei==1 && typej==0)) {



    // if the inputs are out of bounds set them back to a point in bounds



    if (Nij<piCHdom[0][0] || Nij>piCHdom[0][1] ||

        Nji<piCHdom[1][0] || Nji>piCHdom[1][1] ||

        Nijconj<piCHdom[2][0] || Nijconj>piCHdom[2][1]) {

    if (Nij <  piCHdom[0][0]) Nij=piCHdom[0][0];

      if (Nij>piCHdom[0][1]) Nij=piCHdom[0][1];

    if (Nij >= piCHdom[0][1]) Nij=piCHdom[0][1]-TOL;

    if (Nji <  piCHdom[1][0]) Nji=piCHdom[1][0];

      if (Nji>piCHdom[1][1]) Nji=piCHdom[1][1];

    if (Nji >= piCHdom[1][1]) Nji=piCHdom[1][1]-TOL;

    if (Nijconj <  piCHdom[2][0]) Nijconj=piCHdom[2][0];

      if (Nijconj>piCHdom[2][1]) Nijconj=piCHdom[2][1];

    }



    if (fabs(Nij-floor(Nij))<TOL && fabs(Nji-floor(Nji))<TOL &&

        fabs(Nijconj-floor(Nijconj))<TOL) {

      piRC=piCHf[(int) Nij][(int) Nji][(int) Nijconj];

      dN3[0]=piCHdfdx[(int) Nij][(int) Nji][(int) Nijconj];

      dN3[1]=piCHdfdy[(int) Nij][(int) Nji][(int) Nijconj];

      dN3[2]=piCHdfdz[(int) Nij][(int) Nji][(int) Nijconj];

      done=1;

    }



    if (done==0) {

      for (i=0; i<piCHdom[0][1]; i++)

        if (Nij>=i && Nij<=i+1) x=i;

      for (i=0; i<piCHdom[1][1]; i++)

        if (Nji>=i && Nji<=i+1) y=i;

      for (i=0; i<piCHdom[2][1]; i++)

        if (Nijconj>=i && Nijconj<=i+1) z=i;



      for (i=0; i<64; i++) coeffs[i]=piCH[x][y][z][i];

      piRC=Sptricubic(Nij,Nji,Nijconj,coeffs,dN3);

    }

  }

    if (Nijconj >= piCHdom[2][1]) Nijconj=piCHdom[2][1]-TOL;

    x = (int) floor(Nij);

    y = (int) floor(Nji);

    z = (int) floor(Nijconj);



    if (fabs(Nij-floor(Nij)) < TOL && fabs(Nji-floor(Nji)) < TOL

        && fabs(Nijconj-floor(Nijconj)) < TOL) {

      piRC=piCHf[x][y][z];

      dN3[0]=piCHdfdx[x][y][z];

      dN3[1]=piCHdfdy[x][y][z];

      dN3[2]=piCHdfdz[x][y][z];

    } else {

      piRC=Sptricubic(Nij,Nji,Nijconj,piCH[x][y][z],dN3);

    }

  } else if (typei==1 && typej==1) {

    if (Nij <  piHHdom[0][0]) Nij=piHHdom[0][0];

    if (Nij >= piHHdom[0][1]) Nij=piHHdom[0][1]-TOL;

    if (Nji <  piHHdom[1][0]) Nji=piHHdom[1][0];

    if (Nji >= piHHdom[1][1]) Nji=piHHdom[1][1]-TOL;

    if (Nijconj <  piHHdom[2][0]) Nijconj=piHHdom[2][0];

    if (Nijconj >= piHHdom[2][1]) Nijconj=piHHdom[2][1]-TOL;

    x = (int) floor(Nij);

    y = (int) floor(Nji);

    z = (int) floor(Nijconj);



  if (typei==1 && typej==1) {

    if (Nij<piHHdom[0][0] || Nij>piHHdom[0][1] ||

        Nji<piHHdom[1][0] || Nji>piHHdom[1][1] ||

        Nijconj<piHHdom[2][0] || Nijconj>piHHdom[2][1]) {

      Nij=0.0;

      Nji=0.0;

      Nijconj=0.0;

    }

    if (fabs(Nij-floor(Nij))<TOL && fabs(Nji-floor(Nji))<TOL &&

        fabs(Nijconj-floor(Nijconj))<TOL) {

      piRC=piHHf[(int) Nij][(int) Nji][(int) Nijconj];

      dN3[0]=piHHdfdx[(int) Nij][(int) Nji][(int) Nijconj];

      dN3[1]=piHHdfdy[(int) Nij][(int) Nji][(int) Nijconj];

      dN3[2]=piHHdfdz[(int) Nij][(int) Nji][(int) Nijconj];

      done=1;

    }

    if (done==0) {

      for (i=0; i<piHHdom[0][1]; i++)

        if (Nij>=i && Nij<=i+1) x=i;

      for (i=0; i<piHHdom[1][1]; i++)

        if (Nji>=i && Nji<=i+1) y=i;

      for (i=0; i<piHHdom[2][1]; i++)

        if (Nijconj>=i && Nijconj<=i+1) z=i;



      for (i=0; i<64; i++) coeffs[i]=piHH[x][y][z][i];

      piRC=Sptricubic(Nij,Nji,Nijconj,coeffs,dN3);

      piRC=piHHf[x][y][z];

      dN3[0]=piHHdfdx[x][y][z];

      dN3[1]=piHHdfdy[x][y][z];

      dN3[2]=piHHdfdz[x][y][z];

    } else {

      piRC=Sptricubic(Nij,Nji,Nijconj,piHH[x][y][z],dN3);

    }

  }
@@ -3293,45 +3244,37 @@ double PairAIREBO::piRCSpline(double Nij, double Nji, double Nijconj, 
double PairAIREBO::TijSpline(double Nij, double Nji,

                             double Nijconj, double dN3[3])

{

  int x,y,z,i,done;

  double Tijf,coeffs[64];

  int x,y,z;

  double Tijf;



  x=0;

  y=0;

  z=0;

  i=0;

  Tijf=0.0;

  done=0;

  for (i=0; i<64; i++) coeffs[i]=0.0;

  dN3[0]=0.0;

  dN3[1]=0.0;

  dN3[2]=0.0;



  //if the inputs are out of bounds set them back to a point in bounds



  if (Nij <  Tijdom[0][0]) Nij=Tijdom[0][0];

  if (Nij>Tijdom[0][1]) Nij=Tijdom[0][1];

  if (Nij >= Tijdom[0][1]) Nij=Tijdom[0][1]-TOL;

  if (Nji <  Tijdom[1][0]) Nji=Tijdom[1][0];

  if (Nji>Tijdom[1][1]) Nji=Tijdom[1][1];

  if (Nji >= Tijdom[1][1]) Nji=Tijdom[1][1]-TOL;

  if (Nijconj <  Tijdom[2][0]) Nijconj=Tijdom[2][0];

  if (Nijconj>Tijdom[2][1]) Nijconj=Tijdom[2][1];



  if (fabs(Nij-floor(Nij))<TOL && fabs(Nji-floor(Nji))<TOL &&

      fabs(Nijconj-floor(Nijconj))<TOL) {

    Tijf=Tf[(int) Nij][(int) Nji][(int) Nijconj];

    dN3[0]=Tdfdx[(int) Nij][(int) Nji][(int) Nijconj];

    dN3[1]=Tdfdy[(int) Nij][(int) Nji][(int) Nijconj];

    dN3[2]=Tdfdz[(int) Nij][(int) Nji][(int) Nijconj];

    done=1;

  }



  if (done==0) {

    for (i=0; i<Tijdom[0][1]; i++)

      if (Nij>=i && Nij<=i+1) x=i;

    for (i=0; i<Tijdom[1][1]; i++)

      if (Nji>=i && Nji<=i+1) y=i;

    for (i=0; i<Tijdom[2][1]; i++)

      if (Nijconj>=i && Nijconj<=i+1) z=i;



    for (i=0; i<64; i++) coeffs[i]=Tijc[x][y][z][i];

    Tijf=Sptricubic(Nij,Nji,Nijconj,coeffs,dN3);

  if (Nijconj >= Tijdom[2][1]) Nijconj=Tijdom[2][1]-TOL;

  x = (int) floor(Nij);

  y = (int) floor(Nji);

  z = (int) floor(Nijconj);



  if (fabs(Nij-floor(Nij)) < TOL && fabs(Nji-floor(Nji)) < TOL

      && fabs(Nijconj-floor(Nijconj)) < TOL) {

    Tijf=Tf[x][y][z];

    dN3[0]=Tdfdx[x][y][z];

    dN3[1]=Tdfdy[x][y][z];

    dN3[2]=Tdfdz[x][y][z];

  } else {

    Tijf=Sptricubic(Nij,Nji,Nijconj,Tijc[x][y][z],dN3);

  }



  return Tijf;

CRA Git | Maintained and supported by SUSTech CRA and CCSE