diff --git a/doc/fix_indent.html b/doc/fix_indent.html index 89f61fe843..1481bff711 100644 --- a/doc/fix_indent.html +++ b/doc/fix_indent.html @@ -13,13 +13,13 @@

Syntax:

-
fix ID group-ID indent k keyword values ... 
+
fix ID group-ID indent K keyword values ... 
 

Examples:

-
fix 1 all indent 10.0 sphere 0.0 0.0 15.0 3.0 vel 0.0 0.0 -1.0
+
fix 1 all indent 10.0 sphere 0.0 0.0 15.0 3.0
+fix 1 all indent 10.0 sphere v_x v_y 0.0 v_radius side in
 fix 2 flow indent 10.0 cylinder z 0.0 0.0 10.0 units box 
 

Description: @@ -67,21 +66,26 @@ must set one of those 3 keywords.

A spherical indenter exerts a force of magnitude

-
F(r) = - k (r - R)^2 
+
F(r) = - K (r - R)^2 
 
-

on each atom where k is the specified force constant, r is the +

on each atom where K is the specified force constant, r is the distance from the atom to the center of the indenter, and R is the radius of the indenter. The force is repulsive and F(r) = 0 for r > -R. The calculation of distance to the indenter center accounts -for periodic boundaries, which means the indenter can effectively -straddle one or more periodic boundaries. +R.

A cylindrical indenter exerts the same force, except that r is the distance from the atom to the center axis of the cylinder. The -cylinder extends infinitely along its axis. The calculation of -distance to the indenter axis accounts for periodic boundaries, which -means the indenter can effectively straddle one or more periodic -boundaries. +cylinder extends infinitely along its axis. +

+

Spherical and cylindrical indenters account for periodic boundaries in +two ways. First, the center point of a spherical indenter (x,y,z) or +axis of a cylindrical indenter (c1,c2) is remapped back into the +simulation box, if the box is periodic in a particular dimension. +This occurs every timestep if the indenter geometry is specified with +a variable (see below), e.g. it is moving over time. Second, the +calculation of distance to the indenter center or axis accounts for +periodic boundaries. Both of these mean that an indenter can +effectively move through and straddle one or more periodic boundaries.

A planar indenter is really an axis-aligned infinite-extent wall exerting the same force on atoms in the system, where R is the @@ -93,24 +97,46 @@ towards the hi end of the box and atoms with a coordinate higher than the plane's current position will feel no force. Vice versa if side is specified as hi.

-

If the vel keyword is specified, the center (or axis or position) of -the spherical (or cylindrical or planar) indenter will move during the -simulation, based on its initial position (x,y,z), the specified -(vx,vy,vz), and the time elapsed since the beginning of the -simulation. For periodic systems and spherical or cylindrical -indenters, the new position of the center or axis is wrapped back into -the periodic box as needed. See the note below about making the -indenter move continuously across multiple runs. +

Any of the 4 quantities defining a spherical indenter's geometry can +be specified as an equal-style variable, namely x, y, +z, or R. Similarly, for a cylindrical indenter, any of c1, +c2, or R, can be a variable. For a planar indenter, pos can be +a variable. If the value is a variable, it should be specified as +v_ID, where ID is the variable ID. In this case, the variable will be +evaluated each timestep, and its value used to define the indenter +geometry.

-

If the rstart keyword is specified, then the radius of the indenter -is a time-dependent quantity. This only applies to spherical or -cylindrical indenters. R0 is the value assigned at the start of the -run; R is the value at the end. At intermediate times, the radius is -linearly interpolated between these two values. This option can be -used, for example, to grow/shrink a void within the simulation box. -See the note below about making the radius change continuously across -multiple runs. +

Note that equal-style variables can specify formulas with various +mathematical functions, and include thermo_style +command keywords for the simulation box parameters and timestep and +elapsed time. Thus it is easy to specify indenter properties that +change as a function of time or span consecutive runs in a continuous +fashion. For the latter, see the start and stop keywords of the +run command and the elaplong keyword of thermo_style +custom for details.

+

For example, if a spherical indenter's x-position is specfied as v_x, +then this variable definition will keep it's center at a relative +position in the simulation box, 1/4 of the way from the left edge to +the right edge, even if the box size changes: +

+
variable x equal "xlo + 0.25*lx" 
+
+

Similarly, these variable definitions will move the indenter at a +constant velocity: +

+
variable x0 equal 2.5
+variable vx equal 5.0
+variable x equal "v_x0 + step*dt*v_vx" 
+
+

If a spherical indenter's radius is specified as v_r, then these +variable definitions will grow the size of the indenter at a specfied +rate. +

+
variable r0 equal 0.0
+variable rate equal 1.0
+variable r equal "v_r0 + step*dt*v_rate" 
+

If the side keyword is specified as out, which is the default, then particles outside the indenter are pushded away from its outer surface, as described above. This only applies to spherical or @@ -121,15 +147,27 @@ is now a containing wall that traps the particles inside it. If the radius shrinks over time, it will squeeze the particles.

The units keyword determines the meaning of the distance units used -to define the indenter. A box value selects standard distance units -as defined by the units command, e.g. Angstroms for units -= real or metal. A lattice value means the distance units are in -lattice spacings. The lattice command must have been -previously used to define the lattice spacing. Note that the units -choice affects not only the indenter's physical geometry, but also its -velocity and force constant since they are defined in terms of -distance as well. +to define the indenter geometry. A box value selects standard +distance units as defined by the units command, +e.g. Angstroms for units = real or metal. A lattice value means the +distance units are in lattice spacings. The lattice +command must have been previously used to define the lattice spacing. +Note that the units keyword only affects indenter geometry parameters +specified directly with numbers, not those specified as variables. In +the latter case, you should use the xlat, ylat, zlat keywords of +the thermo_style command if you want to include +lattice spacings in a variable formula.

+

The force constant K is not affected by the units keyword. It is +always in force/distance^2 units where force and distance are defined +by the units command. If you wish K to be scaled by the +lattice spacing, you can define K with a variable whose formula +contains xlat, ylat, zlat keywords of the +thermo_style command, e.g. +

+
variable k equal 100.0/xlat/xlat
+fix 1 all indent $k sphere ... 
+

Restart, fix_modify, output, run start/stop, minimize info:

No information about this fix is written to binary restart @@ -146,17 +184,13 @@ indenter), which can be accessed by various commands. The scalar and vector values calculated by this fix are "extensive".

-

This fix can adjust the indenter position and radius over multiple -runs, using the start and stop keywords of the run -command. See the run command for details of how to do -this. If you do not do this, the indenter position and readius will -be reset to their specified initial values at the beginning of each -run. -

The forces due to this fix are imposed during an energy minimization, -invoked by the minimize command. The rstart keyword -does not change the indenter radius during an energy minimization; the -indenter always has a radius of its final value R in that case. +invoked by the minimize command. Note that if you +define the indenter geometry with a variable using a time-dependent +formula, LAMMPS uses the iteration count in the minimizer as the +timestep. But it is almost certainly a bad idea to have the indenter +change its position or size during a minimization. LAMMPS does not +check if you have done this.

IMPORTANT NOTE: If you want the atom/indenter interaction energy to be included in the total potential energy of the system (the quantity @@ -169,6 +203,6 @@ being minimized), you must enable the fix_modify

Default:

-

The option defaults are vel = 0,0,0, side = out, and units = lattice. +

The option defaults are side = out and units = lattice.

diff --git a/doc/fix_indent.txt b/doc/fix_indent.txt index 0460bdb916..e936a3f1e3 100644 --- a/doc/fix_indent.txt +++ b/doc/fix_indent.txt @@ -10,29 +10,27 @@ fix indent command :h3 [Syntax:] -fix ID group-ID indent k keyword values ... :pre +fix ID group-ID indent K keyword values ... :pre ID, group-ID are documented in "fix"_fix.html command :ulb,l indent = style name of this fix command :l -k = force constant for indenter surface (force/distance^2 units) :l +K = force constant for indenter surface (force/distance^2 units) :l one or more keyword/value pairs may be appended :l keyword = {sphere} or {cylinder} or {plane} or {vel} or {rstart} or {side} or {units} :l {sphere} args = x y z R x,y,z = initial position of center of indenter (distance units) R = sphere radius of indenter (distance units) + any of x,y,z,R can be a variable (see below) {cylinder} args = dim c1 c2 R dim = {x} or {y} or {z} = axis of cylinder c1,c2 = coords of cylinder axis in other 2 dimensions (distance units) R = cylinder radius of indenter (distance units) + any of c1,c2,R can be a variable (see below) {plane} args = dim pos side dim = {x} or {y} or {z} = plane perpendicular to this dimension pos = position of plane in dimension x, y, or z (distance units) + pos can be a variable (see below) side = {lo} or {hi} - {vel} args = vx vy vz - vx,vy,vz = velocity of center of indenter (velocity units) - {rstart} value = R0 - R0 = sphere or cylinder radius at start of run (distance units) - R is value at end of run, so indenter expands/contracts over time {side} value = {in} or {out} {in} = the indenter acts on particles inside the sphere or cylinder {out} = the indenter acts on particles outside the sphere or cylinder @@ -43,7 +41,8 @@ keyword = {sphere} or {cylinder} or {plane} or {vel} or {rstart} or {side} or {u [Examples:] -fix 1 all indent 10.0 sphere 0.0 0.0 15.0 3.0 vel 0.0 0.0 -1.0 +fix 1 all indent 10.0 sphere 0.0 0.0 15.0 3.0 +fix 1 all indent 10.0 sphere v_x v_y 0.0 v_radius side in fix 2 flow indent 10.0 cylinder z 0.0 0.0 10.0 units box :pre [Description:] @@ -58,21 +57,26 @@ must set one of those 3 keywords. A spherical indenter exerts a force of magnitude -F(r) = - k (r - R)^2 :pre +F(r) = - K (r - R)^2 :pre -on each atom where {k} is the specified force constant, {r} is the +on each atom where {K} is the specified force constant, {r} is the distance from the atom to the center of the indenter, and {R} is the radius of the indenter. The force is repulsive and F(r) = 0 for {r} > -{R}. The calculation of distance to the indenter center accounts -for periodic boundaries, which means the indenter can effectively -straddle one or more periodic boundaries. +{R}. A cylindrical indenter exerts the same force, except that {r} is the distance from the atom to the center axis of the cylinder. The -cylinder extends infinitely along its axis. The calculation of -distance to the indenter axis accounts for periodic boundaries, which -means the indenter can effectively straddle one or more periodic -boundaries. +cylinder extends infinitely along its axis. + +Spherical and cylindrical indenters account for periodic boundaries in +two ways. First, the center point of a spherical indenter (x,y,z) or +axis of a cylindrical indenter (c1,c2) is remapped back into the +simulation box, if the box is periodic in a particular dimension. +This occurs every timestep if the indenter geometry is specified with +a variable (see below), e.g. it is moving over time. Second, the +calculation of distance to the indenter center or axis accounts for +periodic boundaries. Both of these mean that an indenter can +effectively move through and straddle one or more periodic boundaries. A planar indenter is really an axis-aligned infinite-extent wall exerting the same force on atoms in the system, where {R} is the @@ -84,23 +88,45 @@ towards the hi end of the box and atoms with a coordinate higher than the plane's current position will feel no force. Vice versa if {side} is specified as {hi}. -If the {vel} keyword is specified, the center (or axis or position) of -the spherical (or cylindrical or planar) indenter will move during the -simulation, based on its initial position (x,y,z), the specified -(vx,vy,vz), and the time elapsed since the beginning of the -simulation. For periodic systems and spherical or cylindrical -indenters, the new position of the center or axis is wrapped back into -the periodic box as needed. See the note below about making the -indenter move continuously across multiple runs. +Any of the 4 quantities defining a spherical indenter's geometry can +be specified as an equal-style "variable"_variable, namely {x}, {y}, +{z}, or {R}. Similarly, for a cylindrical indenter, any of {c1}, +{c2}, or {R}, can be a variable. For a planar indenter, {pos} can be +a variable. If the value is a variable, it should be specified as +v_ID, where ID is the variable ID. In this case, the variable will be +evaluated each timestep, and its value used to define the indenter +geometry. -If the {rstart} keyword is specified, then the radius of the indenter -is a time-dependent quantity. This only applies to spherical or -cylindrical indenters. R0 is the value assigned at the start of the -run; R is the value at the end. At intermediate times, the radius is -linearly interpolated between these two values. This option can be -used, for example, to grow/shrink a void within the simulation box. -See the note below about making the radius change continuously across -multiple runs. +Note that equal-style variables can specify formulas with various +mathematical functions, and include "thermo_style"_thermo_style.html +command keywords for the simulation box parameters and timestep and +elapsed time. Thus it is easy to specify indenter properties that +change as a function of time or span consecutive runs in a continuous +fashion. For the latter, see the {start} and {stop} keywords of the +"run"_run.html command and the {elaplong} keyword of "thermo_style +custom"_thermo_style.html for details. + +For example, if a spherical indenter's x-position is specfied as v_x, +then this variable definition will keep it's center at a relative +position in the simulation box, 1/4 of the way from the left edge to +the right edge, even if the box size changes: + +variable x equal "xlo + 0.25*lx" :pre + +Similarly, these variable definitions will move the indenter at a +constant velocity: + +variable x0 equal 2.5 +variable vx equal 5.0 +variable x equal "v_x0 + step*dt*v_vx" :pre + +If a spherical indenter's radius is specified as v_r, then these +variable definitions will grow the size of the indenter at a specfied +rate. + +variable r0 equal 0.0 +variable rate equal 1.0 +variable r equal "v_r0 + step*dt*v_rate" :pre If the {side} keyword is specified as {out}, which is the default, then particles outside the indenter are pushded away from its outer @@ -112,14 +138,26 @@ is now a containing wall that traps the particles inside it. If the radius shrinks over time, it will squeeze the particles. The {units} keyword determines the meaning of the distance units used -to define the indenter. A {box} value selects standard distance units -as defined by the "units"_units.html command, e.g. Angstroms for units -= real or metal. A {lattice} value means the distance units are in -lattice spacings. The "lattice"_lattice.html command must have been -previously used to define the lattice spacing. Note that the units -choice affects not only the indenter's physical geometry, but also its -velocity and force constant since they are defined in terms of -distance as well. +to define the indenter geometry. A {box} value selects standard +distance units as defined by the "units"_units.html command, +e.g. Angstroms for units = real or metal. A {lattice} value means the +distance units are in lattice spacings. The "lattice"_lattice.html +command must have been previously used to define the lattice spacing. +Note that the units keyword only affects indenter geometry parameters +specified directly with numbers, not those specified as variables. In +the latter case, you should use the {xlat}, {ylat}, {zlat} keywords of +the "thermo_style"_thermo_style.html command if you want to include +lattice spacings in a variable formula. + +The force constant {K} is not affected by the {units} keyword. It is +always in force/distance^2 units where force and distance are defined +by the "units"_units.html command. If you wish K to be scaled by the +lattice spacing, you can define K with a variable whose formula +contains {xlat}, {ylat}, {zlat} keywords of the +"thermo_style"_thermo_style.html command, e.g. + +variable k equal 100.0/xlat/xlat +fix 1 all indent $k sphere ... :pre [Restart, fix_modify, output, run start/stop, minimize info:] @@ -137,17 +175,13 @@ indenter), which can be accessed by various "output commands"_Section_howto.html#4_15. The scalar and vector values calculated by this fix are "extensive". -This fix can adjust the indenter position and radius over multiple -runs, using the {start} and {stop} keywords of the "run"_run.html -command. See the "run"_run.html command for details of how to do -this. If you do not do this, the indenter position and readius will -be reset to their specified initial values at the beginning of each -run. - The forces due to this fix are imposed during an energy minimization, -invoked by the "minimize"_minimize.html command. The {rstart} keyword -does not change the indenter radius during an energy minimization; the -indenter always has a radius of its final value R in that case. +invoked by the "minimize"_minimize.html command. Note that if you +define the indenter geometry with a variable using a time-dependent +formula, LAMMPS uses the iteration count in the minimizer as the +timestep. But it is almost certainly a bad idea to have the indenter +change its position or size during a minimization. LAMMPS does not +check if you have done this. IMPORTANT NOTE: If you want the atom/indenter interaction energy to be included in the total potential energy of the system (the quantity @@ -160,4 +194,4 @@ being minimized), you must enable the "fix_modify"_fix_modify.html [Default:] -The option defaults are vel = 0,0,0, side = out, and units = lattice. +The option defaults are side = out and units = lattice. diff --git a/doc/thermo_style.html b/doc/thermo_style.html index 887b7b6021..ae7929e475 100644 --- a/doc/thermo_style.html +++ b/doc/thermo_style.html @@ -22,8 +22,8 @@
  one args = none
   multi args = none
   custom args = list of attributes
-    possible attributes = step, atoms, cpu, temp, press,
-                          pe, ke, etotal, enthalpy,
+    possible attributes = step, elapsed, elaplong, dt, cpu,
+                          atoms, temp, press, pe, ke, etotal, enthalpy,
                           evdwl, ecoul, epair, ebond, eangle, edihed, eimp,
                           emol, elong, etail,
                           vol, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi,
@@ -33,8 +33,11 @@
                           f_ID, f_ID[I], f_ID[I][J],
                           v_name
       step = timestep
-      atoms = # of atoms
+      elapsed = timesteps since start of this run
+      elaplong = timesteps since start of initial run in a series of runs
+      dt = timestep size
       cpu = elapsed CPU time
+      atoms = # of atoms
       temp = temperature
       press = pressure
       pe = total potential energy
@@ -192,6 +195,16 @@ etc.
 


+

The step, elapsed, and elaplong keywords refer to timestep +count. Step is the current timestep, or iteration count when a +minimization is being performed. Elapsed is the +number of timesteps elapsed since the beginning of this run. +Elaplong is the number of timesteps elapsed since the beginning of +an initial run in a series of runs. See the start and stop +keywords for the run for info on how to invoke a series of +runs that keep track of an initial starting time. If these keywords +are not used, then elapsed and elaplong are the same value. +

The c_ID and c_ID[I] and c_ID[I][J] keywords allow global values calculated by a compute to be output. As discussed on the compute doc page, computes can calculate global, diff --git a/doc/thermo_style.txt b/doc/thermo_style.txt index 09051ab15d..4013a0441e 100644 --- a/doc/thermo_style.txt +++ b/doc/thermo_style.txt @@ -17,8 +17,8 @@ args = list of arguments for a particular style :l {one} args = none {multi} args = none {custom} args = list of attributes - possible attributes = step, atoms, cpu, temp, press, - pe, ke, etotal, enthalpy, + possible attributes = step, elapsed, elaplong, dt, cpu, + atoms, temp, press, pe, ke, etotal, enthalpy, evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail, vol, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, @@ -28,8 +28,11 @@ args = list of arguments for a particular style :l f_ID, f_ID\[I\], f_ID\[I\]\[J\], v_name step = timestep - atoms = # of atoms + elapsed = timesteps since start of this run + elaplong = timesteps since start of initial run in a series of runs + dt = timestep size cpu = elapsed CPU time + atoms = # of atoms temp = temperature press = pressure pe = total potential energy @@ -186,6 +189,16 @@ etc. :line +The {step}, {elapsed}, and {elaplong} keywords refer to timestep +count. {Step} is the current timestep, or iteration count when a +"minimization"_minimize.html is being performed. {Elapsed} is the +number of timesteps elapsed since the beginning of this run. +{Elaplong} is the number of timesteps elapsed since the beginning of +an initial run in a series of runs. See the {start} and {stop} +keywords for the "run"_run.html for info on how to invoke a series of +runs that keep track of an initial starting time. If these keywords +are not used, then {elapsed} and {elaplong} are the same value. + The {c_ID} and {c_ID\[I\]} and {c_ID\[I\]\[J\]} keywords allow global values calculated by a compute to be output. As discussed on the "compute"_compute.html doc page, computes can calculate global, diff --git a/examples/indent/in.indent b/examples/indent/in.indent index d65df80498..5eac9aafb7 100644 --- a/examples/indent/in.indent +++ b/examples/indent/in.indent @@ -37,21 +37,22 @@ fix 1 all nve fix 2 lower setforce 0.0 0.0 0.0 fix 3 all temp/rescale 100 0.1 0.1 0.01 1.0 -# indenter - -fix 4 all indent 1000.0 sphere 10 13 0 5.0 vel 0.0 -0.02 0.0 -fix 5 all enforce2d - -# Run with indenter +# run with indenter timestep 0.003 +variable k equal 1000.0/xlat +variable y equal "13.0*ylat - step*dt*0.02*ylat" + +fix 4 all indent $k sphere 10 v_y 0 5.0 +fix 5 all enforce2d + thermo 1000 thermo_modify temp new dump 1 all atom 250 dump.indent run 30000 -# Run without indenter +# run without indenter unfix 4 run 30000 diff --git a/examples/indent/in.indent.min b/examples/indent/in.indent.min index d354ba3b09..a09f3dafc5 100644 --- a/examples/indent/in.indent.min +++ b/examples/indent/in.indent.min @@ -45,26 +45,29 @@ dump_modify 1 scale no minimize 1.0e-6 1.0e-6 1000 1000 -fix 4 all indent 5000.0 sphere 10 13.0 0 6.0 +variable k equal 5000.0/xlat +variable k1 equal 1000.0/xlat + +fix 4 all indent $k sphere 10 13.0 0 6.0 fix_modify 4 energy yes minimize 1.0e-6 1.0e-6 1000 1000 -fix 4 all indent 1000.0 sphere 10 12.5 0 6.0 +fix 4 all indent ${k1} sphere 10 12.5 0 6.0 fix_modify 4 energy yes minimize 1.0e-6 1.0e-6 1000 1000 -fix 4 all indent 1000.0 sphere 10 12.0 0 6.0 +fix 4 all indent ${k1} sphere 10 12.0 0 6.0 fix_modify 4 energy yes minimize 1.0e-6 1.0e-6 1000 1000 -fix 4 all indent 1000.0 sphere 10 11.4 0 6.0 +fix 4 all indent ${k1} sphere 10 11.4 0 6.0 fix_modify 4 energy yes minimize 1.0e-6 1.0e-6 1000 1000 -fix 4 all indent 1000.0 sphere 10 11.2 0 6.0 +fix 4 all indent ${k1} sphere 10 11.2 0 6.0 fix_modify 4 energy yes minimize 1.0e-6 1.0e-6 1000 1000 -fix 4 all indent 1000.0 sphere 10 11.0 0 6.0 +fix 4 all indent ${k1} sphere 10 11.0 0 6.0 fix_modify 4 energy yes minimize 1.0e-6 1.0e-6 1000 1000 diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index b440726499..cf88c7ca76 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -20,6 +20,8 @@ #include "stdlib.h" #include "fix_indent.h" #include "atom.h" +#include "input.h" +#include "variable.h" #include "domain.h" #include "lattice.h" #include "update.h" @@ -47,6 +49,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : extvector = 1; k = atof(arg[3]); + k3 = k/3.0; // read options from end of input line @@ -65,41 +68,17 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : } else xscale = yscale = zscale = 1.0; - // apply scaling factors to force constant, velocity, and geometry + // apply scaling factors to geometry - k /= xscale; - k3 = k/3.0; - vx *= xscale; - vy *= yscale; - vz *= zscale; - - if (istyle == SPHERE) { - x0 *= xscale; - y0 *= yscale; - z0 *= zscale; - r0_stop *= xscale; - r0_start *= xscale; - } else if (istyle == CYLINDER) { - if (cdim == 0) { - c1 *= yscale; - c2 *= zscale; - r0_stop *= xscale; - r0_start *= xscale; - } else if (cdim == 1) { - c1 *= xscale; - c2 *= zscale; - r0_stop *= yscale; - r0_start *= yscale; - } else if (cdim == 2) { - c1 *= xscale; - c2 *= yscale; - r0_stop *= zscale; - r0_start *= zscale; - } + if (istyle == SPHERE || istyle == CYLINDER) { + if (!xstr) xvalue *= xscale; + if (!ystr) yvalue *= yscale; + if (!zstr) zvalue *= zscale; + if (!rstr) rvalue *= xscale; } else if (istyle == PLANE) { - if (cdim == 0) planepos *= xscale; - else if (cdim == 1) planepos *= yscale; - else if (cdim == 2) planepos *= zscale; + if (cdim == 0 && !pstr) pvalue *= xscale; + else if (cdim == 1 && !pstr) pvalue *= yscale; + else if (cdim == 2 && !pstr) pvalue *= zscale; } else error->all("Illegal fix indent command"); indenter_flag = 0; @@ -108,6 +87,17 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ +FixIndent::~FixIndent() +{ + delete [] xstr; + delete [] ystr; + delete [] zstr; + delete [] rstr; + delete [] pstr; +} + +/* ---------------------------------------------------------------------- */ + int FixIndent::setmask() { int mask = 0; @@ -122,6 +112,37 @@ int FixIndent::setmask() void FixIndent::init() { + if (xstr) { + xvar = input->variable->find(xstr); + if (xvar < 0) error->all("Variable name for fix indent does not exist"); + if (!input->variable->equalstyle(xvar)) + error->all("Variable for fix indent is not equal style"); + } + if (ystr) { + yvar = input->variable->find(ystr); + if (yvar < 0) error->all("Variable name for fix indent does not exist"); + if (!input->variable->equalstyle(yvar)) + error->all("Variable for fix indent is not equal style"); + } + if (zstr) { + zvar = input->variable->find(zstr); + if (zvar < 0) error->all("Variable name for fix indent does not exist"); + if (!input->variable->equalstyle(zvar)) + error->all("Variable for fix indent is not equal style"); + } + if (rstr) { + rvar = input->variable->find(rstr); + if (rvar < 0) error->all("Variable name for fix indent does not exist"); + if (!input->variable->equalstyle(rvar)) + error->all("Variable for fix indent is not equal style"); + } + if (pstr) { + pvar = input->variable->find(pstr); + if (pvar < 0) error->all("Variable name for fix indent does not exist"); + if (!input->variable->equalstyle(pvar)) + error->all("Variable for fix indent is not equal style"); + } + if (strcmp(update->integrate_style,"respa") == 0) nlevels_respa = ((Respa *) update->integrate)->nlevels; } @@ -150,17 +171,6 @@ void FixIndent::min_setup(int vflag) void FixIndent::post_force(int vflag) { - // set current r0 - // for minimization, always set to r0_stop - - double r0; - if (!radflag || update->whichflag == 2) r0 = r0_stop; - else { - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - r0 = r0_start + delta * (r0_stop-r0_start); - } - // indenter values, 0 = energy, 1-3 = force components indenter_flag = 0; @@ -170,16 +180,22 @@ void FixIndent::post_force(int vflag) if (istyle == SPHERE) { - // ctr = current indenter center from original x0,y0,z0 + // ctr = current indenter center // remap into periodic box double ctr[3]; - double delta = (update->ntimestep - update->beginstep) * update->dt; - ctr[0] = x0 + delta*vx; - ctr[1] = y0 + delta*vy; - ctr[2] = z0 + delta*vz; + if (xstr) ctr[0] = input->variable->compute_equal(xvar); + else ctr[0] = xvalue; + if (ystr) ctr[1] = input->variable->compute_equal(yvar); + else ctr[1] = yvalue; + if (zstr) ctr[2] = input->variable->compute_equal(zvar); + else ctr[2] = zvalue; domain->remap(ctr); + double radius; + if (rstr) radius = input->variable->compute_equal(zvar); + else radius = rvalue; + double **x = atom->x; double **f = atom->f; int *mask = atom->mask; @@ -195,10 +211,10 @@ void FixIndent::post_force(int vflag) domain->minimum_image(delx,dely,delz); r = sqrt(delx*delx + dely*dely + delz*delz); if (side == OUTSIDE) { - dr = r - r0; + dr = r - radius; fmag = k*dr*dr; } else { - dr = r0 - r; + dr = radius - r; fmag = -k*dr*dr; } if (dr >= 0.0) continue; @@ -218,27 +234,36 @@ void FixIndent::post_force(int vflag) } else if (istyle == CYLINDER) { - // ctr = current indenter axis from original c1,c2 + // ctr = current indenter axis // remap into periodic box // 3rd coord is just near box for remap(), since isn't used double ctr[3]; - double delta = (update->ntimestep - update->beginstep) * update->dt; if (cdim == 0) { ctr[0] = domain->boxlo[0]; - ctr[1] = c1 + delta*vy; - ctr[2] = c2 + delta*vz; + if (ystr) ctr[1] = input->variable->compute_equal(yvar); + else ctr[1] = yvalue; + if (zstr) ctr[2] = input->variable->compute_equal(zvar); + else ctr[2] = zvalue; } else if (cdim == 1) { - ctr[0] = c1 + delta*vx; + if (xstr) ctr[0] = input->variable->compute_equal(xvar); + else ctr[0] = xvalue; ctr[1] = domain->boxlo[1]; - ctr[2] = c2 + delta*vz; + if (zstr) ctr[2] = input->variable->compute_equal(zvar); + else ctr[2] = zvalue; } else { - ctr[0] = c1 + delta*vx; - ctr[1] = c2 + delta*vy; + if (xstr) ctr[0] = input->variable->compute_equal(xvar); + else ctr[0] = xvalue; + if (ystr) ctr[1] = input->variable->compute_equal(yvar); + else ctr[1] = yvalue; ctr[2] = domain->boxlo[2]; } domain->remap(ctr); + double radius; + if (rstr) radius = input->variable->compute_equal(zvar); + else radius = rvalue; + double **x = atom->x; double **f = atom->f; int *mask = atom->mask; @@ -264,10 +289,10 @@ void FixIndent::post_force(int vflag) domain->minimum_image(delx,dely,delz); r = sqrt(delx*delx + dely*dely + delz*delz); if (side == OUTSIDE) { - dr = r - r0; + dr = r - radius; fmag = k*dr*dr; } else { - dr = r0 - r; + dr = radius - r; fmag = -k*dr*dr; } if (dr >= 0.0) continue; @@ -287,13 +312,11 @@ void FixIndent::post_force(int vflag) } else { - // posnew = current coord of plane from original planepos + // plane = current plane position - double delta = (update->ntimestep - update->beginstep) * update->dt; - double posnew; - if (cdim == 0) posnew = planepos + delta*vx; - else if (cdim == 1) posnew = planepos + delta*vy; - else if (cdim == 2) posnew = planepos + delta*vz; + double plane; + if (pstr) plane = input->variable->compute_equal(pvar); + else plane = pvalue; double **x = atom->x; double **f = atom->f; @@ -304,7 +327,7 @@ void FixIndent::post_force(int vflag) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - dr = planeside * (posnew - x[i][cdim]); + dr = planeside * (plane - x[i][cdim]); if (dr >= 0.0) continue; fatom = -planeside * k*dr*dr; f[i][cdim] += fatom; @@ -367,9 +390,8 @@ void FixIndent::options(int narg, char **arg) if (narg < 0) error->all("Illegal fix indent command"); istyle = NONE; - vx = vy = vz = 0.0; - radflag = 0; - r0_start = 0.0; + xstr = ystr = zstr = rstr = pstr = NULL; + xvalue = yvalue = zvalue = rvalue = pvalue = 0.0; scaleflag = 1; side = OUTSIDE; @@ -377,52 +399,107 @@ void FixIndent::options(int narg, char **arg) while (iarg < narg) { if (strcmp(arg[iarg],"sphere") == 0) { if (iarg+5 > narg) error->all("Illegal fix indent command"); - x0 = atof(arg[iarg+1]); - y0 = atof(arg[iarg+2]); - z0 = atof(arg[iarg+3]); - r0_stop = atof(arg[iarg+4]); + + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { + int n = strlen(&arg[iarg+1][2]) + 1; + xstr = new char[n]; + strcpy(xstr,&arg[iarg+1][2]); + } else xvalue = atof(arg[iarg+1]); + if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { + int n = strlen(&arg[iarg+2][2]) + 1; + ystr = new char[n]; + strcpy(ystr,&arg[iarg+2][2]); + } else yvalue = atof(arg[iarg+2]); + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { + int n = strlen(&arg[iarg+3][2]) + 1; + zstr = new char[n]; + strcpy(zstr,&arg[iarg+3][2]); + } else zvalue = atof(arg[iarg+3]); + if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) { + int n = strlen(&arg[iarg+4][2]) + 1; + rstr = new char[n]; + strcpy(rstr,&arg[iarg+4][2]); + } else rvalue = atof(arg[iarg+4]); + istyle = SPHERE; iarg += 5; + } else if (strcmp(arg[iarg],"cylinder") == 0) { if (iarg+5 > narg) error->all("Illegal fix indent command"); - if (strcmp(arg[iarg+1],"x") == 0) cdim = 0; - else if (strcmp(arg[iarg+1],"y") == 0) cdim = 1; - else if (strcmp(arg[iarg+1],"z") == 0) cdim = 2; - else error->all("Illegal fix indent command"); - c1 = atof(arg[iarg+2]); - c2 = atof(arg[iarg+3]); - r0_stop = atof(arg[iarg+4]); + + if (strcmp(arg[iarg+1],"x") == 0) { + cdim = 0; + if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { + int n = strlen(&arg[iarg+2][2]) + 1; + ystr = new char[n]; + strcpy(ystr,&arg[iarg+2][2]); + } else yvalue = atof(arg[iarg+2]); + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { + int n = strlen(&arg[iarg+3][2]) + 1; + zstr = new char[n]; + strcpy(zstr,&arg[iarg+3][2]); + } else zvalue = atof(arg[iarg+3]); + } else if (strcmp(arg[iarg+1],"y") == 0) { + cdim = 1; + if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { + int n = strlen(&arg[iarg+2][2]) + 1; + xstr = new char[n]; + strcpy(xstr,&arg[iarg+2][2]); + } else xvalue = atof(arg[iarg+2]); + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { + int n = strlen(&arg[iarg+3][2]) + 1; + zstr = new char[n]; + strcpy(zstr,&arg[iarg+3][2]); + } else zvalue = atof(arg[iarg+3]); + } else if (strcmp(arg[iarg+1],"z") == 0) { + cdim = 2; + if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { + int n = strlen(&arg[iarg+2][2]) + 1; + xstr = new char[n]; + strcpy(xstr,&arg[iarg+2][2]); + } else xvalue = atof(arg[iarg+2]); + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { + int n = strlen(&arg[iarg+3][2]) + 1; + ystr = new char[n]; + strcpy(ystr,&arg[iarg+3][2]); + } else yvalue = atof(arg[iarg+3]); + } else error->all("Illegal fix indent command"); + + if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) { + int n = strlen(&arg[iarg+4][2]) + 1; + rstr = new char[n]; + strcpy(rstr,&arg[iarg+4][2]); + } else rvalue = atof(arg[iarg+4]); + istyle = CYLINDER; iarg += 5; + } else if (strcmp(arg[iarg],"plane") == 0) { if (iarg+4 > narg) error->all("Illegal fix indent command"); if (strcmp(arg[iarg+1],"x") == 0) cdim = 0; else if (strcmp(arg[iarg+1],"y") == 0) cdim = 1; else if (strcmp(arg[iarg+1],"z") == 0) cdim = 2; else error->all("Illegal fix indent command"); - planepos = atof(arg[iarg+2]); + + if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { + int n = strlen(&arg[iarg+2][2]) + 1; + pstr = new char[n]; + strcpy(pstr,&arg[iarg+2][2]); + } else pvalue = atof(arg[iarg+2]); + if (strcmp(arg[iarg+3],"lo") == 0) planeside = -1; else if (strcmp(arg[iarg+3],"hi") == 0) planeside = 1; else error->all("Illegal fix indent command"); istyle = PLANE; iarg += 4; - } else if (strcmp(arg[iarg],"vel") == 0) { - if (iarg+4 > narg) error->all("Illegal fix indent command"); - vx = atof(arg[iarg+1]); - vy = atof(arg[iarg+2]); - vz = atof(arg[iarg+3]); - iarg += 4; - } else if (strcmp(arg[iarg],"rstart") == 0) { - if (iarg+2 > narg) error->all("Illegal fix indent command"); - radflag = 1; - r0_start = atof(arg[iarg+1]); - iarg += 2; + } else if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all("Illegal fix indent command"); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; else error->all("Illegal fix indent command"); iarg += 2; + } else if (strcmp(arg[iarg],"side") == 0) { if (iarg+2 > narg) error->all("Illegal fix indent command"); if (strcmp(arg[iarg+1],"in") == 0) side = INSIDE; diff --git a/src/fix_indent.h b/src/fix_indent.h index 937c8679b2..d9faf392f0 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -27,6 +27,7 @@ namespace LAMMPS_NS { class FixIndent : public Fix { public: FixIndent(class LAMMPS *, int, char **); + ~FixIndent(); int setmask(); void init(); void setup(int); @@ -38,14 +39,14 @@ class FixIndent : public Fix { double compute_vector(int); private: - int istyle,scaleflag,radflag,thermo_flag,eflag_enable,side; + int istyle,scaleflag,thermo_flag,eflag_enable,side; double k,k3; - double x0,y0,z0,r0_stop,r0_start,planepos; + char *xstr,*ystr,*zstr,*rstr,*pstr; + int xvar,yvar,zvar,rvar,pvar; + double xvalue,yvalue,zvalue,rvalue,pvalue; int indenter_flag,planeside; double indenter[4],indenter_all[4]; int cdim; - double c1,c2; - double vx,vy,vz; int nlevels_respa; void options(int, char **); diff --git a/src/thermo.cpp b/src/thermo.cpp index bedd8459f6..1e83376864 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -41,7 +41,8 @@ using namespace LAMMPS_NS; // customize a new keyword by adding to this list: -// step, atoms, cpu, dt, temp, press, pe, ke, etotal, enthalpy +// step, elapsed, elaplong, dt, cpu +// atoms, temp, press, pe, ke, etotal, enthalpy // evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail // vol, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz, xlat, ylat, zlat // pxx, pyy, pzz, pxy, pxz, pyz @@ -60,7 +61,7 @@ enum{SCALAR,VECTOR,ARRAY}; #define INVOKED_VECTOR 2 #define INVOKED_ARRAY 4 -#define MAXLINE 1024 +#define MAXLINE 8192 // make this 4x longer than Input::MAXLINE #define DELTA 8 #define MIN(A,B) ((A) < (B)) ? (A) : (B) @@ -583,13 +584,17 @@ void Thermo::parse_fields(char *str) if (strcmp(word,"step") == 0) { addfield("Step",&Thermo::compute_step,INT); - } else if (strcmp(word,"atoms") == 0) { - addfield("Atoms",&Thermo::compute_atoms,INT); - } else if (strcmp(word,"cpu") == 0) { - addfield("CPU",&Thermo::compute_cpu,FLOAT); + } else if (strcmp(word,"elapsed") == 0) { + addfield("Elapsed",&Thermo::compute_elapsed,INT); + } else if (strcmp(word,"elaplong") == 0) { + addfield("Elaplong",&Thermo::compute_elapsed_long,INT); } else if (strcmp(word,"dt") == 0) { addfield("Dt",&Thermo::compute_dt,FLOAT); + } else if (strcmp(word,"cpu") == 0) { + addfield("CPU",&Thermo::compute_cpu,FLOAT); + } else if (strcmp(word,"atoms") == 0) { + addfield("Atoms",&Thermo::compute_atoms,INT); } else if (strcmp(word,"temp") == 0) { addfield("Temp",&Thermo::compute_temp,FLOAT); index_temp = add_compute(id_temp,SCALAR); @@ -884,16 +889,28 @@ int Thermo::evaluate_keyword(char *word, double *answer) compute_step(); dvalue = ivalue; - } else if (strcmp(word,"atoms") == 0) { - compute_atoms(); + } else if (strcmp(word,"elapsed") == 0) { + if (update->whichflag == 0) + error->all("This variable thermo keyword cannot be used between runs"); + compute_elapsed(); dvalue = ivalue; + } else if (strcmp(word,"elaplong") == 0) { + if (update->whichflag == 0) + error->all("This variable thermo keyword cannot be used between runs"); + compute_elapsed_long(); + dvalue = ivalue; + + } else if (strcmp(word,"dt") == 0) { + compute_dt(); + } else if (strcmp(word,"cpu") == 0) { if (update->whichflag == 0) firststep = 0; compute_cpu(); - } else if (strcmp(word,"dt") == 0) { - compute_dt(); + } else if (strcmp(word,"atoms") == 0) { + compute_atoms(); + dvalue = ivalue; } else if (strcmp(word,"temp") == 0) { if (!temperature) @@ -1282,9 +1299,23 @@ void Thermo::compute_step() /* ---------------------------------------------------------------------- */ -void Thermo::compute_atoms() +void Thermo::compute_elapsed() { - ivalue = static_cast (natoms); + ivalue = update->ntimestep - update->firststep; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_elapsed_long() +{ + ivalue = update->ntimestep - update->beginstep; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_dt() +{ + dvalue = update->dt; } /* ---------------------------------------------------------------------- */ @@ -1297,9 +1328,9 @@ void Thermo::compute_cpu() /* ---------------------------------------------------------------------- */ -void Thermo::compute_dt() +void Thermo::compute_atoms() { - dvalue = update->dt; + ivalue = static_cast (natoms); } /* ---------------------------------------------------------------------- */ diff --git a/src/thermo.h b/src/thermo.h index cf3febee74..0f9d06732e 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -106,11 +106,13 @@ class Thermo : protected Pointers { // functions that compute a single value // customize a new keyword by adding a method prototype - void compute_step(); - void compute_atoms(); - void compute_cpu(); + void compute_step(); + void compute_elapsed(); + void compute_elapsed_long(); void compute_dt(); + void compute_cpu(); + void compute_atoms(); void compute_temp(); void compute_press(); void compute_pe();