diff --git a/src/COLLOID/fix_wall_colloid.cpp b/src/COLLOID/fix_wall_colloid.cpp index e4500ca230..764ca7fe12 100644 --- a/src/COLLOID/fix_wall_colloid.cpp +++ b/src/COLLOID/fix_wall_colloid.cpp @@ -16,14 +16,11 @@ ------------------------------------------------------------------------- */ #include "math.h" -#include "stdlib.h" #include "string.h" #include "fix_wall_colloid.h" #include "atom.h" #include "atom_vec.h" -#include "domain.h" #include "update.h" -#include "output.h" #include "respa.h" #include "error.h" @@ -31,75 +28,8 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -FixWallColloid::FixWallColloid(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg != 8) error->all("Illegal fix wall/colloid command"); - - scalar_flag = 1; - vector_flag = 1; - size_vector = 3; - scalar_vector_freq = 1; - extscalar = 1; - extvector = 1; - - if (strcmp(arg[3],"xlo") == 0) { - dim = 0; - side = -1; - } else if (strcmp(arg[3],"xhi") == 0) { - dim = 0; - side = 1; - } else if (strcmp(arg[3],"ylo") == 0) { - dim = 1; - side = -1; - } else if (strcmp(arg[3],"yhi") == 0) { - dim = 1; - side = 1; - } else if (strcmp(arg[3],"zlo") == 0) { - dim = 2; - side = -1; - } else if (strcmp(arg[3],"zhi") == 0) { - dim = 2; - side = 1; - } else error->all("Illegal fix wall/colloid command"); - - coord = atof(arg[4]); - hamaker = atof(arg[5]); - sigma = atof(arg[6]); - cutoff = atof(arg[7]); - - coeff1 = -4.0/315.0 * hamaker * pow(sigma,6.0); - coeff2 = -2.0/3.0 * hamaker; - coeff3 = hamaker * pow(sigma,6.0)/7560.0; - coeff4 = hamaker/6.0; - - double rinv = 1.0/cutoff; - double r2inv = rinv*rinv; - double r4inv = r2inv*r2inv; - offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv; - - if (dim == 0 && domain->xperiodic) - error->all("Cannot use wall in periodic dimension"); - if (dim == 1 && domain->yperiodic) - error->all("Cannot use wall in periodic dimension"); - if (dim == 2 && domain->zperiodic) - error->all("Cannot use wall in periodic dimension"); - - wall_flag = 0; - wall[0] = wall[1] = wall[2] = wall[3] = 0.0; -} - -/* ---------------------------------------------------------------------- */ - -int FixWallColloid::setmask() -{ - int mask = 0; - mask |= POST_FORCE; - mask |= THERMO_ENERGY; - mask |= POST_FORCE_RESPA; - mask |= MIN_POST_FORCE; - return mask; -} +FixWallColloid::FixWallColloid(LAMMPS *lmp, int narg, char **arg) : + FixWall(lmp, narg, arg) {} /* ---------------------------------------------------------------------- */ @@ -135,34 +65,33 @@ void FixWallColloid::init() MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall) error->all("Fix wall/colloid requires extended particles"); - if (strcmp(update->integrate_style,"respa") == 0) - nlevels_respa = ((Respa *) update->integrate)->nlevels; + FixWall::init(); } /* ---------------------------------------------------------------------- */ -void FixWallColloid::setup(int vflag) +void FixWallColloid::precompute(int m) { - if (strcmp(update->integrate_style,"verlet") == 0) - post_force(vflag); - else { - ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(vflag,nlevels_respa-1,0); - ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); - } + coeff1[m] = -4.0/315.0 * epsilon[m] * pow(sigma[m],6.0); + coeff2[m] = -2.0/3.0 * epsilon[m]; + coeff3[m] = epsilon[m] * pow(sigma[m],6.0)/7560.0; + coeff4[m] = epsilon[m]/6.0; + + double rinv = 1.0/cutoff[m]; + double r2inv = rinv*rinv; + double r4inv = r2inv*r2inv; + offset[m] = coeff3[m]*r4inv*r4inv*rinv - coeff4[m]*r2inv*rinv; } /* ---------------------------------------------------------------------- */ -void FixWallColloid::min_setup(int vflag) +void FixWallColloid::wall_particle(int m, double coord) { - post_force(vflag); -} + double delta,delta2,rinv,r2inv,r4inv,r8inv,fwall; + double r2,rinv2,r2inv2,r4inv2,r6inv2; + double r3,rinv3,r2inv3,r4inv3,r6inv3; + double rad,rad2,rad4,rad8,diam,new_coeff2; -/* ---------------------------------------------------------------------- */ - -void FixWallColloid::post_force(int vflag) -{ double **x = atom->x; double **f = atom->f; double **shape = atom->shape; @@ -170,21 +99,17 @@ void FixWallColloid::post_force(int vflag) int *type = atom->type; int nlocal = atom->nlocal; - double delta,delta2,rinv,r2inv,r4inv,r8inv,fwall; - double r2,rinv2,r2inv2,r4inv2,r6inv2; - double r3,rinv3,r2inv3,r4inv3,r6inv3; - double rad,rad2,rad4,rad8,diam,new_coeff2; - wall[0] = wall[1] = wall[2] = wall[3] = 0.0; - wall_flag = 0; + int dim = m/2; + int side = m % 2; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (side == -1) delta = x[i][dim] - coord; + if (side == 0) delta = x[i][dim] - coord; else delta = coord - x[i][dim]; if (delta <= 0.0) continue; - if (delta > cutoff) continue; + if (delta > cutoff[m]) continue; rad = shape[type[i]][0]; - new_coeff2 = coeff2*rad*rad*rad; + new_coeff2 = coeff2[m]*rad*rad*rad; diam = 2.0*rad; rad2 = rad*rad; rad4 = rad2*rad2; @@ -194,9 +119,9 @@ void FixWallColloid::post_force(int vflag) r2inv = rinv*rinv; r4inv = r2inv*r2inv; r8inv = r4inv*r4inv; - fwall = (coeff1*(rad8*rad + 27.0*rad4*rad2*rad*pow(delta,2.0) - + 63.0*rad4*rad*pow(delta,4.0) - + 21.0*rad2*rad*pow(delta,6.0))*r8inv - + fwall = (coeff1[m]*(rad8*rad + 27.0*rad4*rad2*rad*pow(delta,2.0) + + 63.0*rad4*rad*pow(delta,4.0) + + 21.0*rad2*rad*pow(delta,6.0))*r8inv - new_coeff2*r2inv) * side; f[i][dim] -= fwall; r2 = 0.5*diam - delta; @@ -204,59 +129,15 @@ void FixWallColloid::post_force(int vflag) r2inv2 = rinv2*rinv2; r4inv2 = r2inv2*r2inv2; r6inv2 = r4inv2*r2inv2; - r3 = delta+0.5*diam; + r3 = delta + 0.5*diam; rinv3 = 1.0/r3; r2inv3 = rinv3*rinv3; r4inv3 = r2inv3*r2inv3; r6inv3 = r4inv3*r2inv3; - wall[0] += coeff3*((-3.5*diam+delta)*r4inv2*r2inv2*rinv2 - + (3.5*diam+delta)*r4inv3*r2inv3*rinv3) - - coeff4*((-diam*delta+r2*r3*(log(-r2)-log(r3)))* - (-rinv2)*rinv3) - offset; - wall[dim+1] += fwall; + ewall[0] += coeff3[m]*((-3.5*diam+delta)*r4inv2*r2inv2*rinv2 + + (3.5*diam+delta)*r4inv3*r2inv3*rinv3) - + coeff4[m]*((-diam*delta+r2*r3*(log(-r2)-log(r3)))* + (-rinv2)*rinv3) - offset[m]; + ewall[m+1] += fwall; } } - -/* ---------------------------------------------------------------------- */ - -void FixWallColloid::post_force_respa(int vflag, int ilevel, int iloop) -{ - if (ilevel == nlevels_respa-1) post_force(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixWallColloid::min_post_force(int vflag) -{ - post_force(vflag); -} - -/* ---------------------------------------------------------------------- - energy of wall interaction -------------------------------------------------------------------------- */ - -double FixWallColloid::compute_scalar() -{ - // only sum across procs one time - - if (wall_flag == 0) { - MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world); - wall_flag = 1; - } - return wall_all[0]; -} - -/* ---------------------------------------------------------------------- - components of force on wall -------------------------------------------------------------------------- */ - -double FixWallColloid::compute_vector(int n) -{ - // only sum across procs one time - - if (wall_flag == 0) { - MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world); - wall_flag = 1; - } - return wall_all[n+1]; -} diff --git a/src/COLLOID/fix_wall_colloid.h b/src/COLLOID/fix_wall_colloid.h index 291338b705..215970aa29 100644 --- a/src/COLLOID/fix_wall_colloid.h +++ b/src/COLLOID/fix_wall_colloid.h @@ -14,31 +14,19 @@ #ifndef FIX_WALL_COLLOID_H #define FIX_WALL_COLLOID_H -#include "fix.h" +#include "fix_wall.h" namespace LAMMPS_NS { -class FixWallColloid : public Fix { +class FixWallColloid : public FixWall { public: FixWallColloid(class LAMMPS *, int, char **); - ~FixWallColloid() {} - int setmask(); void init(); - void setup(int); - void min_setup(int); - void post_force(int); - void post_force_respa(int, int, int); - void min_post_force(int); - double compute_scalar(); - double compute_vector(int); + void precompute(int); + void wall_particle(int, double); private: - int dim,side,thermo_flag,eflag_enable; - double coord,hamaker,sigma,cutoff; - double coeff1,coeff2,coeff3,coeff4,offset; - double wall[4],wall_all[4]; - int wall_flag; - int nlevels_respa; + double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6]; }; } diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp new file mode 100644 index 0000000000..d8cdc18162 --- /dev/null +++ b/src/fix_wall.cpp @@ -0,0 +1,222 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "fix_wall.h" +#include "atom.h" +#include "domain.h" +#include "update.h" +#include "output.h" +#include "respa.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{XLO,XHI,YLO,YHI,ZLO,ZHI}; + +/* ---------------------------------------------------------------------- */ + +FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + scalar_flag = 1; + vector_flag = 1; + size_vector = 6; + scalar_vector_freq = 1; + extscalar = 1; + extvector = 1; + + // parse args + + for (int m = 0; m < 6; m++) wflag[m] = 0; + vflag = 0; + aflag = 0; + + int iarg = 3; + if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) || + (strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) || + (strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) { + if (iarg+5 > narg) error->all("Illegal fix wall command"); + int m; + if (strcmp(arg[iarg],"xlo") == 0) m = XLO; + else if (strcmp(arg[iarg],"xhi") == 0) m = XHI; + else if (strcmp(arg[iarg],"ylo") == 0) m = YLO; + else if (strcmp(arg[iarg],"yhi") == 0) m = YHI; + else if (strcmp(arg[iarg],"zlo") == 0) m = ZLO; + else if (strcmp(arg[iarg],"zhi") == 0) m = ZHI; + wflag[m] = 1; + coord0[m] = atof(arg[iarg+1]); + epsilon[m] = atof(arg[iarg+2]); + sigma[m] = atof(arg[iarg+3]); + cutoff[m] = atof(arg[iarg+4]); + iarg += 5; + } else if (strcmp(arg[iarg],"vel") == 0) { + if (iarg+2 > narg) error->all("Illegal fix wall command"); + vflag = 1; + vel = atof(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"wiggle") == 0) { + if (iarg+3 > narg) error->all("Illegal fix wall command"); + aflag = 1; + amplitude = atof(arg[iarg+1]); + period = atof(arg[iarg+2]); + iarg += 3; + } else error->all("Illegal fix wall command"); + + // error check + + int flag = 0; + for (int m = 0; m < 6; m++) if (wflag[m]) flag; + if (!flag) error->all("Illegal fix wall command"); + + if (vflag && aflag) + error->all("Cannot set both vel and wiggle in fix wall command"); + + if ((wflag[XLO] || wflag[XHI]) && domain->xperiodic) + error->all("Cannot use fix wall in periodic dimension"); + if ((wflag[YLO] || wflag[YHI]) && domain->yperiodic) + error->all("Cannot use fix wall in periodic dimension"); + if ((wflag[ZLO] || wflag[ZHI]) && domain->xperiodic) + error->all("Cannot use fix wall in periodic zimension"); + + // setup oscillations + + if (aflag) { + double PI = 4.0 * atan(1.0); + omega = 2.0*PI / period; + } + + eflag = 0; + for (int m = 0; m < 7; m++) ewall[m] = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +int FixWall::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= THERMO_ENERGY; + mask |= POST_FORCE_RESPA; + mask |= MIN_POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixWall::init() +{ + dt = update->dt; + + // setup coefficients + + for (int m = 0; m < 6; m++) + if (wflag[m]) precompute(m); + + if (strcmp(update->integrate_style,"respa") == 0) + nlevels_respa = ((Respa *) update->integrate)->nlevels; +} + +/* ---------------------------------------------------------------------- */ + +void FixWall::setup(int vflag) +{ + if (strcmp(update->integrate_style,"verlet") == 0) + post_force(vflag); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(vflag,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixWall::min_setup(int vflag) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixWall::post_force(int vflag) +{ + eflag = 0; + for (int m = 0; m < 7; m++) ewall[m] = 0.0; + + // coord0 = initial position of wall + // coord = current position of wall + + double delta = (update->ntimestep - update->beginstep) * dt; + double coord; + + for (int m = 0; m < 6; m++) { + if (wflag[m] == 0) continue; + + if (vflag) { + if (m/2 == 0) coord = coord0[m] + delta*vel; + else coord = coord0[m] - delta*vel; + } else if (aflag) { + if (m/2 == 0) coord = coord0[m] + amplitude*sin(omega*delta); + else coord = coord0[m] - amplitude*sin(omega*delta); + } else coord = coord0[m]; + + wall_particle(m,coord); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixWall::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixWall::min_post_force(int vflag) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- + energy of wall interaction +------------------------------------------------------------------------- */ + +double FixWall::compute_scalar() +{ + // only sum across procs one time + + if (eflag == 0) { + MPI_Allreduce(ewall,ewall_all,7,MPI_DOUBLE,MPI_SUM,world); + eflag = 1; + } + return ewall_all[0]; +} + +/* ---------------------------------------------------------------------- + components of force on wall +------------------------------------------------------------------------- */ + +double FixWall::compute_vector(int n) +{ + // only sum across procs one time + + if (eflag == 0) { + MPI_Allreduce(ewall,ewall_all,7,MPI_DOUBLE,MPI_SUM,world); + eflag = 1; + } + return ewall_all[n+1]; +} diff --git a/src/fix_wall.h b/src/fix_wall.h new file mode 100644 index 0000000000..7e51ac46a4 --- /dev/null +++ b/src/fix_wall.h @@ -0,0 +1,51 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef FIX_WALL_H +#define FIX_WALL_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixWall : public Fix { + public: + FixWall(class LAMMPS *, int, char **); + virtual ~FixWall() {} + int setmask(); + virtual void init(); + void setup(int); + void min_setup(int); + void post_force(int); + void post_force_respa(int, int, int); + void min_post_force(int); + double compute_scalar(); + double compute_vector(int); + + virtual void precompute(int) = 0; + virtual void wall_particle(int, double) = 0; + + protected: + int wflag[6]; + double coord0[6],epsilon[6],sigma[6],cutoff[6]; + int vflag,aflag; + double vel,amplitude,period,omega; + int eflag; + double ewall[7],ewall_all[7]; + int nlevels_respa; + double dt; +}; + +} + +#endif diff --git a/src/fix_wall_lj126.cpp b/src/fix_wall_lj126.cpp index c857bae82c..89359dd252 100644 --- a/src/fix_wall_lj126.cpp +++ b/src/fix_wall_lj126.cpp @@ -11,212 +11,57 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- - Contributing author: Mark Stevens (SNL) -------------------------------------------------------------------------- */ - #include "math.h" -#include "stdlib.h" -#include "string.h" #include "fix_wall_lj126.h" #include "atom.h" -#include "domain.h" -#include "update.h" -#include "output.h" -#include "respa.h" -#include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -FixWallLJ126::FixWallLJ126(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) +FixWallLJ126::FixWallLJ126(LAMMPS *lmp, int narg, char **arg) : + FixWall(lmp, narg, arg) {} + +/* ---------------------------------------------------------------------- */ + +void FixWallLJ126::precompute(int m) { - if (narg < 8) error->all("Illegal fix wall/lj126 command"); + coeff1[m] = 48.0 * epsilon[m] * pow(sigma[m],12.0); + coeff2[m] = 24.0 * epsilon[m] * pow(sigma[m],6.0); + coeff3[m] = 4.0 * epsilon[m] * pow(sigma[m],12.0); + coeff4[m] = 4.0 * epsilon[m] * pow(sigma[m],6.0); - scalar_flag = 1; - vector_flag = 1; - size_vector = 3; - scalar_vector_freq = 1; - extscalar = 1; - extvector = 1; - - // set defaults - - if (strcmp(arg[3],"xlo") == 0) { - dim = 0; - side = -1; - } else if (strcmp(arg[3],"xhi") == 0) { - dim = 0; - side = 1; - } else if (strcmp(arg[3],"ylo") == 0) { - dim = 1; - side = -1; - } else if (strcmp(arg[3],"yhi") == 0) { - dim = 1; - side = 1; - } else if (strcmp(arg[3],"zlo") == 0) { - dim = 2; - side = -1; - } else if (strcmp(arg[3],"zhi") == 0) { - dim = 2; - side = 1; - } else error->all("Illegal fix wall/lj126 command"); - - coord0 = atof(arg[4]); - epsilon = atof(arg[5]); - sigma = atof(arg[6]); - cutoff = atof(arg[7]); - - // read options - - vel = 0.0; - - int iarg = 8; - while (iarg < narg) { - if (strcmp(arg[iarg],"vel") == 0) { - if (iarg+2 > narg) error->all("Illegal fix wall/lj126 command"); - vel = atof(arg[iarg+1]); - iarg += 2; - } else error->all("Illegal fix wall/lj126 command"); - } - - coeff1 = 48.0 * epsilon * pow(sigma,12.0); - coeff2 = 24.0 * epsilon * pow(sigma,6.0); - coeff3 = 4.0 * epsilon * pow(sigma,12.0); - coeff4 = 4.0 * epsilon * pow(sigma,6.0); - - double r2inv = 1.0/(cutoff*cutoff); + double r2inv = 1.0/(cutoff[m]*cutoff[m]); double r6inv = r2inv*r2inv*r2inv; - offset = r6inv*(coeff3*r6inv - coeff4); - - if (dim == 0 && domain->xperiodic) - error->all("Cannot use wall in periodic dimension"); - if (dim == 1 && domain->yperiodic) - error->all("Cannot use wall in periodic dimension"); - if (dim == 2 && domain->zperiodic) - error->all("Cannot use wall in periodic dimension"); - - wall_flag = 0; - wall[0] = wall[1] = wall[2] = wall[3] = 0.0; + offset[m] = r6inv*(coeff3[m]*r6inv - coeff4[m]); } /* ---------------------------------------------------------------------- */ -int FixWallLJ126::setmask() +void FixWallLJ126::wall_particle(int m, double coord) { - int mask = 0; - mask |= POST_FORCE; - mask |= THERMO_ENERGY; - mask |= POST_FORCE_RESPA; - mask |= MIN_POST_FORCE; - return mask; -} + double delta,rinv,r2inv,r6inv,fwall; -/* ---------------------------------------------------------------------- */ - -void FixWallLJ126::init() -{ - if (strcmp(update->integrate_style,"respa") == 0) - nlevels_respa = ((Respa *) update->integrate)->nlevels; -} - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ126::setup(int vflag) -{ - if (strcmp(update->integrate_style,"verlet") == 0) - post_force(vflag); - else { - ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(vflag,nlevels_respa-1,0); - ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ126::min_setup(int vflag) -{ - post_force(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ126::post_force(int vflag) -{ double **x = atom->x; double **f = atom->f; int *mask = atom->mask; int nlocal = atom->nlocal; - double delta,rinv,r2inv,r6inv,fwall; - wall[0] = wall[1] = wall[2] = wall[3] = 0.0; - wall_flag = 0; - - // coord = current position of wall - // coord0 = initial position of wall - - double delt = (update->ntimestep - update->beginstep) * update->dt; - double coord = coord0 + delt*vel; + int dim = m/2; + int side = m % 2; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (side == -1) delta = x[i][dim] - coord; + if (side == 0) delta = x[i][dim] - coord; else delta = coord - x[i][dim]; if (delta <= 0.0) continue; - if (delta > cutoff) continue; + if (delta > cutoff[m]) continue; rinv = 1.0/delta; r2inv = rinv*rinv; r6inv = r2inv*r2inv*r2inv; - fwall = side * r6inv*(coeff1*r6inv - coeff2) * rinv; + fwall = side * r6inv*(coeff1[m]*r6inv - coeff2[m]) * rinv; f[i][dim] -= fwall; - wall[0] += r6inv*(coeff3*r6inv - coeff4) - offset; - wall[dim+1] += fwall; + ewall[0] += r6inv*(coeff3[m]*r6inv - coeff4[m]) - offset[m]; + ewall[m+1] += fwall; } } - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ126::post_force_respa(int vflag, int ilevel, int iloop) -{ - if (ilevel == nlevels_respa-1) post_force(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ126::min_post_force(int vflag) -{ - post_force(vflag); -} - -/* ---------------------------------------------------------------------- - energy of wall interaction -------------------------------------------------------------------------- */ - -double FixWallLJ126::compute_scalar() -{ - // only sum across procs one time - - if (wall_flag == 0) { - MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world); - wall_flag = 1; - } - return wall_all[0]; -} - -/* ---------------------------------------------------------------------- - components of force on wall -------------------------------------------------------------------------- */ - -double FixWallLJ126::compute_vector(int n) -{ - // only sum across procs one time - - if (wall_flag == 0) { - MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world); - wall_flag = 1; - } - return wall_all[n+1]; -} diff --git a/src/fix_wall_lj126.h b/src/fix_wall_lj126.h index b31a8089f4..2c2ddebacb 100644 --- a/src/fix_wall_lj126.h +++ b/src/fix_wall_lj126.h @@ -14,34 +14,19 @@ #ifndef FIX_WALL_LJ126_H #define FIX_WALL_LJ126_H -#include "fix.h" +#include "fix_wall.h" namespace LAMMPS_NS { -class FixWallLJ126 : public Fix { +class FixWallLJ126 : public FixWall { public: FixWallLJ126(class LAMMPS *, int, char **); - int setmask(); - void init(); - void setup(int); - void min_setup(int); - void post_force(int); - void post_force_respa(int, int, int); - void min_post_force(int); - double compute_scalar(); - double compute_vector(int); + void precompute(int); + void wall_particle(int, double); private: - int dim,side,thermo_flag,eflag_enable; - double coord,epsilon,sigma,cutoff; - double coeff1,coeff2,coeff3,coeff4,offset; - double wall[4],wall_all[4]; - double vel,coord0; - int wall_flag; - int nlevels_respa; + double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6]; }; #endif } - - diff --git a/src/fix_wall_lj93.cpp b/src/fix_wall_lj93.cpp index 833785f958..496719eba1 100644 --- a/src/fix_wall_lj93.cpp +++ b/src/fix_wall_lj93.cpp @@ -12,209 +12,59 @@ ------------------------------------------------------------------------- */ #include "math.h" -#include "stdlib.h" -#include "string.h" #include "fix_wall_lj93.h" #include "atom.h" -#include "domain.h" -#include "update.h" -#include "output.h" -#include "respa.h" -#include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -FixWallLJ93::FixWallLJ93(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) +FixWallLJ93::FixWallLJ93(LAMMPS *lmp, int narg, char **arg) : + FixWall(lmp, narg, arg) {} + +/* ---------------------------------------------------------------------- */ + +void FixWallLJ93::precompute(int m) { - if (narg < 8) error->all("Illegal fix wall/lj93 command"); + coeff1[m] = 6.0/5.0 * epsilon[m] * pow(sigma[m],9.0); + coeff2[m] = 3.0 * epsilon[m] * pow(sigma[m],3.0); + coeff3[m] = 2.0/15.0 * epsilon[m] * pow(sigma[m],9.0); + coeff4[m] = epsilon[m] * pow(sigma[m],3.0); - scalar_flag = 1; - vector_flag = 1; - size_vector = 3; - scalar_vector_freq = 1; - extscalar = 1; - extvector = 1; - - // set defaults - - if (strcmp(arg[3],"xlo") == 0) { - dim = 0; - side = -1; - } else if (strcmp(arg[3],"xhi") == 0) { - dim = 0; - side = 1; - } else if (strcmp(arg[3],"ylo") == 0) { - dim = 1; - side = -1; - } else if (strcmp(arg[3],"yhi") == 0) { - dim = 1; - side = 1; - } else if (strcmp(arg[3],"zlo") == 0) { - dim = 2; - side = -1; - } else if (strcmp(arg[3],"zhi") == 0) { - dim = 2; - side = 1; - } else error->all("Illegal fix wall/lj93 command"); - - coord0 = atof(arg[4]); - epsilon = atof(arg[5]); - sigma = atof(arg[6]); - cutoff = atof(arg[7]); - - // read options - - vel = 0.0; - - int iarg = 8; - while (iarg < narg) { - if (strcmp(arg[iarg],"vel") == 0) { - if (iarg+2 > narg) error->all("Illegal fix wall/lj93 command"); - vel = atof(arg[iarg+1]); - iarg += 2; - } else error->all("Illegal fix wall/lj93 command"); - } - - coeff1 = 6.0/5.0 * epsilon * pow(sigma,9.0); - coeff2 = 3.0 * epsilon * pow(sigma,3.0); - coeff3 = 2.0/15.0 * epsilon * pow(sigma,9.0); - coeff4 = epsilon * pow(sigma,3.0); - - double rinv = 1.0/cutoff; + double rinv = 1.0/cutoff[m]; double r2inv = rinv*rinv; double r4inv = r2inv*r2inv; - offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv; - - if (dim == 0 && domain->xperiodic) - error->all("Cannot use wall in periodic dimension"); - if (dim == 1 && domain->yperiodic) - error->all("Cannot use wall in periodic dimension"); - if (dim == 2 && domain->zperiodic) - error->all("Cannot use wall in periodic dimension"); - - wall_flag = 0; - wall[0] = wall[1] = wall[2] = wall[3] = 0.0; + offset[m] = coeff3[6]*r4inv*r4inv*rinv - coeff4[6]*r2inv*rinv; } /* ---------------------------------------------------------------------- */ -int FixWallLJ93::setmask() +void FixWallLJ93::wall_particle(int m, double coord) { - int mask = 0; - mask |= POST_FORCE; - mask |= THERMO_ENERGY; - mask |= POST_FORCE_RESPA; - mask |= MIN_POST_FORCE; - return mask; -} + double delta,rinv,r2inv,r4inv,r10inv,fwall; -/* ---------------------------------------------------------------------- */ - -void FixWallLJ93::init() -{ - if (strcmp(update->integrate_style,"respa") == 0) - nlevels_respa = ((Respa *) update->integrate)->nlevels; -} - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ93::setup(int vflag) -{ - if (strcmp(update->integrate_style,"verlet") == 0) - post_force(vflag); - else { - ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(vflag,nlevels_respa-1,0); - ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ93::min_setup(int vflag) -{ - post_force(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ93::post_force(int vflag) -{ double **x = atom->x; double **f = atom->f; int *mask = atom->mask; int nlocal = atom->nlocal; - double delta,rinv,r2inv,r4inv,r10inv,fwall; - wall[0] = wall[1] = wall[2] = wall[3] = 0.0; - wall_flag = 0; - - // coord = current position of wall - // coord0 = initial position of wall - - double delt = (update->ntimestep - update->beginstep) * update->dt; - double coord = coord0 + delt*vel; + int dim = m/2; + int side = m % 2; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (side == -1) delta = x[i][dim] - coord; + if (side == 0) delta = x[i][dim] - coord; else delta = coord - x[i][dim]; if (delta <= 0.0) continue; - if (delta > cutoff) continue; + if (delta > cutoff[m]) continue; rinv = 1.0/delta; r2inv = rinv*rinv; r4inv = r2inv*r2inv; r10inv = r4inv*r4inv*r2inv; - fwall = (coeff1*r10inv - coeff2*r4inv) * side; + fwall = (coeff1[m]*r10inv - coeff2[m]*r4inv) * side; f[i][dim] -= fwall; - wall[0] += coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv - offset; - wall[dim+1] += fwall; + ewall[0] += coeff3[m]*r4inv*r4inv*rinv - + coeff4[m]*r2inv*rinv - offset[m]; + ewall[m+1] += fwall; } } - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ93::post_force_respa(int vflag, int ilevel, int iloop) -{ - if (ilevel == nlevels_respa-1) post_force(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixWallLJ93::min_post_force(int vflag) -{ - post_force(vflag); -} - -/* ---------------------------------------------------------------------- - energy of wall interaction -------------------------------------------------------------------------- */ - -double FixWallLJ93::compute_scalar() -{ - // only sum across procs one time - - if (wall_flag == 0) { - MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world); - wall_flag = 1; - } - return wall_all[0]; -} - -/* ---------------------------------------------------------------------- - components of force on wall -------------------------------------------------------------------------- */ - -double FixWallLJ93::compute_vector(int n) -{ - // only sum across procs one time - - if (wall_flag == 0) { - MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world); - wall_flag = 1; - } - return wall_all[n+1]; -} diff --git a/src/fix_wall_lj93.h b/src/fix_wall_lj93.h index 1836a7f1ad..9020731b9b 100644 --- a/src/fix_wall_lj93.h +++ b/src/fix_wall_lj93.h @@ -14,32 +14,18 @@ #ifndef FIX_WALL_LJ93_H #define FIX_WALL_LJ93_H -#include "fix.h" +#include "fix_wall.h" namespace LAMMPS_NS { -class FixWallLJ93 : public Fix { +class FixWallLJ93 : public FixWall { public: FixWallLJ93(class LAMMPS *, int, char **); - ~FixWallLJ93() {} - int setmask(); - void init(); - void setup(int); - void min_setup(int); - void post_force(int); - void post_force_respa(int, int, int); - void min_post_force(int); - double compute_scalar(); - double compute_vector(int); + void precompute(int); + void wall_particle(int, double); private: - int dim,side,thermo_flag,eflag_enable; - double coord,epsilon,sigma,cutoff; - double coeff1,coeff2,coeff3,coeff4,offset; - double wall[4],wall_all[4]; - double vel,coord0; - int wall_flag; - int nlevels_respa; + double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6]; }; }