Commit 925481c3 authored by Tim Mattox's avatar Tim Mattox
Browse files

USER-DPD: Fix hard-wall force interaction bug, and ensure fraction is >= 0 

pair_exp6_rx.cpp patch by Jim Larentzos.  Applied by Tim Mattox.
parent f509f133

src/USER-DPD/pair_exp6_rx.cpp

+36 −10
Original line number Diff line number Diff line
@@ -310,13 +310,13 @@ void PairExp6rx::compute(int eflag, int vflag) 


            uin1 = buck1*(6.0*rin1exp - alphaOld12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);



            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;

            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;



            aRep = -1.0*win1*powint(rin1,nRep)/nRep;

            aRep = win1*powint(rin1,nRep)/nRep;



            uin1rep = aRep/powint(rin1,nRep);



            forceExp6 = -double(nRep)*aRep/powint(r,nRep);

            forceExp6 = double(nRep)*aRep/powint(r,nRep);

            fpairOldEXP6_12 = factor_lj*forceExp6*r2inv;



            evdwlOldEXP6_12 = uin1 - uin1rep + aRep/powint(r,nRep);
@@ -350,13 +350,13 @@ void PairExp6rx::compute(int eflag, int vflag) 


            uin1 = buck1*(6.0*rin1exp - alphaOld21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);



            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;

            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;



            aRep = -1.0*win1*powint(rin1,nRep)/nRep;

            aRep = win1*powint(rin1,nRep)/nRep;



            uin1rep = aRep/powint(rin1,nRep);



            forceExp6 = -double(nRep)*aRep/powint(r,nRep);

            forceExp6 = double(nRep)*aRep/powint(r,nRep);

            fpairOldEXP6_21 = factor_lj*forceExp6*r2inv;



            evdwlOldEXP6_21 = uin1 - uin1rep + aRep/powint(r,nRep);
@@ -405,9 +405,9 @@ void PairExp6rx::compute(int eflag, int vflag) 


            uin1 = buck1*(6.0*rin1exp - alpha12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);



            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;

            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;



            aRep = -1.0*win1*powint(rin1,nRep)/nRep;

            aRep = win1*powint(rin1,nRep)/nRep;



            uin1rep = aRep/powint(rin1,nRep);
@@ -437,9 +437,9 @@ void PairExp6rx::compute(int eflag, int vflag) 


            uin1 = buck1*(6.0*rin1exp - alpha21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);



            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;

            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;



            aRep = -1.0*win1*powint(rin1,nRep)/nRep;

            aRep = win1*powint(rin1,nRep)/nRep;



            uin1rep = aRep/powint(rin1,nRep);
@@ -1102,6 +1102,32 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm 
      rm2_old *= pow(nTotalOFA_old,fuchslinR);

    }

  }



  // Check that no fractions are less than zero

  if(fraction1 < 0.0){

    if(fraction1 < -1.0e-10){

      error->all(FLERR,"Computed fraction less than -1.0e-10");

    }

    fraction1 = 0.0;

  }

  if(fraction2 < 0.0){

    if(fraction2 < -1.0e-10){

      error->all(FLERR,"Computed fraction less than -1.0e-10");

    }

    fraction2 = 0.0;

  }

  if(fraction1_old < 0.0){

    if(fraction1_old < -1.0e-10){

      error->all(FLERR,"Computed fraction less than -1.0e-10");

    }

    fraction1_old = 0.0;

  }

  if(fraction2_old < 0.0){

    if(fraction2_old < -1.0e-10){

      error->all(FLERR,"Computed fraction less than -1.0e-10");

    }

    fraction2_old = 0.0;

  }

}



/* ---------------------------------------------------------------------- */

CRA Git | Maintained and supported by SUSTech CRA and CCSE