The dynamic keyword determines whether the number of atoms N in the compute group is re-computed each time a temperature is computed. diff -Naur lammps-4Jul09/doc/compute_modify.txt lammps-5Jul09/doc/compute_modify.txt --- lammps-4Jul09/doc/compute_modify.txt 2009-04-01 16:25:09.000000000 -0600 +++ lammps-5Jul09/doc/compute_modify.txt 2009-07-02 08:39:37.000000000 -0600 @@ -38,7 +38,10 @@ temperature computation. Only computes that compute a temperature use this option. The default is 2 or 3 for "2d or 3d systems"_dimension.html which is a correction factor for an ensemble -of velocities with zero total linear momentum. +of velocities with zero total linear momentum. You can use a negative +number for the {extra} parameter if you need to add +degrees-of-freedom. See the "compute +temp/asphere"_compute_temp_aspher.html command for an example. The {dynamic} keyword determines whether the number of atoms N in the compute group is re-computed each time a temperature is computed. diff -Naur lammps-4Jul09/doc/compute_temp_asphere.html lammps-5Jul09/doc/compute_temp_asphere.html --- lammps-4Jul09/doc/compute_temp_asphere.html 2009-06-26 12:41:31.000000000 -0600 +++ lammps-5Jul09/doc/compute_temp_asphere.html 2009-07-02 08:39:37.000000000 -0600 @@ -32,16 +32,18 @@ usual compute temp command, which assumes point particles with only translational kinetic energy.
-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 +
Only finite-size particles (aspherical or spherical) can be included +in the group. For 3d finite-size particles, each has 6 degrees of +freedom (3 translational, 3 rotational). For 2d finite-size +particles, each has 3 degrees of freedom (2 translational, 1 +rotational). +
+IMPORTANT NOTE: This choice for degrees of freedom (dof) assumes that +all finite-size aspherical or 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.
diff -Naur lammps-4Jul09/doc/compute_temp_asphere.txt lammps-5Jul09/doc/compute_temp_asphere.txt --- lammps-4Jul09/doc/compute_temp_asphere.txt 2009-06-26 12:41:31.000000000 -0600 +++ lammps-5Jul09/doc/compute_temp_asphere.txt 2009-07-02 08:39:37.000000000 -0600 @@ -29,16 +29,18 @@ usual "compute temp"_compute_temp.html command, which assumes point particles with only translational kinetic energy. -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 +Only finite-size particles (aspherical or spherical) can be included +in the group. For 3d finite-size particles, each has 6 degrees of +freedom (3 translational, 3 rotational). For 2d finite-size +particles, each has 3 degrees of freedom (2 translational, 1 +rotational). + +IMPORTANT NOTE: This choice for degrees of freedom (dof) assumes that +all finite-size aspherical or 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. diff -Naur lammps-4Jul09/doc/compute_temp_sphere.html lammps-5Jul09/doc/compute_temp_sphere.html --- lammps-4Jul09/doc/compute_temp_sphere.html 2009-06-26 12:41:31.000000000 -0600 +++ lammps-5Jul09/doc/compute_temp_sphere.html 2009-07-02 08:39:37.000000000 -0600 @@ -32,13 +32,15 @@ usual compute temp command, which assumes point particles with only translational kinetic energy. -For 3d spherical particles, each has 6 degrees of freedom (3 -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 +
Both point and finite-size particles can be included in the group. +Point particles do not rotate, so they have only translational degrees +of freedom. For 3d spherical particles, each has 6 degrees of freedom +(3 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) assumes that +all finite-size 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 @@ -56,7 +58,7 @@
A 6-component kinetic energy tensor is also calculated by this compute. The formula for the components of the tensor is the same as -the above formula, except that v^2 and w^2 are replaced by vx*vy and +the above formulas, except that v^2 and w^2 are replaced by vx*vy and wx*wy for the xy component.
The number of atoms contributing to the temperature is assumed to be diff -Naur lammps-4Jul09/doc/compute_temp_sphere.txt lammps-5Jul09/doc/compute_temp_sphere.txt --- lammps-4Jul09/doc/compute_temp_sphere.txt 2009-06-26 12:41:31.000000000 -0600 +++ lammps-5Jul09/doc/compute_temp_sphere.txt 2009-07-02 08:39:37.000000000 -0600 @@ -29,13 +29,15 @@ usual "compute temp"_compute_temp.html command, which assumes point particles with only translational kinetic energy. -For 3d spherical particles, each has 6 degrees of freedom (3 -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 +Both point and finite-size particles can be included in the group. +Point particles do not rotate, so they have only translational degrees +of freedom. For 3d spherical particles, each has 6 degrees of freedom +(3 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) assumes that +all finite-size 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 @@ -53,7 +55,7 @@ A 6-component kinetic energy tensor is also calculated by this compute. The formula for the components of the tensor is the same as -the above formula, except that v^2 and w^2 are replaced by vx*vy and +the above formulas, except that v^2 and w^2 are replaced by vx*vy and wx*wy for the xy component. The number of atoms contributing to the temperature is assumed to be diff -Naur lammps-4Jul09/doc/fix_rigid.html lammps-5Jul09/doc/fix_rigid.html --- lammps-4Jul09/doc/fix_rigid.html 2009-06-30 14:58:06.000000000 -0600 +++ lammps-5Jul09/doc/fix_rigid.html 2009-07-02 09:09:14.000000000 -0600 @@ -156,31 +156,30 @@ includes all the desired rigid bodies. LAMMPS will allow multiple rigid fixes to be defined, but it is more expensive.
-This fix uses constant-energy integration, so you may need to impose -additional constraints to control the temperature of an ensemble of -rigid bodies. You can use fix langevin for this -purpose to treat the system as effectively immersed in an implicit -solvent, e.g. a Brownian dynamics model. Or you can thermostat only -the non-rigid atoms that surround one or more rigid bodies -(i.e. explicit solvent) by appropriate choice of groups in the compute -and fix commands for temperature and thermostatting. +
This fix uses constant-energy NVE-style integration, so you may need +to impose additional constraints to control the temperature of an +ensemble of rigid bodies. You can use fix +langevin for this purpose to treat the system as +effectively immersed in an implicit solvent, e.g. a Brownian dynamics +model. Or you can thermostat only the non-rigid atoms that surround +one or more rigid bodies (i.e. explicit solvent) by appropriate choice +of groups in the compute and fix commands for temperature and +thermostatting.
-If you calculate a temperature for the rigid bodies, the +
If you calculate a temperature for particles in the rigid bodies, the degrees-of-freedom removed by each rigid body are accounted for in the temperature (and pressure) computation, but only if the temperature -group includes the entire rigid body. Rigid bodies in 3d have 6 -degrees of freedom (3 translational, 3 rotational), except for dimers -which only have 5. Rigid bodies in 2d have 3 degrees of freedom. -
-IMPORTANT NOTE: Linear rigid bodies are ones consisting of point -particles in a straight line. Linear rigid bodies in 3d with three or -more atoms also have 5 degrees of freedom instead of 6, but LAMMPS -will not detect this. Thus you should use the -compute_modify command to subtract an additional -degree of freedom per rigid body. You may also wish to explicitly -subtract additional degrees-of-freedom if you use the force and -torque keywords to eliminate certain motions of the rigid body, as -LAMMPS does not do this automatically. +group includes all the particles in a particular rigid body. +
+A 3d rigid body has 6 degrees of freedom (3 translational, 3 +rotational), except for a collection of point particles lying on a +straight line, which has only 5, e.g a dimer. A 2d rigid body has 3 +degrees of freedom (2 translational, 1 rotational). +
+IMPORTANT NOTE: You may wish to explicitly subtract additional +degrees-of-freedom if you use the force and torque keywords to +eliminate certain motions of one or more rigid bodies, as LAMMPS does +not do this automatically.
The rigid body contribution to the pressure of the system (virial) is also accounted for by this fix. diff -Naur lammps-4Jul09/doc/fix_rigid.txt lammps-5Jul09/doc/fix_rigid.txt --- lammps-4Jul09/doc/fix_rigid.txt 2009-06-30 14:58:06.000000000 -0600 +++ lammps-5Jul09/doc/fix_rigid.txt 2009-07-02 09:09:14.000000000 -0600 @@ -147,31 +147,30 @@ includes all the desired rigid bodies. LAMMPS will allow multiple rigid fixes to be defined, but it is more expensive. -This fix uses constant-energy integration, so you may need to impose -additional constraints to control the temperature of an ensemble of -rigid bodies. You can use "fix langevin"_fix_langevin.html for this -purpose to treat the system as effectively immersed in an implicit -solvent, e.g. a Brownian dynamics model. Or you can thermostat only -the non-rigid atoms that surround one or more rigid bodies -(i.e. explicit solvent) by appropriate choice of groups in the compute -and fix commands for temperature and thermostatting. +This fix uses constant-energy NVE-style integration, so you may need +to impose additional constraints to control the temperature of an +ensemble of rigid bodies. You can use "fix +langevin"_fix_langevin.html for this purpose to treat the system as +effectively immersed in an implicit solvent, e.g. a Brownian dynamics +model. Or you can thermostat only the non-rigid atoms that surround +one or more rigid bodies (i.e. explicit solvent) by appropriate choice +of groups in the compute and fix commands for temperature and +thermostatting. -If you calculate a temperature for the rigid bodies, the +If you calculate a temperature for particles in the rigid bodies, the degrees-of-freedom removed by each rigid body are accounted for in the temperature (and pressure) computation, but only if the temperature -group includes the entire rigid body. Rigid bodies in 3d have 6 -degrees of freedom (3 translational, 3 rotational), except for dimers -which only have 5. Rigid bodies in 2d have 3 degrees of freedom. - -IMPORTANT NOTE: Linear rigid bodies are ones consisting of point -particles in a straight line. Linear rigid bodies in 3d with three or -more atoms also have 5 degrees of freedom instead of 6, but LAMMPS -will not detect this. Thus you should use the -"compute_modify"_compute_modify.html command to subtract an additional -degree of freedom per rigid body. You may also wish to explicitly -subtract additional degrees-of-freedom if you use the {force} and -{torque} keywords to eliminate certain motions of the rigid body, as -LAMMPS does not do this automatically. +group includes all the particles in a particular rigid body. + +A 3d rigid body has 6 degrees of freedom (3 translational, 3 +rotational), except for a collection of point particles lying on a +straight line, which has only 5, e.g a dimer. A 2d rigid body has 3 +degrees of freedom (2 translational, 1 rotational). + +IMPORTANT NOTE: You may wish to explicitly subtract additional +degrees-of-freedom if you use the {force} and {torque} keywords to +eliminate certain motions of one or more rigid bodies, as LAMMPS does +not do this automatically. The rigid body contribution to the pressure of the system (virial) is also accounted for by this fix. diff -Naur lammps-4Jul09/doc/fix_viscous.html lammps-5Jul09/doc/fix_viscous.html --- lammps-4Jul09/doc/fix_viscous.html 2009-06-30 16:51:36.000000000 -0600 +++ lammps-5Jul09/doc/fix_viscous.html 2009-07-02 07:41:57.000000000 -0600 @@ -62,12 +62,12 @@ LAMMPS units options like "real" or "metal" that are not self-consistent.
-In a Brownian dynamics context, gamma = kT / D, where k = Boltzmann's -constant, T = temperature, and D = particle diffusion coefficient. D -can be written as kT / (3 pi eta d), where eta = dynamic viscosity of -the frictional fluid and d = diameter of particle. This means gamma = -3 pi eta d, and thus is proportional to the viscosity of the fluid and -the particle diameter. +
In a Brownian dynamics context, gamma = Kb T / D, where Kb = +Boltzmann's constant, T = temperature, and D = particle diffusion +coefficient. D can be written as Kb T / (3 pi eta d), where eta = +dynamic viscosity of the frictional fluid and d = diameter of +particle. This means gamma = 3 pi eta d, and thus is proportional to +the viscosity of the fluid and the particle diameter.
In the current implementation, rather than have the user specify a viscosity, gamma is specified directly in force/velocity units. If diff -Naur lammps-4Jul09/doc/fix_viscous.txt lammps-5Jul09/doc/fix_viscous.txt --- lammps-4Jul09/doc/fix_viscous.txt 2009-06-30 16:51:36.000000000 -0600 +++ lammps-5Jul09/doc/fix_viscous.txt 2009-07-02 07:41:57.000000000 -0600 @@ -52,12 +52,12 @@ LAMMPS "units"_units.html options like "real" or "metal" that are not self-consistent. -In a Brownian dynamics context, gamma = kT / D, where k = Boltzmann's -constant, T = temperature, and D = particle diffusion coefficient. D -can be written as kT / (3 pi eta d), where eta = dynamic viscosity of -the frictional fluid and d = diameter of particle. This means gamma = -3 pi eta d, and thus is proportional to the viscosity of the fluid and -the particle diameter. +In a Brownian dynamics context, gamma = Kb T / D, where Kb = +Boltzmann's constant, T = temperature, and D = particle diffusion +coefficient. D can be written as Kb T / (3 pi eta d), where eta = +dynamic viscosity of the frictional fluid and d = diameter of +particle. This means gamma = 3 pi eta d, and thus is proportional to +the viscosity of the fluid and the particle diameter. In the current implementation, rather than have the user specify a viscosity, gamma is specified directly in force/velocity units. If diff -Naur lammps-4Jul09/src/ASPHERE/compute_temp_asphere.cpp lammps-5Jul09/src/ASPHERE/compute_temp_asphere.cpp --- lammps-4Jul09/src/ASPHERE/compute_temp_asphere.cpp 2009-07-01 14:28:37.000000000 -0600 +++ lammps-5Jul09/src/ASPHERE/compute_temp_asphere.cpp 2009-07-02 08:32:20.000000000 -0600 @@ -127,11 +127,12 @@ { // 6 dof for 3d, 3 dof for 2d // assume full rotation of extended particles - // user can correct this via compute_modify if needed + // user should correct this via compute_modify if needed double natoms = group->count(igroup); - if (domain->dimension == 3) dof = 6*natoms; - else dof = 3*natoms; + int nper = 6; + if (domain->dimension == 2) nper = 3; + dof = nper*natoms; // additional adjustments to dof @@ -145,8 +146,7 @@ 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; + dof -= nper*count_all; } dof -= extra_dof + fix_dof; diff -Naur lammps-4Jul09/src/compute_temp_sphere.cpp lammps-5Jul09/src/compute_temp_sphere.cpp --- lammps-4Jul09/src/compute_temp_sphere.cpp 2009-07-01 14:28:37.000000000 -0600 +++ lammps-5Jul09/src/compute_temp_sphere.cpp 2009-07-02 08:39:59.000000000 -0600 @@ -127,7 +127,7 @@ // 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 + // user should correct this via compute_modify if needed int dimension = domain->dimension; diff -Naur lammps-4Jul09/src/fix_ave_time.cpp lammps-5Jul09/src/fix_ave_time.cpp --- lammps-4Jul09/src/fix_ave_time.cpp 2009-01-09 14:31:31.000000000 -0700 +++ lammps-5Jul09/src/fix_ave_time.cpp 2009-07-01 14:44:39.000000000 -0600 @@ -356,7 +356,7 @@ compute->compute_scalar(); compute->invoked_flag |= INVOKED_SCALAR; } - vector[i] += compute->compute_scalar(); + vector[i] += compute->scalar; } else { if (!(compute->invoked_flag & INVOKED_VECTOR)) { compute->compute_vector(); diff -Naur lammps-4Jul09/src/fix_rigid.cpp lammps-5Jul09/src/fix_rigid.cpp --- lammps-4Jul09/src/fix_rigid.cpp 2009-06-30 11:18:31.000000000 -0600 +++ lammps-5Jul09/src/fix_rigid.cpp 2009-07-02 09:17:08.000000000 -0600 @@ -418,6 +418,8 @@ // extended = 1 if any particle in a rigid body is finite size + extended = dorientflag = qorientflag = 0; + double *radius = atom->radius; double *rmass = atom->rmass; double *mass = atom->mass; @@ -1357,51 +1359,68 @@ { int groupbit = group->bitmask[igroup]; - // nall = # of atoms in each rigid body that are also in temperature group + // nall = # of point particles in each rigid body + // mall = # of finite-size particles in each rigid body + // particles must also be in temperature group int *mask = atom->mask; int nlocal = atom->nlocal; int *ncount = new int[nbody]; - for (int ibody = 0; ibody < nbody; ibody++) ncount[ibody] = 0; + int *mcount = new int[nbody]; + for (int ibody = 0; ibody < nbody; ibody++) + ncount[ibody] = mcount[ibody] = 0; for (int i = 0; i < nlocal; i++) - if (body[i] >= 0 && mask[i] & groupbit) ncount[body[i]]++; + if (body[i] >= 0 && mask[i] & groupbit) { + if (extended && eflags[i]) mcount[body[i]]++; + else ncount[body[i]]++; + } int *nall = new int[nbody]; + int *mall = new int[nbody]; MPI_Allreduce(ncount,nall,nbody,MPI_INT,MPI_SUM,world); + MPI_Allreduce(mcount,mall,nbody,MPI_INT,MPI_SUM,world); - // warn if nall != nrigid for any body included in temperature group + // warn if nall+mall != nrigid for any body included in temperature group int flag = 0; for (int ibody = 0; ibody < nbody; ibody++) { - if (nall[ibody] && nall[ibody] != nrigid[ibody]) flag = 1; + if (nall[ibody]+mall[ibody] > 0 && + nall[ibody]+mall[ibody] != nrigid[ibody]) flag = 1; } if (flag && comm->me == 0) error->warning("Computing temperature of portions of rigid bodies"); // remove appropriate DOFs for each rigid body wholly in temperature group - // each has 3N (3d) or 2N (2d) to start with - // 3d N>2 rigid body should have 6 dof, so remove 3N-6 dof - // 3d N=2 rigid dimer should have 5 dof, so remove 1 dof - // 2d N>=2 rigid body should have 3 dof, so remove 2N-3 dof + // N = # of point particles in body + // M = # of finite-size particles in body + // 3d body has 3N + 6M dof to start with + // 2d body has 2N + 3M dof to start with + // 3d point-particle body with all non-zero I should have 6 dof, remove 3N-6 + // 3d point-particle body (linear) with a 0 I should have 5 dof, remove 3N-5 + // 2d point-particle body should have 3 dof, remove 2N-3 + // 3d body with any finite-size M should have 6 dof, remove (3N+6M) - 6 + // 2d body with any finite-size M should have 3 dof, remove (2N+3M) - 3 int n = 0; if (domain->dimension == 3) { for (int ibody = 0; ibody < nbody; ibody++) - if (nall[ibody] == nrigid[ibody]) { - n += 3*nrigid[ibody] - 6; - if (nrigid[ibody] == 2) n++; + if (nall[ibody]+mall[ibody] == nrigid[ibody]) { + n += 3*nall[ibody] + 6*mall[ibody] - 6; + if (inertia[ibody][0] == 0.0 || inertia[ibody][1] == 0.0 || + inertia[ibody][2] == 0.0) n++; } } else if (domain->dimension == 2) { for (int ibody = 0; ibody < nbody; ibody++) - if (nall[ibody] == nrigid[ibody]) { - n += 2*nrigid[ibody] - 3; - } + if (nall[ibody]+mall[ibody] == nrigid[ibody]) + n += 2*nall[ibody] + 3*mall[ibody] - 3; } delete [] ncount; + delete [] mcount; delete [] nall; + delete [] mall; return n; } diff -Naur lammps-4Jul09/src/fix_shake.cpp lammps-5Jul09/src/fix_shake.cpp --- lammps-4Jul09/src/fix_shake.cpp 2009-02-12 14:36:52.000000000 -0700 +++ lammps-5Jul09/src/fix_shake.cpp 2009-07-02 09:17:08.000000000 -0600 @@ -532,8 +532,8 @@ int *tag = atom->tag; int nlocal = atom->nlocal; - // count dof in a cluster if and only if the central atom is in group - // and this atom i is the central atom + // count dof in a cluster if and only if + // the central atom is in group and atom i is the central atom int n = 0; for (int i = 0; i < nlocal; i++) { diff -Naur lammps-4Jul09/src/modify.cpp lammps-5Jul09/src/modify.cpp --- lammps-4Jul09/src/modify.cpp 2009-03-12 17:42:57.000000000 -0600 +++ lammps-5Jul09/src/modify.cpp 2009-07-02 10:19:33.000000000 -0600 @@ -174,6 +174,20 @@ list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force); list_init(MIN_ENERGY,n_min_energy,list_min_energy); + // init each fix + // needs to come before compute init + // this is b/c some computes call fix->dof() + // FixRigid::dof() depends on its own init having been called + + comm->maxforward_fix = comm->maxreverse_fix = 0; + for (i = 0; i < nfix; i++) fix[i]->init(); + + // set global flag if any fix has its restart_pbc flag set + + restart_pbc_any = 0; + for (i = 0; i < nfix; i++) + if (fix[i]->restart_pbc) restart_pbc_any = 1; + // create list of computes that store invocation times list_init_compute(); @@ -187,17 +201,6 @@ for (i = 0; i < ncompute; i++) compute[i]->init(); addstep_compute_all(update->ntimestep); - // init each fix - - comm->maxforward_fix = comm->maxreverse_fix = 0; - for (i = 0; i < nfix; i++) fix[i]->init(); - - // set global flag if any fix has its restart_pbc flag set - - restart_pbc_any = 0; - for (i = 0; i < nfix; i++) - if (fix[i]->restart_pbc) restart_pbc_any = 1; - // warn if any particle is time integrated more than once int nlocal = atom->nlocal;