diff --git a/doc/src/fix_cmap.txt b/doc/src/fix_cmap.txt index 2b14a20c1d..b126382b80 100644 --- a/doc/src/fix_cmap.txt +++ b/doc/src/fix_cmap.txt @@ -97,6 +97,11 @@ The "fix_modify"_fix_modify.html {energy} option is supported by this fix to add the potential "energy" of the CMAP interactions system's potential energy as part of "thermodynamic output"_thermo_style.html. +The "fix_modify"_fix_modify.html {virial} option is supported by this +fix to add the contribution due to the interaction between +atoms and each wall to the system's virial as part of "thermodynamic +output"_thermo_style.html. The default is {virial yes} + This fix computes a global scalar which can be accessed by various "output commands"_Section_howto.html#howto_15. The scalar is the potential energy discussed above. The scalar value calculated by this diff --git a/doc/src/fix_efield.txt b/doc/src/fix_efield.txt index 9e3d7797d8..6537898315 100644 --- a/doc/src/fix_efield.txt +++ b/doc/src/fix_efield.txt @@ -124,6 +124,11 @@ can include the forces added by this fix in a consistent manner. I.e. there is a decrease in potential energy when atoms move in the direction of the added force due to the electric field. +The "fix_modify"_fix_modify.html {virial} option is supported by this +fix to add the contribution due to the interaction between +atoms and each wall to the system's virial as part of "thermodynamic +output"_thermo_style.html. The default is {virial no} + The "fix_modify"_fix_modify.html {respa} option is supported by this fix. This allows to set at which level of the "r-RESPA"_run_style.html integrator the fix adding its forces. Default is the outermost level. diff --git a/doc/src/fix_external.txt b/doc/src/fix_external.txt index 25158be0db..c77a54d800 100644 --- a/doc/src/fix_external.txt +++ b/doc/src/fix_external.txt @@ -131,6 +131,11 @@ forces added by this fix in a consistent manner. I.e. there is a decrease in potential energy when atoms move in the direction of the added force. +The "fix_modify"_fix_modify.html {virial} option is supported by this +fix to add the contribution due to the interaction between +atoms and each wall to the system's virial as part of "thermodynamic +output"_thermo_style.html. The default is {virial yes} + This fix computes a global scalar which can be accessed by various "output commands"_Section_howto.html#howto_15. The scalar is the potential energy discussed above. The scalar stored by this fix diff --git a/doc/src/fix_modify.txt b/doc/src/fix_modify.txt index 4e31092888..1c7dcf77a5 100644 --- a/doc/src/fix_modify.txt +++ b/doc/src/fix_modify.txt @@ -14,10 +14,11 @@ fix_modify fix-ID keyword value ... :pre fix-ID = ID of the fix to modify :ulb,l one or more keyword/value pairs may be appended :l -keyword = {temp} or {press} or {energy} or {respa} or {dynamic/dof} :l +keyword = {temp} or {press} or {energy} or {virial} or {respa} or {dynamic/dof} :l {temp} value = compute ID that calculates a temperature {press} value = compute ID that calculates a pressure {energy} value = {yes} or {no} + {virial} value = {yes} or {no} {respa} value = {1} to {max respa level} or {0} (for outermost level) {dynamic/dof} value = {yes} or {no} yes/no = do or do not recompute the number of degrees of freedom (DOF) contributing to the temperature :pre @@ -52,11 +53,10 @@ define their own compute by default, as described in their documentation. Thus this option allows the user to override the default method for computing P. -For fixes that calculate a contribution to the potential energy of the -system, the {energy} keyword will include that contribution in -thermodynamic output of potential energy. This is because the {energy -yes} setting must be specified to include the fix's global or per-atom -energy in the calculation performed by the "compute +The {energy} keyword can be used with fixes that support it. +{energy yes} adds a contribution to the potential energy of the +system. The fix's global and per-atom +energy is included in the calculation performed by the "compute pe"_compute_pe.html or "compute pe/atom"_compute_pe_atom.html commands. See the "thermo_style"_thermo_style.html command for info on how potential energy is output. For fixes that tally a global @@ -69,6 +69,25 @@ are using it when performing an "energy minimization"_minimize.html and if you want the energy and forces it produces to be part of the optimization criteria. +The {virial} keyword can be used with fixes that support it. +{virial yes} adds a contribution to the virial of the +system. The fix's global and per-atom +virial is included in the calculation performed by the "compute +pressure"_compute_pressure.html or +"compute stress/atom"_compute_stress_atom.html +commands. See the "thermo_style"_thermo_style.html command for info +on how pressure is output. + +NOTE: You must specify the {virial yes} setting for a fix if you +are doing "box relaxation"_fix_box_relax.html and +if you want virial contribution of the fix to be part of the +relaxation criteria, although this seems unlikely. + +NOTE: This option is only supported by fixes that explicitly say +so. For some of these (e.g. the +"fix shake"_fix_shake.html command) the default setting is +{virial yes}, for others it is {virial no}. + For fixes that set or modify forces, it may be possible to select at which "r-RESPA"_run_style.html level the fix operates via the {respa} keyword. The RESPA level at which the fix is active can be selected. @@ -111,4 +130,4 @@ pressure"_compute_pressure.html, "thermo_style"_thermo_style.html [Default:] The option defaults are temp = ID defined by fix, press = ID defined -by fix, energy = no, respa = 0. +by fix, energy = no, virial = different for each fix style, respa = 0. diff --git a/doc/src/fix_rigid.txt b/doc/src/fix_rigid.txt index 62969112f7..af00882c5b 100644 --- a/doc/src/fix_rigid.txt +++ b/doc/src/fix_rigid.txt @@ -703,6 +703,11 @@ NVT, NPT, NPH rigid styles to add the energy change induced by the thermostatting to the system's potential energy as part of "thermodynamic output"_thermo_style.html. +The "fix_modify"_fix_modify.html {virial} option is supported by this +fix to add the contribution due to the interaction between +atoms and each wall to the system's virial as part of "thermodynamic +output"_thermo_style.html. The default is {virial yes} + The "fix_modify"_fix_modify.html {temp} and {press} options are supported by the 4 NPT and NPH rigid styles to change the computes used to calculate the instantaneous pressure tensor. Note that the 2 diff --git a/doc/src/fix_shake.txt b/doc/src/fix_shake.txt index c187b17c6c..c39ccdd07e 100644 --- a/doc/src/fix_shake.txt +++ b/doc/src/fix_shake.txt @@ -186,6 +186,11 @@ to 1 and recompiling LAMMPS. [Restart, fix_modify, output, run start/stop, minimize info:] +The "fix_modify"_fix_modify.html {virial} option is supported by this +fix to add the contribution due to the interaction between +atoms and each wall to the system's virial as part of "thermodynamic +output"_thermo_style.html. The default is {virial yes} + No information about these fixes is written to "binary restart files"_restart.html. None of the "fix_modify"_fix_modify.html options are relevant to these fixes. No global or per-atom quantities are diff --git a/doc/src/fix_wall.txt b/doc/src/fix_wall.txt index 6bbfccf9db..e814c89a07 100644 --- a/doc/src/fix_wall.txt +++ b/doc/src/fix_wall.txt @@ -252,6 +252,11 @@ fix to add the energy of interaction between atoms and each wall to the system's potential energy as part of "thermodynamic output"_thermo_style.html. +The "fix_modify"_fix_modify.html {virial} option is supported by this +fix to add the contribution due to the interaction between +atoms and each wall to the system's virial as part of "thermodynamic +output"_thermo_style.html. The default is {virial no} + The "fix_modify"_fix_modify.html {respa} option is supported by this fix. This allows to set at which level of the "r-RESPA"_run_style.html integrator the fix is adding its forces. Default is the outermost level. diff --git a/doc/src/fix_wall_region.txt b/doc/src/fix_wall_region.txt index ca5335e3fb..9700545dc9 100644 --- a/doc/src/fix_wall_region.txt +++ b/doc/src/fix_wall_region.txt @@ -15,7 +15,7 @@ fix ID group-ID wall/region region-ID style epsilon sigma cutoff :pre ID, group-ID are documented in "fix"_fix.html command wall/region = style name of this fix command region-ID = region whose boundary will act as wall -style = {lj93} or {lj126} or {colloid} or {harmonic} +style = {lj93} or {lj126} or {lj1043} or {colloid} or {harmonic} epsilon = strength factor for wall-particle interaction (energy or energy/distance^2 units) sigma = size factor for wall-particle interaction (distance units) cutoff = distance from wall at which wall-particle interaction is cut off (distance units) :ul @@ -112,6 +112,10 @@ For style {lj126}, the energy E is given by the 12/6 potential: :c,image(Eqs/pair_lj.jpg) +For style {wall/lj1043}, the energy E is given by the 10/4/3 potential: + +:c,image(Eqs/fix_wall_lj1043.jpg) + For style {colloid}, the energy E is given by an integrated form of the "pair_style colloid"_pair_colloid.html potential: @@ -128,49 +132,8 @@ surface no longer interact. The energy of the wall potential is shifted so that the wall-particle interaction energy is 0.0 at the cutoff distance. -For the {lj93} and {lj126} styles, {epsilon} and {sigma} are the usual -Lennard-Jones parameters, which determine the strength and size of the -particle as it interacts with the wall. Epsilon has energy units. -Note that this {epsilon} and {sigma} may be different than any -{epsilon} or {sigma} values defined for a pair style that computes -particle-particle interactions. - -The {lj93} interaction is derived by integrating over a 3d -half-lattice of Lennard-Jones 12/6 particles. The {lj126} interaction -is effectively a harder, more repulsive wall interaction. - -For the {colloid} style, {epsilon} is effectively a Hamaker constant -with energy units for the colloid-wall interaction, {R} is the radius -of the colloid particle, {D} is the distance from the surface of the -colloid particle to the wall (r-R), and {sigma} is the size of a -constituent LJ particle inside the colloid particle. Note that the -cutoff distance Rc in this case is the distance from the colloid -particle center to the wall. - -The {colloid} interaction is derived by integrating over constituent -LJ particles of size {sigma} within the colloid particle and a 3d -half-lattice of Lennard-Jones 12/6 particles of size {sigma} in the -wall. - -For the {wall/harmonic} style, {epsilon} is effectively the spring -constant K, and has units (energy/distance^2). The input parameter -{sigma} is ignored. The minimum energy position of the harmonic -spring is at the {cutoff}. This is a repulsive-only spring since the -interaction is truncated at the {cutoff} - -NOTE: For all of the styles, you must insure that r is always > 0 for -all particles in the group, or LAMMPS will generate an error. This -means you cannot start your simulation with particles on the region -surface (r = 0) or with particles on the wrong side of the region -surface (r < 0). For the {wall/lj93} and {wall/lj126} styles, the -energy of the wall/particle interaction (and hence the force on the -particle) blows up as r -> 0. The {wall/colloid} style is even more -restrictive, since the energy blows up as D = r-R -> 0. This means -the finite-size particles of radius R must be a distance larger than R -from the region surface. The {harmonic} style is a softer potential -and does not blow up as r -> 0, but you must use a large enough -{epsilon} that particles always reamin on the correct side of the -region surface (r > 0). +For a full description of these wall styles, see fix_style +"wall"_fix_wall.html [Restart, fix_modify, output, run start/stop, minimize info:] @@ -182,6 +145,11 @@ fix to add the energy of interaction between atoms and the wall to the system's potential energy as part of "thermodynamic output"_thermo_style.html. +The "fix_modify"_fix_modify.html {virial} option is supported by this +fix to add the contribution due to the interaction between +atoms and each wall to the system's virial as part of "thermodynamic +output"_thermo_style.html. The default is {virial no} + The "fix_modify"_fix_modify.html {respa} option is supported by this fix. This allows to set at which level of the "r-RESPA"_run_style.html integrator the fix is adding its forces. Default is the outermost level. diff --git a/src/COLLOID/fix_wall_colloid.cpp b/src/COLLOID/fix_wall_colloid.cpp index 7c0203ac5f..651e0da4f8 100644 --- a/src/COLLOID/fix_wall_colloid.cpp +++ b/src/COLLOID/fix_wall_colloid.cpp @@ -81,6 +81,7 @@ void FixWallColloid::wall_particle(int m, int which, double coord) double r3,rinv3,r2inv3,r4inv3; double rad,rad2,rad4,rad8,diam,new_coeff2; double eoffset; + double vn; double **x = atom->x; double **f = atom->f; @@ -151,6 +152,12 @@ void FixWallColloid::wall_particle(int m, int which, double coord) ewall[0] -= eoffset; ewall[m+1] += fwall; + + if (evflag) { + if (side < 0) vn = -fwall*delta; + else vn = fwall*delta; + v_tally(dim, i, vn); + } } if (onflag) error->one(FLERR,"Particle on or inside fix wall surface"); diff --git a/src/KOKKOS/fix_wall_lj93_kokkos.cpp b/src/KOKKOS/fix_wall_lj93_kokkos.cpp index b0f7e0bda4..720b844601 100644 --- a/src/KOKKOS/fix_wall_lj93_kokkos.cpp +++ b/src/KOKKOS/fix_wall_lj93_kokkos.cpp @@ -31,6 +31,7 @@ FixWallLJ93Kokkos::FixWallLJ93Kokkos(LAMMPS *lmp, int narg, char **a execution_space = ExecutionSpaceFromDevice::space; datamask_read = EMPTY_MASK; datamask_modify = EMPTY_MASK; + virial_flag = 0 } /* ---------------------------------------------------------------------- diff --git a/src/MISC/fix_efield.cpp b/src/MISC/fix_efield.cpp index 8a6b7937b5..c0384adc41 100644 --- a/src/MISC/fix_efield.cpp +++ b/src/MISC/fix_efield.cpp @@ -55,6 +55,7 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) : extscalar = 1; respa_level_support = 1; ilevel_respa = 0; + virial_flag = 1; qe2f = force->qe2f; xstr = ystr = zstr = NULL; @@ -257,6 +258,11 @@ void FixEfield::post_force(int vflag) imageint *image = atom->image; int nlocal = atom->nlocal; + // energy and virial setup + + if (vflag) v_setup(vflag); + else evflag = 0; + // reallocate efield array if necessary if (varflag == ATOM && atom->nmax > maxatom) { @@ -281,6 +287,7 @@ void FixEfield::post_force(int vflag) double **x = atom->x; double fx,fy,fz; + double v[6]; // constant efield @@ -306,6 +313,15 @@ void FixEfield::post_force(int vflag) fsum[1] += fx; fsum[2] += fy; fsum[3] += fz; + if (evflag) { + v[0] = fx*unwrap[0]; + v[1] = fy*unwrap[1]; + v[2] = fz*unwrap[2]; + v[3] = fx*unwrap[1]; + v[4] = fx*unwrap[2]; + v[5] = fy*unwrap[2]; + v_tally(i, v); + } } } diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp index 551e049238..56106cc3e0 100644 --- a/src/MOLECULE/fix_cmap.cpp +++ b/src/MOLECULE/fix_cmap.cpp @@ -76,6 +76,7 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) : restart_peratom = 1; peatom_flag = 1; virial_flag = 1; + thermo_virial = 1; peratom_freq = 1; scalar_flag = 1; global_freq = 1; diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 9d46968273..d2a770cc47 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -74,6 +74,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : time_integrate = 1; rigid_flag = 1; virial_flag = 1; + thermo_virial = 1; create_attribute = 1; dof_flag = 1; enforce2d_flag = 1; diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 60afeecbf0..1404c3bf58 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -78,6 +78,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : time_integrate = 1; rigid_flag = 1; virial_flag = 1; + thermo_virial = 1; create_attribute = 1; dof_flag = 1; enforce2d_flag = 1; diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index f08228f3d3..c624fe2fe8 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -60,6 +60,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : MPI_Comm_size(world,&nprocs); virial_flag = 1; + thermo_virial = 1; create_attribute = 1; dof_flag = 1; diff --git a/src/fix_external.cpp b/src/fix_external.cpp index 85523b912c..ab942f874c 100644 --- a/src/fix_external.cpp +++ b/src/fix_external.cpp @@ -37,6 +37,7 @@ FixExternal::FixExternal(LAMMPS *lmp, int narg, char **arg) : scalar_flag = 1; global_freq = 1; virial_flag = 1; + thermo_virial = 1; extscalar = 1; if (strcmp(arg[3],"pf/callback") == 0) { diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index 8b569cafc6..2c3aa4afa0 100644 --- a/src/fix_wall.cpp +++ b/src/fix_wall.cpp @@ -45,6 +45,7 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : extvector = 1; respa_level_support = 1; ilevel_respa = 0; + virial_flag = 1; // parse args @@ -298,7 +299,12 @@ void FixWall::pre_force(int vflag) void FixWall::post_force(int vflag) { + + // energy and virial setup + eflag = 0; + if (vflag) v_setup(vflag); + else evflag = 0; for (int m = 0; m <= nwall; m++) ewall[m] = 0.0; // coord = current position of wall diff --git a/src/fix_wall_harmonic.cpp b/src/fix_wall_harmonic.cpp index 41d630b7fe..f91295d469 100644 --- a/src/fix_wall_harmonic.cpp +++ b/src/fix_wall_harmonic.cpp @@ -34,6 +34,7 @@ FixWallHarmonic::FixWallHarmonic(LAMMPS *lmp, int narg, char **arg) : void FixWallHarmonic::wall_particle(int m, int which, double coord) { double delta,dr,fwall; + double vn; double **x = atom->x; double **f = atom->f; @@ -60,6 +61,12 @@ void FixWallHarmonic::wall_particle(int m, int which, double coord) f[i][dim] -= fwall; ewall[0] += epsilon[m]*dr*dr; ewall[m+1] += fwall; + + if (evflag) { + if (side < 0) vn = -fwall*delta; + else vn = fwall*delta; + v_tally(dim, i, vn); + } } if (onflag) error->one(FLERR,"Particle on or inside fix wall surface"); diff --git a/src/fix_wall_lj1043.cpp b/src/fix_wall_lj1043.cpp index 72d71fbb3c..e8e706ed5e 100644 --- a/src/fix_wall_lj1043.cpp +++ b/src/fix_wall_lj1043.cpp @@ -52,6 +52,7 @@ void FixWallLJ1043::precompute(int m) void FixWallLJ1043::wall_particle(int m, int which, double coord) { double delta,rinv,r2inv,r4inv,r10inv,fwall; + double vn; double **x = atom->x; double **f = atom->f; @@ -79,5 +80,11 @@ void FixWallLJ1043::wall_particle(int m, int which, double coord) ewall[0] += coeff1[m]*r10inv - coeff2[m]*r4inv - coeff3[m]*pow(delta+coeff4[m],-3.0) - offset[m]; ewall[m+1] += fwall; + + if (evflag) { + if (side < 0) vn = -fwall*delta; + else vn = fwall*delta; + v_tally(dim, i, vn); + } } } diff --git a/src/fix_wall_lj126.cpp b/src/fix_wall_lj126.cpp index f6d8654eea..22199fed5b 100644 --- a/src/fix_wall_lj126.cpp +++ b/src/fix_wall_lj126.cpp @@ -48,6 +48,7 @@ void FixWallLJ126::precompute(int m) void FixWallLJ126::wall_particle(int m, int which, double coord) { double delta,rinv,r2inv,r6inv,fwall; + double vn; double **x = atom->x; double **f = atom->f; @@ -76,6 +77,12 @@ void FixWallLJ126::wall_particle(int m, int which, double coord) f[i][dim] -= fwall; ewall[0] += r6inv*(coeff3[m]*r6inv - coeff4[m]) - offset[m]; ewall[m+1] += fwall; + + if (evflag) { + if (side < 0) vn = -fwall*delta; + else vn = fwall*delta; + v_tally(dim, i, vn); + } } if (onflag) error->one(FLERR,"Particle on or inside fix wall surface"); diff --git a/src/fix_wall_lj93.cpp b/src/fix_wall_lj93.cpp index f7a3aaf14c..281a1fafc4 100644 --- a/src/fix_wall_lj93.cpp +++ b/src/fix_wall_lj93.cpp @@ -49,6 +49,7 @@ void FixWallLJ93::precompute(int m) void FixWallLJ93::wall_particle(int m, int which, double coord) { double delta,rinv,r2inv,r4inv,r10inv,fwall; + double vn; double **x = atom->x; double **f = atom->f; @@ -79,6 +80,12 @@ void FixWallLJ93::wall_particle(int m, int which, double coord) ewall[0] += coeff3[m]*r4inv*r4inv*rinv - coeff4[m]*r2inv*rinv - offset[m]; ewall[m+1] += fwall; + + if (evflag) { + if (side < 0) vn = -fwall*delta; + else vn = fwall*delta; + v_tally(dim, i, vn); + } } if (onflag) error->one(FLERR,"Particle on or inside fix wall surface"); diff --git a/src/fix_wall_region.cpp b/src/fix_wall_region.cpp index f543c12171..1d22e6141b 100644 --- a/src/fix_wall_region.cpp +++ b/src/fix_wall_region.cpp @@ -25,11 +25,13 @@ #include "output.h" #include "respa.h" #include "error.h" +#include "math_const.h" using namespace LAMMPS_NS; using namespace FixConst; +using namespace MathConst; -enum{LJ93,LJ126,COLLOID,HARMONIC}; +enum{LJ93,LJ126,LJ1043,COLLOID,HARMONIC}; /* ---------------------------------------------------------------------- */ @@ -47,6 +49,7 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) : extvector = 1; respa_level_support = 1; ilevel_respa = 0; + virial_flag = 1; // parse args @@ -59,6 +62,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],"colloid") == 0) style = COLLOID; else if (strcmp(arg[4],"harmonic") == 0) style = HARMONIC; else error->all(FLERR,"Illegal fix wall/region command"); @@ -143,6 +147,20 @@ void FixWallRegion::init() double r2inv = 1.0/(cutoff*cutoff); double r6inv = r2inv*r2inv*r2inv; offset = r6inv*(coeff3*r6inv - coeff4); + } else if (style == LJ1043) { + coeff1 = MY_2PI * 2.0/5.0 * epsilon * pow(sigma,10.0); + coeff2 = MY_2PI * epsilon * pow(sigma,4.0); + coeff3 = MY_2PI * pow(2.0,1/2.0) / 3 * epsilon * pow(sigma,3.0); + coeff4 = 0.61 / pow(2.0,1/2.0) * sigma; + 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 == COLLOID) { coeff1 = -4.0/315.0 * epsilon * pow(sigma,6.0); coeff2 = -2.0/3.0 * epsilon; @@ -186,6 +204,7 @@ void FixWallRegion::post_force(int vflag) { int i,m,n; double rinv,fx,fy,fz,tooclose; + double delx, dely, delz, v[6]; double **x = atom->x; double **f = atom->f; @@ -198,13 +217,18 @@ void FixWallRegion::post_force(int vflag) int onflag = 0; + // energy and virial setup + + eflag = 0; + if (vflag) v_setup(vflag); + else evflag = 0; + // region->match() insures particle is in region or on surface, else error // if returned contact dist r = 0, is on surface, also an error // in COLLOID case, r <= radius is an error // initilize ewall after region->prematch(), // so a dynamic region can access last timestep values - eflag = 0; ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0; for (i = 0; i < nlocal; i++) @@ -226,19 +250,32 @@ 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 == COLLOID) colloid(region->contact[m].r,radius[i]); else harmonic(region->contact[m].r); - ewall[0] += eng; - fx = fwall * region->contact[m].delx * rinv; - fy = fwall * region->contact[m].dely * rinv; - fz = fwall * region->contact[m].delz * rinv; + delx = region->contact[m].delx; + dely = region->contact[m].dely; + delz = region->contact[m].delz; + fx = fwall * delx * rinv; + fy = fwall * dely * rinv; + fz = fwall * delz * rinv; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; ewall[1] -= fx; ewall[2] -= fy; ewall[3] -= fz; + ewall[0] += eng; + if (evflag) { + v[0] = fx*delx; + v[1] = fy*dely; + v[2] = fz*delz; + v[3] = fx*dely; + v[4] = fx*delz; + v[5] = fy*delz; + v_tally(i, v); + } } } @@ -319,6 +356,23 @@ void FixWallRegion::lj126(double r) eng = r6inv*(coeff3*r6inv - coeff4) - offset; } +/* ---------------------------------------------------------------------- + LJ 10/4/3 interaction for particle with wall + compute eng and fwall = magnitude of wall force +------------------------------------------------------------------------- */ + +void FixWallRegion::lj1043(double r) +{ + double rinv = 1.0/r; + double r2inv = rinv*rinv; + double r4inv = r2inv*r2inv; + double r10inv = r4inv*r4inv*r2inv; + fwall = coeff5*r10inv*rinv - coeff6*r4inv*rinv - + coeff7*pow(r+coeff4,-4.0); + eng = coeff1*r10inv - coeff2*r4inv - + coeff3*pow(r+coeff4,-3.0) - offset; +} + /* ---------------------------------------------------------------------- colloid interaction for finite-size particle of rad with wall compute eng and fwall = magnitude of wall force diff --git a/src/fix_wall_region.h b/src/fix_wall_region.h index cecd7faa33..e3688c99ee 100644 --- a/src/fix_wall_region.h +++ b/src/fix_wall_region.h @@ -47,10 +47,12 @@ class FixWallRegion : public Fix { char *idregion; double coeff1,coeff2,coeff3,coeff4,offset; + double coeff5,coeff6,coeff7; double eng,fwall; void lj93(double); void lj126(double); + void lj1043(double); void colloid(double, double); void harmonic(double); };