jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #649 from akohlmey/fix-virial-aidan 

Add support for selected fixes to optionally contribute to the virial
This commit is contained in:
Steve Plimpton 2017-09-15 15:25:09 -06:00 committed by GitHub
parent 27f99534ec 99e5dc7104
commit d101fe3e79
30 changed files with 361 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Section_packages.txt
View File

@ -2191,7 +2191,7 @@ src/USER-MESO/README
"pair_style edpd"_pair_meso.html
"pair_style mdpd"_pair_meso.html
"pair_style tdpd"_pair_meso.html
"fix mvv/dpd"_fix_mvv.html
"fix mvv/dpd"_fix_mvv_dpd.html
examples/USER/meso
http://lammps.sandia.gov/movies.html#mesodpd :ul

5
doc/src/fix_addforce.txt
View File

@ -139,6 +139,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 added forces on atoms 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

5
doc/src/fix_cmap.txt
View File

@ -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 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

5
doc/src/fix_efield.txt
View File

@ -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 added forces on atoms 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.

5
doc/src/fix_external.txt
View File

@ -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 interactions computed by the
external program 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

33
doc/src/fix_modify.txt
View File

@ -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.

5
doc/src/fix_rigid.txt
View File

@ -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 keeping the objects rigid 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

5
doc/src/fix_shake.txt
View File

@ -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 keeping the constraints 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

5
doc/src/fix_smd.txt
View File

@ -101,6 +101,11 @@ See the "read_restart"_read_restart.html command for info on how to
re-specify a fix in an input script that reads a restart file, so that
the operation of the fix continues in an uninterrupted fashion.
The "fix_modify"_fix_modify.html {virial} option is supported by this
fix to add the contribution due to the added forces on atoms 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.

5
doc/src/fix_wall.txt
View File

@ -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.

56
doc/src/fix_wall_region.txt
View File

@ -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.

7
src/COLLOID/fix_wall_colloid.cpp
View File

@ -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");

1
src/KOKKOS/fix_wall_lj93_kokkos.cpp
View File

@ -31,6 +31,7 @@ FixWallLJ93Kokkos<DeviceType>::FixWallLJ93Kokkos(LAMMPS *lmp, int narg, char **a
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
virial_flag = 0
}
/* ----------------------------------------------------------------------

16
src/MISC/fix_efield.cpp
View File

@ -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);
}
}
}

1
src/MOLECULE/fix_cmap.cpp
View File

@ -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;

1
src/RIGID/fix_rigid.cpp
View File

@ -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;

1
src/RIGID/fix_rigid_small.cpp
View File

@ -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;

1
src/RIGID/fix_shake.cpp
View File

@ -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;

29
src/USER-MISC/fix_smd.cpp
View File

@ -58,6 +58,7 @@ FixSMD::FixSMD(LAMMPS *lmp, int narg, char **arg) :
extvector = 1;
respa_level_support = 1;
ilevel_respa = 0;
virial_flag = 1;
int argoffs=3;
if (strcmp(arg[argoffs],"cvel") == 0) {
@ -181,6 +182,11 @@ void FixSMD::setup(int vflag)
void FixSMD::post_force(int vflag)
{
// energy and virial setup
if (vflag) v_setup(vflag);
else evflag = 0;
if (styleflag & SMD_TETHER) smd_tether();
else smd_couple();
@ -238,12 +244,15 @@ void FixSMD::smd_tether()
// apply restoring force to atoms in group
// f = -k*(r-r0)*mass/masstotal
double **x = atom->x;
double **f = atom->f;
imageint *image = atom->image;
int *mask = atom->mask;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
double massfrac;
double unwrap[3],v[6];
int nlocal = atom->nlocal;
ftotal[0] = ftotal[1] = ftotal[2] = 0.0;
@ -259,6 +268,16 @@ void FixSMD::smd_tether()
ftotal[0] -= fx*massfrac;
ftotal[1] -= fy*massfrac;
ftotal[2] -= fz*massfrac;
if (evflag) {
domain->unmap(x[i],image[i],unwrap);
v[0] = fx*massfrac*unwrap[0];
v[1] = fy*massfrac*unwrap[1];
v[2] = fz*massfrac*unwrap[2];
v[3] = fx*massfrac*unwrap[1];
v[4] = fx*massfrac*unwrap[2];
v[5] = fy*massfrac*unwrap[2];
v_tally(i, v);
}
}
} else {
for (int i = 0; i < nlocal; i++)
@ -270,6 +289,16 @@ void FixSMD::smd_tether()
ftotal[0] -= fx*massfrac;
ftotal[1] -= fy*massfrac;
ftotal[2] -= fz*massfrac;
if (evflag) {
domain->unmap(x[i],image[i],unwrap);
v[0] = fx*massfrac*unwrap[0];
v[1] = fy*massfrac*unwrap[1];
v[2] = fz*massfrac*unwrap[2];
v[3] = fx*massfrac*unwrap[1];
v[4] = fx*massfrac*unwrap[2];
v[5] = fy*massfrac*unwrap[2];
v_tally(i, v);
}
}
}
}

87
src/fix.cpp
View File

@ -31,7 +31,7 @@ int Fix::instance_total = 0;
/* ---------------------------------------------------------------------- */
Fix::Fix(LAMMPS *lmp, int narg, char **arg) :
Fix::Fix(LAMMPS *lmp, int narg, char **arg) :
Pointers(lmp),
id(NULL), style(NULL), extlist(NULL), vector_atom(NULL), array_atom(NULL),
vector_local(NULL), array_local(NULL), eatom(NULL), vatom(NULL)
@ -61,6 +61,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) :
force_reneighbor = 0;
box_change_size = box_change_shape = box_change_domain = 0;
thermo_energy = 0;
thermo_virial = 0;
rigid_flag = 0;
peatom_flag = 0;
virial_flag = 0;
@ -140,8 +141,20 @@ void Fix::modify_params(int narg, char **arg)
} else if (strcmp(arg[iarg],"energy") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
if (strcmp(arg[iarg+1],"no") == 0) thermo_energy = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) thermo_energy = 1;
else error->all(FLERR,"Illegal fix_modify command");
else if (strcmp(arg[iarg+1],"yes") == 0) {
if (!(THERMO_ENERGY & setmask()))
error->all(FLERR,"Illegal fix_modify command");
thermo_energy = 1;
} else error->all(FLERR,"Illegal fix_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"virial") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
if (strcmp(arg[iarg+1],"no") == 0) thermo_virial = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) {
if (virial_flag == 0)
error->all(FLERR,"Illegal fix_modify command");
thermo_virial = 1;
} else error->all(FLERR,"Illegal fix_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"respa") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
@ -214,15 +227,22 @@ void Fix::ev_setup(int eflag, int vflag)
}
/* ----------------------------------------------------------------------
setup for virial computation
see integrate::ev_set() for values of vflag (0-6)
fixes call this if use v_tally()
if thermo_virial is on:
setup for virial computation
see integrate::ev_set() for values of vflag (0-6)
fixes call this if use v_tally()
else: set evflag=0
------------------------------------------------------------------------- */
void Fix::v_setup(int vflag)
{
int i,n;
if (!thermo_virial) {
evflag = 0;
return;
}
evflag = 1;
vflag_global = vflag % 4;
@ -316,3 +336,58 @@ void Fix::v_tally(int n, int *list, double total, double *v)
}
}
}
/* ----------------------------------------------------------------------
tally virial into global and per-atom accumulators
i = local index of atom
v = total virial for the interaction
increment global virial by v
increment per-atom virial by v
this method can be used when fix computes forces in post_force()
and the force depends on a distance to some external object
e.g. fix wall/lj93: compute virial only on owned atoms
------------------------------------------------------------------------- */
void Fix::v_tally(int i, double *v)
{
if (vflag_global) {
virial[0] += v[0];
virial[1] += v[1];
virial[2] += v[2];
virial[3] += v[3];
virial[4] += v[4];
virial[5] += v[5];
}
if (vflag_atom) {
vatom[i][0] += v[0];
vatom[i][1] += v[1];
vatom[i][2] += v[2];
vatom[i][3] += v[3];
vatom[i][4] += v[4];
vatom[i][5] += v[5];
}
}
/* ----------------------------------------------------------------------
tally virial component into global and per-atom accumulators
n = index of virial component (0-5)
i = local index of atom
vn = nth component of virial for the interaction
increment nth component of global virial by vn
increment nth component of per-atom virial by vn
this method can be used when fix computes forces in post_force()
and the force depends on a distance to some external object
e.g. fix wall/lj93: compute virial only on owned atoms
------------------------------------------------------------------------- */
void Fix::v_tally(int n, int i, double vn)
{
if (vflag_global)
virial[n] += vn;
if (vflag_atom)
vatom[i][n] += vn;
}

3
src/fix.h
View File

@ -36,6 +36,7 @@ class Fix : protected Pointers {
bigint next_reneighbor; // next timestep to force a reneighboring
int thermo_energy; // 1 if fix_modify enabled ThEng, 0 if not
int thermo_virial; // 1 if fix_modify enabled ThVir, 0 if not
int nevery; // how often to call an end_of_step fix
int rigid_flag; // 1 if Fix integrates rigid bodies, 0 if not
int peatom_flag; // 1 if Fix contributes per-atom eng, 0 if not
@ -223,6 +224,8 @@ class Fix : protected Pointers {
void ev_tally(int, int *, double, double, double *);
void v_setup(int);
void v_tally(int, int *, double, double *);
void v_tally(int, double *);
void v_tally(int, int, double);
// union data struct for packing 32-bit and 64-bit ints into double bufs
// see atom_vec.h for documentation

52
src/fix_addforce.cpp
View File

@ -51,6 +51,7 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
extvector = 1;
respa_level_support = 1;
ilevel_respa = 0;
virial_flag = 1;
xstr = ystr = zstr = NULL;
@ -237,10 +238,16 @@ void FixAddForce::post_force(int vflag)
double **f = atom->f;
int *mask = atom->mask;
imageint *image = atom->image;
double v[6];
int nlocal = atom->nlocal;
if (update->ntimestep % nevery) return;
// energy and virial setup
if (vflag) v_setup(vflag);
else evflag = 0;
if (lmp->kokkos)
atom->sync_modify(Host, (unsigned int) (F_MASK | MASK_MASK),
(unsigned int) F_MASK);
@ -283,6 +290,15 @@ void FixAddForce::post_force(int vflag)
f[i][0] += xvalue;
f[i][1] += yvalue;
f[i][2] += zvalue;
if (evflag) {
v[0] = xvalue * unwrap[0];
v[1] = yvalue * unwrap[1];
v[2] = zvalue * unwrap[2];
v[3] = xvalue * unwrap[1];
v[4] = xvalue * unwrap[2];
v[5] = yvalue * unwrap[2];
v_tally(i,v);
}
}
// variable force, wrap with clear/add
@ -290,6 +306,7 @@ void FixAddForce::post_force(int vflag)
// wrap with clear/add
} else {
double unwrap[3];
modify->clearstep_compute();
@ -307,20 +324,39 @@ void FixAddForce::post_force(int vflag)
modify->addstep_compute(update->ntimestep + 1);
for (int i = 0; i < nlocal; i++)
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;
if (estyle == ATOM) foriginal[0] += sforce[i][3];
domain->unmap(x[i],image[i],unwrap);
if (xstyle == ATOM) xvalue = sforce[i][0];
if (ystyle == ATOM) yvalue = sforce[i][1];
if (zstyle == ATOM) zvalue = sforce[i][2];
if (estyle == ATOM) {
foriginal[0] += sforce[i][3];
} else {
if (xstyle) foriginal[0] -= xvalue*unwrap[0];
if (ystyle) foriginal[0] -= yvalue*unwrap[1];
if (zstyle) foriginal[0] -= zvalue*unwrap[2];
}
foriginal[1] += f[i][0];
foriginal[2] += f[i][1];
foriginal[3] += f[i][2];
if (xstyle == ATOM) f[i][0] += sforce[i][0];
else if (xstyle) f[i][0] += xvalue;
if (ystyle == ATOM) f[i][1] += sforce[i][1];
else if (ystyle) f[i][1] += yvalue;
if (zstyle == ATOM) f[i][2] += sforce[i][2];
else if (zstyle) f[i][2] += zvalue;
if (xstyle) f[i][0] += xvalue;
if (ystyle) f[i][1] += yvalue;
if (zstyle) f[i][2] += zvalue;
if (evflag) {
v[0] = xstyle ? xvalue*unwrap[0] : 0.0;
v[1] = ystyle ? yvalue*unwrap[1] : 0.0;
v[2] = zstyle ? zvalue*unwrap[2] : 0.0;
v[3] = xstyle ? xvalue*unwrap[1] : 0.0;
v[4] = xstyle ? xvalue*unwrap[2] : 0.0;
v[5] = ystyle ? yvalue*unwrap[2] : 0.0;
v_tally(i, v);
}
}
}
}
}

1
src/fix_external.cpp
View File

@ -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) {

6
src/fix_wall.cpp
View File

@ -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

7
src/fix_wall_harmonic.cpp
View File

@ -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");

7
src/fix_wall_lj1043.cpp
View File

@ -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);
}
}
}

7
src/fix_wall_lj126.cpp
View File

@ -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");

7
src/fix_wall_lj93.cpp
View File

@ -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");

66
src/fix_wall_region.cpp
View File

@ -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

2
src/fix_wall_region.h
View File

@ -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);
};