diff -Naur lammps-27Jun09/doc/atom_style.html lammps-28Jun09/doc/atom_style.html --- lammps-27Jun09/doc/atom_style.html 2009-03-04 08:25:43.000000000 -0700 +++ lammps-28Jun09/doc/atom_style.html 2009-06-24 10:08:00.000000000 -0600 @@ -92,7 +92,7 @@ "molecular" package. The granular style is part of the "granular" package. The dpd style is part of the "dpd" package. The dipole style is part of the "dipole" package. The ellipsoid style is part -of the "ellipsoid" package. The peri style is part of the "peri" +of the "asphere" package. The peri style is part of the "peri" package for Peridynamics. They are only enabled if LAMMPS was built with that package. See the Making LAMMPS section for more info. diff -Naur lammps-27Jun09/doc/atom_style.txt lammps-28Jun09/doc/atom_style.txt --- lammps-27Jun09/doc/atom_style.txt 2009-03-04 08:25:43.000000000 -0700 +++ lammps-28Jun09/doc/atom_style.txt 2009-06-24 10:08:00.000000000 -0600 @@ -88,7 +88,7 @@ "molecular" package. The {granular} style is part of the "granular" package. The {dpd} style is part of the "dpd" package. The {dipole} style is part of the "dipole" package. The {ellipsoid} style is part -of the "ellipsoid" package. The {peri} style is part of the "peri" +of the "asphere" package. The {peri} style is part of the "peri" package for Peridynamics. They are only enabled if LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#2_3 section for more info. diff -Naur lammps-27Jun09/doc/compute_erotate_asphere.html lammps-28Jun09/doc/compute_erotate_asphere.html --- lammps-27Jun09/doc/compute_erotate_asphere.html 2008-03-19 11:04:31.000000000 -0600 +++ lammps-28Jun09/doc/compute_erotate_asphere.html 2009-06-26 12:22:33.000000000 -0600 @@ -47,8 +47,19 @@ angular momentum and a shape which is determined by the shape command.

This compute requires that atoms store angular momentum and a +quaternion to represent their orientation, as defined by the +atom_style. It also require they store a per-type +shape. The particles cannot store a per-particle +diameter or per-particle mass. +

All particles in the group must be finite-size. They cannot be point +particles, but they can be aspherical or spherical. +

Related commands: none

compute erotate/sphere +

Default: none

diff -Naur lammps-27Jun09/doc/compute_erotate_asphere.txt lammps-28Jun09/doc/compute_erotate_asphere.txt --- lammps-27Jun09/doc/compute_erotate_asphere.txt 2008-03-19 11:04:31.000000000 -0600 +++ lammps-28Jun09/doc/compute_erotate_asphere.txt 2009-06-26 12:22:33.000000000 -0600 @@ -44,6 +44,17 @@ angular momentum and a shape which is determined by the "shape"_shape.html command. +This compute requires that atoms store angular momentum and a +quaternion to represent their orientation, as defined by the +"atom_style"_atom_style.html. It also require they store a per-type +"shape"_shape.html. The particles cannot store a per-particle +diameter or per-particle mass. + +All particles in the group must be finite-size. They cannot be point +particles, but they can be aspherical or spherical. + [Related commands:] none +"compute erotate/sphere"_compute_rotate_sphere.html + [Default:] none diff -Naur lammps-27Jun09/doc/compute_erotate_sphere.html lammps-28Jun09/doc/compute_erotate_sphere.html --- lammps-27Jun09/doc/compute_erotate_sphere.html 2008-03-18 15:02:30.000000000 -0600 +++ lammps-28Jun09/doc/compute_erotate_sphere.html 2009-06-26 12:22:33.000000000 -0600 @@ -41,13 +41,17 @@

Restrictions:

This compute requires that particles be represented as extended -spheres and not point particles. This means they will have an angular -velocity and a diameter which is determined either by the -shape command or by each particle being assigned an -individual radius, e.g. for atom_style granular. +

This compute requires that atoms store angular velocity (omega) as +defined by the atom_style. It also require they +store either a per-particle diameter or per-type shape.

Related commands: none +

All particles in the group must be finite-size spheres or point +particles. They cannot be aspherical. Point particles will not +contribute to the rotational energy. +

Related commands: +

compute erotate/asphere

Default: none

diff -Naur lammps-27Jun09/doc/compute_erotate_sphere.txt lammps-28Jun09/doc/compute_erotate_sphere.txt --- lammps-27Jun09/doc/compute_erotate_sphere.txt 2008-03-18 15:02:30.000000000 -0600 +++ lammps-28Jun09/doc/compute_erotate_sphere.txt 2009-06-26 12:22:33.000000000 -0600 @@ -38,12 +38,16 @@ [Restrictions:] -This compute requires that particles be represented as extended -spheres and not point particles. This means they will have an angular -velocity and a diameter which is determined either by the -"shape"_shape.html command or by each particle being assigned an -individual radius, e.g. for "atom_style granular"_atom_style.html. +This compute requires that atoms store angular velocity (omega) as +defined by the "atom_style"_atom_style.html. It also require they +store either a per-particle diameter or per-type "shape"_shape.html. -[Related commands:] none +All particles in the group must be finite-size spheres or point +particles. They cannot be aspherical. Point particles will not +contribute to the rotational energy. + +[Related commands:] + +"compute erotate/asphere"_compute_rotate_asphere.html [Default:] none diff -Naur lammps-27Jun09/doc/compute_temp_asphere.html lammps-28Jun09/doc/compute_temp_asphere.html --- lammps-27Jun09/doc/compute_temp_asphere.html 2009-04-29 10:54:14.000000000 -0600 +++ lammps-28Jun09/doc/compute_temp_asphere.html 2009-06-26 12:41:31.000000000 -0600 @@ -32,23 +32,26 @@ usual compute temp command, which assumes point particles with only translational kinetic energy.

For 3d aspherical particles, each has 3, 5, or 6 degrees of freedom (3 -translational, remainder rotational), depending on whether the -particle is spherical, uniaxial, or biaxial. This is determined by -the shape command. Uniaxial means two of its three shape -parameters are equal. Biaxial means all 3 shape parameters are -different. -

For 2d aspherical particles, each has 3 or 4 degrees of freedom (3 -translational, remainder rotational), depending on whether the -particle is spherical, or biaxial. Biaxial means the x,y shape -parameters are unequal. -

IMPORTANT NOTE: These degrees of freedom assume that the interaction -potential between degenerate aspherical particles does not impart -rotational motion to the extra degrees of freedom. E.g. the GayBerne -pair potential does not impart torque to spherical -particles, so they do not rotate. +

For 3d aspherical particles, each has 6 degrees of freedom (3 +translational, 3 rotational). For 2d aspherical particles, each has 3 +degrees of freedom (2 translational, 1 rotational). +

IMPORTANT NOTE: This choice for degrees of freedom (dof) makes the +assumption that all aspherical particles in your model will freely +rotate, sampling all their rotational dof. It is possible to use a +combination of interaction potentials and fixes that induce no torque +or otherwise constrain some of all of your particles so that this is +not the case. Then there are less dof and you should use the +compute_modify extra command to adjust the dof +accordingly. +

For example, an aspherical particle with all three of its +shape parameters the same is a sphere. If it does not +rotate, then it should have 3 dof instead of 6 in 3d (or 2 instead of +3 in 2d). A uniaxial aspherical particle has two of its three shape +parameters the same. If it does not rotate around the axis +perpendicular to its circular cross section, then it should have 5 dof +instead of 6 in 3d.

The translational kinetic energy is computed the same as is described by the compute temp command. The rotational diff -Naur lammps-27Jun09/doc/compute_temp_asphere.txt lammps-28Jun09/doc/compute_temp_asphere.txt --- lammps-27Jun09/doc/compute_temp_asphere.txt 2009-04-29 10:54:14.000000000 -0600 +++ lammps-28Jun09/doc/compute_temp_asphere.txt 2009-06-26 12:41:31.000000000 -0600 @@ -29,23 +29,26 @@ usual "compute temp"_compute_temp.html command, which assumes point particles with only translational kinetic energy. -For 3d aspherical particles, each has 3, 5, or 6 degrees of freedom (3 -translational, remainder rotational), depending on whether the -particle is spherical, uniaxial, or biaxial. This is determined by -the "shape"_shape.html command. Uniaxial means two of its three shape -parameters are equal. Biaxial means all 3 shape parameters are -different. - -For 2d aspherical particles, each has 3 or 4 degrees of freedom (3 -translational, remainder rotational), depending on whether the -particle is spherical, or biaxial. Biaxial means the x,y shape -parameters are unequal. - -IMPORTANT NOTE: These degrees of freedom assume that the interaction -potential between degenerate aspherical particles does not impart -rotational motion to the extra degrees of freedom. E.g. the "GayBerne -pair potential"_pair_gayberne.html does not impart torque to spherical -particles, so they do not rotate. +For 3d aspherical particles, each has 6 degrees of freedom (3 +translational, 3 rotational). For 2d aspherical particles, each has 3 +degrees of freedom (2 translational, 1 rotational). + +IMPORTANT NOTE: This choice for degrees of freedom (dof) makes the +assumption that all aspherical particles in your model will freely +rotate, sampling all their rotational dof. It is possible to use a +combination of interaction potentials and fixes that induce no torque +or otherwise constrain some of all of your particles so that this is +not the case. Then there are less dof and you should use the +"compute_modify extra"_compute_modify.html command to adjust the dof +accordingly. + +For example, an aspherical particle with all three of its +"shape"_shape.html parameters the same is a sphere. If it does not +rotate, then it should have 3 dof instead of 6 in 3d (or 2 instead of +3 in 2d). A uniaxial aspherical particle has two of its three shape +parameters the same. If it does not rotate around the axis +perpendicular to its circular cross section, then it should have 5 dof +instead of 6 in 3d. The translational kinetic energy is computed the same as is described by the "compute temp"_compute_temp.html command. The rotational diff -Naur lammps-27Jun09/doc/compute_temp_sphere.html lammps-28Jun09/doc/compute_temp_sphere.html --- lammps-27Jun09/doc/compute_temp_sphere.html 2009-04-29 10:54:14.000000000 -0600 +++ lammps-28Jun09/doc/compute_temp_sphere.html 2009-06-26 12:41:31.000000000 -0600 @@ -36,6 +36,15 @@ translational, 3 rotational). For 2d spherical particles, each has 3 degrees of freedom (2 translational, 1 rotational).

IMPORTANT NOTE: This choice for degrees of freedom (dof) makes the +assumption that all spherical particles in your model will freely +rotate, sampling all their rotational dof. It is possible to use a +combination of interaction potentials and fixes that induce no torque +or otherwise constrain some of all of your particles so that this is +not the case. Then there are less dof and you should use the +compute_modify extra command to adjust the dof +accordingly. +

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 diff -Naur lammps-27Jun09/doc/compute_temp_sphere.txt lammps-28Jun09/doc/compute_temp_sphere.txt --- lammps-27Jun09/doc/compute_temp_sphere.txt 2009-04-29 10:54:14.000000000 -0600 +++ lammps-28Jun09/doc/compute_temp_sphere.txt 2009-06-26 12:41:31.000000000 -0600 @@ -33,6 +33,15 @@ translational, 3 rotational). For 2d spherical particles, each has 3 degrees of freedom (2 translational, 1 rotational). +IMPORTANT NOTE: This choice for degrees of freedom (dof) makes the +assumption that all spherical particles in your model will freely +rotate, sampling all their rotational dof. It is possible to use a +combination of interaction potentials and fixes that induce no torque +or otherwise constrain some of all of your particles so that this is +not the case. Then there are less dof and you should use the +"compute_modify extra"_compute_modify.html command to adjust the dof +accordingly. + 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 diff -Naur lammps-27Jun09/doc/fix_npt_asphere.html lammps-28Jun09/doc/fix_npt_asphere.html --- lammps-27Jun09/doc/fix_npt_asphere.html 2009-05-12 07:33:47.000000000 -0600 +++ lammps-28Jun09/doc/fix_npt_asphere.html 2009-06-26 12:22:33.000000000 -0600 @@ -217,10 +217,14 @@ LAMMPS was built with that package. See the Making LAMMPS section for more info.

This fix requires that particles be represented as extended ellipsoids -and not point particles. This means they will have an angular -momentum and a shape which is determined by the shape -command. +

This fix requires that atoms store torque and angular momentum and a +quaternion to represent their orientation, as defined by the +atom_style. It also require they store a per-type +shape. The particles cannot store a per-particle +diameter or per-particle mass. +

All particles in the group must be finite-size. They cannot be point +particles, but they can be aspherical or spherical.

Any dimension being adjusted by this fix must be periodic. A dimension whose target pressures are specified as NULL can be diff -Naur lammps-27Jun09/doc/fix_npt_asphere.txt lammps-28Jun09/doc/fix_npt_asphere.txt --- lammps-27Jun09/doc/fix_npt_asphere.txt 2009-05-12 07:33:47.000000000 -0600 +++ lammps-28Jun09/doc/fix_npt_asphere.txt 2009-06-26 12:22:33.000000000 -0600 @@ -206,10 +206,14 @@ LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#2_3 section for more info. -This fix requires that particles be represented as extended ellipsoids -and not point particles. This means they will have an angular -momentum and a shape which is determined by the "shape"_shape.html -command. +This fix requires that atoms store torque and angular momentum and a +quaternion to represent their orientation, as defined by the +"atom_style"_atom_style.html. It also require they store a per-type +"shape"_shape.html. The particles cannot store a per-particle +diameter or per-particle mass. + +All particles in the group must be finite-size. They cannot be point +particles, but they can be aspherical or spherical. Any dimension being adjusted by this fix must be periodic. A dimension whose target pressures are specified as NULL can be diff -Naur lammps-27Jun09/doc/fix_npt_sphere.html lammps-28Jun09/doc/fix_npt_sphere.html --- lammps-27Jun09/doc/fix_npt_sphere.html 2009-05-12 07:33:47.000000000 -0600 +++ lammps-28Jun09/doc/fix_npt_sphere.html 2009-06-26 12:22:33.000000000 -0600 @@ -215,10 +215,12 @@

Restrictions:

This fix requires that particles be represented as extended spheres -and not point particles. This means they will have an angular -velocity and a diameter which is determined by the shape -command. +

This fix requires that atoms store torque and angular velocity (omega) +as defined by the atom_style. It also require they +store either a per-particle diameter or per-type shape. +

All particles in the group must be finite-size spheres. They cannot +be point particles, nor can they be aspherical.

Any dimension being adjusted by this fix must be periodic. A dimension whose target pressures are specified as NULL can be diff -Naur lammps-27Jun09/doc/fix_npt_sphere.txt lammps-28Jun09/doc/fix_npt_sphere.txt --- lammps-27Jun09/doc/fix_npt_sphere.txt 2009-05-12 07:33:47.000000000 -0600 +++ lammps-28Jun09/doc/fix_npt_sphere.txt 2009-06-26 12:22:33.000000000 -0600 @@ -204,10 +204,12 @@ [Restrictions:] -This fix requires that particles be represented as extended spheres -and not point particles. This means they will have an angular -velocity and a diameter which is determined by the "shape"_shape.html -command. +This fix requires that atoms store torque and angular velocity (omega) +as defined by the "atom_style"_atom_style.html. It also require they +store either a per-particle diameter or per-type "shape"_shape.html. + +All particles in the group must be finite-size spheres. They cannot +be point particles, nor can they be aspherical. Any dimension being adjusted by this fix must be periodic. A dimension whose target pressures are specified as NULL can be diff -Naur lammps-27Jun09/doc/fix_nve_asphere.html lammps-28Jun09/doc/fix_nve_asphere.html --- lammps-27Jun09/doc/fix_nve_asphere.html 2008-03-19 16:51:12.000000000 -0600 +++ lammps-28Jun09/doc/fix_nve_asphere.html 2009-06-26 12:22:33.000000000 -0600 @@ -48,10 +48,14 @@ LAMMPS was built with that package. See the Making LAMMPS section for more info.

All particles in the group must be finite-size. They cannot be point +particles, but they can be aspherical or spherical.

Related commands:

diff -Naur lammps-27Jun09/doc/fix_nve_asphere.txt lammps-28Jun09/doc/fix_nve_asphere.txt --- lammps-27Jun09/doc/fix_nve_asphere.txt 2008-03-19 16:51:12.000000000 -0600 +++ lammps-28Jun09/doc/fix_nve_asphere.txt 2009-06-26 12:22:33.000000000 -0600 @@ -45,10 +45,14 @@ LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#2_3 section for more info. -This fix requires that particles be represented as extended ellipsoids -and not point particles. This means they will have an angular -momentum and a shape which is determined by the "shape"_shape.html -command. +This fix requires that atoms store torque and angular momentum and a +quaternion to represent their orientation, as defined by the +"atom_style"_atom_style.html. It also require they store a per-type +"shape"_shape.html. The particles cannot store a per-particle +diameter or per-particle mass. + +All particles in the group must be finite-size. They cannot be point +particles, but they can be aspherical or spherical. [Related commands:] diff -Naur lammps-27Jun09/doc/fix_nve_sphere.html lammps-28Jun09/doc/fix_nve_sphere.html --- lammps-27Jun09/doc/fix_nve_sphere.html 2008-03-19 16:51:12.000000000 -0600 +++ lammps-28Jun09/doc/fix_nve_sphere.html 2009-06-26 12:22:33.000000000 -0600 @@ -61,11 +61,13 @@

Restrictions:

This fix requires that particles be represented as extended spheres -and not point particles. This means they will have an angular -velocity and a diameter which is determined either by the -shape command or by each particle being assigned an -individual radius, e.g. for atom_style granular. +

This fix requires that atoms store torque and angular velocity (omega) +as defined by the atom_style. It also require they +store either a per-particle diameter or per-type shape. If +the dipole keyword is used, then they must store a dipole moment. +

All particles in the group must be finite-size spheres. They cannot +be point particles, nor can they be aspherical.

Related commands:

diff -Naur lammps-27Jun09/doc/fix_nve_sphere.txt lammps-28Jun09/doc/fix_nve_sphere.txt --- lammps-27Jun09/doc/fix_nve_sphere.txt 2008-03-19 16:51:12.000000000 -0600 +++ lammps-28Jun09/doc/fix_nve_sphere.txt 2009-06-26 12:22:33.000000000 -0600 @@ -53,11 +53,13 @@ [Restrictions:] -This fix requires that particles be represented as extended spheres -and not point particles. This means they will have an angular -velocity and a diameter which is determined either by the -"shape"_shape.html command or by each particle being assigned an -individual radius, e.g. for "atom_style granular"_atom_style.html. +This fix requires that atoms store torque and angular velocity (omega) +as defined by the "atom_style"_atom_style.html. It also require they +store either a per-particle diameter or per-type "shape"_shape.html. If +the {dipole} keyword is used, then they must store a dipole moment. + +All particles in the group must be finite-size spheres. They cannot +be point particles, nor can they be aspherical. [Related commands:] diff -Naur lammps-27Jun09/doc/fix_nvt_asphere.html lammps-28Jun09/doc/fix_nvt_asphere.html --- lammps-27Jun09/doc/fix_nvt_asphere.html 2009-05-12 07:33:47.000000000 -0600 +++ lammps-28Jun09/doc/fix_nvt_asphere.html 2009-06-26 12:22:33.000000000 -0600 @@ -140,10 +140,14 @@ LAMMPS was built with that package. See the Making LAMMPS section for more info.

All particles in the group must be finite-size. They cannot be point +particles, but they can be aspherical or spherical.

The final Tstop cannot be 0.0 since it would make the target T = 0.0 at some timestep during the simulation which is not allowed in diff -Naur lammps-27Jun09/doc/fix_nvt_asphere.txt lammps-28Jun09/doc/fix_nvt_asphere.txt --- lammps-27Jun09/doc/fix_nvt_asphere.txt 2009-05-12 07:33:47.000000000 -0600 +++ lammps-28Jun09/doc/fix_nvt_asphere.txt 2009-06-26 12:22:33.000000000 -0600 @@ -131,10 +131,14 @@ LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#2_3 section for more info. -This fix requires that particles be represented as extended ellipsoids -and not point particles. This means they will have an angular -momentum and a shape which is determined by the "shape"_shape.html -command. +This fix requires that atoms store torque and angular momentum and a +quaternion to represent their orientation, as defined by the +"atom_style"_atom_style.html. It also require they store a per-type +"shape"_shape.html. The particles cannot store a per-particle +diameter or per-particle mass. + +All particles in the group must be finite-size. They cannot be point +particles, but they can be aspherical or spherical. The final Tstop cannot be 0.0 since it would make the target T = 0.0 at some timestep during the simulation which is not allowed in diff -Naur lammps-27Jun09/doc/fix_nvt_sphere.html lammps-28Jun09/doc/fix_nvt_sphere.html --- lammps-27Jun09/doc/fix_nvt_sphere.html 2009-05-12 07:33:47.000000000 -0600 +++ lammps-28Jun09/doc/fix_nvt_sphere.html 2009-06-26 12:22:33.000000000 -0600 @@ -139,10 +139,12 @@

Restrictions:

This fix requires that atoms store torque and angular velocity (omega) +as defined by the atom_style. It also require they +store either a per-particle radius or per-type shape. +

All particles in the group must be finite-size spheres. They cannot +be point particles, nor can they be aspherical.

The final Tstop cannot be 0.0 since it would make the target T = 0.0 at some timestep during the simulation which is not allowed in diff -Naur lammps-27Jun09/doc/fix_nvt_sphere.txt lammps-28Jun09/doc/fix_nvt_sphere.txt --- lammps-27Jun09/doc/fix_nvt_sphere.txt 2009-05-12 07:33:47.000000000 -0600 +++ lammps-28Jun09/doc/fix_nvt_sphere.txt 2009-06-26 12:22:33.000000000 -0600 @@ -130,10 +130,12 @@ [Restrictions:] -This fix requires that particles be represented as extended spheres -and not point particles. This means they will have an angular -velocity and a diameter which is determined by the "shape"_shape.html -command. +This fix requires that atoms store torque and angular velocity (omega) +as defined by the "atom_style"_atom_style.html. It also require they +store either a per-particle radius or per-type "shape"_shape.html. + +All particles in the group must be finite-size spheres. They cannot +be point particles, nor can they be aspherical. The final Tstop cannot be 0.0 since it would make the target T = 0.0 at some timestep during the simulation which is not allowed in diff -Naur lammps-27Jun09/doc/mass.html lammps-28Jun09/doc/mass.html --- lammps-27Jun09/doc/mass.html 2009-05-19 08:57:46.000000000 -0600 +++ lammps-28Jun09/doc/mass.html 2009-06-26 12:22:33.000000000 -0600 @@ -66,8 +66,8 @@

If you define a hybrid atom style which includes one (or more) sub-styles which require per-type mass and one (or more) sub-styles which require per-atom mass, then you must define both. In -this case the per-type mass will be ignored; only the per-atom mass is -used by LAMMPS. +this case the per-type mass will be ignored; only the per-atom mass +will be used by LAMMPS.

Restrictions:

diff -Naur lammps-27Jun09/doc/mass.txt lammps-28Jun09/doc/mass.txt --- lammps-27Jun09/doc/mass.txt 2009-05-19 08:57:46.000000000 -0600 +++ lammps-28Jun09/doc/mass.txt 2009-06-26 12:22:33.000000000 -0600 @@ -63,8 +63,8 @@ If you define a "hybrid atom style"_atom_style.html which includes one (or more) sub-styles which require per-type mass and one (or more) sub-styles which require per-atom mass, then you must define both. In -this case the per-type mass will be ignored; only the per-atom mass is -used by LAMMPS. +this case the per-type mass will be ignored; only the per-atom mass +will be used by LAMMPS. [Restrictions:] diff -Naur lammps-27Jun09/doc/pair_gayberne.html lammps-28Jun09/doc/pair_gayberne.html --- lammps-27Jun09/doc/pair_gayberne.html 2009-05-19 08:43:06.000000000 -0600 +++ lammps-28Jun09/doc/pair_gayberne.html 2009-06-26 12:22:33.000000000 -0600 @@ -157,10 +157,16 @@

This style is part of the "asphere" package. It is only enabled if LAMMPS was built with that package. See the Making -LAMMPS section for more info. You must also -define a size and shape for each particle type via the -shape command which requires atom_style -ellipsoid. +LAMMPS section for more info. +

This pair style requires that atoms store torque and a quaternion to +represent their orientation, as defined by the +atom_style. It also require they store a per-type +shape. The particles cannot store a per-particle +diameter. +

Particles acted on by the potential can be extended aspherical or +spherical particles, or point particles.

The Gay-Berne potential does not become isotropic as r increases (Everaers). The distance-of-closest-approach diff -Naur lammps-27Jun09/doc/pair_gayberne.txt lammps-28Jun09/doc/pair_gayberne.txt --- lammps-27Jun09/doc/pair_gayberne.txt 2009-05-19 08:43:06.000000000 -0600 +++ lammps-28Jun09/doc/pair_gayberne.txt 2009-06-26 12:22:33.000000000 -0600 @@ -154,10 +154,16 @@ This style is part of the "asphere" package. It is only enabled if LAMMPS was built with that package. See the "Making -LAMMPS"_Section_start.html#2_3 section for more info. You must also -define a size and shape for each particle type via the -"shape"_shape.html command which requires "atom_style -ellipsoid"_atom_style.html. +LAMMPS"_Section_start.html#2_3 section for more info. + +This pair style requires that atoms store torque and a quaternion to +represent their orientation, as defined by the +"atom_style"_atom_style.html. It also require they store a per-type +"shape"_shape.html. The particles cannot store a per-particle +diameter. + +Particles acted on by the potential can be extended aspherical or +spherical particles, or point particles. The Gay-Berne potential does not become isotropic as r increases "(Everaers)"_#Everaers. The distance-of-closest-approach diff -Naur lammps-27Jun09/doc/pair_gran.html lammps-28Jun09/doc/pair_gran.html --- lammps-27Jun09/doc/pair_gran.html 2009-05-12 07:34:00.000000000 -0600 +++ lammps-28Jun09/doc/pair_gran.html 2009-06-26 12:22:33.000000000 -0600 @@ -189,6 +189,12 @@ is only enabled if LAMMPS was built with that package. See the Making LAMMPS section for more info.

These pair styles require that atoms store torque and angular velocity +(omega) as defined by the atom_style. They also +require they store a per-particle radius and that velocities are +communicated with ghost atoms. The granular and dpd atom styles +are the ones that currently do this. +

Related commands:

pair_coeff diff -Naur lammps-27Jun09/doc/pair_gran.txt lammps-28Jun09/doc/pair_gran.txt --- lammps-27Jun09/doc/pair_gran.txt 2009-05-12 07:34:00.000000000 -0600 +++ lammps-28Jun09/doc/pair_gran.txt 2009-06-26 12:22:33.000000000 -0600 @@ -179,6 +179,12 @@ is only enabled if LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#2_3 section for more info. +These pair styles require that atoms store torque and angular velocity +(omega) as defined by the "atom_style"_atom_style.html. They also +require they store a per-particle radius and that velocities are +communicated with ghost atoms. The {granular} and {dpd} atom styles +are the ones that currently do this. + [Related commands:] "pair_coeff"_pair_coeff.html diff -Naur lammps-27Jun09/doc/pair_lubricate.html lammps-28Jun09/doc/pair_lubricate.html --- lammps-27Jun09/doc/pair_lubricate.html 2009-04-16 08:29:04.000000000 -0600 +++ lammps-28Jun09/doc/pair_lubricate.html 2009-06-26 12:22:33.000000000 -0600 @@ -118,12 +118,16 @@ LAMMPS was built with that package. See the Making LAMMPS section for more info.

Because this potential computes forces and torques on particles, the -atom style must support particles whose size is set via the -shape command. This is atom_style -ellipsoid and dipole. Since only spherical mono-disperse particles -are currently allowed for pair_style lubricate, this means the 3 shape -radii for all particle types must be the same. +

All the shape settings must be for finite-size spheres. They cannot +be point particles, nor can they be aspherical. Additionally all the +shape types must specify particles of the same size, i.e. a +monodisperse system.

Related commands:

diff -Naur lammps-27Jun09/doc/pair_lubricate.txt lammps-28Jun09/doc/pair_lubricate.txt --- lammps-27Jun09/doc/pair_lubricate.txt 2009-04-16 08:29:04.000000000 -0600 +++ lammps-28Jun09/doc/pair_lubricate.txt 2009-06-26 12:22:33.000000000 -0600 @@ -115,12 +115,16 @@ LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#2_3 section for more info. -Because this potential computes forces and torques on particles, the -atom style must support particles whose size is set via the -"shape"_shape.html command. This is "atom_style"_atom_style.html -ellipsoid and dipole. Since only spherical mono-disperse particles -are currently allowed for pair_style lubricate, this means the 3 shape -radii for all particle types must be the same. +This pair style requires that atoms store torque and a quaternion to +represent their orientation, as defined by the +"atom_style"_atom_style.html. It also require they store a per-type +"shape"_shape.html. The particles cannot store a per-particle +diameter or per-particle mass. + +All the shape settings must be for finite-size spheres. They cannot +be point particles, nor can they be aspherical. Additionally all the +shape types must specify particles of the same size, i.e. a +monodisperse system. [Related commands:] diff -Naur lammps-27Jun09/doc/pair_resquared.html lammps-28Jun09/doc/pair_resquared.html --- lammps-27Jun09/doc/pair_resquared.html 2009-05-19 08:43:06.000000000 -0600 +++ lammps-28Jun09/doc/pair_resquared.html 2009-06-26 12:22:33.000000000 -0600 @@ -181,10 +181,16 @@

Particles acted on by the potential can be extended aspherical or +spherical particles, or point particles.

The distance-of-closest-approach approximation used by LAMMPS becomes less accurate when high-aspect ratio ellipsoids are used. diff -Naur lammps-27Jun09/doc/pair_resquared.txt lammps-28Jun09/doc/pair_resquared.txt --- lammps-27Jun09/doc/pair_resquared.txt 2009-05-19 08:43:06.000000000 -0600 +++ lammps-28Jun09/doc/pair_resquared.txt 2009-06-26 12:22:33.000000000 -0600 @@ -178,10 +178,16 @@ This style is part of the "asphere" package. It is only enabled if LAMMPS was built with that package. See the "Making -LAMMPS"_Section_start.html#2_3 section for more info. You must also -define a size and shape for each particle type via the -"shape"_shape.html command which requires "atom_style -ellipsoid"_atom_style.html. +LAMMPS"_Section_start.html#2_3 section for more info. + +This pair style requires that atoms store torque and a quaternion to +represent their orientation, as defined by the +"atom_style"_atom_style.html. It also require they store a per-type +"shape"_shape.html. The particles cannot store a per-particle +diameter. + +Particles acted on by the potential can be extended aspherical or +spherical particles, or point particles. The distance-of-closest-approach approximation used by LAMMPS becomes less accurate when high-aspect ratio ellipsoids are used. diff -Naur lammps-27Jun09/doc/read_data.html lammps-28Jun09/doc/read_data.html --- lammps-27Jun09/doc/read_data.html 2009-04-03 10:19:32.000000000 -0600 +++ lammps-28Jun09/doc/read_data.html 2009-06-26 12:22:33.000000000 -0600 @@ -314,6 +314,19 @@ to. It can be 0 if it is an unbonded atom or if you don't care to keep track of molecule assignments.

The diameter specifies the size of a finite size particle, analagous +to the shape command which sets the size on a per-type +basis. A diameter can be set to 0.0, which means that atom is a point +particle and not a finite-size particles. Some pair styles and fixes +and computes that operate on finite-size particles allow for a mixture +of finite-size and point particles. See the doc pages of individual +commands for details. +

The density is used in conjunction with the diameter to set the mass +of a particle as mass = density * volume. If the diameter and volume +are 0.0 meaning a point particle, then the mass is not 0.0 but is set +as mass = density. +

The values quatw, quati, quatj, and quatk set the orientation of the atom as a quaternion (4-vector). Note that the shape command or "Shapes" section of the data file diff -Naur lammps-27Jun09/doc/read_data.txt lammps-28Jun09/doc/read_data.txt --- lammps-27Jun09/doc/read_data.txt 2009-04-03 10:19:32.000000000 -0600 +++ lammps-28Jun09/doc/read_data.txt 2009-06-26 12:22:33.000000000 -0600 @@ -291,6 +291,19 @@ to. It can be 0 if it is an unbonded atom or if you don't care to keep track of molecule assignments. +The diameter specifies the size of a finite size particle, analagous +to the "shape"_shape.html command which sets the size on a per-type +basis. A diameter can be set to 0.0, which means that atom is a point +particle and not a finite-size particles. Some pair styles and fixes +and computes that operate on finite-size particles allow for a mixture +of finite-size and point particles. See the doc pages of individual +commands for details. + +The density is used in conjunction with the diameter to set the mass +of a particle as mass = density * volume. If the diameter and volume +are 0.0 meaning a point particle, then the mass is not 0.0 but is set +as mass = density. + The values {quatw}, {quati}, {quatj}, and {quatk} set the orientation of the atom as a quaternion (4-vector). Note that the "shape"_shape.html command or "Shapes" section of the data file diff -Naur lammps-27Jun09/doc/shape.html lammps-28Jun09/doc/shape.html --- lammps-27Jun09/doc/shape.html 2009-01-19 10:17:01.000000000 -0700 +++ lammps-28Jun09/doc/shape.html 2009-06-26 12:22:33.000000000 -0600 @@ -30,7 +30,9 @@

Set the shape for all atoms of one or more atom types. In LAMMPS, particles that have a finite size are said to have a "shape", as -opposed to being a point mass. Shape values can also be set in the +opposed to being a point mass. The shape can be spherical or +aspherical, depending on whether the 3 shape values are the same or +different. Shape values can also be set in the read_data data file using the "Shapes" keyword. See the units command for what distance units to use.

@@ -52,10 +54,11 @@

1 1.0 1.0 1.0

The shape values can be set to 0.0, which means that atoms of that -type are point masses and not finite-size particles. Pair styles and -fixes that rely on particles having a finite size should not be used -for such particles. +

The shape values can be set to all 0.0, which means that atoms of that +type are point particles and not finite-size particles. Some pair +styles and fixes and computes that operate on finite-size particles +allow for a mixture of finite-size and point particles. See the doc +pages of individual commands for details.

Note that the shape command can only be used if the atom style requires per-type atom shape to be set. @@ -81,10 +84,10 @@

If you define a hybrid atom style which includes one (or more) sub-styles which require per-type shape and one (or more) -sub-styles which require per-atom shape, then you must define both. -If the per-atom diameter is set to 0.0, the per-type shape is used. -If the per-atom diameter is non-zero, then the per-type shape is -ignored. +sub-styles which require per-atom diameter, then you must define both. +In this case the per-type shape will be ignored; only the per-atom +diameter will be used by LAMMPS. Note that this means you can not +currently mix aspherical particles with per-atom diameter particles.

Restrictions:

diff -Naur lammps-27Jun09/doc/shape.txt lammps-28Jun09/doc/shape.txt --- lammps-27Jun09/doc/shape.txt 2009-01-19 10:17:01.000000000 -0700 +++ lammps-28Jun09/doc/shape.txt 2009-06-26 12:22:33.000000000 -0600 @@ -27,7 +27,9 @@ Set the shape for all atoms of one or more atom types. In LAMMPS, particles that have a finite size are said to have a "shape", as -opposed to being a point mass. Shape values can also be set in the +opposed to being a point mass. The shape can be spherical or +aspherical, depending on whether the 3 shape values are the same or +different. Shape values can also be set in the "read_data"_read_data.html data file using the "Shapes" keyword. See the "units"_units.html command for what distance units to use. @@ -49,10 +51,11 @@ 1 1.0 1.0 1.0 :pre -The shape values can be set to 0.0, which means that atoms of that -type are point masses and not finite-size particles. Pair styles and -fixes that rely on particles having a finite size should not be used -for such particles. +The shape values can be set to all 0.0, which means that atoms of that +type are point particles and not finite-size particles. Some pair +styles and fixes and computes that operate on finite-size particles +allow for a mixture of finite-size and point particles. See the doc +pages of individual commands for details. Note that the shape command can only be used if the "atom style"_atom_style.html requires per-type atom shape to be set. @@ -78,10 +81,10 @@ If you define a "hybrid atom style"_atom_style.html which includes one (or more) sub-styles which require per-type shape and one (or more) -sub-styles which require per-atom shape, then you must define both. -If the per-atom diameter is set to 0.0, the per-type shape is used. -If the per-atom diameter is non-zero, then the per-type shape is -ignored. +sub-styles which require per-atom diameter, then you must define both. +In this case the per-type shape will be ignored; only the per-atom +diameter will be used by LAMMPS. Note that this means you can not +currently mix aspherical particles with per-atom diameter particles. [Restrictions:] diff -Naur lammps-27Jun09/src/ASPHERE/compute_erotate_asphere.cpp lammps-28Jun09/src/ASPHERE/compute_erotate_asphere.cpp --- lammps-27Jun09/src/ASPHERE/compute_erotate_asphere.cpp 2009-01-05 15:26:08.000000000 -0700 +++ lammps-28Jun09/src/ASPHERE/compute_erotate_asphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -15,6 +15,7 @@ #include "compute_erotate_asphere.h" #include "math_extra.h" #include "atom.h" +#include "atom_vec.h" #include "update.h" #include "force.h" #include "memory.h" @@ -24,19 +25,27 @@ /* ---------------------------------------------------------------------- */ -ComputeERotateAsphere::ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) : +ComputeERotateAsphere:: +ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all("Illegal compute erotate/asphere command"); - if (!atom->angmom_flag || !atom->quat_flag) - error->all("Compute erotate/asphere requires atom attributes angmom, quat"); - scalar_flag = 1; extscalar = 1; inertia = memory->create_2d_double_array(atom->ntypes+1,3,"fix_temp_sphere:inertia"); + + // error checks + + if (!atom->angmom_flag || !atom->quat_flag || + !atom->avec->shape_type) + error->all("Compute erotate/asphere requires atom attributes " + "angmom, quat, shape"); + if (atom->radius_flag || atom->rmass_flag) + error->all("Compute erotate/asphere cannot be used with atom attributes " + "diameter or rmass"); } /* ---------------------------------------------------------------------- */ @@ -50,11 +59,21 @@ void ComputeERotateAsphere::init() { - pfactor = 0.5 * force->mvv2e; + // check that all particles are finite-size + // no point particles allowed, spherical is OK + + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (!atom->shape) - error->all("Compute erotate/asphere requires atom attribute shape"); + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (shape[type[i]][0] == 0.0) + error->one("Compue erotate/asphere requires extended particles"); + pfactor = 0.5 * force->mvv2e; calculate_inertia(); } @@ -62,6 +81,8 @@ double ComputeERotateAsphere::compute_scalar() { + int i,itype; + invoked_scalar = update->ntimestep; double **quat = atom->quat; @@ -70,12 +91,14 @@ int *type = atom->type; int nlocal = atom->nlocal; - int itype; + // sum rotational energy for each particle + // no point particles since divide by inertia + double wbody[3]; double rot[3][3]; double erotate = 0.0; - for (int i = 0; i < nlocal; i++) + for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { itype = type[i]; diff -Naur lammps-27Jun09/src/ASPHERE/compute_temp_asphere.cpp lammps-28Jun09/src/ASPHERE/compute_temp_asphere.cpp --- lammps-27Jun09/src/ASPHERE/compute_temp_asphere.cpp 2009-01-05 15:26:08.000000000 -0700 +++ lammps-28Jun09/src/ASPHERE/compute_temp_asphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -19,6 +19,7 @@ #include "compute_temp_asphere.h" #include "math_extra.h" #include "atom.h" +#include "atom_vec.h" #include "update.h" #include "force.h" #include "domain.h" @@ -38,9 +39,6 @@ if (narg != 3 && narg != 4) error->all("Illegal compute temp/asphere command"); - if (!atom->quat_flag || !atom->angmom_flag) - error->all("Compute temp/asphere requires atom attributes quat, angmom"); - scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; @@ -59,6 +57,16 @@ vector = new double[6]; inertia = memory->create_2d_double_array(atom->ntypes+1,3,"fix_temp_sphere:inertia"); + + // error checks + + if (!atom->angmom_flag || !atom->quat_flag || + !atom->avec->shape_type) + error->all("Compute temp/asphere requires atom attributes " + "angmom, quat, shape"); + if (atom->radius_flag || atom->rmass_flag) + error->all("Compute temp/asphere cannot be used with atom attributes " + "diameter or rmass"); } /* ---------------------------------------------------------------------- */ @@ -74,6 +82,20 @@ void ComputeTempAsphere::init() { + // check that all particles are finite-size + // no point particles allowed, spherical is OK + + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (shape[type[i]][0] == 0.0) + error->one("Compue temp/asphere requires extended particles"); + if (tempbias) { int i = modify->find_compute(id_bias); if (i < 0) error->all("Could not find compute ID for temperature bias"); @@ -101,44 +123,30 @@ void ComputeTempAsphere::dof_compute() { - double natoms = group->count(igroup); - int dimension = domain->dimension; - dof = dimension * natoms; - - if (tempbias == 1) dof -= tbias->dof_remove(-1) * natoms; - - // rotational degrees of freedom - // 0 for sphere, 2 for uniaxial, 3 for biaxial + // 6 dof for 3d, 3 dof for 2d + // assume full rotation of extended particles + // user can correct this via compute_modify if needed - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; + double natoms = group->count(igroup); + if (domain->dimension == 3) dof = 6*natoms; + else dof = 3*natoms; - int itype; - int count = 0; + // additional adjustments to dof - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (tempbias == 2 && tbias->dof_remove(i)) continue; - itype = type[i]; - if (dimension == 2) { - if (shape[itype][0] == shape[itype][1]) continue; - else count++; - } else { - if (shape[itype][0] == shape[itype][1] && - shape[itype][1] == shape[itype][2]) continue; - else if (shape[itype][0] == shape[itype][1] || - shape[itype][1] == shape[itype][2] || - shape[itype][0] == shape[itype][2]) count += 2; - else count += 3; - } - } + if (tempbias == 1) dof -= tbias->dof_remove(-1) * natoms; + else if (tempbias == 2) { + int *mask = atom->mask; + int nlocal = atom->nlocal; + int count = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (tbias->dof_remove(i)) count++; + int count_all; + MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); + if (domain->dimension == 3) dof -= 6*count_all; + else dof -= 3*count_all; + } - int count_all; - MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); - dof += count_all; - dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; @@ -169,31 +177,26 @@ double rot[3][3]; double t = 0.0; + // sum translationals and rotational energy for each particle + // no point particles since divide by inertia + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - // translational kinetic energy - itype = type[i]; t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[itype]; // wbody = angular velocity in body frame + + MathExtra::quat_to_mat(quat[i],rot); + MathExtra::transpose_times_column3(rot,angmom[i],wbody); + wbody[0] /= inertia[itype][0]; + wbody[1] /= inertia[itype][1]; + wbody[2] /= inertia[itype][2]; - if (!(shape[itype][0] == shape[itype][1] && - shape[itype][1] == shape[itype][2])) { - - MathExtra::quat_to_mat(quat[i],rot); - MathExtra::transpose_times_column3(rot,angmom[i],wbody); - wbody[0] /= inertia[itype][0]; - wbody[1] /= inertia[itype][1]; - wbody[2] /= inertia[itype][2]; - - // rotational kinetic energy - - t += inertia[itype][0]*wbody[0]*wbody[0]+ - inertia[itype][1]*wbody[1]*wbody[1]+ - inertia[itype][2]*wbody[2]*wbody[2]; - } + t += inertia[itype][0]*wbody[0]*wbody[0] + + inertia[itype][1]*wbody[1]*wbody[1] + + inertia[itype][2]*wbody[2]*wbody[2]; } if (tempbias) tbias->restore_bias_all(); diff -Naur lammps-27Jun09/src/ASPHERE/fix_npt_asphere.cpp lammps-28Jun09/src/ASPHERE/fix_npt_asphere.cpp --- lammps-27Jun09/src/ASPHERE/fix_npt_asphere.cpp 2008-12-16 12:46:30.000000000 -0700 +++ lammps-28Jun09/src/ASPHERE/fix_npt_asphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -27,6 +27,7 @@ #include "kspace.h" #include "update.h" #include "domain.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -38,23 +39,56 @@ FixNPTAsphere::FixNPTAsphere(LAMMPS *lmp, int narg, char **arg) : FixNPT(lmp, narg, arg) { + inertia = + memory->create_2d_double_array(atom->ntypes+1,3,"fix_npt_asphere:inertia"); + + // error checks + if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag || - !atom->shape) + !atom->avec->shape_type) error->all("Fix npt/asphere requires atom attributes " "quat, angmom, torque, shape"); + if (atom->radius_flag || atom->rmass_flag) + error->all("Fix npt/asphere cannot be used with atom attributes " + "diameter or rmass"); } -/* ---------------------------------------------------------------------- - 1st half of Verlet update -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ + +FixNPTAsphere::~FixNPTAsphere() +{ + memory->destroy_2d_double_array(inertia); +} + +/* ---------------------------------------------------------------------- */ + +void FixNPTAsphere::init() +{ + // check that all particles are finite-size + // no point particles allowed, spherical is OK + + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (shape[type[i]][0] == 0.0) + error->one("Fix nvt/asphere requires extended particles"); + + FixNPT::init(); + calculate_inertia(); +} + +/* ---------------------------------------------------------------------- */ void FixNPTAsphere::initial_integrate(int vflag) { int i; double dtfm; - dtq = 0.5 * dtv; - double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; @@ -85,8 +119,6 @@ } factor_rotate = exp(-dthalf*eta_dot); - // update v of only atoms in group - double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -135,7 +167,11 @@ } } - // update angular momentum by 1/2 step + // set timestep here since dt may have changed or come via rRESPA + + dtq = 0.5 * dtv; + + // update angular momentum by 1/2 step for all particles // update quaternion a full step via Richardson iteration // returns new normalized quaternion @@ -145,9 +181,7 @@ angmom[i][1] = angmom[i][1]*factor_rotate + dtf*torque[i][1]; angmom[i][2] = angmom[i][2]*factor_rotate + dtf*torque[i][2]; - double inertia[3]; - calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia); - richardson(quat[i],angmom[i],inertia); + richardson(quat[i],angmom[i],inertia[type[i]]); } } @@ -158,9 +192,7 @@ if (kspace_flag) force->kspace->setup(); } -/* ---------------------------------------------------------------------- - 2nd half of Verlet update -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ void FixNPTAsphere::final_integrate() { @@ -320,15 +352,22 @@ MathExtra::times_column3(rot,wbody,w); } + /* ---------------------------------------------------------------------- - calculate the moment of inertia for an ellipsoid, from mass and radii - shape = x,y,z radii in body frame + principal moments of inertia for ellipsoids ------------------------------------------------------------------------- */ -void FixNPTAsphere::calculate_inertia(double mass, double *shape, - double *inertia) +void FixNPTAsphere::calculate_inertia() { - inertia[0] = 0.2*mass*(shape[1]*shape[1]+shape[2]*shape[2]); - inertia[1] = 0.2*mass*(shape[0]*shape[0]+shape[2]*shape[2]); - inertia[2] = 0.2*mass*(shape[0]*shape[0]+shape[1]*shape[1]); + double *mass = atom->mass; + double **shape = atom->shape; + + for (int i = 1; i <= atom->ntypes; i++) { + inertia[i][0] = 0.2*mass[i] * + (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]); + inertia[i][1] = 0.2*mass[i] * + (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]); + inertia[i][2] = 0.2*mass[i] * + (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]); + } } diff -Naur lammps-27Jun09/src/ASPHERE/fix_npt_asphere.h lammps-28Jun09/src/ASPHERE/fix_npt_asphere.h --- lammps-27Jun09/src/ASPHERE/fix_npt_asphere.h 2008-04-08 16:55:40.000000000 -0600 +++ lammps-28Jun09/src/ASPHERE/fix_npt_asphere.h 2009-06-26 12:23:16.000000000 -0600 @@ -21,16 +21,19 @@ class FixNPTAsphere : public FixNPT { public: FixNPTAsphere(class LAMMPS *, int, char **); + ~FixNPTAsphere(); + void init(); void initial_integrate(int); void final_integrate(); private: double dtq; double factor_rotate; + double **inertia; void richardson(double *, double *, double *); void omega_from_mq(double *, double *, double *, double *); - void calculate_inertia(double mass, double *shape, double *inertia); + void calculate_inertia(); }; } diff -Naur lammps-27Jun09/src/ASPHERE/fix_nve_asphere.cpp lammps-28Jun09/src/ASPHERE/fix_nve_asphere.cpp --- lammps-27Jun09/src/ASPHERE/fix_nve_asphere.cpp 2008-12-16 12:46:30.000000000 -0700 +++ lammps-28Jun09/src/ASPHERE/fix_nve_asphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -37,12 +37,18 @@ FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) : FixNVE(lmp, narg, arg) { - if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag || - !atom->shape) - error->all("Fix nve/asphere requires atom attributes " - "quat, angmom, torque, shape"); inertia = - memory->create_2d_double_array(atom->ntypes+1,3,"fix_temp_sphere:inertia"); + memory->create_2d_double_array(atom->ntypes+1,3,"fix_nve_asphere:inertia"); + + // error checks + + if (!atom->angmom_flag || !atom->quat_flag || !atom->torque_flag || + !atom->avec->shape_type) + error->all("Fix nve/asphere requires atom attributes " + "angmom, quat, torque, shape"); + if (atom->radius_flag || atom->rmass_flag) + error->all("Fix nve/asphere cannot be used with atom attributes " + "diameter or rmass"); } /* ---------------------------------------------------------------------- */ @@ -56,6 +62,20 @@ void FixNVEAsphere::init() { + // check that all particles are finite-size + // no point particles allowed, spherical is OK + + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (shape[type[i]][0] == 0.0) + error->one("Fix nve/asphere requires extended particles"); + FixNVE::init(); calculate_inertia(); } @@ -66,8 +86,6 @@ { double dtfm; - dtq = 0.5 * dtv; - double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -80,6 +98,10 @@ int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; + // set timestep here since dt may have changed or come via rRESPA + + dtq = 0.5 * dtv; + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { dtfm = dtf / mass[type[i]]; @@ -216,7 +238,7 @@ { double *mass = atom->mass; double **shape = atom->shape; - + for (int i = 1; i <= atom->ntypes; i++) { inertia[i][0] = 0.2*mass[i] * (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]); diff -Naur lammps-27Jun09/src/ASPHERE/fix_nvt_asphere.cpp lammps-28Jun09/src/ASPHERE/fix_nvt_asphere.cpp --- lammps-27Jun09/src/ASPHERE/fix_nvt_asphere.cpp 2008-12-16 12:46:30.000000000 -0700 +++ lammps-28Jun09/src/ASPHERE/fix_nvt_asphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -28,6 +28,7 @@ #include "update.h" #include "modify.h" #include "compute.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -39,10 +40,47 @@ FixNVTAsphere::FixNVTAsphere(LAMMPS *lmp, int narg, char **arg) : FixNVT(lmp, narg, arg) { + inertia = + memory->create_2d_double_array(atom->ntypes+1,3,"fix_nvt_asphere:inertia"); + + // error checks + if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag || - !atom->shape) + !atom->avec->shape_type) error->all("Fix nvt/asphere requires atom attributes " "quat, angmom, torque, shape"); + if (atom->radius_flag || atom->rmass_flag) + error->all("Fix nvt/asphere cannot be used with atom attributes " + "diameter or rmass"); +} + +/* ---------------------------------------------------------------------- */ + +FixNVTAsphere::~FixNVTAsphere() +{ + memory->destroy_2d_double_array(inertia); +} + +/* ---------------------------------------------------------------------- */ + +void FixNVTAsphere::init() +{ + // check that all particles are finite-size + // no point particles allowed, spherical is OK + + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (shape[type[i]][0] == 0.0) + error->one("Fix nvt/asphere requires extended particles"); + + FixNVT::init(); + calculate_inertia(); } /* ---------------------------------------------------------------------- */ @@ -51,8 +89,6 @@ { double dtfm; - dtq = 0.5 * dtv; - double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; t_target = t_start + delta * (t_stop-t_start); @@ -65,8 +101,6 @@ eta += dtv*eta_dot; factor = exp(-dthalf*eta_dot); - // update v and x of only atoms in group - double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -79,6 +113,10 @@ int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; + // set timestep here since dt may have changed or come via rRESPA + + dtq = 0.5 * dtv; + if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -98,9 +136,7 @@ angmom[i][1] = angmom[i][1]*factor + dtf*torque[i][1]; angmom[i][2] = angmom[i][2]*factor + dtf*torque[i][2]; - double inertia[3]; - calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia); - richardson(quat[i],angmom[i],inertia); + richardson(quat[i],angmom[i],inertia[type[i]]); } } @@ -125,9 +161,7 @@ angmom[i][1] = angmom[i][1]*factor + dtf*torque[i][1]; angmom[i][2] = angmom[i][2]*factor + dtf*torque[i][2]; - double inertia[3]; - calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia); - richardson(quat[i],angmom[i],inertia); + richardson(quat[i],angmom[i],inertia[type[i]]); } } } @@ -139,8 +173,6 @@ { double dtfm; - // update v of only atoms in group - double **v = atom->v; double **f = atom->f; double **angmom = atom->angmom; @@ -269,14 +301,20 @@ } /* ---------------------------------------------------------------------- - calculate the moment of inertia for an ELLIPSOID, from mass and radii - shape = x,y,z radii in body frame + principal moments of inertia for ellipsoids ------------------------------------------------------------------------- */ -void FixNVTAsphere::calculate_inertia(double mass, double *shape, - double *inertia) +void FixNVTAsphere::calculate_inertia() { - inertia[0] = 0.2*mass*(shape[1]*shape[1]+shape[2]*shape[2]); - inertia[1] = 0.2*mass*(shape[0]*shape[0]+shape[2]*shape[2]); - inertia[2] = 0.2*mass*(shape[0]*shape[0]+shape[1]*shape[1]); + double *mass = atom->mass; + double **shape = atom->shape; + + for (int i = 1; i <= atom->ntypes; i++) { + inertia[i][0] = 0.2*mass[i] * + (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]); + inertia[i][1] = 0.2*mass[i] * + (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]); + inertia[i][2] = 0.2*mass[i] * + (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]); + } } diff -Naur lammps-27Jun09/src/ASPHERE/fix_nvt_asphere.h lammps-28Jun09/src/ASPHERE/fix_nvt_asphere.h --- lammps-27Jun09/src/ASPHERE/fix_nvt_asphere.h 2008-04-08 16:55:40.000000000 -0600 +++ lammps-28Jun09/src/ASPHERE/fix_nvt_asphere.h 2009-06-26 12:23:16.000000000 -0600 @@ -21,15 +21,18 @@ class FixNVTAsphere : public FixNVT { public: FixNVTAsphere(class LAMMPS *, int, char **); + ~FixNVTAsphere(); + void init(); void initial_integrate(int); void final_integrate(); private: double dtq; + double **inertia; void richardson(double *, double *, double *); void omega_from_mq(double *, double *, double *, double *); - void calculate_inertia(double mass, double *shape, double *inertia); + void calculate_inertia(); }; } diff -Naur lammps-27Jun09/src/ASPHERE/pair_gayberne.cpp lammps-28Jun09/src/ASPHERE/pair_gayberne.cpp --- lammps-27Jun09/src/ASPHERE/pair_gayberne.cpp 2008-08-08 08:41:33.000000000 -0600 +++ lammps-28Jun09/src/ASPHERE/pair_gayberne.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -22,6 +22,7 @@ #include "pair_gayberne.h" #include "math_extra.h" #include "atom.h" +#include "atom_vec.h" #include "comm.h" #include "force.h" #include "neighbor.h" @@ -334,8 +335,10 @@ void PairGayBerne::init_style() { - if (!atom->quat_flag || !atom->torque_flag) - error->all("Pair gayberne requires atom attributes quat, torque"); + if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type) + error->all("Pair gayberne requires atom attributes quat, torque, shape"); + if (atom->radius_flag) + error->all("Pair gayberne cannot be used with atom attribute diameter"); int irequest = neighbor->request(this); diff -Naur lammps-27Jun09/src/ASPHERE/pair_resquared.cpp lammps-28Jun09/src/ASPHERE/pair_resquared.cpp --- lammps-27Jun09/src/ASPHERE/pair_resquared.cpp 2009-01-27 08:33:31.000000000 -0700 +++ lammps-28Jun09/src/ASPHERE/pair_resquared.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -22,6 +22,7 @@ #include "pair_resquared.h" #include "math_extra.h" #include "atom.h" +#include "atom_vec.h" #include "comm.h" #include "force.h" #include "neighbor.h" @@ -323,8 +324,10 @@ void PairRESquared::init_style() { - if (!atom->quat_flag || !atom->torque_flag) - error->all("Pair resquared requires atom attributes quat, torque"); + if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type) + error->all("Pair resquared requires atom attributes quat, torque, shape"); + if (atom->radius_flag) + error->all("Pair resquared cannot be used with atom attribute diameter"); int irequest = neighbor->request(this); @@ -353,9 +356,11 @@ error->all("Pair resquared epsilon a,b,c coeffs are not all set"); int ishape = 0; - if (shape[i][0] != 0 && shape[i][1] != 0 && shape[i][2] != 0) ishape = 1; + if (shape[i][0] != 0.0 && shape[i][1] != 0.0 && shape[i][2] != 0.0) + ishape = 1; int jshape = 0; - if (shape[j][0] != 0 && shape[j][1] != 0 && shape[j][2] != 0) jshape = 1; + if (shape[j][0] != 0.0 && shape[j][1] != 0.0 && shape[j][2] != 0.0) + jshape = 1; if (ishape == 0 && jshape == 0) { form[i][j] = SPHERE_SPHERE; diff -Naur lammps-27Jun09/src/COLLOID/pair_lubricate.cpp lammps-28Jun09/src/COLLOID/pair_lubricate.cpp --- lammps-27Jun09/src/COLLOID/pair_lubricate.cpp 2009-01-09 17:09:53.000000000 -0700 +++ lammps-28Jun09/src/COLLOID/pair_lubricate.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -21,6 +21,7 @@ #include "string.h" #include "pair_lubricate.h" #include "atom.h" +#include "atom_vec.h" #include "comm.h" #include "force.h" #include "neighbor.h" @@ -56,6 +57,7 @@ memory->destroy_2d_double_array(cut); memory->destroy_2d_double_array(cut_inner); } + delete random; } @@ -99,10 +101,10 @@ numneigh = list->numneigh; firstneigh = list->firstneigh; - a_squeeze = a_shear = a_pump = a_twist = 0.0; - // loop over neighbors of my atoms + a_squeeze = a_shear = a_pump = a_twist = 0.0; + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; @@ -388,17 +390,21 @@ void PairLubricate::init_style() { - if (!atom->torque_flag || atom->shape == NULL) - error->all("Pair lubricate requires atom attributes torque, shape"); + if (!atom->quat_flag || !atom->torque_flag || !atom->avec->shape_type) + error->all("Pair lubricate requires atom attributes quat, torque, shape"); + if (atom->radius_flag || atom->rmass_flag) + error->all("Pair lubricate cannot be used with atom attributes" + "diameter or rmass"); - // insure all particle shapes are spherical + // insure all particle shapes are finite-size, spherical, and monodisperse + double value = atom->shape[1][0]; + if (value == 0.0) error->all("Pair lubricate requires extended particles"); for (int i = 1; i <= atom->ntypes; i++) - if ((atom->shape[i][0] != atom->shape[i][1]) || - (atom->shape[i][0] != atom->shape[i][2]) || - (atom->shape[i][1] != atom->shape[i][2]) ) - error->all("Pair lubricate requires spherical particles"); - + if (atom->shape[i][0] != value || atom->shape[i][0] != value || + atom->shape[i][0] != value) + error->all("Pair lubricate requires spherical, mono-disperse particles"); + int irequest = neighbor->request(this); } diff -Naur lammps-27Jun09/src/GRANULAR/atom_vec_granular.cpp lammps-28Jun09/src/GRANULAR/atom_vec_granular.cpp --- lammps-27Jun09/src/GRANULAR/atom_vec_granular.cpp 2009-02-12 14:36:52.000000000 -0700 +++ lammps-28Jun09/src/GRANULAR/atom_vec_granular.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -533,8 +533,11 @@ radius[nlocal] = buf[m++]; density[nlocal] = buf[m++]; - rmass[nlocal] = 4.0*PI/3.0 * - radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; + + if (radius[nlocal] == 0.0) rmass[nlocal] = density[nlocal]; + else + rmass[nlocal] = 4.0*PI/3.0 * + radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; omega[nlocal][0] = buf[m++]; omega[nlocal][1] = buf[m++]; @@ -601,10 +604,17 @@ error->one("Invalid atom type in Atoms section of data file"); radius[nlocal] = 0.5 * atof(values[2]); + if (radius[nlocal] < 0.0) + error->one("Invalid radius in Atoms section of data file"); + density[nlocal] = atof(values[3]); - rmass[nlocal] = 4.0*PI/3.0 * - radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; - if (rmass[nlocal] <= 0.0) error->one("Invalid mass value"); + if (density[nlocal] <= 0.0) + error->one("Invalid density in Atoms section of data file"); + + if (radius[nlocal] == 0.0) rmass[nlocal] = density[nlocal]; + else + rmass[nlocal] = 4.0*PI/3.0 * + radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; x[nlocal][0] = coord[0]; x[nlocal][1] = coord[1]; @@ -631,10 +641,17 @@ int AtomVecGranular::data_atom_hybrid(int nlocal, char **values) { radius[nlocal] = 0.5 * atof(values[0]); + if (radius[nlocal] < 0.0) + error->one("Invalid radius in Atoms section of data file"); + density[nlocal] = atof(values[1]); - rmass[nlocal] = 4.0*PI/3.0 * - radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; - if (rmass[nlocal] <= 0.0) error->one("Invalid mass value"); + if (density[nlocal] <= 0.0) + error->one("Invalid density in Atoms section of data file"); + + if (radius[nlocal] == 0.0) rmass[nlocal] = density[nlocal]; + else + rmass[nlocal] = 4.0*PI/3.0 * + radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; diff -Naur lammps-27Jun09/src/GRANULAR/pair_gran_hertz_history.cpp lammps-28Jun09/src/GRANULAR/pair_gran_hertz_history.cpp --- lammps-27Jun09/src/GRANULAR/pair_gran_hertz_history.cpp 2009-01-08 14:46:20.000000000 -0700 +++ lammps-28Jun09/src/GRANULAR/pair_gran_hertz_history.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -43,7 +43,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum; + int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; double radi,radj,radsum,rsq,r,rinv,rsqinv; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; @@ -66,6 +66,8 @@ double **torque = atom->torque; double *radius = atom->radius; double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -141,9 +143,18 @@ // normal force = Hertzian contact + normal velocity damping - meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); - if (mask[i] & freeze_group_bit) meff = rmass[j]; - if (mask[j] & freeze_group_bit) meff = rmass[i]; + if (rmass) { + meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); + if (mask[i] & freeze_group_bit) meff = rmass[j]; + if (mask[j] & freeze_group_bit) meff = rmass[i]; + } else { + itype = type[i]; + jtype = type[j]; + meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]); + if (mask[i] & freeze_group_bit) meff = mass[jtype]; + if (mask[j] & freeze_group_bit) meff = mass[itype]; + } + damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); diff -Naur lammps-27Jun09/src/GRANULAR/pair_gran_hooke.cpp lammps-28Jun09/src/GRANULAR/pair_gran_hooke.cpp --- lammps-27Jun09/src/GRANULAR/pair_gran_hooke.cpp 2009-01-08 14:46:20.000000000 -0700 +++ lammps-28Jun09/src/GRANULAR/pair_gran_hooke.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -39,7 +39,7 @@ void PairGranHooke::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum; + int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; double radi,radj,radsum,rsq,r,rinv,rsqinv; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; @@ -59,6 +59,8 @@ double **torque = atom->torque; double *radius = atom->radius; double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; int newton_pair = force->newton_pair; @@ -121,9 +123,18 @@ // normal forces = Hookian contact + normal velocity damping - meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); - if (mask[i] & freeze_group_bit) meff = rmass[j]; - if (mask[j] & freeze_group_bit) meff = rmass[i]; + if (rmass) { + meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); + if (mask[i] & freeze_group_bit) meff = rmass[j]; + if (mask[j] & freeze_group_bit) meff = rmass[i]; + } else { + itype = type[i]; + jtype = type[j]; + meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]); + if (mask[i] & freeze_group_bit) meff = mass[jtype]; + if (mask[j] & freeze_group_bit) meff = mass[itype]; + } + damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; diff -Naur lammps-27Jun09/src/GRANULAR/pair_gran_hooke_history.cpp lammps-28Jun09/src/GRANULAR/pair_gran_hooke_history.cpp --- lammps-27Jun09/src/GRANULAR/pair_gran_hooke_history.cpp 2009-01-08 14:46:20.000000000 -0700 +++ lammps-28Jun09/src/GRANULAR/pair_gran_hooke_history.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -67,7 +67,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum; + int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; double radi,radj,radsum,rsq,r,rinv,rsqinv; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; @@ -90,6 +90,8 @@ double **torque = atom->torque; double *radius = atom->radius; double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -165,9 +167,18 @@ // normal forces = Hookian contact + normal velocity damping - meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); - if (mask[i] & freeze_group_bit) meff = rmass[j]; - if (mask[j] & freeze_group_bit) meff = rmass[i]; + if (rmass) { + meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); + if (mask[i] & freeze_group_bit) meff = rmass[j]; + if (mask[j] & freeze_group_bit) meff = rmass[i]; + } else { + itype = type[i]; + jtype = type[j]; + meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]); + if (mask[i] & freeze_group_bit) meff = mass[jtype]; + if (mask[j] & freeze_group_bit) meff = mass[itype]; + } + damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; @@ -328,6 +339,8 @@ { int i; + // error and warning checks + if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag) error->all("Pair granular requires atom attributes radius, omega, torque"); if (atom->avec->ghost_velocity == 0) diff -Naur lammps-27Jun09/src/MOLECULE/fix_bond_swap.cpp lammps-28Jun09/src/MOLECULE/fix_bond_swap.cpp --- lammps-27Jun09/src/MOLECULE/fix_bond_swap.cpp 2009-01-06 10:16:16.000000000 -0700 +++ lammps-28Jun09/src/MOLECULE/fix_bond_swap.cpp 2009-06-24 10:10:08.000000000 -0600 @@ -124,6 +124,9 @@ if (force->pair == NULL || force->bond == NULL) error->all("Fix bond/swap requires pair and bond styles"); + if (force->pair->single_enable == 0) + error->all("Pair style does not support fix bond/swap"); + if (force->angle == NULL && atom->nangles > 0 && comm->me == 0) error->warning("Fix bond/swap will ignore defined angles"); diff -Naur lammps-27Jun09/src/atom.cpp lammps-28Jun09/src/atom.cpp --- lammps-27Jun09/src/atom.cpp 2009-05-19 09:01:02.000000000 -0600 +++ lammps-28Jun09/src/atom.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -1124,17 +1124,15 @@ shape[itype][2] = 0.5*c; shape_setflag[itype] = 1; - if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || shape[itype][2] < 0.0) - error->all("Invalid shape value"); - if (shape[itype][0] == 0.0 && - (shape[itype][1] != 0.0 || shape[itype][2] != 0.0)) - error->all("Invalid shape value"); - if (shape[itype][1] == 0.0 && - (shape[itype][0] != 0.0 || shape[itype][2] != 0.0)) - error->all("Invalid shape value"); - if (shape[itype][2] == 0.0 && - (shape[itype][0] != 0.0 || shape[itype][1] != 0.0)) + if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || + shape[itype][2] < 0.0) error->all("Invalid shape value"); + if (shape[itype][0] > 0.0 || shape[itype][1] > 0.0 || + shape[itype][2] > 0.0) { + if (shape[itype][0] == 0.0 || shape[itype][1] == 0.0 || + shape[itype][2] == 0.0) + error->all("Invalid shape value"); + } } /* ---------------------------------------------------------------------- @@ -1159,18 +1157,15 @@ shape[itype][2] = 0.5*atof(arg[3]); shape_setflag[itype] = 1; - if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || + if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || shape[itype][2] < 0.0) error->all("Invalid shape value"); - if (shape[itype][0] == 0.0 && - (shape[itype][1] != 0.0 || shape[itype][2] != 0.0)) - error->all("Invalid shape value"); - if (shape[itype][1] == 0.0 && - (shape[itype][0] != 0.0 || shape[itype][2] != 0.0)) - error->all("Invalid shape value"); - if (shape[itype][2] == 0.0 && - (shape[itype][0] != 0.0 || shape[itype][1] != 0.0)) - error->all("Invalid shape value"); + if (shape[itype][0] > 0.0 || shape[itype][1] > 0.0 || + shape[itype][2] > 0.0) { + if (shape[itype][0] == 0.0 || shape[itype][1] == 0.0 || + shape[itype][2] == 0.0) + error->all("Invalid shape value"); + } } } diff -Naur lammps-27Jun09/src/compute_erotate_sphere.cpp lammps-28Jun09/src/compute_erotate_sphere.cpp --- lammps-27Jun09/src/compute_erotate_sphere.cpp 2009-02-12 14:36:52.000000000 -0700 +++ lammps-28Jun09/src/compute_erotate_sphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -14,6 +14,7 @@ #include "mpi.h" #include "compute_erotate_sphere.h" #include "atom.h" +#include "atom_vec.h" #include "update.h" #include "force.h" #include "domain.h" @@ -31,71 +32,100 @@ { if (narg != 3) error->all("Illegal compute erotate/sphere command"); - if (!atom->omega_flag) - error->all("Compute erotate/sphere requires atom attribute omega"); - scalar_flag = 1; extscalar = 1; - inertia = new double[atom->ntypes+1]; -} + // error checks -/* ---------------------------------------------------------------------- */ - -ComputeERotateSphere::~ComputeERotateSphere() -{ - delete [] inertia; + if (!atom->omega_flag) + error->all("Compute erotate/sphere requires atom attribute omega"); + if (!atom->radius_flag && !atom->avec->shape_type) + error->all("Compute erotate/sphere requires atom attribute " + "radius or shape"); } /* ---------------------------------------------------------------------- */ void ComputeERotateSphere::init() { - pfactor = 0.5 * force->mvv2e * INERTIA; + int i,itype; - if (atom->mass && !atom->shape) - error->all("Compute erotate/sphere requires atom attribute shape"); - if (atom->rmass && !atom->radius_flag) - error->all("Compute erotate/sphere requires atom attribute radius"); + // if shape used, check that all particles are spherical + // point particles are allowed - if (atom->mass) { - double *mass = atom->mass; + if (atom->radius == NULL) { double **shape = atom->shape; - - for (int i = 1; i <= atom->ntypes; i++) { - if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) - error->all("Compute erotate/sphere requires spherical particle shapes"); - inertia[i] = shape[i][0]*shape[i][0] * mass[i]; - } + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + if (shape[itype][0] != shape[itype][1] || + shape[itype][0] != shape[itype][2]) + error->one("Compute erotate/sphere requires " + "spherical particle shapes"); + } } + + pfactor = 0.5 * force->mvv2e * INERTIA; } /* ---------------------------------------------------------------------- */ double ComputeERotateSphere::compute_scalar() { + int i,itype; + invoked_scalar = update->ntimestep; double **omega = atom->omega; double *radius = atom->radius; double *rmass = atom->rmass; double *mass = atom->mass; + double **shape = atom->shape; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; + // sum rotational energy for each particle + // point particles will not contribute due to radius or shape = 0 + double erotate = 0.0; - if (rmass) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i]; + if (radius) { + if (rmass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i]; + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * + radius[i]*radius[i]*mass[itype]; + } + } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * inertia[type[i]]; + if (rmass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * + shape[itype][0]*shape[itype][0]*rmass[i]; + } + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * + shape[itype][0]*shape[itype][0]*mass[itype]; + } + } } MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world); diff -Naur lammps-27Jun09/src/compute_erotate_sphere.h lammps-28Jun09/src/compute_erotate_sphere.h --- lammps-27Jun09/src/compute_erotate_sphere.h 2008-03-17 17:37:19.000000000 -0600 +++ lammps-28Jun09/src/compute_erotate_sphere.h 2009-06-26 12:23:16.000000000 -0600 @@ -21,13 +21,12 @@ class ComputeERotateSphere : public Compute { public: ComputeERotateSphere(class LAMMPS *, int, char **); - ~ComputeERotateSphere(); + ~ComputeERotateSphere() {} void init(); double compute_scalar(); private: double pfactor; - double *inertia; }; } diff -Naur lammps-27Jun09/src/compute_group_group.cpp lammps-28Jun09/src/compute_group_group.cpp --- lammps-27Jun09/src/compute_group_group.cpp 2009-03-17 09:40:57.000000000 -0600 +++ lammps-28Jun09/src/compute_group_group.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -12,7 +12,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing authors: Naveen Michaud-Agrawal (Johns Hopkins U) + Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U) ------------------------------------------------------------------------- */ #include "mpi.h" @@ -65,8 +65,9 @@ void ComputeGroupGroup::init() { - if (force->pair == NULL) - error->all("Compute group/group requires a pair style be defined"); + if (force->pair == NULL || force->pair->single_enable == 0) + error->all("Pair style does not support compute group/group"); + pair = force->pair; cutsq = force->pair->cutsq; diff -Naur lammps-27Jun09/src/compute_temp_sphere.cpp lammps-28Jun09/src/compute_temp_sphere.cpp --- lammps-27Jun09/src/compute_temp_sphere.cpp 2009-02-12 14:36:52.000000000 -0700 +++ lammps-28Jun09/src/compute_temp_sphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -15,6 +15,7 @@ #include "string.h" #include "compute_temp_sphere.h" #include "atom.h" +#include "atom_vec.h" #include "update.h" #include "force.h" #include "domain.h" @@ -35,9 +36,6 @@ if (narg != 3 && narg != 4) error->all("Illegal compute temp/sphere command"); - if (!atom->omega_flag) - error->all("Compute temp/sphere requires atom attribute omega"); - scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; @@ -54,23 +52,51 @@ } vector = new double[6]; - inertia = new double[atom->ntypes+1]; + + // error checks + + if (!atom->omega_flag) + error->all("Compute temp/sphere requires atom attribute omega"); + if (!atom->radius_flag && !atom->avec->shape_type) + error->all("Compute temp/sphere requires atom attribute " + "radius or shape"); } /* ---------------------------------------------------------------------- */ ComputeTempSphere::~ComputeTempSphere() { + delete [] id_bias; delete [] vector; - delete [] inertia; } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::init() { + int i,itype; + + // if shape used, check that all particles are spherical + // point particles are allowed + + if (atom->radius == NULL) { + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + if (shape[itype][0] != shape[itype][1] || + shape[itype][0] != shape[itype][2]) + error->one("Compute temp/sphere requires " + "spherical particle shapes"); + } + } + if (tempbias) { - int i = modify->find_compute(id_bias); + i = modify->find_compute(id_bias); if (i < 0) error->all("Could not find compute ID for temperature bias"); tbias = modify->compute[i]; if (tbias->tempflag == 0) @@ -85,44 +111,119 @@ } fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) + for (i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); dof_compute(); - - if (atom->mass) { - double *mass = atom->mass; - double **shape = atom->shape; - - for (int i = 1; i <= atom->ntypes; i++) { - if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) - error->all("Compute temp/sphere requires spherical particle shapes"); - inertia[i] = INERTIA * shape[i][0]*shape[i][0] * mass[i]; - } - } } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::dof_compute() { - double natoms = group->count(igroup); - int nper = 6; - if (domain->dimension == 2) nper = 3; - dof = nper * natoms; + int count,count_all; - if (tempbias) { - if (tempbias == 1) dof -= tbias->dof_remove(-1) * natoms; - else { - int *mask = atom->mask; - int nlocal = atom->nlocal; - int count = 0; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - if (tbias->dof_remove(i)) count++; - int count_all; - MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); - dof -= nper * count_all; + // 6 or 3 dof for extended/point particles for 3d + // 3 or 2 dof for extended/point particles for 2d + // assume full rotation of extended particles + // user can correct this via compute_modify if needed + + int dimension = domain->dimension; + + double *radius = atom->radius; + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + count = 0; + if (dimension == 3) { + if (radius) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (radius[i] == 0.0) count += 3; + else count += 6; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (shape[type[i]][0] == 0.0) count += 3; + else count += 6; + } + } } + } else { + if (radius) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (radius[i] == 0.0) count += 2; + else count += 3; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (shape[type[i]][0] == 0.0) count += 2; + else count += 3; + } + } + } + } + + MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); + dof = count_all; + + // additional adjustments to dof + + if (tempbias == 1) { + double natoms = group->count(igroup); + dof -= tbias->dof_remove(-1) * natoms; + + } else if (tempbias == 2) { + int *mask = atom->mask; + int nlocal = atom->nlocal; + + count = 0; + if (dimension == 3) { + if (radius) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (tbias->dof_remove(i)) { + if (radius[i] == 0.0) count += 3; + else count += 6; + } + } + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (tbias->dof_remove(i)) { + if (shape[type[i]][0] == 0.0) count += 3; + else count += 6; + } + } + } + } else { + if (radius) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (tbias->dof_remove(i)) { + if (radius[i] == 0.0) count += 2; + else count += 3; + } + } + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (tbias->dof_remove(i)) { + if (shape[type[i]][0] == 0.0) count += 2; + else count += 3; + } + } + } + } + + MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); + dof -= count_all; } dof -= extra_dof + fix_dof; @@ -134,6 +235,8 @@ double ComputeTempSphere::compute_scalar() { + int i,itype; + invoked_scalar = update->ntimestep; if (tempbias) { @@ -143,30 +246,65 @@ double **v = atom->v; double **omega = atom->omega; - double *mass = atom->mass; - double *rmass = atom->rmass; double *radius = atom->radius; + double *rmass = atom->rmass; + double *mass = atom->mass; + double **shape = atom->shape; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + // 4 cases depending on radius vs shape and rmass vs mass + // point particles will not contribute rotation due to radius or shape = 0 + double t = 0.0; - if (rmass) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; - t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i]; - } + if (radius) { + if (rmass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + rmass[i]; + t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * + INERTIA*radius[i]*radius[i]*rmass[i]; + } + + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + mass[itype]; + t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * + INERTIA*radius[i]*radius[i]*mass[itype]; + } + } + } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - mass[type[i]]; - t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + - omega[i][2]*omega[i][2]) * inertia[type[i]]; - } + if (rmass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + rmass[i]; + t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * + INERTIA*shape[itype][0]*shape[itype][0]*rmass[i]; + } + + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + mass[itype]; + t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * + INERTIA*shape[itype][0]*shape[itype][0]*mass[itype]; + } + } } if (tempbias) tbias->restore_bias_all(); @@ -181,7 +319,7 @@ void ComputeTempSphere::compute_vector() { - int i; + int i,itype; invoked_vector = update->ntimestep; @@ -195,51 +333,103 @@ double *mass = atom->mass; double *rmass = atom->rmass; double *radius = atom->radius; + double **shape = atom->shape; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + // 4 cases depending on radius vs shape and rmass vs mass + // point particles will not contribute rotation due to radius or shape = 0 + double massone,inertiaone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; - if (rmass) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - massone = rmass[i]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; - - inertiaone = INERTIA*radius[i]*radius[i]*rmass[i]; - t[0] += inertiaone * omega[i][0]*omega[i][0]; - t[1] += inertiaone * omega[i][1]*omega[i][1]; - t[2] += inertiaone * omega[i][2]*omega[i][2]; - t[3] += inertiaone * omega[i][0]*omega[i][1]; - t[4] += inertiaone * omega[i][0]*omega[i][2]; - t[5] += inertiaone * omega[i][1]*omega[i][2]; - } + if (radius) { + if (rmass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + massone = rmass[i]; + t[0] += massone * v[i][0]*v[i][0]; + t[1] += massone * v[i][1]*v[i][1]; + t[2] += massone * v[i][2]*v[i][2]; + t[3] += massone * v[i][0]*v[i][1]; + t[4] += massone * v[i][0]*v[i][2]; + t[5] += massone * v[i][1]*v[i][2]; + + inertiaone = INERTIA*radius[i]*radius[i]*rmass[i]; + t[0] += inertiaone * omega[i][0]*omega[i][0]; + t[1] += inertiaone * omega[i][1]*omega[i][1]; + t[2] += inertiaone * omega[i][2]*omega[i][2]; + t[3] += inertiaone * omega[i][0]*omega[i][1]; + t[4] += inertiaone * omega[i][0]*omega[i][2]; + t[5] += inertiaone * omega[i][1]*omega[i][2]; + } + + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + massone = mass[itype]; + t[0] += massone * v[i][0]*v[i][0]; + t[1] += massone * v[i][1]*v[i][1]; + t[2] += massone * v[i][2]*v[i][2]; + t[3] += massone * v[i][0]*v[i][1]; + t[4] += massone * v[i][0]*v[i][2]; + t[5] += massone * v[i][1]*v[i][2]; + + inertiaone = INERTIA*radius[i]*radius[i]*mass[itype]; + t[0] += inertiaone * omega[i][0]*omega[i][0]; + t[1] += inertiaone * omega[i][1]*omega[i][1]; + t[2] += inertiaone * omega[i][2]*omega[i][2]; + t[3] += inertiaone * omega[i][0]*omega[i][1]; + t[4] += inertiaone * omega[i][0]*omega[i][2]; + t[5] += inertiaone * omega[i][1]*omega[i][2]; + } + } + } else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - massone = mass[type[i]]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; - - inertiaone = inertia[type[i]]; - t[0] += inertiaone * omega[i][0]*omega[i][0]; - t[1] += inertiaone * omega[i][1]*omega[i][1]; - t[2] += inertiaone * omega[i][2]*omega[i][2]; - t[3] += inertiaone * omega[i][0]*omega[i][1]; - t[4] += inertiaone * omega[i][0]*omega[i][2]; - t[5] += inertiaone * omega[i][1]*omega[i][2]; - } + if (rmass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + massone = rmass[i]; + t[0] += massone * v[i][0]*v[i][0]; + t[1] += massone * v[i][1]*v[i][1]; + t[2] += massone * v[i][2]*v[i][2]; + t[3] += massone * v[i][0]*v[i][1]; + t[4] += massone * v[i][0]*v[i][2]; + t[5] += massone * v[i][1]*v[i][2]; + + inertiaone = INERTIA*shape[itype][0]*shape[itype][0]*rmass[i]; + t[0] += inertiaone * omega[i][0]*omega[i][0]; + t[1] += inertiaone * omega[i][1]*omega[i][1]; + t[2] += inertiaone * omega[i][2]*omega[i][2]; + t[3] += inertiaone * omega[i][0]*omega[i][1]; + t[4] += inertiaone * omega[i][0]*omega[i][2]; + t[5] += inertiaone * omega[i][1]*omega[i][2]; + } + + } else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + massone = mass[itype]; + t[0] += massone * v[i][0]*v[i][0]; + t[1] += massone * v[i][1]*v[i][1]; + t[2] += massone * v[i][2]*v[i][2]; + t[3] += massone * v[i][0]*v[i][1]; + t[4] += massone * v[i][0]*v[i][2]; + t[5] += massone * v[i][1]*v[i][2]; + + inertiaone = INERTIA*shape[itype][0]*shape[itype][0]*mass[itype]; + t[0] += inertiaone * omega[i][0]*omega[i][0]; + t[1] += inertiaone * omega[i][1]*omega[i][1]; + t[2] += inertiaone * omega[i][2]*omega[i][2]; + t[3] += inertiaone * omega[i][0]*omega[i][1]; + t[4] += inertiaone * omega[i][0]*omega[i][2]; + t[5] += inertiaone * omega[i][1]*omega[i][2]; + } + } } if (tempbias) tbias->restore_bias_all(); diff -Naur lammps-27Jun09/src/fix_nph.cpp lammps-28Jun09/src/fix_nph.cpp --- lammps-27Jun09/src/fix_nph.cpp 2008-05-16 11:20:43.000000000 -0600 +++ lammps-28Jun09/src/fix_nph.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -250,8 +250,6 @@ void FixNPH::init() { if (domain->triclinic) error->all("Cannot use fix nph with triclinic box"); - if (atom->mass == NULL) - error->all("Cannot use fix nph without per-type mass defined"); for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { @@ -353,6 +351,7 @@ void FixNPH::initial_integrate(int vflag) { int i; + double dtfm; double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; @@ -380,18 +379,29 @@ double **x = atom->x; double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - double dtfm; - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + if (rmass) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } + } + } else { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } } } @@ -423,23 +433,35 @@ void FixNPH::final_integrate() { int i; + double dtfm; // v update only for atoms in NPH group double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - double dtfm; - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; - v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; - v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + if (rmass) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + } + } + } else { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + } } } @@ -499,6 +521,7 @@ double **x = atom->x; double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; @@ -533,12 +556,23 @@ // v update only for atoms in NPH group - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } } } @@ -550,12 +584,23 @@ // v update only for atoms in NPH group - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } } } } @@ -595,17 +640,29 @@ double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } } } } diff -Naur lammps-27Jun09/src/fix_npt.cpp lammps-28Jun09/src/fix_npt.cpp --- lammps-27Jun09/src/fix_npt.cpp 2009-02-27 14:23:34.000000000 -0700 +++ lammps-28Jun09/src/fix_npt.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -262,8 +262,6 @@ void FixNPT::init() { if (domain->triclinic) error->all("Cannot use fix npt with triclinic box"); - if (atom->mass == NULL) - error->all("Cannot use fix npt without per-type mass defined"); for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { @@ -359,6 +357,7 @@ void FixNPT::initial_integrate(int vflag) { int i; + double dtfm; double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; @@ -389,37 +388,59 @@ dilation[i] = exp(dthalf*omega_dot[i]); } - // v update only for atoms in group - double **x = atom->x; double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - double dtfm; - - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + if (rmass) { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } - } else if (which == BIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); + + } else { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } } @@ -452,37 +473,60 @@ void FixNPT::final_integrate() { int i; - - // v update only for atoms in group + double dtfm; double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - double dtfm; - - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; - v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; - v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + if (rmass) { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + temperature->restore_bias(i,v[i]); + } } } - } else if (which == BIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]]; - v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; - v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; - v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; - temperature->restore_bias(i,v[i]); + + } else { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[type[i]]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + temperature->restore_bias(i,v[i]); + } } } } @@ -549,6 +593,7 @@ double **x = atom->x; double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; @@ -592,12 +637,23 @@ // v update only for atoms in group - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } } } @@ -609,24 +665,49 @@ // v update only for atoms in group - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; + if (rmass) { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); + + } else { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } } @@ -668,30 +749,56 @@ double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; + if (rmass) { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); + + } else { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } } diff -Naur lammps-27Jun09/src/fix_npt_sphere.cpp lammps-28Jun09/src/fix_npt_sphere.cpp --- lammps-27Jun09/src/fix_npt_sphere.cpp 2009-01-09 11:18:56.000000000 -0700 +++ lammps-28Jun09/src/fix_npt_sphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -35,41 +35,61 @@ FixNPTSphere::FixNPTSphere(LAMMPS *lmp, int narg, char **arg) : FixNPT(lmp, narg, arg) { + // error checks + if (!atom->omega_flag || !atom->torque_flag) error->all("Fix npt/sphere requires atom attributes omega, torque"); - - dttype = new double[atom->ntypes+1]; + if (!atom->radius_flag && !atom->avec->shape_type) + error->all("Fix npt/sphere requires atom attribute radius or shape"); } /* ---------------------------------------------------------------------- */ -FixNPTSphere::~FixNPTSphere() +void FixNPTSphere::init() { - delete [] dttype; -} + int i,itype; -/* ---------------------------------------------------------------------- */ + // check that all particles are finite-size and spherical + // no point particles allowed -void FixNPTSphere::init() -{ - FixNPT::init(); + if (atom->radius_flag) { + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (!atom->shape) - error->all("Fix npt/sphere requires atom attribute shape"); + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (radius[i] == 0.0) + error->one("Fix nvt/sphere requires extended particles"); + } - double **shape = atom->shape; - for (int i = 1; i <= atom->ntypes; i++) - if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) - error->all("Fix npt/sphere requires spherical particle shapes"); + } else { + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + if (shape[itype][0] == 0.0) + error->one("Fix nvt/sphere requires extended particles"); + if (shape[itype][0] != shape[itype][1] || + shape[itype][0] != shape[itype][2]) + error->one("Fix nvt/sphere requires spherical particle shapes"); + } + } + + FixNPT::init(); } -/* ---------------------------------------------------------------------- - 1st half of Verlet update -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ void FixNPTSphere::initial_integrate(int vflag) { - int i; + int i,itype; double dtfm,dtirotate; double delta = update->ntimestep - update->beginstep; @@ -102,37 +122,63 @@ } factor_rotate = exp(-dthalf*eta_dot); - // v update only for atoms in group - double **x = atom->x; double **v = atom->v; double **f = atom->f; double **omega = atom->omega; double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; double *mass = atom->mass; + double **shape = atom->shape; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + if (rmass) { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } - } else if (which == BIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); + + } else { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } } @@ -151,24 +197,57 @@ } } - // recompute timesteps since dt may have changed or come via rRESPA + // set timestep here since dt may have changed or come via rRESPA double dtfrotate = dtf / INERTIA; - int ntypes = atom->ntypes; - double **shape = atom->shape; - for (int i = 1; i <= ntypes; i++) - dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]); - // update angular momentum by 1/2 step - // update quaternion a full step via Richardson iteration - // returns new normalized quaternion + // update omega for all particles + // d_omega/dt = torque / inertia + // 4 cases depending on radius vs shape and rmass vs mass + + if (radius) { + if (rmass) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2]; + } + } + } else { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[type[i]]); + omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2]; + } + } + } - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtirotate = dttype[type[i]]; - omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2]; + } else { + if (rmass) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2]; + } + } + } else { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2]; + } + } } } @@ -179,68 +258,199 @@ if (kspace_flag) force->kspace->setup(); } -/* ---------------------------------------------------------------------- - 2nd half of Verlet update -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ void FixNPTSphere::final_integrate() { int i,itype; double dtfm,dtirotate; - // v update only for atoms in group - double **v = atom->v; double **f = atom->f; double **omega = atom->omega; double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; double *mass = atom->mass; + double **shape = atom->shape; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - // recompute timesteps since dt may have changed or come via rRESPA + // set timestep here since dt may have changed or come via rRESPA double dtfrotate = dtf / INERTIA; - int ntypes = atom->ntypes; - double **shape = atom->shape; - for (int i = 1; i <= ntypes; i++) - dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]); - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - - dtfm = dtf / mass[itype]; - v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; - v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; - v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; - - dtirotate = dttype[itype]; - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor_rotate; + // update v,omega for all particles + // d_omega/dt = torque / inertia + // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias + + if (radius) { + if (rmass) { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * + factor_rotate; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * + factor_rotate; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * + factor_rotate; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + temperature->restore_bias(i,v[i]); + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * + factor_rotate; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * + factor_rotate; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * + factor_rotate; + } + } + } + + } else { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * + factor_rotate; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * + factor_rotate; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * + factor_rotate; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[itype]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + temperature->restore_bias(i,v[i]); + + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * + factor_rotate; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * + factor_rotate; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * + factor_rotate; + } + } } } - } else if (which == BIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[itype]; - v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; - v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; - v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; - temperature->restore_bias(i,v[i]); - - dtirotate = dttype[itype]; - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor_rotate; + } else { + if (rmass) { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / rmass[i]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * + factor_rotate; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * + factor_rotate; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * + factor_rotate; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + temperature->restore_bias(i,v[i]); + + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * + factor_rotate; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * + factor_rotate; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * + factor_rotate; + } + } + } + + } else { + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * + factor_rotate; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * + factor_rotate; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * + factor_rotate; + } + } + } else if (which == BIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[itype]; + v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0]; + v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1]; + v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; + temperature->restore_bias(i,v[i]); + + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * + factor_rotate; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * + factor_rotate; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * + factor_rotate; + } + } } } } diff -Naur lammps-27Jun09/src/fix_npt_sphere.h lammps-28Jun09/src/fix_npt_sphere.h --- lammps-27Jun09/src/fix_npt_sphere.h 2008-04-08 16:55:40.000000000 -0600 +++ lammps-28Jun09/src/fix_npt_sphere.h 2009-06-26 12:23:16.000000000 -0600 @@ -21,13 +21,12 @@ class FixNPTSphere : public FixNPT { public: FixNPTSphere(class LAMMPS *, int, char **); - ~FixNPTSphere(); + ~FixNPTSphere() {} void init(); void initial_integrate(int); void final_integrate(); private: - double *dttype; double factor_rotate; }; diff -Naur lammps-27Jun09/src/fix_nve_sphere.cpp lammps-28Jun09/src/fix_nve_sphere.cpp --- lammps-27Jun09/src/fix_nve_sphere.cpp 2009-02-12 14:36:52.000000000 -0700 +++ lammps-28Jun09/src/fix_nve_sphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -51,21 +51,14 @@ } else error->all("Illegal fix nve/sphere command"); } - // error check + // error checks if (!atom->omega_flag || !atom->torque_flag) error->all("Fix nve/sphere requires atom attributes omega, torque"); + if (!atom->radius_flag && !atom->avec->shape_type) + error->all("Fix nve/sphere requires atom attribute radius or shape"); if (extra == DIPOLE && !atom->mu_flag) error->all("Fix nve/sphere requires atom attribute mu"); - - dttype = new double[atom->ntypes+1]; -} - -/* ---------------------------------------------------------------------- */ - -FixNVESphere::~FixNVESphere() -{ - delete [] dttype; } /* ---------------------------------------------------------------------- */ @@ -84,23 +77,42 @@ void FixNVESphere::init() { - dtv = update->dt; - dtf = 0.5 * update->dt * force->ftm2v; + int i,itype; - if (strcmp(update->integrate_style,"respa") == 0) - step_respa = ((Respa *) update->integrate)->step; + // check that all particles are finite-size and spherical + // no point particles allowed - if (atom->mass && !atom->shape) - error->all("Fix nve/sphere requires atom attribute shape"); - if (atom->rmass && !atom->radius_flag) - error->all("Fix nve/sphere requires atom attribute radius"); + if (atom->radius_flag) { + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (atom->mass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (radius[i] == 0.0) + error->one("Fix nve/sphere requires extended particles"); + } + + } else { double **shape = atom->shape; - for (int i = 1; i <= atom->ntypes; i++) - if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) - error->all("Fix nve/sphere requires spherical particle shapes"); + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + if (shape[itype][0] == 0.0) + error->one("Fix nve/sphere requires extended particles"); + if (shape[itype][0] != shape[itype][1] || + shape[itype][0] != shape[itype][2]) + error->one("Fix nve/sphere requires spherical particle shapes"); + } } + + FixNVE::init(); } /* ---------------------------------------------------------------------- */ @@ -119,58 +131,97 @@ double *radius = atom->radius; double *rmass = atom->rmass; double *mass = atom->mass; - int *mask = atom->mask; + double **shape = atom->shape; int *type = atom->type; + int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - // recompute timesteps since dt may have changed or come via rRESPA + // set timestep here since dt may have changed or come via rRESPA - dtfrotate = dtf / INERTIA; - if (mass) { - double **shape = atom->shape; - int ntypes = atom->ntypes; - for (int i = 1; i <= ntypes; i++) - dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]); - } + double dtfrotate = dtf / INERTIA; // update v,x,omega for all particles // d_omega/dt = torque / inertia + // 4 cases depending on radius vs shape and rmass vs mass - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; + if (radius) { + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } + } + + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } } } } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dttype[itype]; - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } + } + + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } } } } @@ -213,50 +264,85 @@ double *mass = atom->mass; double *rmass = atom->rmass; double *radius = atom->radius; + double **shape = atom->shape; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - // recompute timesteps since dt may have changed or come via rRESPA + // set timestep here since dt may have changed or come via rRESPA - dtfrotate = dtf / INERTIA; - if (mass) { - double **shape = atom->shape; - int ntypes = atom->ntypes; - for (int i = 1; i <= ntypes; i++) - dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]); - } + double dtfrotate = dtf / INERTIA; - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; + // update v,omega for all particles + // d_omega/dt = torque / inertia + // 4 cases depending on radius vs shape and rmass vs mass + + if (radius) { + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } + } + + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } } } - + } else { - dtfrotate = dtf / INERTIA; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - - dtirotate = dttype[itype]; - omega[i][0] += dtirotate * torque[i][0]; - omega[i][1] += dtirotate * torque[i][1]; - omega[i][2] += dtirotate * torque[i][2]; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } + } + + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } } } } diff -Naur lammps-27Jun09/src/fix_nve_sphere.h lammps-28Jun09/src/fix_nve_sphere.h --- lammps-27Jun09/src/fix_nve_sphere.h 2008-04-08 16:55:40.000000000 -0600 +++ lammps-28Jun09/src/fix_nve_sphere.h 2009-06-26 12:23:16.000000000 -0600 @@ -21,7 +21,7 @@ class FixNVESphere : public FixNVE { public: FixNVESphere(class LAMMPS *, int, char **); - ~FixNVESphere(); + ~FixNVESphere() {} int setmask(); void init(); void initial_integrate(int); @@ -29,8 +29,6 @@ private: int extra; - double dtfrotate; - double *dttype; }; } diff -Naur lammps-27Jun09/src/fix_nvt.cpp lammps-28Jun09/src/fix_nvt.cpp --- lammps-27Jun09/src/fix_nvt.cpp 2009-01-05 15:26:08.000000000 -0700 +++ lammps-28Jun09/src/fix_nvt.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -117,9 +117,6 @@ void FixNVT::init() { - if (atom->mass == NULL) - error->all("Cannot use fix nvt without per-type mass defined"); - int icompute = modify->find_compute(id_temp); if (icompute < 0) error->all("Temp ID for fix nvt does not exist"); temperature = modify->compute[icompute]; @@ -132,7 +129,6 @@ dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; dthalf = 0.5 * update->dt; - drag_factor = 1.0 - (update->dt * t_freq * drag); if (strcmp(update->integrate_style,"respa") == 0) { @@ -172,37 +168,70 @@ double **x = atom->x; double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; + if (rmass) { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } + + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } } } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; + } else { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } + + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } } } } @@ -218,31 +247,58 @@ double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]] * factor; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + if (rmass) { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + } + } + + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]] * factor; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); + } else { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + } + } + + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[type[i]] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } } @@ -276,6 +332,7 @@ double **x = atom->x; double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; @@ -300,25 +357,51 @@ factor = exp(-dthalf*eta_dot); } else factor = 1.0; - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + if (rmass) { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + } + } + + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); + } else { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + } + } + + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + } } } } @@ -351,18 +434,31 @@ else { double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } } } } @@ -440,7 +536,6 @@ dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; dthalf = 0.5 * update->dt; - drag_factor = 1.0 - (update->dt * t_freq * drag); } diff -Naur lammps-27Jun09/src/fix_nvt_sphere.cpp lammps-28Jun09/src/fix_nvt_sphere.cpp --- lammps-27Jun09/src/fix_nvt_sphere.cpp 2008-10-22 09:52:45.000000000 -0600 +++ lammps-28Jun09/src/fix_nvt_sphere.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -36,32 +36,54 @@ FixNVTSphere::FixNVTSphere(LAMMPS *lmp, int narg, char **arg) : FixNVT(lmp, narg, arg) { + // error checks + if (!atom->omega_flag || !atom->torque_flag) error->all("Fix nvt/sphere requires atom attributes omega, torque"); - - dttype = new double[atom->ntypes+1]; + if (!atom->radius_flag && !atom->avec->shape_type) + error->all("Fix nvt/sphere requires atom attribute radius or shape"); } /* ---------------------------------------------------------------------- */ -FixNVTSphere::~FixNVTSphere() +void FixNVTSphere::init() { - delete [] dttype; -} + int i,itype; -/* ---------------------------------------------------------------------- */ + // check that all particles are finite-size and spherical + // no point particles allowed -void FixNVTSphere::init() -{ - FixNVT::init(); + if (atom->radius_flag) { + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (radius[i] == 0.0) + error->one("Fix nvt/sphere requires extended particles"); + } - if (!atom->shape) - error->all("Fix nvt/sphere requires atom attribute shape"); + } else { + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; - double **shape = atom->shape; - for (int i = 1; i <= atom->ntypes; i++) - if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) - error->all("Fix nvt/sphere requires spherical particle shapes"); + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + if (shape[itype][0] == 0.0) + error->one("Fix nvt/sphere requires extended particles"); + if (shape[itype][0] != shape[itype][1] || + shape[itype][0] != shape[itype][2]) + error->one("Fix nvt/sphere requires spherical particle shapes"); + } + } + + FixNVT::init(); } /* ---------------------------------------------------------------------- */ @@ -90,59 +112,188 @@ double **f = atom->f; double **omega = atom->omega; double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; double *mass = atom->mass; + double **shape = atom->shape; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - // recompute timesteps since dt may have changed or come via rRESPA + // set timestep here since dt may have changed or come via rRESPA double dtfrotate = dtf / INERTIA; - int ntypes = atom->ntypes; - double **shape = atom->shape; - for (int i = 1; i <= ntypes; i++) - dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]); - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; + // update v,x,omega for all particles + // d_omega/dt = torque / inertia + // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias + + if (radius) { + if (rmass) { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } + } + } - dtfm = dtf / mass[itype]; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dttype[itype]; - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } else { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); + omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[itype]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); + omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } + } } } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; + } else { + if (rmass) { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } + } + } - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[itype]; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dttype[itype]; - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } else { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[itype]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; + omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; + omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; + } + } } } } @@ -161,53 +312,164 @@ double **f = atom->f; double **omega = atom->omega; double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; double *mass = atom->mass; + double **shape = atom->shape; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - // recompute timesteps since dt may have changed or come via rRESPA + // set timestep here since dt may have changed or come via rRESPA double dtfrotate = dtf / INERTIA; - int ntypes = atom->ntypes; - double **shape = atom->shape; - for (int i = 1; i <= ntypes; i++) - dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]); - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; + // update v,omega for all particles + // d_omega/dt = torque / inertia + // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias + + if (radius) { + if (rmass) { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } + } + } - dtfm = dtf / mass[itype] * factor; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; - - dtirotate = dttype[itype]; - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } else { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[itype] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } + } } } } else { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; + if (rmass) { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / rmass[i] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + temperature->remove_bias(i,v[i]); + dtfm = dtf / rmass[i] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } + } + } - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[itype] * factor; - v[i][0] = v[i][0]*factor + dtfm*f[i][0]; - v[i][1] = v[i][1]*factor + dtfm*f[i][1]; - v[i][2] = v[i][2]*factor + dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); - - dtirotate = dttype[itype]; - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } else { + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtfm = dtf / mass[itype] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + temperature->remove_bias(i,v[i]); + dtfm = dtf / mass[itype] * factor; + v[i][0] = v[i][0]*factor + dtfm*f[i][0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2]; + temperature->restore_bias(i,v[i]); + + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; + omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; + omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; + } + } } } } diff -Naur lammps-27Jun09/src/fix_nvt_sphere.h lammps-28Jun09/src/fix_nvt_sphere.h --- lammps-27Jun09/src/fix_nvt_sphere.h 2008-04-08 16:55:40.000000000 -0600 +++ lammps-28Jun09/src/fix_nvt_sphere.h 2009-06-26 12:23:16.000000000 -0600 @@ -21,13 +21,10 @@ class FixNVTSphere : public FixNVT { public: FixNVTSphere(class LAMMPS *, int, char **); - ~FixNVTSphere(); + ~FixNVTSphere() {} void init(); void initial_integrate(int); void final_integrate(); - - private: - double *dttype; }; } diff -Naur lammps-27Jun09/src/fix_press_berendsen.cpp lammps-28Jun09/src/fix_press_berendsen.cpp --- lammps-27Jun09/src/fix_press_berendsen.cpp 2008-07-08 13:45:01.000000000 -0600 +++ lammps-28Jun09/src/fix_press_berendsen.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -234,15 +234,14 @@ { if (domain->triclinic) error->all("Cannot use fix press/berendsen with triclinic box"); - if (atom->mass == NULL) - error->all("Cannot use fix press/berendsen without per-type mass defined"); for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { int *dimflag = ((FixDeform *) modify->fix[i])->dimflag; if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || (p_flag[2] && dimflag[2])) - error->all("Cannot use fix press/berendsen and fix deform on same dimension"); + error->all("Cannot use fix press/berendsen and " + "fix deform on same dimension"); } // set temperature and pressure ptrs diff -Naur lammps-27Jun09/src/fix_rigid.cpp lammps-28Jun09/src/fix_rigid.cpp --- lammps-27Jun09/src/fix_rigid.cpp 2009-04-10 17:25:32.000000000 -0600 +++ lammps-28Jun09/src/fix_rigid.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -368,9 +368,6 @@ triclinic = domain->triclinic; - if (atom->mass == NULL) - error->all("Cannot use fix rigid without per-type mass defined"); - // warn if more than one rigid fix int count = 0; @@ -401,10 +398,11 @@ // compute masstotal & center-of-mass of each rigid body - int *type = atom->type; - int *image = atom->image; - double *mass = atom->mass; double **x = atom->x; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *image = atom->image; + int *type = atom->type; int nlocal = atom->nlocal; double xprd = domain->xprd; @@ -426,7 +424,8 @@ xbox = (image[i] & 1023) - 512; ybox = (image[i] >> 10 & 1023) - 512; zbox = (image[i] >> 20) - 512; - massone = mass[type[i]]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; if (triclinic == 0) { xunwrap = x[i][0] + xbox*xprd; @@ -486,7 +485,8 @@ dx = xunwrap - xcm[ibody][0]; dy = yunwrap - xcm[ibody][1]; dz = zunwrap - xcm[ibody][2]; - massone = mass[type[i]]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; sum[ibody][0] += massone * (dy*dy + dz*dz); sum[ibody][1] += massone * (dx*dx + dz*dz); @@ -613,7 +613,8 @@ for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; - massone = mass[type[i]]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; sum[ibody][0] += massone * (displace[i][1]*displace[i][1] + displace[i][2]*displace[i][2]); @@ -669,6 +670,7 @@ int *type = atom->type; double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int nlocal = atom->nlocal; @@ -679,7 +681,9 @@ for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; - massone = mass[type[i]]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + sum[ibody][0] += v[i][0] * massone; sum[ibody][1] += v[i][1] * massone; sum[ibody][2] += v[i][2] * massone; @@ -739,7 +743,9 @@ dy = yunwrap - xcm[ibody][1]; dz = zunwrap - xcm[ibody][2]; - massone = mass[type[i]]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + sum[ibody][0] += dy * massone*v[i][2] - dz * massone*v[i][1]; sum[ibody][1] += dz * massone*v[i][0] - dx * massone*v[i][2]; sum[ibody][2] += dx * massone*v[i][1] - dy * massone*v[i][0]; @@ -1413,6 +1419,7 @@ double **x = atom->x; double **v = atom->v; double **f = atom->f; + double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; int nlocal = atom->nlocal; @@ -1495,7 +1502,9 @@ // assume per-atom contribution is due to constraint force on that atom if (evflag) { - massone = mass[type[i]]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; @@ -1526,10 +1535,11 @@ double xy,xz,yz; double vr[6]; - double *mass = atom->mass; double **f = atom->f; double **v = atom->v; double **x = atom->x; + double *rmass = atom->rmass; + double *mass = atom->mass; int *type = atom->type; int *image = atom->image; int nlocal = atom->nlocal; @@ -1578,7 +1588,9 @@ // assume per-atom contribution is due to constraint force on that atom if (evflag) { - massone = mass[type[i]]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; diff -Naur lammps-27Jun09/src/fix_temp_berendsen.cpp lammps-28Jun09/src/fix_temp_berendsen.cpp --- lammps-27Jun09/src/fix_temp_berendsen.cpp 2008-03-20 09:32:33.000000000 -0600 +++ lammps-28Jun09/src/fix_temp_berendsen.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -87,9 +87,6 @@ void FixTempBerendsen::init() { - if (atom->mass == NULL) - error->all("Cannot use fix temp/berendsen without per-type mass defined"); - int icompute = modify->find_compute(id_temp); if (icompute < 0) error->all("Temp ID for fix temp/berendsen does not exist"); @@ -127,7 +124,6 @@ v[i][2] *= lamda; } } - } else if (which == BIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { diff -Naur lammps-27Jun09/src/variable.cpp lammps-28Jun09/src/variable.cpp --- lammps-27Jun09/src/variable.cpp 2009-05-06 07:07:08.000000000 -0600 +++ lammps-28Jun09/src/variable.cpp 2009-06-26 12:23:16.000000000 -0600 @@ -1424,13 +1424,18 @@ double *argstack, int &nargstack) { // parse contents for arg1,arg2,arg3 separated by commas - - char *ptr1 = strchr(contents,','); - if (ptr1) *ptr1 = '\0'; - char *ptr2 = strchr(ptr1+1,','); - if (ptr2) *ptr2 = '\0'; + // ptr1,ptr2 = location of 1st and 2nd comma, NULL if none char *arg1,*arg2,*arg3; + char *ptr1,*ptr2; + + ptr1 = strchr(contents,','); + if (ptr1) { + *ptr1 = '\0'; + ptr2 = strchr(ptr1+1,','); + if (ptr2) *ptr2 = '\0'; + } else ptr2 = NULL; + int n = strlen(contents) + 1; arg1 = new char[n]; strcpy(arg1,contents);