Added a Morse potential option to 'fix wall/region'

This commit is contained in:
Eugen Rožić 2019-02-21 01:49:04 +01:00
parent 52d3b9f325
commit 101948ce1e
2 changed files with 32 additions and 4 deletions

View File

@ -31,7 +31,7 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{LJ93,LJ126,LJ1043,COLLOID,HARMONIC};
enum{LJ93,LJ126,LJ1043,MORSE,COLLOID,HARMONIC};
/* ---------------------------------------------------------------------- */
@ -39,7 +39,7 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
idregion(NULL)
{
if (narg != 8) error->all(FLERR,"Illegal fix wall/region command");
if (narg < 8) error->all(FLERR,"Illegal fix wall/region command");
scalar_flag = 1;
vector_flag = 1;
@ -63,6 +63,7 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[4],"lj93") == 0) style = LJ93;
else if (strcmp(arg[4],"lj126") == 0) style = LJ126;
else if (strcmp(arg[4],"lj1043") == 0) style = LJ1043;
else if (strcmp(arg[4],"morse") == 0) style = MORSE;
else if (strcmp(arg[4],"colloid") == 0) style = COLLOID;
else if (strcmp(arg[4],"harmonic") == 0) style = HARMONIC;
else error->all(FLERR,"Illegal fix wall/region command");
@ -73,6 +74,14 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
sigma = force->numeric(FLERR,arg[6]);
cutoff = force->numeric(FLERR,arg[7]);
if (style == MORSE) {
if (narg != 9)
error->all(FLERR,"Illegal fix wall/region command");
else
alpha = force->numeric(FLERR,arg[8]);
} else if (narg != 8)
error->all(FLERR,"Illegal fix wall/region command");
if (cutoff <= 0.0) error->all(FLERR,"Fix wall/region cutoff <= 0.0");
eflag = 0;
@ -157,12 +166,15 @@ void FixWallRegion::init()
coeff5 = coeff1 * 10.0;
coeff6 = coeff2 * 4.0;
coeff7 = coeff3 * 3.0;
double rinv = 1.0/cutoff;
double r2inv = rinv*rinv;
double r4inv = r2inv*r2inv;
offset = coeff1*r4inv*r4inv*r2inv - coeff2*r4inv -
coeff3*pow(cutoff+coeff4,-3.0);
} else if (style == MORSE) {
coeff1 = 2 * epsilon * alpha;
double alpha_dr = -alpha * (cutoff - sigma);
offset = epsilon * (exp(2.0*alpha_dr) - 2.0*exp(alpha_dr));
} else if (style == COLLOID) {
coeff1 = -4.0/315.0 * epsilon * pow(sigma,6.0);
coeff2 = -2.0/3.0 * epsilon;
@ -253,6 +265,7 @@ void FixWallRegion::post_force(int vflag)
if (style == LJ93) lj93(region->contact[m].r);
else if (style == LJ126) lj126(region->contact[m].r);
else if (style == LJ1043) lj1043(region->contact[m].r);
else if (style == MORSE) morse(region->contact[m].r);
else if (style == COLLOID) colloid(region->contact[m].r,radius[i]);
else harmonic(region->contact[m].r);
@ -287,7 +300,7 @@ void FixWallRegion::post_force(int vflag)
/* ---------------------------------------------------------------------- */
void FixWallRegion::post_force_respa(int vflag, int ilevel, int /*iloop*/)
void FixWallRegion::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == ilevel_respa) post_force(vflag);
}
@ -375,6 +388,19 @@ void FixWallRegion::lj1043(double r)
coeff3*pow(r+coeff4,-3.0) - offset;
}
/* ----------------------------------------------------------------------
Morse interaction for particle with wall
compute eng and fwall = magnitude of wall force
------------------------------------------------------------------------- */
void FixWallRegion::morse(double r)
{
double dr = r - sigma;
double dexp = exp(-alpha * dr);
fwall = coeff1 * (dexp*dexp - dexp) / r;
eng = epsilon * (dexp*dexp - 2.0*dexp) - offset;
}
/* ----------------------------------------------------------------------
colloid interaction for finite-size particle of rad with wall
compute eng and fwall = magnitude of wall force

View File

@ -41,6 +41,7 @@ class FixWallRegion : public Fix {
private:
int style,iregion;
double epsilon,sigma,cutoff;
double alpha;
int eflag;
double ewall[4],ewall_all[4];
int ilevel_respa;
@ -53,6 +54,7 @@ class FixWallRegion : public Fix {
void lj93(double);
void lj126(double);
void lj1043(double);
void morse(double);
void colloid(double, double);
void harmonic(double);
};