From 964577314cc1a2bc99b9001273af5f6ba1b9428b Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 10 Nov 2006 23:52:37 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@144 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/fix_wall_lj126.cpp | 179 +++++++++++++++++++++++++++++++++++++++++ src/fix_wall_lj126.h | 43 ++++++++++ src/fix_wall_lj93.cpp | 8 +- src/fix_wall_lj93.h | 2 +- src/style.h | 2 + 5 files changed, 232 insertions(+), 2 deletions(-) create mode 100644 src/fix_wall_lj126.cpp create mode 100644 src/fix_wall_lj126.h diff --git a/src/fix_wall_lj126.cpp b/src/fix_wall_lj126.cpp new file mode 100644 index 0000000000..2b7f906859 --- /dev/null +++ b/src/fix_wall_lj126.cpp @@ -0,0 +1,179 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mark Stevens (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "fix_wall_lj126.h" +#include "atom.h" +#include "update.h" +#include "output.h" +#include "respa.h" +#include "error.h" + +/* ---------------------------------------------------------------------- */ + +FixWallLJ126::FixWallLJ126(int narg, char **arg) : Fix(narg, arg) +{ + if (narg != 8) error->all("Illegal fix wall/lj126 command"); + + 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"); + + coord = atof(arg[4]); + epsilon = atof(arg[5]); + sigma = atof(arg[6]); + cutoff = atof(arg[7]); + + 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 r6inv = r2inv*r2inv*r2inv; + offset = r6inv*(coeff3*r6inv - coeff4); +} + +/* ---------------------------------------------------------------------- */ + +int FixWallLJ126::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= THERMO; + mask |= POST_FORCE_RESPA; + mask |= MIN_POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixWallLJ126::init() +{ + if (strcmp(update->integrate_style,"respa") == 0) + nlevels_respa = ((Respa *) update->integrate)->nlevels; + + if (thermo_print || thermo_energy) thermo_flag = 1; + else thermo_flag = 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixWallLJ126::setup() +{ + eflag_on = 1; + if (strcmp(update->integrate_style,"verlet") == 0) + post_force(1); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(1,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } + eflag_on = 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixWallLJ126::min_setup() +{ + eflag_on = 1; + post_force(1); +} + +/* ---------------------------------------------------------------------- */ + +void FixWallLJ126::post_force(int vflag) +{ + bool eflag = false; + if (thermo_flag) { + if (eflag_on) eflag = true; + else if (output->next_thermo == update->ntimestep) eflag = true; + } + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double delta,rinv,r2inv,r6inv; + if (eflag) eng = 0.0; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (side == -1) delta = x[i][dim] - coord; + else delta = coord - x[i][dim]; + if (delta <= 0.0) continue; + if (delta > cutoff) continue; + rinv = 1.0/delta; + r2inv = rinv*rinv; + r6inv = r2inv*r2inv*r2inv; + f[i][dim] -= r6inv*(coeff1*r6inv - coeff2) * side; + if (eflag) eng += r6inv*(coeff3*r6inv - coeff4) - offset; + } + + if (eflag) MPI_Allreduce(&eng,&etotal,1,MPI_DOUBLE,MPI_SUM,world); +} + +/* ---------------------------------------------------------------------- */ + +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); +} + +/* ---------------------------------------------------------------------- */ + +int FixWallLJ126::thermo_fields(int n, int *flags, char **keywords) +{ + if (n == 0) return 1; + flags[0] = 3; + strcpy(keywords[0],"Wall"); + return 1; +} + +/* ---------------------------------------------------------------------- */ + +int FixWallLJ126::thermo_compute(double *values) +{ + values[0] = etotal; + return 1; +} diff --git a/src/fix_wall_lj126.h b/src/fix_wall_lj126.h new file mode 100644 index 0000000000..490fd0a966 --- /dev/null +++ b/src/fix_wall_lj126.h @@ -0,0 +1,43 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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_LJ126_H +#define FIX_WALL_LJ126_H + +#include "fix.h" + +class FixWallLJ126 : public Fix { + public: + FixWallLJ126(int, char **); + ~FixWallLJ126() {} + int setmask(); + void init(); + void setup(); + void min_setup(); + void post_force(int); + void post_force_respa(int, int, int); + void min_post_force(int); + + int thermo_fields(int, int *, char **); + int thermo_compute(double *); + + private: + int dim,side,thermo_flag,eflag_on; + double coord,epsilon,sigma,cutoff; + double coeff1,coeff2,coeff3,coeff4,offset; + double eng,etotal; + int nlevels_respa; +}; + +#endif + diff --git a/src/fix_wall_lj93.cpp b/src/fix_wall_lj93.cpp index 24b4b8a474..7677403c0e 100644 --- a/src/fix_wall_lj93.cpp +++ b/src/fix_wall_lj93.cpp @@ -56,6 +56,12 @@ FixWallLJ93::FixWallLJ93(int narg, char **arg) : Fix(narg, arg) 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 r2inv = rinv*rinv; + double r4inv = r2inv*r2inv; + double r10inv = r4inv*r4inv*r2inv; + offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv; } /* ---------------------------------------------------------------------- */ @@ -133,7 +139,7 @@ void FixWallLJ93::post_force(int vflag) r4inv = r2inv*r2inv; r10inv = r4inv*r4inv*r2inv; f[i][dim] -= (coeff1*r10inv - coeff2*r4inv) * side; - if (eflag) eng += coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv; + if (eflag) eng += coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv - offset; } if (eflag) MPI_Allreduce(&eng,&etotal,1,MPI_DOUBLE,MPI_SUM,world); diff --git a/src/fix_wall_lj93.h b/src/fix_wall_lj93.h index c8cc09af00..2febca5d0c 100644 --- a/src/fix_wall_lj93.h +++ b/src/fix_wall_lj93.h @@ -33,7 +33,7 @@ class FixWallLJ93 : public Fix { private: int dim,side,thermo_flag,eflag_on; double coord,epsilon,sigma,cutoff; - double coeff1,coeff2,coeff3,coeff4; + double coeff1,coeff2,coeff3,coeff4,offset; double eng,etotal; int nlevels_respa; }; diff --git a/src/style.h b/src/style.h index 5a479f89ff..d75feff684 100644 --- a/src/style.h +++ b/src/style.h @@ -130,6 +130,7 @@ DumpStyle(xyz,DumpXYZ) #include "fix_uniaxial.h" #include "fix_viscous.h" #include "fix_volume_rescale.h" +#include "fix_wall_lj126.h" #include "fix_wall_lj93.h" #include "fix_wall_reflect.h" #include "fix_wiggle.h" @@ -174,6 +175,7 @@ FixStyle(tmd,FixTMD) FixStyle(uniaxial,FixUniaxial) FixStyle(viscous,FixViscous) FixStyle(volume/rescale,FixVolRescale) +FixStyle(wall/lj126,FixWallLJ126) FixStyle(wall/lj93,FixWallLJ93) FixStyle(wall/reflect,FixWallReflect) FixStyle(wiggle,FixWiggle)