These are compute styles contributed by users, which can be used if diff -Naur lammps-28Apr09/doc/Section_commands.txt lammps-30Apr09/doc/Section_commands.txt --- lammps-28Apr09/doc/Section_commands.txt 2009-04-28 08:56:27.000000000 -0600 +++ lammps-30Apr09/doc/Section_commands.txt 2009-04-29 10:54:14.000000000 -0600 @@ -474,6 +474,7 @@ "temp/com"_compute_temp_com.html, "temp/deform"_compute_temp_deform.html, "temp/partial"_compute_temp_partial.html, +"temp/profile"_compute_temp_profile.html, "temp/ramp"_compute_temp_ramp.html, "temp/region"_compute_temp_region.html, "temp/sphere"_compute_temp_sphere.html :tb(c=6,ea=c) diff -Naur lammps-28Apr09/doc/Section_howto.html lammps-30Apr09/doc/Section_howto.html --- lammps-28Apr09/doc/Section_howto.html 2009-01-26 13:14:19.000000000 -0700 +++ lammps-30Apr09/doc/Section_howto.html 2009-04-29 10:54:14.000000000 -0600 @@ -1178,6 +1178,7 @@
Thermostatting in LAMMPS is performed by fixes. Four thermostatting fixes are currently available: Nose-Hoover (nvt), @@ -1203,7 +1205,7 @@
Fix nvt only thermostats the translational velocity of
-particles. Fix nvt/sllod does as well, except
+particles. Fix nvt/sllod also does this, except
that it subtracts out a velocity bias due to a deforming box and
integrates the SLLOD equations of motion. See the NEMD
simulations section of this page for further details. Fix
@@ -1225,8 +1227,9 @@
IMPORTANT NOTE: Only the nvt fixes perform time integration, meaning
they update the velocities and positions of particles due to forces
and velocities respectively. The other thermostat fixes only adjust
-velocities; they do NOT perform time integration. Thus they should be
-used in conjunction with a constant NVE integration fix such as these:
+velocities; they do NOT perform time integration updates. Thus they
+should be used in conjunction with a constant NVE integration fix such
+as these:
All of the barostatting fixes use the compute pressure compute to calculate a current -pressure. The barostatting fixes can also use temperature computes -that remove bias for the purpose of computing the current temperature -which contributes to the current pressure. See the doc pages for the -individual fixes and for the fix_modify command for -instructions on how to assign a temperature or pressure compute to a -barostatting fix. +pressure. By default, this compute is created with a simple compute +temp (see the last argument of the compute +pressure command), which is used to calculated +the kinetic componenet of the pressure. The barostatting fixes can +also use temperature computes that remove bias for the purpose of +computing the kinetic componenet which contributes to the current +pressure. See the doc pages for the individual fixes and for the +fix_modify command for instructions on how to assign +a temperature or pressure compute to a barostatting fix.
IMPORTANT NOTE: As with the thermostats, the Nose/Hoover methods (fix npt and fix nph) perform time diff -Naur lammps-28Apr09/doc/Section_howto.txt lammps-30Apr09/doc/Section_howto.txt --- lammps-28Apr09/doc/Section_howto.txt 2009-01-26 13:14:19.000000000 -0700 +++ lammps-30Apr09/doc/Section_howto.txt 2009-04-29 10:54:14.000000000 -0600 @@ -1169,6 +1169,7 @@ "compute temp/com"_compute_temp_com.html "compute temp/deform"_compute_temp_deform.html "compute temp/partial"_compute_temp_partial.html +"compute temp/profile"_compute_temp_profile.html "compute temp/ramp"_compute_temp_ramp.html "compute temp/region"_compute_temp_region.html :ul @@ -1177,9 +1178,10 @@ "Fix temp/sphere"_fix_temp_sphere.html and "fix temp/asphere"_fix_temp_asphere.html compute kinetic energy for extended particles that includes rotational degrees of freedom. They -both allow, as an extra argument, another temperature compute that -subtracts a velocity bias, so the translational velocity of extended -spherical or aspherical particles can be adjusted in prescribed ways. +both allow, as an extra argument, which is another temperature compute +that subtracts a velocity bias. This allows the translational +velocity of extended spherical or aspherical particles to be adjusted +in prescribed ways. Thermostatting in LAMMPS is performed by "fixes"_fix.html. Four thermostatting fixes are currently available: Nose-Hoover (nvt), @@ -1194,7 +1196,7 @@ "fix temp/rescale"_fix_temp_rescale.html :ul "Fix nvt"_fix_nvt.html only thermostats the translational velocity of -particles. "Fix nvt/sllod"_fix_nvt_sllod.html does as well, except +particles. "Fix nvt/sllod"_fix_nvt_sllod.html also does this, except that it subtracts out a velocity bias due to a deforming box and integrates the SLLOD equations of motion. See the "NEMD simulations"_#4_13 section of this page for further details. "Fix @@ -1216,8 +1218,9 @@ IMPORTANT NOTE: Only the nvt fixes perform time integration, meaning they update the velocities and positions of particles due to forces and velocities respectively. The other thermostat fixes only adjust -velocities; they do NOT perform time integration. Thus they should be -used in conjunction with a constant NVE integration fix such as these: +velocities; they do NOT perform time integration updates. Thus they +should be used in conjunction with a constant NVE integration fix such +as these: "fix nve"_fix_nve.html "fix nve/sphere"_fix_nve_sphere.html @@ -1249,12 +1252,15 @@ All of the barostatting fixes use the "compute pressure"_compute_pressure.html compute to calculate a current -pressure. The barostatting fixes can also use temperature computes -that remove bias for the purpose of computing the current temperature -which contributes to the current pressure. See the doc pages for the -individual fixes and for the "fix_modify"_fix_modify.html command for -instructions on how to assign a temperature or pressure compute to a -barostatting fix. +pressure. By default, this compute is created with a simple "compute +temp"_compute_temp.html (see the last argument of the "compute +pressure"_compute_pressure.html command), which is used to calculated +the kinetic componenet of the pressure. The barostatting fixes can +also use temperature computes that remove bias for the purpose of +computing the kinetic componenet which contributes to the current +pressure. See the doc pages for the individual fixes and for the +"fix_modify"_fix_modify.html command for instructions on how to assign +a temperature or pressure compute to a barostatting fix. IMPORTANT NOTE: As with the thermostats, the Nose/Hoover methods ("fix npt"_fix_npt.html and "fix nph"_fix_nph.html) perform time diff -Naur lammps-28Apr09/doc/compute.html lammps-30Apr09/doc/compute.html --- lammps-28Apr09/doc/compute.html 2009-04-28 08:56:27.000000000 -0600 +++ lammps-30Apr09/doc/compute.html 2009-04-29 10:54:14.000000000 -0600 @@ -130,6 +130,7 @@
The rotational kinetic energy is computed as 1/2 I w^2, where I is the -inertia tensor for the aspherical particle and w is its angular -velocity, which is computed from its angular momentum. +
The translational kinetic energy is computed the same as is described +by the compute temp command. The rotational +kinetic energy is computed as 1/2 I w^2, where I is the inertia tensor +for the aspherical particle and w is its angular velocity, which is +computed from its angular momentum.
IMPORTANT NOTE: For 2d models, particles are treated as ellipsoids, not ellipses, meaning their moments of inertia will be diff -Naur lammps-28Apr09/doc/compute_temp_asphere.txt lammps-30Apr09/doc/compute_temp_asphere.txt --- lammps-28Apr09/doc/compute_temp_asphere.txt 2008-04-07 09:58:34.000000000 -0600 +++ lammps-30Apr09/doc/compute_temp_asphere.txt 2009-04-29 10:54:14.000000000 -0600 @@ -47,9 +47,11 @@ pair potential"_pair_gayberne.html does not impart torque to spherical particles, so they do not rotate. -The rotational kinetic energy is computed as 1/2 I w^2, where I is the -inertia tensor for the aspherical particle and w is its angular -velocity, which is computed from its angular momentum. +The translational kinetic energy is computed the same as is described +by the "compute temp"_compute_temp.html command. The rotational +kinetic energy is computed as 1/2 I w^2, where I is the inertia tensor +for the aspherical particle and w is its angular velocity, which is +computed from its angular momentum. IMPORTANT NOTE: For "2d models"_dimension.html, particles are treated as ellipsoids, not ellipses, meaning their moments of inertia will be diff -Naur lammps-28Apr09/doc/compute_temp_com.html lammps-30Apr09/doc/compute_temp_com.html --- lammps-28Apr09/doc/compute_temp_com.html 2008-04-04 15:15:50.000000000 -0600 +++ lammps-30Apr09/doc/compute_temp_com.html 2009-04-29 10:54:14.000000000 -0600 @@ -33,7 +33,8 @@ e.g. thermo_modify, fix temp/rescale, fix npt, etc.
-The temperature is calculated by the formula KE = dim/2 N k T, where +
After the center-of-mass velocity has been subtracted from each atom, +the temperature is calculated by the formula KE = dim/2 N k T, where KE = total kinetic energy of the group of atoms (sum of 1/2 m v^2), dim = 2 or 3 = dimensionality of the simulation, N = number of atoms in the group, k = Boltzmann constant, and T = temperature. diff -Naur lammps-28Apr09/doc/compute_temp_com.txt lammps-30Apr09/doc/compute_temp_com.txt --- lammps-28Apr09/doc/compute_temp_com.txt 2008-04-04 15:15:50.000000000 -0600 +++ lammps-30Apr09/doc/compute_temp_com.txt 2009-04-29 10:54:14.000000000 -0600 @@ -30,7 +30,8 @@ e.g. "thermo_modify"_thermo_modify.html, "fix temp/rescale"_fix_temp_rescale.html, "fix npt"_fix_npt.html, etc. -The temperature is calculated by the formula KE = dim/2 N k T, where +After the center-of-mass velocity has been subtracted from each atom, +the temperature is calculated by the formula KE = dim/2 N k T, where KE = total kinetic energy of the group of atoms (sum of 1/2 m v^2), dim = 2 or 3 = dimensionality of the simulation, N = number of atoms in the group, k = Boltzmann constant, and T = temperature. diff -Naur lammps-28Apr09/doc/compute_temp_deform.html lammps-30Apr09/doc/compute_temp_deform.html --- lammps-28Apr09/doc/compute_temp_deform.html 2008-04-22 08:54:46.000000000 -0600 +++ lammps-30Apr09/doc/compute_temp_deform.html 2009-04-29 10:54:14.000000000 -0600 @@ -36,10 +36,10 @@ fix temp/rescale, fix npt, etc.
The deformation fix changes the box size and/or shape over time, so -each point in the simulation box can be thought of as having a +each atom in the simulation box can be thought of as having a "streaming" velocity. For example, if the box is being sheared in x, -relative to y, then points at the bottom of the box (low y) have a -small x velocity, while points at the top of the box (hi y) have a +relative to y, then atoms at the bottom of the box (low y) have a +small x velocity, while atoms at the top of the box (hi y) have a large x velocity. This position-dependent streaming velocity is subtracted from each atom's actual velocity to yield a thermal velocity which is used to compute the temperature. @@ -56,11 +56,12 @@ by this compute. LAMMPS will warn you if fix deform is defined and its remap setting is not consistent with this compute.
-The temperature is calculated by the formula KE = dim/2 N k T, where -KE = total kinetic energy of the group of atoms (sum of 1/2 m v^2), -dim = 2 or 3 = dimensionality of the simulation, N = number of atoms -in the group, k = Boltzmann constant, and T = temperature. Note that -v in the kinetic energy formula is the atom's thermal velocity. +
After the streaming velocity has been subtracted from each atom, the +temperature is calculated by the formula KE = dim/2 N k T, where KE = +total kinetic energy of the group of atoms (sum of 1/2 m v^2), dim = 2 +or 3 = dimensionality of the simulation, N = number of atoms in the +group, k = Boltzmann constant, and T = temperature. Note that v in +the kinetic energy formula is the atom's thermal velocity.
A 6-component kinetic energy tensor is also calculated by this compute for use in the computation of a pressure tensor. The formula for the @@ -104,8 +105,9 @@
Related commands:
-compute temp/ramp, fix -deform, fix nvt/sllod +
compute temp/ramp, compute +temp/profile, fix deform, +fix nvt/sllod
Default: none
diff -Naur lammps-28Apr09/doc/compute_temp_deform.txt lammps-30Apr09/doc/compute_temp_deform.txt --- lammps-28Apr09/doc/compute_temp_deform.txt 2008-04-22 08:54:46.000000000 -0600 +++ lammps-30Apr09/doc/compute_temp_deform.txt 2009-04-29 10:54:14.000000000 -0600 @@ -33,10 +33,10 @@ "fix temp/rescale"_fix_temp_rescale.html, "fix npt"_fix_npt.html, etc. The deformation fix changes the box size and/or shape over time, so -each point in the simulation box can be thought of as having a +each atom in the simulation box can be thought of as having a "streaming" velocity. For example, if the box is being sheared in x, -relative to y, then points at the bottom of the box (low y) have a -small x velocity, while points at the top of the box (hi y) have a +relative to y, then atoms at the bottom of the box (low y) have a +small x velocity, while atoms at the top of the box (hi y) have a large x velocity. This position-dependent streaming velocity is subtracted from each atom's actual velocity to yield a thermal velocity which is used to compute the temperature. @@ -53,11 +53,12 @@ by this compute. LAMMPS will warn you if fix deform is defined and its remap setting is not consistent with this compute. -The temperature is calculated by the formula KE = dim/2 N k T, where -KE = total kinetic energy of the group of atoms (sum of 1/2 m v^2), -dim = 2 or 3 = dimensionality of the simulation, N = number of atoms -in the group, k = Boltzmann constant, and T = temperature. Note that -v in the kinetic energy formula is the atom's thermal velocity. +After the streaming velocity has been subtracted from each atom, the +temperature is calculated by the formula KE = dim/2 N k T, where KE = +total kinetic energy of the group of atoms (sum of 1/2 m v^2), dim = 2 +or 3 = dimensionality of the simulation, N = number of atoms in the +group, k = Boltzmann constant, and T = temperature. Note that v in +the kinetic energy formula is the atom's thermal velocity. A 6-component kinetic energy tensor is also calculated by this compute for use in the computation of a pressure tensor. The formula for the @@ -101,7 +102,8 @@ [Related commands:] -"compute temp/ramp"_compute_temp_ramp.html, "fix -deform"_fix_deform.html, "fix nvt/sllod"_fix_nvt_sllod.html +"compute temp/ramp"_compute_temp_ramp.html, "compute +temp/profile"_compute_temp_profile.html, "fix deform"_fix_deform.html, +"fix nvt/sllod"_fix_nvt_sllod.html [Default:] none diff -Naur lammps-28Apr09/doc/compute_temp_profile.html lammps-30Apr09/doc/compute_temp_profile.html --- lammps-28Apr09/doc/compute_temp_profile.html 1969-12-31 17:00:00.000000000 -0700 +++ lammps-30Apr09/doc/compute_temp_profile.html 2009-04-29 11:29:52.000000000 -0600 @@ -0,0 +1,144 @@ + +Syntax: +
+compute ID group-ID temp/profile xflag yflag zflag binstyle args ++
x arg = Nx + y arg = Ny + z arg = Nz + xy args = Nx Ny + yz args = Ny Nz + xz args = Nx Nz + xyz args = Nx Ny Nz + Nx,Ny,Nz = number of velocity bins in x,y,z dimensions ++ +
Examples: +
+compute myTemp flow temp/profile 1 1 1 x 10 +compute myTemp flow temp/profile 0 1 1 xyz 20 20 20 ++
Description: +
+Define a computation that calculates the temperature of a group of +atoms, after subtracting out a spatially-averaged velocity field, +before computing the kinetic energy. This can be useful for +thermostatting a collection of atoms undergoing a complex flow, +e.g. via a profile-unbiased thermostat (PUT) as described in +(Evans). A compute of this style can be used by any command +that computes a temperature, e.g. thermo_modify, +fix temp/rescale, fix npt, etc. +
+The xflag, yflag, zflag settings determine which components of +velocity are subtracted out. +
+The binstyle setting and its Nx, Ny, Nz arguments determine +how bins are setup to perform spatial averaging. "Bins" can be 1d +slabs, 2d pencils, or 3d bricks depending on which binstyle is used. +The simulation box is partitioned conceptually into Nx by Ny by +Nz bins. Depending on the binstyle, you may only specify one or +two of these values; the others are effectively set to 1 (no binning +in that dimension). For non-orthogonal (triclinic) simulation boxes, +the bins are "tilted" slabs or pencils or bricks that are parallel to +the tilted faces of the box. See the region prism +command for a discussion of the geometry of tilted boxes in LAMMPS. +
+When a temperature is computed, the velocity for the set of atoms that +are both in the compute group and in the same spatial bin is summed to +compute an average velocity for the bin. This bias velocity is then +subtracted from the velocities of individual atoms in the bin to yield +a thermal velocity. +
+After the spatially-averaged velocity field has been subtracted from +each atom, the temperature is calculated by the formula KE = dim/2 N k +T, where KE = total kinetic energy of the group of atoms (sum of 1/2 m +v^2), dim = 2 or 3 = dimensionality of the simulation, N = number of +atoms in the group, k = Boltzmann constant, and T = temperature. +
+A 6-component kinetic energy tensor is also calculated by this compute +for use in the calculation of a pressure tensor. The formula for the +components of the tensor is the same as the above formula, except that +v^2 is replaced by vx * vy for the xy component, etc. +
+The number of atoms contributing to the temperature is assumed to be +constant for the duration of the run; use the dynamic option of the +compute_modify command if this is not the case. +
+The removal of the spatially-averaged velocity field by this fix is +essentially computing the temperature after a "bias" has been removed +from the velocity of the atoms. If this compute is used with a fix +command that performs thermostatting then this bias will be subtracted +from each atom, thermostatting of the remaining thermal velocity will +be performed, and the bias will be added back in. Thermostatting +fixes that work in this way include fix nvt, fix +temp/rescale, fix +temp/berendsen, and fix +langevin. +
+This compute subtracts out degrees-of-freedom due to fixes that +constrain molecular motion, such as fix shake and +fix rigid. This means the temperature of groups of +atoms that include these constraints will be computed correctly. If +needed, the subtracted degrees-of-freedom can be altered using the +extra option of the compute_modify command. +
+See this howto section of the manual for a +discussion of different ways to compute temperature and perform +thermostatting. Using this compute in conjunction with a +thermostatting fix will effectively implement a PUT, as described in +(Evans). +
+Output info: +
+The scalar value calculated by this compute is "intensive", meaning it +is independent of the number of atoms in the simulation. The vector +values are "extensive", meaning they scale with the number of atoms in +the simulation. +
+Restrictions: +
+You should not use too large a velocity-binning grid, especially in +3d. In the current implementation, the binned velocity averages are +summed across all processors, so this will be inefficient if the grid +is too large, and the operation is performed every timestep, as it +will be for most thermostats. +
+Related commands: +
+compute temp, compute +temp/ramp, compute +temp/deform, compute +pressure +
+Default: +
+The option default is units = lattice. +
+(Evans) Evans and Morriss, Phys Rev Lett, 56, 2172-2175 (1986). +
+ diff -Naur lammps-28Apr09/doc/compute_temp_profile.txt lammps-30Apr09/doc/compute_temp_profile.txt --- lammps-28Apr09/doc/compute_temp_profile.txt 1969-12-31 17:00:00.000000000 -0700 +++ lammps-30Apr09/doc/compute_temp_profile.txt 2009-04-29 11:29:52.000000000 -0600 @@ -0,0 +1,133 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +compute temp/profile command :h3 + +[Syntax:] + +compute ID group-ID temp/profile xflag yflag zflag binstyle args :pre + +ID, group-ID are documented in "compute"_compute.html command :ulb,l +temp/profile = style name of this compute command :l +xflag,yflag,zflag = 0/1 for whether to exclude/include this dimension :l +binstyle = {x} or {y} or {z} or {xy} or {yz} or {xz} or {xyz} :l + {x} arg = Nx + {y} arg = Ny + {z} arg = Nz + {xy} args = Nx Ny + {yz} args = Ny Nz + {xz} args = Nx Nz + {xyz} args = Nx Ny Nz + Nx,Ny,Nz = number of velocity bins in x,y,z dimensions :pre +:ule + +[Examples:] + +compute myTemp flow temp/profile 1 1 1 x 10 +compute myTemp flow temp/profile 0 1 1 xyz 20 20 20 :pre + +[Description:] + +Define a computation that calculates the temperature of a group of +atoms, after subtracting out a spatially-averaged velocity field, +before computing the kinetic energy. This can be useful for +thermostatting a collection of atoms undergoing a complex flow, +e.g. via a profile-unbiased thermostat (PUT) as described in +"(Evans)"_#Evans. A compute of this style can be used by any command +that computes a temperature, e.g. "thermo_modify"_thermo_modify.html, +"fix temp/rescale"_fix_temp_rescale.html, "fix npt"_fix_npt.html, etc. + +The {xflag}, {yflag}, {zflag} settings determine which components of +velocity are subtracted out. + +The {binstyle} setting and its {Nx}, {Ny}, {Nz} arguments determine +how bins are setup to perform spatial averaging. "Bins" can be 1d +slabs, 2d pencils, or 3d bricks depending on which {binstyle} is used. +The simulation box is partitioned conceptually into {Nx} by {Ny} by +{Nz} bins. Depending on the {binstyle}, you may only specify one or +two of these values; the others are effectively set to 1 (no binning +in that dimension). For non-orthogonal (triclinic) simulation boxes, +the bins are "tilted" slabs or pencils or bricks that are parallel to +the tilted faces of the box. See the "region prism"_region.html +command for a discussion of the geometry of tilted boxes in LAMMPS. + +When a temperature is computed, the velocity for the set of atoms that +are both in the compute group and in the same spatial bin is summed to +compute an average velocity for the bin. This bias velocity is then +subtracted from the velocities of individual atoms in the bin to yield +a thermal velocity. + +After the spatially-averaged velocity field has been subtracted from +each atom, the temperature is calculated by the formula KE = dim/2 N k +T, where KE = total kinetic energy of the group of atoms (sum of 1/2 m +v^2), dim = 2 or 3 = dimensionality of the simulation, N = number of +atoms in the group, k = Boltzmann constant, and T = temperature. + +A 6-component kinetic energy tensor is also calculated by this compute +for use in the calculation of a pressure tensor. The formula for the +components of the tensor is the same as the above formula, except that +v^2 is replaced by vx * vy for the xy component, etc. + +The number of atoms contributing to the temperature is assumed to be +constant for the duration of the run; use the {dynamic} option of the +"compute_modify"_compute_modify.html command if this is not the case. + +The removal of the spatially-averaged velocity field by this fix is +essentially computing the temperature after a "bias" has been removed +from the velocity of the atoms. If this compute is used with a fix +command that performs thermostatting then this bias will be subtracted +from each atom, thermostatting of the remaining thermal velocity will +be performed, and the bias will be added back in. Thermostatting +fixes that work in this way include "fix nvt"_fix_nvt.html, "fix +temp/rescale"_fix_temp_rescale.html, "fix +temp/berendsen"_fix_temp_berendsen, and "fix +langevin"_fix_langevin.html. + +This compute subtracts out degrees-of-freedom due to fixes that +constrain molecular motion, such as "fix shake"_fix_shake.html and +"fix rigid"_fix_rigid.html. This means the temperature of groups of +atoms that include these constraints will be computed correctly. If +needed, the subtracted degrees-of-freedom can be altered using the +{extra} option of the "compute_modify"_compute_modify.html command. + +See "this howto section"_Section_howto.html#4_16 of the manual for a +discussion of different ways to compute temperature and perform +thermostatting. Using this compute in conjunction with a +thermostatting fix will effectively implement a PUT, as described in +"(Evans)"_#Evans. + +[Output info:] + +The scalar value calculated by this compute is "intensive", meaning it +is independent of the number of atoms in the simulation. The vector +values are "extensive", meaning they scale with the number of atoms in +the simulation. + +[Restrictions:] + +You should not use too large a velocity-binning grid, especially in +3d. In the current implementation, the binned velocity averages are +summed across all processors, so this will be inefficient if the grid +is too large, and the operation is performed every timestep, as it +will be for most thermostats. + +[Related commands:] + +"compute temp"_compute_temp.html, "compute +temp/ramp"_compute_temp_ramp.html, "compute +temp/deform"_compute_temp_deform.html, "compute +pressure"_compute_pressure.html + +[Default:] + +The option default is units = lattice. + +:line + +:link(Evans) +[(Evans)] Evans and Morriss, Phys Rev Lett, 56, 2172-2175 (1986). diff -Naur lammps-28Apr09/doc/compute_temp_ramp.html lammps-30Apr09/doc/compute_temp_ramp.html --- lammps-28Apr09/doc/compute_temp_ramp.html 2008-04-04 15:15:50.000000000 -0600 +++ lammps-30Apr09/doc/compute_temp_ramp.html 2009-04-29 10:54:14.000000000 -0600 @@ -28,20 +28,27 @@Examples:
-temperature 2nd middle ramp vx 0 8 y 2 12 units lattice +compute 2nd middle temp/ramp vx 0 8 y 2 12 units latticeDescription:
Define a computation that calculates the temperature of a group of -atoms, after subtracting out an imposed velocity on the system before +atoms, after subtracting out an ramped velocity profile before computing the kinetic energy. A compute of this style can be used by any command that computes a temperature, e.g. thermo_modify, fix temp/rescale, fix npt, etc.
-The meaning of the arguments for this command is the same as for the -velocity ramp command which was presumably used to -impose the velocity. +
The meaning of the arguments for this command which define the +velocity ramp are the same as for the velocity ramp +command which was presumably used to impose the velocity. +
+After the ramp velocity has been subtracted from the specified +dimension for each atom, the temperature is calculated by the formula +KE = dim/2 N k T, where KE = total kinetic energy of the group of +atoms (sum of 1/2 m v^2), dim = 2 or 3 = dimensionality of the +simulation, N = number of atoms in the group, k = Boltzmann constant, +and T = temperature.
The units keyword determines the meaning of the distance units used for coordinates (c1,c2) and velocities (vlo,vhi). A box value @@ -93,8 +100,8 @@
Related commands:
-compute temp, compute -temp/region, compute +
-compute temp, compute +temp/profie, compute temp/deform, compute pressure
diff -Naur lammps-28Apr09/doc/compute_temp_ramp.txt lammps-30Apr09/doc/compute_temp_ramp.txt --- lammps-28Apr09/doc/compute_temp_ramp.txt 2008-04-04 15:15:50.000000000 -0600 +++ lammps-30Apr09/doc/compute_temp_ramp.txt 2009-04-29 10:54:14.000000000 -0600 @@ -24,20 +24,27 @@ [Examples:] -temperature 2nd middle ramp vx 0 8 y 2 12 units lattice :pre +compute 2nd middle temp/ramp vx 0 8 y 2 12 units lattice :pre [Description:] Define a computation that calculates the temperature of a group of -atoms, after subtracting out an imposed velocity on the system before +atoms, after subtracting out an ramped velocity profile before computing the kinetic energy. A compute of this style can be used by any command that computes a temperature, e.g. "thermo_modify"_thermo_modify.html, "fix temp/rescale"_fix_temp_rescale.html, "fix npt"_fix_npt.html, etc. -The meaning of the arguments for this command is the same as for the -"velocity ramp"_velocity.html command which was presumably used to -impose the velocity. +The meaning of the arguments for this command which define the +velocity ramp are the same as for the "velocity ramp"_velocity.html +command which was presumably used to impose the velocity. + +After the ramp velocity has been subtracted from the specified +dimension for each atom, the temperature is calculated by the formula +KE = dim/2 N k T, where KE = total kinetic energy of the group of +atoms (sum of 1/2 m v^2), dim = 2 or 3 = dimensionality of the +simulation, N = number of atoms in the group, k = Boltzmann constant, +and T = temperature. The {units} keyword determines the meaning of the distance units used for coordinates (c1,c2) and velocities (vlo,vhi). A {box} value @@ -90,11 +97,10 @@ [Related commands:] "compute temp"_compute_temp.html, "compute -temp/region"_compute_temp_region.html, "compute +temp/profie"_compute_temp_profile.html, "compute temp/deform"_compute_temp_deform.html, "compute pressure"_compute_pressure.html [Default:] The option default is units = lattice. - diff -Naur lammps-28Apr09/doc/compute_temp_sphere.html lammps-30Apr09/doc/compute_temp_sphere.html --- lammps-28Apr09/doc/compute_temp_sphere.html 2008-04-04 15:15:50.000000000 -0600 +++ lammps-30Apr09/doc/compute_temp_sphere.html 2009-04-29 10:54:14.000000000 -0600 @@ -36,9 +36,10 @@ translational, 3 rotational). For 2d spherical particles, each has 3 degrees of freedom (2 translational, 1 rotational).The rotational kinetic energy is computed as 1/2 I w^2, where I is the -moment of inertia for a sphere and w is the particle's angular -velocity. +
The translational kinetic energy is computed the same as is described +by the compute temp command. The rotational +kinetic energy is computed as 1/2 I w^2, where I is the moment of +inertia for a sphere and w is the particle's angular velocity.
IMPORTANT NOTE: For 2d models, particles are treated as spheres, not disks, meaning their moment of inertia will be the diff -Naur lammps-28Apr09/doc/compute_temp_sphere.txt lammps-30Apr09/doc/compute_temp_sphere.txt --- lammps-28Apr09/doc/compute_temp_sphere.txt 2008-04-04 15:15:50.000000000 -0600 +++ lammps-30Apr09/doc/compute_temp_sphere.txt 2009-04-29 10:54:14.000000000 -0600 @@ -33,9 +33,10 @@ translational, 3 rotational). For 2d spherical particles, each has 3 degrees of freedom (2 translational, 1 rotational). -The rotational kinetic energy is computed as 1/2 I w^2, where I is the -moment of inertia for a sphere and w is the particle's angular -velocity. +The translational kinetic energy is computed the same as is described +by the "compute temp"_compute_temp.html command. The rotational +kinetic energy is computed as 1/2 I w^2, where I is the moment of +inertia for a sphere and w is the particle's angular velocity. IMPORTANT NOTE: For "2d models"_dimension.html, particles are treated as spheres, not disks, meaning their moment of inertia will be the diff -Naur lammps-28Apr09/doc/fix_langevin.html lammps-30Apr09/doc/fix_langevin.html --- lammps-28Apr09/doc/fix_langevin.html 2008-12-04 09:15:03.000000000 -0700 +++ lammps-30Apr09/doc/fix_langevin.html 2009-04-29 10:54:14.000000000 -0600 @@ -97,19 +97,17 @@ run from Tstart to Tstop.
Like other fixes that perform thermostatting, this fix can be used -with compute commands that calculate a temperature -after removing a "bias" from the atom velocities. E.g. removing the -center-of-mass velocity from a group of atoms or only calculating -temperature on the x-component of velocity or only calculating -temperature for atoms in a geometric region. This is not done by -default, but only if the fix_modify command is used -to assign a temperature compute to this fix that includes such a bias -term. See the doc pages for individual compute -commands to determine which ones include a bias. In -this case, the thermostat works in the following manner: the current -temperature is calculated taking the bias into account, bias is -removed from each atom, thermostatting is performed on the remaining -thermal degrees of freedom, and the bias is added back in. +with compute commands that remove a "bias" from the +atom velocities. E.g. removing the center-of-mass velocity from a +group of atoms or removing the x-component of velocity from the +calculation. This is not done by default, but only if the +fix_modify command is used to assign a temperature +compute to this fix that includes such a bias term. See the doc pages +for individual compute commands to determine which ones +include a bias. In this case, the thermostat works in the following +manner: bias is removed from each atom, thermostatting is performed on +the remaining thermal degrees of freedom, and the bias is added back +in.
The damp parameter is specified in time units and determines how rapidly the temperature is relaxed. For example, a value of 100.0 diff -Naur lammps-28Apr09/doc/fix_langevin.txt lammps-30Apr09/doc/fix_langevin.txt --- lammps-28Apr09/doc/fix_langevin.txt 2008-12-04 09:15:03.000000000 -0700 +++ lammps-30Apr09/doc/fix_langevin.txt 2009-04-29 10:54:14.000000000 -0600 @@ -88,19 +88,17 @@ run from {Tstart} to {Tstop}. Like other fixes that perform thermostatting, this fix can be used -with "compute commands"_compute.html that calculate a temperature -after removing a "bias" from the atom velocities. E.g. removing the -center-of-mass velocity from a group of atoms or only calculating -temperature on the x-component of velocity or only calculating -temperature for atoms in a geometric region. This is not done by -default, but only if the "fix_modify"_fix_modify.html command is used -to assign a temperature compute to this fix that includes such a bias -term. See the doc pages for individual "compute -commands"_compute.html to determine which ones include a bias. In -this case, the thermostat works in the following manner: the current -temperature is calculated taking the bias into account, bias is -removed from each atom, thermostatting is performed on the remaining -thermal degrees of freedom, and the bias is added back in. +with "compute commands"_compute.html that remove a "bias" from the +atom velocities. E.g. removing the center-of-mass velocity from a +group of atoms or removing the x-component of velocity from the +calculation. This is not done by default, but only if the +"fix_modify"_fix_modify.html command is used to assign a temperature +compute to this fix that includes such a bias term. See the doc pages +for individual "compute commands"_compute.html to determine which ones +include a bias. In this case, the thermostat works in the following +manner: bias is removed from each atom, thermostatting is performed on +the remaining thermal degrees of freedom, and the bias is added back +in. The {damp} parameter is specified in time units and determines how rapidly the temperature is relaxed. For example, a value of 100.0 diff -Naur lammps-28Apr09/src/compute.h lammps-30Apr09/src/compute.h --- lammps-28Apr09/src/compute.h 2009-01-05 15:26:08.000000000 -0700 +++ lammps-30Apr09/src/compute.h 2009-04-29 10:54:06.000000000 -0600 @@ -96,6 +96,10 @@ int extra_dof; // extra DOF for temperature computes int dynamic; // recount atoms for temperature computes int thermoflag; // 1 if include fix PE for PE computes + + double vbias[3]; // stored velocity bias for one atom + double **vbiasall; // stored velocity bias for all atoms + int maxbias; // size of vbiasall array }; } diff -Naur lammps-28Apr09/src/compute_pressure.cpp lammps-30Apr09/src/compute_pressure.cpp --- lammps-28Apr09/src/compute_pressure.cpp 2009-03-05 13:55:21.000000000 -0700 +++ lammps-30Apr09/src/compute_pressure.cpp 2009-04-29 10:54:06.000000000 -0600 @@ -110,8 +110,8 @@ nktv2p = force->nktv2p; dimension = domain->dimension; - // set temperature used by pressure - // must be done in init() since user can change it via modify command + // set temperature compute, must be done in init() + // fixes could have changed or compute_modify could have changed it int icompute = modify->find_compute(id_pre); if (icompute < 0) error->all("Could not find compute pressure temp ID"); diff -Naur lammps-28Apr09/src/compute_temp_partial.cpp lammps-30Apr09/src/compute_temp_partial.cpp --- lammps-28Apr09/src/compute_temp_partial.cpp 2009-02-12 14:36:52.000000000 -0700 +++ lammps-30Apr09/src/compute_temp_partial.cpp 2009-04-29 10:54:06.000000000 -0600 @@ -33,10 +33,6 @@ { if (narg != 6) error->all("Illegal compute temp/partial command"); - xflag = atoi(arg[3]); - yflag = atoi(arg[4]); - zflag = atoi(arg[5]); - scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; @@ -44,6 +40,10 @@ tempflag = 1; tempbias = 1; + xflag = atoi(arg[3]); + yflag = atoi(arg[4]); + zflag = atoi(arg[5]); + maxbias = 0; vbiasall = NULL; vector = new double[6]; @@ -193,21 +193,27 @@ "compute/temp:vbiasall"); } - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (!xflag) { + if (!xflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { vbiasall[i][0] = v[i][0]; v[i][0] = 0.0; } - if (!yflag) { + } + if (!yflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { vbiasall[i][1] = v[i][1]; v[i][1] = 0.0; } - if (!zflag) { + } + if (!zflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { vbiasall[i][2] = v[i][2]; v[i][2] = 0.0; } - } + } } /* ---------------------------------------------------------------------- @@ -233,12 +239,21 @@ int *mask = atom->mask; int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (!xflag) v[i][0] += vbiasall[i][0]; - if (!yflag) v[i][1] += vbiasall[i][1]; - if (!zflag) v[i][2] += vbiasall[i][2]; - } + if (!xflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][0] += vbiasall[i][0]; + } + if (!yflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][1] += vbiasall[i][1]; + } + if (!zflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][2] += vbiasall[i][2]; + } } /* ---------------------------------------------------------------------- */ diff -Naur lammps-28Apr09/src/compute_temp_partial.h lammps-30Apr09/src/compute_temp_partial.h --- lammps-28Apr09/src/compute_temp_partial.h 2008-03-20 18:14:38.000000000 -0600 +++ lammps-30Apr09/src/compute_temp_partial.h 2009-04-29 10:54:06.000000000 -0600 @@ -37,9 +37,6 @@ int xflag,yflag,zflag; int fix_dof; double tfactor; - double vbias[3]; // stored velocity bias for one atom - double **vbiasall; // stored velocity bias for all atoms - int maxbias; // size of vbiasall array void dof_compute(); }; diff -Naur lammps-28Apr09/src/compute_temp_profile.cpp lammps-30Apr09/src/compute_temp_profile.cpp --- lammps-28Apr09/src/compute_temp_profile.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-30Apr09/src/compute_temp_profile.cpp 2009-04-29 11:22:30.000000000 -0600 @@ -0,0 +1,494 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "stdlib.h" +#include "string.h" +#include "compute_temp_profile.h" +#include "atom.h" +#include "update.h" +#include "force.h" +#include "group.h" +#include "modify.h" +#include "fix.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +/* ---------------------------------------------------------------------- */ + +ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 7) error->all("Illegal compute temp/profile command"); + + scalar_flag = vector_flag = 1; + size_vector = 6; + extscalar = 0; + extvector = 1; + tempflag = 1; + tempbias = 1; + + xflag = atoi(arg[3]); + yflag = atoi(arg[4]); + zflag = atoi(arg[5]); + + nbinx = nbiny = nbinz = 1; + + if (strcmp(arg[6],"x") == 0) { + if (narg != 8) error->all("Illegal compute temp/profile command"); + nbinx = atoi(arg[7]); + ivx = 0; + ncount = 1; + } else if (strcmp(arg[6],"y") == 0) { + if (narg != 8) error->all("Illegal compute temp/profile command"); + nbiny = atoi(arg[7]); + ivy = 0; + ncount = 1; + } else if (strcmp(arg[6],"z") == 0) { + if (narg != 8) error->all("Illegal compute temp/profile command"); + if (domain->dimension == 2) + error->all("Compute temp/profile command cannot bin z for 2d systems"); + nbinz = atoi(arg[7]); + ivz = 0; + ncount = 1; + } else if (strcmp(arg[6],"xy") == 0) { + if (narg != 9) error->all("Illegal compute temp/profile command"); + nbinx = atoi(arg[7]); + nbiny = atoi(arg[8]); + ivx = 0; + ivy = 1; + ncount = 2; + } else if (strcmp(arg[6],"yz") == 0) { + if (narg != 9) error->all("Illegal compute temp/profile command"); + if (domain->dimension == 2) + error->all("Compute temp/profile command cannot bin z for 2d systems"); + nbiny = atoi(arg[7]); + nbinz = atoi(arg[8]); + ivy = 0; + ivz = 1; + ncount = 2; + } else if (strcmp(arg[6],"xz") == 0) { + if (narg != 9) error->all("Illegal compute temp/profile command"); + if (domain->dimension == 2) + error->all("Compute temp/profile command cannot bin z for 2d systems"); + nbinx = atoi(arg[7]); + nbinz = atoi(arg[8]); + ivx = 0; + ivz = 1; + ncount = 2; + } else if (strcmp(arg[6],"xyz") == 0) { + if (narg != 10) error->all("Illegal compute temp/profile command"); + if (domain->dimension == 2) + error->all("Compute temp/profile command cannot bin z for 2d systems"); + nbinx = atoi(arg[7]); + nbiny = atoi(arg[8]); + nbinz = atoi(arg[9]); + ivx = 0; + ivy = 1; + ivz = 2; + ncount = 3; + } else error->all("Illegal compute temp/profile command"); + + nbins = nbinx*nbiny*nbinz; + if (nbins <= 0) error->all("Illegal compute temp/profile command"); + + vbin = memory->create_2d_double_array(nbins,ncount+1,"temp/profile:vbin"); + binave = memory->create_2d_double_array(nbins,ncount+1, + "temp/profile:binave"); + + maxatom = 0; + bin = NULL; + + maxbias = 0; + vbiasall = NULL; + vector = new double[6]; +} + +/* ---------------------------------------------------------------------- */ + +ComputeTempProfile::~ComputeTempProfile() +{ + memory->destroy_2d_double_array(vbin); + memory->destroy_2d_double_array(binave); + memory->sfree(bin); + memory->destroy_2d_double_array(vbiasall); + delete [] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempProfile::init() +{ + fix_dof = 0; + for (int i = 0; i < modify->nfix; i++) + fix_dof += modify->fix[i]->dof(igroup); + dof_compute(); + + // ptrs to domain data + + box_change = domain->box_change; + triclinic = domain->triclinic; + periodicity = domain->periodicity; + + if (triclinic) { + boxlo = domain->boxlo_lamda; + boxhi = domain->boxhi_lamda; + prd = domain->prd_lamda; + } else { + boxlo = domain->boxlo; + boxhi = domain->boxhi; + prd = domain->prd; + } + + if (!box_change) bin_setup(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempProfile::dof_compute() +{ + double natoms = group->count(igroup); + int nper = domain->dimension; + dof = nper * natoms; + dof -= extra_dof + fix_dof; + if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); + else tfactor = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +double ComputeTempProfile::compute_scalar() +{ + int ibin; + double vthermal[3]; + + invoked_scalar = update->ntimestep; + + bin_average(); + + double **x = atom->x; + double **v = atom->v; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double t = 0.0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + ibin = bin[i]; + if (xflag) vthermal[0] = v[i][0] - binave[ibin][ivx]; + else vthermal[0] = v[i][0]; + if (yflag) vthermal[1] = v[i][1] - binave[ibin][ivy]; + else vthermal[1] = v[i][1]; + if (zflag) vthermal[2] = v[i][2] - binave[ibin][ivz]; + else vthermal[2] = v[i][2]; + + if (rmass) + t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + + vthermal[2]*vthermal[2]) * rmass[i]; + else + t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + + vthermal[2]*vthermal[2]) * mass[type[i]]; + } + + MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); + if (dynamic) dof_compute(); + scalar *= tfactor; + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempProfile::compute_vector() +{ + int i,ibin; + double vthermal[3]; + + invoked_vector = update->ntimestep; + + bin_average(); + + double **x = atom->x; + double **v = atom->v; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double massone,t[6]; + for (i = 0; i < 6; i++) t[i] = 0.0; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + ibin = bin[i]; + if (xflag) vthermal[0] = v[i][0] - binave[ibin][ivx]; + else vthermal[0] = v[i][0]; + if (yflag) vthermal[1] = v[i][1] - binave[ibin][ivy]; + else vthermal[1] = v[i][1]; + if (zflag) vthermal[2] = v[i][2] - binave[ibin][ivz]; + else vthermal[2] = v[i][2]; + + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + t[0] += massone * vthermal[0]*vthermal[0]; + t[1] += massone * vthermal[1]*vthermal[1]; + t[2] += massone * vthermal[2]*vthermal[2]; + t[3] += massone * vthermal[0]*vthermal[1]; + t[4] += massone * vthermal[0]*vthermal[2]; + t[5] += massone * vthermal[1]*vthermal[2]; + } + + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; +} + +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempProfile::remove_bias(int i, double *v) +{ + int ibin = bin[i]; + if (xflag) { + vbias[0] = binave[ibin][ivx]; + v[0] -= vbias[0]; + } + if (yflag) { + vbias[1] = binave[ibin][ivy]; + v[1] -= vbias[1]; + } + if (zflag) { + vbias[2] = binave[ibin][ivz]; + v[2] -= vbias[2]; + } +} + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempProfile::remove_bias_all() +{ + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (nlocal > maxbias) { + memory->destroy_2d_double_array(vbiasall); + maxbias = atom->nmax; + vbiasall = memory->create_2d_double_array(maxbias,3, + "compute/temp:vbiasall"); + } + + int ibin; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + ibin = bin[i]; + if (xflag) { + vbiasall[i][0] = binave[ibin][ivx]; + v[i][0] -= vbiasall[i][0]; + } + if (yflag) { + vbiasall[i][1] = binave[ibin][ivy]; + v[i][1] -= vbiasall[i][1]; + } + if (zflag) { + vbiasall[i][2] = binave[ibin][ivz]; + v[i][2] -= vbiasall[i][2]; + } + } +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to atom I removed by remove_bias() + assume remove_bias() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempProfile::restore_bias(int i, double *v) +{ + if (xflag) v[0] += vbias[0]; + if (yflag) v[1] += vbias[1]; + if (zflag) v[2] += vbias[2]; +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to all atoms removed by remove_bias_all() + assume remove_bias_all() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempProfile::restore_bias_all() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (xflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][0] += vbiasall[i][0]; + } + if (yflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][1] += vbiasall[i][1]; + } + if (zflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][2] += vbiasall[i][2]; + } +} + +/* ---------------------------------------------------------------------- + compute average velocity in each bin +------------------------------------------------------------------------- */ + +void ComputeTempProfile::bin_average() +{ + int i,j,ibin; + + if (box_change) bin_setup(); + bin_assign(); + + // clear bins, including particle count + + int nc1 = ncount + 1; + for (i = 0; i < nbins; i++) + for (j = 0; j < nc1; j++) + vbin[i][j] = 0.0; + + // sum each particle's velocity to appropriate bin + + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + ibin = bin[i]; + if (xflag) vbin[ibin][ivx] += v[i][0]; + if (yflag) vbin[ibin][ivy] += v[i][1]; + if (zflag) vbin[ibin][ivz] += v[i][2]; + vbin[ibin][ncount] += 1.0; + } + + // sum bins across processors + + MPI_Allreduce(vbin[0],binave[0],nbins*(ncount+1),MPI_DOUBLE,MPI_SUM,world); + + // compute ave velocity in each bin, checking for no particles + + for (i = 0; i < nbins; i++) + if (binave[i][ncount] > 0.0) + for (j = 0; j < ncount; j++) + binave[i][j] /= binave[i][ncount]; +} + +/* ---------------------------------------------------------------------- + set bin sizes, redo if box size changes +------------------------------------------------------------------------- */ + +void ComputeTempProfile::bin_setup() +{ + invdelta[0] = nbinx / prd[0]; + invdelta[1] = nbiny / prd[1]; + invdelta[2] = nbinz / prd[2]; +} + +/* ---------------------------------------------------------------------- + assign all atoms to bins +------------------------------------------------------------------------- */ + +void ComputeTempProfile::bin_assign() +{ + // reallocate bin array if necessary + + if (atom->nlocal > maxatom) { + maxatom = atom->nmax; + memory->sfree(bin); + bin = (int *) memory->smalloc(maxatom*sizeof(int),"temp/profile:bin"); + } + + // assign each atom to a bin, accounting for PBC + // if triclinic, do this in lamda space + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int ibinx,ibiny,ibinz; + double coord; + + if (triclinic) domain->x2lamda(nlocal); + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (nbinx > 1) { + coord = x[i][0]; + if (periodicity[0]) { + if (coord < boxlo[0]) coord += prd[0]; + if (coord >= boxhi[0]) coord -= prd[0]; + } + ibinx = static_cast
((coord - boxlo[0]) * invdelta[0]); + ibinx = MAX(ibinx,0); + ibinx = MIN(ibinx,nbinx-1); + } else ibinx = 0; + + if (nbiny > 1) { + coord = x[i][1]; + if (periodicity[1]) { + if (coord < boxlo[1]) coord += prd[1]; + if (coord >= boxhi[1]) coord -= prd[1]; + } + ibiny = static_cast ((coord - boxlo[1]) * invdelta[1]); + ibiny = MAX(ibiny,0); + ibiny = MIN(ibiny,nbiny-1); + } else ibiny = 0; + + if (nbinz > 1) { + coord = x[i][2]; + if (periodicity[2]) { + if (coord < boxlo[2]) coord += prd[2]; + if (coord >= boxhi[2]) coord -= prd[2]; + } + ibinz = static_cast ((coord - boxlo[2]) * invdelta[2]); + ibinz = MAX(ibinz,0); + ibinz = MIN(ibinz,nbinz-1); + } else ibinz = 0; + + bin[i] = nbinx*nbiny*ibinz + nbinx*ibiny + ibinx; + } + + if (triclinic) domain->lamda2x(nlocal); +} + +/* ---------------------------------------------------------------------- */ + +double ComputeTempProfile::memory_usage() +{ + double bytes = maxbias * sizeof(double); + bytes += maxatom * sizeof(int); + bytes += nbins*ncount * sizeof(double); + return bytes; +} diff -Naur lammps-28Apr09/src/compute_temp_profile.h lammps-30Apr09/src/compute_temp_profile.h --- lammps-28Apr09/src/compute_temp_profile.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-30Apr09/src/compute_temp_profile.h 2009-04-29 10:54:06.000000000 -0600 @@ -0,0 +1,59 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef COMPUTE_TEMP_PROFILE_H +#define COMPUTE_TEMP_PROFILE_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeTempProfile : public Compute { + public: + ComputeTempProfile(class LAMMPS *, int, char **); + ~ComputeTempProfile(); + void init(); + double compute_scalar(); + void compute_vector(); + + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(int, double *); + void restore_bias_all(); + double memory_usage(); + + private: + int xflag,yflag,zflag,ncount; + int nbinx,nbiny,nbinz,nbins; + int ivx,ivy,ivz; + int fix_dof; + double tfactor; + + int box_change,triclinic; + int *periodicity; + double *boxlo,*boxhi,*prd; + double invdelta[3]; + + int maxatom; + int *bin; + double **vbin,**binave; + + void dof_compute(); + void bin_average(); + void bin_setup(); + void bin_assign(); +}; + +} + +#endif diff -Naur lammps-28Apr09/src/compute_temp_ramp.cpp lammps-30Apr09/src/compute_temp_ramp.cpp --- lammps-28Apr09/src/compute_temp_ramp.cpp 2009-02-12 14:36:52.000000000 -0700 +++ lammps-30Apr09/src/compute_temp_ramp.cpp 2009-04-29 10:54:06.000000000 -0600 @@ -38,6 +38,13 @@ { if (narg < 9) error->all("Illegal compute temp command"); + scalar_flag = vector_flag = 1; + size_vector = 6; + extscalar = 0; + extvector = 1; + tempflag = 1; + tempbias = 1; + // parse optional args scaleflag = 1; @@ -99,15 +106,6 @@ coord_hi = zscale*atof(arg[8]); } - // settings - - scalar_flag = vector_flag = 1; - size_vector = 6; - extscalar = 0; - extvector = 1; - tempflag = 1; - tempbias = 1; - maxbias = 0; vbiasall = NULL; vector = new double[6]; diff -Naur lammps-28Apr09/src/compute_temp_ramp.h lammps-30Apr09/src/compute_temp_ramp.h --- lammps-28Apr09/src/compute_temp_ramp.h 2008-03-20 18:14:38.000000000 -0600 +++ lammps-30Apr09/src/compute_temp_ramp.h 2009-04-29 10:54:06.000000000 -0600 @@ -39,9 +39,6 @@ double v_lo,v_hi; int scaleflag,fix_dof; double tfactor,xscale,yscale,zscale; - double vbias[3]; // stored velocity bias for one atom - double **vbiasall; // stored velocity bias for all atoms - int maxbias; // size of vbiasall array void dof_compute(); }; diff -Naur lammps-28Apr09/src/compute_temp_region.h lammps-30Apr09/src/compute_temp_region.h --- lammps-28Apr09/src/compute_temp_region.h 2008-03-20 18:14:38.000000000 -0600 +++ lammps-30Apr09/src/compute_temp_region.h 2009-04-29 10:54:06.000000000 -0600 @@ -35,9 +35,6 @@ private: int iregion; - double vbias[3]; // stored velocity bias for one atom - double **vbiasall; // stored velocity bias for all atoms - int maxbias; // size of vbiasall array }; } diff -Naur lammps-28Apr09/src/group.cpp lammps-30Apr09/src/group.cpp --- lammps-28Apr09/src/group.cpp 2009-04-28 13:46:52.000000000 -0600 +++ lammps-30Apr09/src/group.cpp 2009-04-29 10:54:06.000000000 -0600 @@ -774,9 +774,11 @@ } MPI_Allreduce(cmone,cm,3,MPI_DOUBLE,MPI_SUM,world); - cm[0] /= masstotal; - cm[1] /= masstotal; - cm[2] /= masstotal; + if (masstotal > 0.0) { + cm[0] /= masstotal; + cm[1] /= masstotal; + cm[2] /= masstotal; + } } /* ---------------------------------------------------------------------- @@ -833,9 +835,11 @@ } MPI_Allreduce(cmone,cm,3,MPI_DOUBLE,MPI_SUM,world); - cm[0] /= masstotal; - cm[1] /= masstotal; - cm[2] /= masstotal; + if (masstotal > 0.0) { + cm[0] /= masstotal; + cm[1] /= masstotal; + cm[2] /= masstotal; + } } /* ---------------------------------------------------------------------- @@ -877,9 +881,11 @@ } MPI_Allreduce(p,cm,3,MPI_DOUBLE,MPI_SUM,world); - cm[0] /= masstotal; - cm[1] /= masstotal; - cm[2] /= masstotal; + if (masstotal > 0.0) { + cm[0] /= masstotal; + cm[1] /= masstotal; + cm[2] /= masstotal; + } } /* ---------------------------------------------------------------------- @@ -923,9 +929,11 @@ } MPI_Allreduce(p,cm,3,MPI_DOUBLE,MPI_SUM,world); - cm[0] /= masstotal; - cm[1] /= masstotal; - cm[2] /= masstotal; + if (masstotal > 0.0) { + cm[0] /= masstotal; + cm[1] /= masstotal; + cm[2] /= masstotal; + } } /* ---------------------------------------------------------------------- @@ -1092,7 +1100,8 @@ double rg_all; MPI_Allreduce(&rg,&rg_all,1,MPI_DOUBLE,MPI_SUM,world); - return sqrt(rg_all/masstotal); + if (masstotal > 0.0) return sqrt(rg_all/masstotal); + return 0.0; } /* ---------------------------------------------------------------------- @@ -1136,7 +1145,8 @@ double rg_all; MPI_Allreduce(&rg,&rg_all,1,MPI_DOUBLE,MPI_SUM,world); - return sqrt(rg_all/masstotal); + if (masstotal > 0.0) return sqrt(rg_all/masstotal); + return 0.0; } /* ---------------------------------------------------------------------- @@ -1264,10 +1274,10 @@ inertia[0][1]*inertia[1][0]*inertia[2][2] - inertia[2][0]*inertia[1][1]*inertia[0][2]; - int i,j; - for (i = 0; i < 3; i++) - for (j = 0; j < 3; j++) - inverse[i][j] /= determinant; + if (determinant > 0.0) + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + inverse[i][j] /= determinant; w[0] = inverse[0][0]*angmom[0] + inverse[0][1]*angmom[1] + inverse[0][2]*angmom[2]; diff -Naur lammps-28Apr09/src/style.h lammps-30Apr09/src/style.h --- lammps-28Apr09/src/style.h 2009-04-28 08:51:48.000000000 -0600 +++ lammps-30Apr09/src/style.h 2009-04-29 10:54:06.000000000 -0600 @@ -94,6 +94,7 @@ #include "compute_temp_com.h" #include "compute_temp_deform.h" #include "compute_temp_partial.h" +#include "compute_temp_profile.h" #include "compute_temp_ramp.h" #include "compute_temp_region.h" #include "compute_temp_sphere.h" @@ -118,6 +119,7 @@ ComputeStyle(temp/com,ComputeTempCOM) ComputeStyle(temp/deform,ComputeTempDeform) ComputeStyle(temp/partial,ComputeTempPartial) +ComputeStyle(temp/profile,ComputeTempProfile) ComputeStyle(temp/ramp,ComputeTempRamp) ComputeStyle(temp/region,ComputeTempRegion) ComputeStyle(temp/sphere,ComputeTempSphere)