diff -Naur lammps-9Nov09/doc/Section_commands.html lammps-10Nov09/doc/Section_commands.html --- lammps-9Nov09/doc/Section_commands.html 2009-11-06 08:25:30.000000000 -0700 +++ lammps-10Nov09/doc/Section_commands.html 2009-11-06 10:05:06.000000000 -0700 @@ -371,18 +371,18 @@ none hybrid hybrid/overlay airebo born/coul/long buck buck/coul/cut buck/coul/long colloid coul/cut coul/debye coul/long -dipole/cut dpd eam eam/opt -eam/alloy eam/alloy/opt eam/fs eam/fs/opt -gayberne gayberne/gpu gran/hertz/history gran/hooke -gran/hooke/history lj/charmm/coul/charmm lj/charmm/coul/charmm/implicit lj/charmm/coul/long -lj/charmm/coul/long/opt lj/class2 lj/class2/coul/cut lj/class2/coul/long -lj/cut lj/cut/gpu lj/cut/opt lj/cut/coul/cut -lj/cut/coul/debye lj/cut/coul/long lj/cut/coul/long/tip4p lj/expand -lj/gromacs lj/gromacs/coul/gromacs lj/smooth lj96/cut -lubricate meam morse morse/opt -peri/pmb reax resquared soft -sw table tersoff tersoff/zbl -yukawa +dipole/cut dpd dsmc eam +eam/opt eam/alloy eam/alloy/opt eam/fs +eam/fs/opt gayberne gayberne/gpu gran/hertz/history +gran/hooke gran/hooke/history lj/charmm/coul/charmm lj/charmm/coul/charmm/implicit +lj/charmm/coul/long lj/charmm/coul/long/opt lj/class2 lj/class2/coul/cut +lj/class2/coul/long lj/cut lj/cut/gpu lj/cut/opt +lj/cut/coul/cut lj/cut/coul/debye lj/cut/coul/long lj/cut/coul/long/tip4p +lj/expand lj/gromacs lj/gromacs/coul/gromacs lj/smooth +lj96/cut lubricate meam morse +morse/opt peri/pmb reax resquared +soft sw table tersoff +tersoff/zbl yukawa

These are pair styles contributed by users, which can be used if diff -Naur lammps-9Nov09/doc/Section_commands.txt lammps-10Nov09/doc/Section_commands.txt --- lammps-9Nov09/doc/Section_commands.txt 2009-11-06 08:25:30.000000000 -0700 +++ lammps-10Nov09/doc/Section_commands.txt 2009-11-06 10:05:06.000000000 -0700 @@ -515,6 +515,7 @@ "coul/long"_pair_coul.html, "dipole/cut"_pair_dipole.html, "dpd"_pair_dpd.html, +"dsmc"_pair_dsmc.html, "eam"_pair_eam.html, "eam/opt"_pair_eam.html, "eam/alloy"_pair_eam.html, diff -Naur lammps-9Nov09/doc/fix_imd.html lammps-10Nov09/doc/fix_imd.html --- lammps-9Nov09/doc/fix_imd.html 2009-11-06 08:46:42.000000000 -0700 +++ lammps-10Nov09/doc/fix_imd.html 2009-11-06 09:15:26.000000000 -0700 @@ -39,10 +39,16 @@

Description:

This fix implements the "Interactive MD" (IMD) protocol which allows -to connect an IMD client, for example the VMD visualization program, -to a running LAMMPS simulation and monitor the progress of the simulation -and interactively apply forces to selected atoms. +

This fix implements the "Interactive MD" (IMD) protocol which allows +to connect an IMD client, for example the VMD visualization +program, to a running LAMMPS simulation and monitor the progress +of the simulation and interactively apply forces to selected atoms. +

The source code for this fix includes code developed by the +Theoretical and Computational Biophysics Group in the Beckman +Institute for Advanced Science and Technology at the University of +Illinois at Urbana-Champaign. We thank them for providing a software +interface that allows codes like LAMMPS to hook to VMD.

Upon initialization of the fix, it will open a communication port on the node with MPI task 0 and wait for an incoming connection. @@ -97,39 +103,35 @@ as if they were real objects. See the VMD IMD Homepage for more details.

If IMD control messages are received, a line of text describing the +message and its effect will be printed to the LAMMPS output screen, if +screen output is active. +

+ +

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

The fix does not write any information to restart files. -

None of the fix_modify options are relevant -to this fix. -

If IMD control messages are received, a line of text describing -the message and its effect will be printed to the screen, if -screen output is active. -

No parameter of this fix can be used with the start/stop keywords of -the run command. This fix is not invoked during energy -minimization. +

No information about this fix is written to binary restart +files. None of the fix_modify options +are relevant to this fix. No global scalar or vector or per-atom +quantities are stored by this fix for access by various output +commands. No parameter of this fix can be +used with the start/stop keywords of the run command. +This fix is not invoked during energy minimization.

Restrictions:

This source code for this fix includes code developed by the -Theoretical and Computational Biophysics Group in the Beckman -Institute for Advanced Science and Technology at the University of -Illinois at Urbana-Champaign. -

This fix is part of the "user-imd" package. It is only enabled if LAMMPS was built with that package. See the Making LAMMPS section for more info.

When used in combination with VMD, a topology or coordinate file -has to be loaded, which matches in number and ordering of atoms -the group the fix is applied to. The fix internally sorts atom -IDs by ascending integer value; in VMD (and thus the IMD protocol) -those will be assigned 0-based consecutive index numbers. +

When used in combination with VMD, a topology or coordinate file has +to be loaded, which matches (in number and ordering of atoms) the +group the fix is applied to. The fix internally sorts atom IDs by +ascending integer value; in VMD (and thus the IMD protocol) those will +be assigned 0-based consecutive index numbers.

When using multiple active IMD connections at the same time, each needs to use a different port number. @@ -138,6 +140,4 @@

Default: none

- diff -Naur lammps-9Nov09/doc/fix_imd.txt lammps-10Nov09/doc/fix_imd.txt --- lammps-9Nov09/doc/fix_imd.txt 2009-11-06 08:46:42.000000000 -0700 +++ lammps-10Nov09/doc/fix_imd.txt 2009-11-06 09:15:26.000000000 -0700 @@ -31,10 +31,16 @@ [Description:] -This fix implements the "Interactive MD" (IMD) protocol which allows -to connect an IMD client, for example the VMD visualization program, -to a running LAMMPS simulation and monitor the progress of the simulation -and interactively apply forces to selected atoms. +This fix implements the "Interactive MD" (IMD) protocol which allows +to connect an IMD client, for example the "VMD visualization +program"_VMD, to a running LAMMPS simulation and monitor the progress +of the simulation and interactively apply forces to selected atoms. + +The source code for this fix includes code developed by the +Theoretical and Computational Biophysics Group in the Beckman +Institute for Advanced Science and Technology at the University of +Illinois at Urbana-Champaign. We thank them for providing a software +interface that allows codes like LAMMPS to hook to "VMD"_VMD. Upon initialization of the fix, it will open a communication port on the node with MPI task 0 and wait for an incoming connection. @@ -89,39 +95,34 @@ as if they were real objects. See the "VMD IMD Homepage"_imdvmd for more details. +If IMD control messages are received, a line of text describing the +message and its effect will be printed to the LAMMPS output screen, if +screen output is active. + +:link(VMD,http://www.ks.uiuc.edu/Research/vmd)x :link(imdvmd,http://www.ks.uiuc.edu/Research/vmd/imd/) [Restart, fix_modify, output, run start/stop, minimize info:] -The fix does not write any information to restart files. - -None of the "fix_modify"_fix_modify.html options are relevant -to this fix. - -If IMD control messages are received, a line of text describing -the message and its effect will be printed to the screen, if -screen output is active. - -No parameter of this fix can be used with the {start/stop} keywords of -the "run"_run.html command. This fix is not invoked during "energy -minimization"_minimize.html. +No information about this fix is written to "binary restart +files"_restart.html. None of the "fix_modify"_fix_modify.html options +are relevant to this fix. No global scalar or vector or per-atom +quantities are stored by this fix for access by various "output +commands"_Section_howto.html#4_15. No parameter of this fix can be +used with the {start/stop} keywords of the "run"_run.html command. +This fix is not invoked during "energy minimization"_minimize.html. [Restrictions:] -This source code for this fix includes code developed by the -Theoretical and Computational Biophysics Group in the Beckman -Institute for Advanced Science and Technology at the University of -Illinois at Urbana-Champaign. - This fix is part of the "user-imd" 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. -When used in combination with VMD, a topology or coordinate file -has to be loaded, which matches in number and ordering of atoms -the group the fix is applied to. The fix internally sorts atom -IDs by ascending integer value; in VMD (and thus the IMD protocol) -those will be assigned 0-based consecutive index numbers. +When used in combination with VMD, a topology or coordinate file has +to be loaded, which matches (in number and ordering of atoms) the +group the fix is applied to. The fix internally sorts atom IDs by +ascending integer value; in VMD (and thus the IMD protocol) those will +be assigned 0-based consecutive index numbers. When using multiple active IMD connections at the same time, each needs to use a different port number. @@ -129,5 +130,3 @@ [Related commands:] none [Default:] none - -:line diff -Naur lammps-9Nov09/doc/pair_coeff.html lammps-10Nov09/doc/pair_coeff.html --- lammps-9Nov09/doc/pair_coeff.html 2009-11-04 15:12:52.000000000 -0700 +++ lammps-10Nov09/doc/pair_coeff.html 2009-11-06 10:05:06.000000000 -0700 @@ -98,6 +98,7 @@

pair_style coul/long - long-range Coulombic potential

pair_style dipole/cut - point dipoles with cutoff

pair_style dpd - dissipative particle dynamics (DPD) +

pair_style dsmc - Direct Simulation Monte Carlo (DSMC)

pair_style eam - embedded atom method (EAM)

pair_style eam/opt - optimized version of EAM

pair_style eam/alloy - alloy EAM diff -Naur lammps-9Nov09/doc/pair_coeff.txt lammps-10Nov09/doc/pair_coeff.txt --- lammps-9Nov09/doc/pair_coeff.txt 2009-11-04 15:12:52.000000000 -0700 +++ lammps-10Nov09/doc/pair_coeff.txt 2009-11-06 10:05:06.000000000 -0700 @@ -95,6 +95,7 @@ "pair_style coul/long"_pair_coul.html - long-range Coulombic potential "pair_style dipole/cut"_pair_dipole.html - point dipoles with cutoff "pair_style dpd"_pair_dpd.html - dissipative particle dynamics (DPD) +"pair_style dsmc"_pair_dsmc.html - Direct Simulation Monte Carlo (DSMC) "pair_style eam"_pair_eam.html - embedded atom method (EAM) "pair_style eam/opt"_pair_eam.html - optimized version of EAM "pair_style eam/alloy"_pair_eam.html - alloy EAM diff -Naur lammps-9Nov09/doc/pair_dsmc.html lammps-10Nov09/doc/pair_dsmc.html --- lammps-9Nov09/doc/pair_dsmc.html 1969-12-31 17:00:00.000000000 -0700 +++ lammps-10Nov09/doc/pair_dsmc.html 2009-11-06 10:05:06.000000000 -0700 @@ -0,0 +1,144 @@ + +LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands + + + + + + + +

+ +

pair_style dsmc command +

Syntax: +

pair_style dsmc max_cell_size seed weighting Tref Nrecompute Nsample 
+

max_cell_size = global maximum cell size for DSMC interactions (distance units) +
seed = random # seed (positive integer) +
weighting = macroparticle weighting +
Tref = reference temperature (temperature units) +
Nrecompute = recompute v*sigma_max every this many timesteps (timesteps) +
Nsample = sample this many times in recomputing v*sigma_max +

Examples: +

pair_style dsmc 2.5 34387 10 1.0 100 20
+pair_coeff * * 1.0
+pair_coeff 1 1 1.0 
+

Description: +

Style dsmc computes collisions between pairs of particles for a +direct simulation Monte Carlo (DSMC) model following the exposition in +(Bird). Each collision resets the velocities of the two +particles involved. The number of pairwise collisions for each pair +or particle types and the length scale within which they occur are +determined by the parameters of the pair_style and pair_coeff +commands. +

Stochastic collisions are performed using the variable hard sphere +(VHS) approach, with the user-defined max_cell_size value used as +the maximum DSMC cell size, and reference cross-sections for +collisions given using the pair_coeff command. +

There is no pairwise energy or virial contributions associated with +this pair style. +

The following coefficient must be defined for each pair of atoms types +via the pair_coeff command as in the examples above, +or in the data file or restart files read by the +read_data or read_restart +commands: +

sigma (area units, i.e. distance-squared) +

The global DSMC max_cell_size determines the maximum cell length +used in the DSMC calculation. A structured mesh is overlayed on the +simulation box such that an integer number of cells are created in +each direction for each processor's sub-domain. Cell lengths are +adjusted up to the user-specified maximum cell size. +

+ +

To perform a DSMC simulation with LAMMPS, several additional options +should be set in your input script, though LAMMPS does not check for +these settings. +

Since this pair style does not compute particle forces, you should use +the "fix nve/noforce" time integration fix for the DSMC particles, +e.g. +

fix 1 all nve/noforce 
+

This pair style assumes that all particles will communicated to +neighboring processors every timestep as they move. This makes it +possible to perform all collisions between pairs of particles that are +on the same processor. To ensure this occurs, you should use +these commands: +

neighbor 0.0 bin
+neigh_modify every 1 delay 0 check no
+communicate single cutoff 0.0 
+

These commands insure that LAMMPS communicates particles to +neighboring processors every timestep and that no ghost atoms are +created. The output statistics for a simulation run should indicate +there are no ghost particles or neighbors. +

+ +

Mixing, shift, table, tail correction, restart, rRESPA info: +

This pair style does not support mixing. Thus, coefficients for all +I,J pairs must be specified explicitly. +

This pair style does not support the pair_modify +shift option for the energy of the pair interaction. +

The pair_modify table option is not relevant +for this pair style. +

This pair style does not support the pair_modify +tail option for adding long-range tail corrections to energy and +pressure. +

This pair style writes its information to binary restart +files, so pair_style and pair_coeff commands do not need +to be specified in an input script that reads a restart file. Note +that the user-specified random number seed is stored in the restart +file, so when a simulation is restarted, each processor will +re-initialize its random number generator the same way it did +initially. This means the random forces will be random, but will not +be the same as they would have been if the original simulation had +continued past the restart time. +

This pair style can only be used via the pair keyword of the +run_style respa command. It does not support the +inner, middle, outer keywords. +

+ +

Restrictions: +

This style is part of the "dsmc" package. It is only enabled if +LAMMPS was built with that package. See the Making +LAMMPS section for more info. +

Related commands: +

pair_coeff, fix nve/noforce, +neigh_modify, neighbor, +communicate +

Default: none +

+ + + +

(Bird) G. A. Bird, "Molecular Gas Dynamics and the Direct Simulation +of Gas Flows" (1994). +

+ diff -Naur lammps-9Nov09/doc/pair_dsmc.txt lammps-10Nov09/doc/pair_dsmc.txt --- lammps-9Nov09/doc/pair_dsmc.txt 1969-12-31 17:00:00.000000000 -0700 +++ lammps-10Nov09/doc/pair_dsmc.txt 2009-11-06 10:05:06.000000000 -0700 @@ -0,0 +1,138 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +pair_style dsmc command :h3 + +[Syntax:] + +pair_style dsmc max_cell_size seed weighting Tref Nrecompute Nsample :pre + +max_cell_size = global maximum cell size for DSMC interactions (distance units) +seed = random # seed (positive integer) +weighting = macroparticle weighting +Tref = reference temperature (temperature units) +Nrecompute = recompute v*sigma_max every this many timesteps (timesteps) +Nsample = sample this many times in recomputing v*sigma_max :ul + +[Examples:] + +pair_style dsmc 2.5 34387 10 1.0 100 20 +pair_coeff * * 1.0 +pair_coeff 1 1 1.0 :pre + +[Description:] + +Style {dsmc} computes collisions between pairs of particles for a +direct simulation Monte Carlo (DSMC) model following the exposition in +"(Bird)"_#Bird. Each collision resets the velocities of the two +particles involved. The number of pairwise collisions for each pair +or particle types and the length scale within which they occur are +determined by the parameters of the pair_style and pair_coeff +commands. + +Stochastic collisions are performed using the variable hard sphere +(VHS) approach, with the user-defined {max_cell_size} value used as +the maximum DSMC cell size, and reference cross-sections for +collisions given using the pair_coeff command. + +There is no pairwise energy or virial contributions associated with +this pair style. + +The following coefficient must be defined for each pair of atoms types +via the "pair_coeff"_pair_coeff.html command as in the examples above, +or in the data file or restart files read by the +"read_data"_read_data.html or "read_restart"_read_restart.html +commands: + +sigma (area units, i.e. distance-squared) :ul + +The global DSMC {max_cell_size} determines the maximum cell length +used in the DSMC calculation. A structured mesh is overlayed on the +simulation box such that an integer number of cells are created in +each direction for each processor's sub-domain. Cell lengths are +adjusted up to the user-specified maximum cell size. + +:line + +To perform a DSMC simulation with LAMMPS, several additional options +should be set in your input script, though LAMMPS does not check for +these settings. + +Since this pair style does not compute particle forces, you should use +the "fix nve/noforce" time integration fix for the DSMC particles, +e.g. + +fix 1 all nve/noforce :pre + +This pair style assumes that all particles will communicated to +neighboring processors every timestep as they move. This makes it +possible to perform all collisions between pairs of particles that are +on the same processor. To ensure this occurs, you should use +these commands: + +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no +communicate single cutoff 0.0 :pre + +These commands insure that LAMMPS communicates particles to +neighboring processors every timestep and that no ghost atoms are +created. The output statistics for a simulation run should indicate +there are no ghost particles or neighbors. + +:line + +[Mixing, shift, table, tail correction, restart, rRESPA info]: + +This pair style does not support mixing. Thus, coefficients for all +I,J pairs must be specified explicitly. + +This pair style does not support the "pair_modify"_pair_modify.html +shift option for the energy of the pair interaction. + +The "pair_modify"_pair_modify.html table option is not relevant +for this pair style. + +This pair style does not support the "pair_modify"_pair_modify.html +tail option for adding long-range tail corrections to energy and +pressure. + +This pair style writes its information to "binary restart +files"_restart.html, so pair_style and pair_coeff commands do not need +to be specified in an input script that reads a restart file. Note +that the user-specified random number seed is stored in the restart +file, so when a simulation is restarted, each processor will +re-initialize its random number generator the same way it did +initially. This means the random forces will be random, but will not +be the same as they would have been if the original simulation had +continued past the restart time. + +This pair style can only be used via the {pair} keyword of the +"run_style respa"_run_style.html command. It does not support the +{inner}, {middle}, {outer} keywords. + +:line + +[Restrictions:] + +This style is part of the "dsmc" 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. + +[Related commands:] + +"pair_coeff"_pair_coeff.html, "fix nve/noforce"_fix_nve_noforce.html, +"neigh_modify"_neigh_modify.html, "neighbor"_neighbor.html, +"communicate"_communicate.html + +[Default:] none + +:line + +:link(Bird) +[(Bird)] G. A. Bird, "Molecular Gas Dynamics and the Direct Simulation +of Gas Flows" (1994). diff -Naur lammps-9Nov09/doc/pair_style.html lammps-10Nov09/doc/pair_style.html --- lammps-9Nov09/doc/pair_style.html 2009-11-04 15:12:52.000000000 -0700 +++ lammps-10Nov09/doc/pair_style.html 2009-11-06 10:05:06.000000000 -0700 @@ -100,6 +100,7 @@

pair_style coul/long - long-range Coulombic potential

pair_style dipole/cut - point dipoles with cutoff

pair_style dpd - dissipative particle dynamics (DPD) +

pair_style dsmc - Direct Simulation Monte Carlo (DSMC)

pair_style eam - embedded atom method (EAM)

pair_style eam/opt - optimized version of EAM

pair_style eam/alloy - alloy EAM diff -Naur lammps-9Nov09/doc/pair_style.txt lammps-10Nov09/doc/pair_style.txt --- lammps-9Nov09/doc/pair_style.txt 2009-11-04 15:12:52.000000000 -0700 +++ lammps-10Nov09/doc/pair_style.txt 2009-11-06 10:05:06.000000000 -0700 @@ -97,6 +97,7 @@ "pair_style coul/long"_pair_coul.html - long-range Coulombic potential "pair_style dipole/cut"_pair_dipole.html - point dipoles with cutoff "pair_style dpd"_pair_dpd.html - dissipative particle dynamics (DPD) +"pair_style dsmc"_pair_dsmc.html - Direct Simulation Monte Carlo (DSMC) "pair_style eam"_pair_eam.html - embedded atom method (EAM) "pair_style eam/opt"_pair_eam.html - optimized version of EAM "pair_style eam/alloy"_pair_eam.html - alloy EAM diff -Naur lammps-9Nov09/src/DSMC/Install.csh lammps-10Nov09/src/DSMC/Install.csh --- lammps-9Nov09/src/DSMC/Install.csh 1969-12-31 17:00:00.000000000 -0700 +++ lammps-10Nov09/src/DSMC/Install.csh 2009-11-06 09:48:55.000000000 -0700 @@ -0,0 +1,20 @@ +# Install/unInstall package classes in LAMMPS + +if ($1 == 1) then + + cp style_dsmc.h .. + + cp pair_dsmc.cpp .. + + cp pair_dsmc.h .. + +else if ($1 == 0) then + + rm ../style_dsmc.h + touch ../style_dsmc.h + + rm ../pair_dsmc.cpp + + rm ../pair_dsmc.h + +endif diff -Naur lammps-9Nov09/src/DSMC/pair_dsmc.cpp lammps-10Nov09/src/DSMC/pair_dsmc.cpp --- lammps-9Nov09/src/DSMC/pair_dsmc.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-10Nov09/src/DSMC/pair_dsmc.cpp 2009-11-06 10:05:48.000000000 -0700 @@ -0,0 +1,529 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_dsmc.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "memory.h" +#include "error.h" +#include "domain.h" +#include "update.h" +#include "random_mars.h" +#include "limits.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairDSMC::PairDSMC(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + + total_number_of_collisions = 0; + max_particles = max_particle_list = 0; + next_particle = NULL; + random = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairDSMC::~PairDSMC() +{ + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + memory->destroy_2d_double_array(sigma); + memory->destroy_2d_double_array(cut); + memory->destroy_2d_double_array(V_sigma_max); + memory->destroy_2d_int_array(particle_list); + memory->destroy_2d_int_array(first); + memory->destroy_2d_int_array(number); + } + + delete [] next_particle; + delete random; +} + +/* ---------------------------------------------------------------------- */ + +void PairDSMC::compute(int eflag, int vflag) +{ + double **x = atom->x; + double *mass = atom->mass; + int *type = atom->type; + int nlocal = atom->nlocal; + + for (int i = 1; i <= atom->ntypes; ++i) + for (int j = 0; j < total_ncells; ++j) { + first[i][j] = -1; + number[i][j] = 0; + } + + if (atom->nmax > max_particles) { + delete [] next_particle; + max_particles = atom->nmax; + next_particle = new int[max_particles]; + } + + // find each particle's cell and sort by type + // assume a constant volume and shape simulation domain + // skip particle if outside processor domain + + for (int i = 0; i < nlocal; ++i) { + int xcell = static_cast((x[i][0] - domain->boxlo[0])/cellx); + int ycell = static_cast((x[i][1] - domain->boxlo[1])/celly); + int zcell = static_cast((x[i][2] - domain->boxlo[2])/cellz); + + if ((xcell < 0) or (xcell > ncellsx-1) or + (ycell < 0) or (ycell > ncellsy-1) or + (zcell < 0) or (zcell > ncellsz-1)) continue; + + int icell = xcell + ycell*ncellsx + zcell*ncellsx*ncellsy; + itype = type[i]; + next_particle[i] = first[itype][icell]; + first[itype][icell] = i; + number[itype][icell]++; + } + + for (int icell = 0; icell < total_ncells; ++icell) { + + for (itype = 1; itype <= atom->ntypes; ++itype) { + number_of_A = number[itype][icell]; + if (number_of_A > max_particle_list) { + max_particle_list = number_of_A; + particle_list = memory->grow_2d_int_array(particle_list,atom->ntypes+1, + max_particle_list, + "pair:particle_list"); + } + + int m = first[itype][icell]; + for (int k = 0; k < number_of_A; k++) { + particle_list[itype][k] = m; + m = next_particle[m]; + } + } + + for (itype = 1; itype <= atom->ntypes; ++itype) { + imass = mass[itype]; + number_of_A = number[itype][icell]; + + for (jtype = itype; jtype <= atom->ntypes; ++jtype) { + jmass = mass[jtype]; + number_of_B = number[jtype][icell]; + + reduced_mass = imass*jmass/(imass + jmass); + total_mass = imass + jmass; + jmass_tmass = jmass/total_mass; + imass_tmass = imass/total_mass; + + // if necessary, recompute V_sigma_max values + + if (recompute_vsigmamax_stride && + (update->ntimestep % recompute_vsigmamax_stride == 0)) + recompute_V_sigma_max(icell); + + // # of collisions to perform for itype-jtype pairs + + double &Vs_max = V_sigma_max[itype][jtype]; + double num_of_collisions_double = 0.5 * number_of_A * number_of_B * + weighting * Vs_max * update->dt / vol; + + if ((itype == jtype) and number_of_B) + num_of_collisions_double *= + double(number_of_B - 1) / double(number_of_B); + + int num_of_collisions = + convert_double_to_equivalent_int(num_of_collisions_double); + + if (num_of_collisions > number_of_A) + error->warning("num_of_collisions > number_of_A"); + if (num_of_collisions > number_of_B) + error->warning("num_of_collisions > number_of_B"); + + // perform collisions on pairs of particles in icell + + for (int k = 0; k < num_of_collisions; k++) { + if ((number_of_A < 1) or (number_of_B < 1)) break; + int ith_A = static_cast(random->uniform()*number_of_A); + int jth_B = static_cast(random->uniform()*number_of_B); + int i = particle_list[itype][ith_A]; + int j = particle_list[jtype][jth_B]; + if (i == j) { + k--; + continue; + } + double probability = V_sigma(i,j)/Vs_max; + if (probability > random->uniform()) scatter_random(i,j,icell); + } + } + } + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairDSMC::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); + sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); + V_sigma_max = memory->create_2d_double_array(n+1,n+1,"pair:V_sigma_max"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairDSMC::settings(int narg, char **arg) +{ + if (narg != 6) error->all("Illegal pair_style command"); + + cut_global = 0.0; + max_cell_size = force->numeric(arg[0]); + seed = force->inumeric(arg[1]); + weighting = force->numeric(arg[2]); + T_ref = force->numeric(arg[3]); + recompute_vsigmamax_stride = force->inumeric(arg[4]); + vsigmamax_samples = force->inumeric(arg[5]); + + // initialize Marsaglia RNG with processor-unique seed + + if (max_cell_size <= 0.0) error->all("Illegal pair_style command"); + if (seed <= 0) error->all("Illegal pair_style command"); + if (random) delete random; + random = new RanMars(lmp,seed + comm->me); + + kT_ref = force->boltz*T_ref; + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairDSMC::coeff(int narg, char **arg) +{ + if (narg < 3 || narg > 4) error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double sigma_one = force->numeric(arg[2]); + + double cut_one = cut_global; + if (narg == 4) cut_one = force->numeric(arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + sigma[i][j] = sigma_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairDSMC::init_style() +{ + ncellsx = ncellsy = ncellsz = 1; + while (((domain->boxhi[0] - domain->boxlo[0])/ncellsx) > max_cell_size) + ncellsx++; + while (((domain->boxhi[1] - domain->boxlo[1])/ncellsy) > max_cell_size) + ncellsy++; + while (((domain->boxhi[2] - domain->boxlo[2])/ncellsz) > max_cell_size) + ncellsz++; + + cellx = (domain->boxhi[0] - domain->boxlo[0])/ncellsx; + celly = (domain->boxhi[1] - domain->boxlo[1])/ncellsy; + cellz = (domain->boxhi[2] - domain->boxlo[2])/ncellsz; + + if (comm->me == 0) { + if (screen) fprintf(screen,"DSMC cell size = %g x %g x %g\n", + cellx,celly,cellz); + if (logfile) fprintf(logfile,"DSMC cell size = %g x %g x %g\n", + cellx,celly,cellz); + } + + total_ncells = ncellsx*ncellsy*ncellsz; + vol = cellx*celly*cellz; + + particle_list = memory->create_2d_int_array(atom->ntypes+1,0, + "pair:particle_list"); + first = memory->create_2d_int_array(atom->ntypes+1,total_ncells, + "pair:first"); + number = memory->create_2d_int_array(atom->ntypes+1,total_ncells, + "pair:number"); + + for (int i = 1; i <= atom->ntypes; i++) + for (int j = 1; j <= atom->ntypes; j++) + V_sigma_max[i][j] = 0.0; + + two_pi = 8.0*atan(1.0); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairDSMC::init_one(int i, int j) +{ + if (setflag[i][j] == 0) cut[i][j] = 0.0; + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairDSMC::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&sigma[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairDSMC::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&sigma[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairDSMC::write_restart_settings(FILE *fp) +{ + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&max_cell_size,sizeof(double),1,fp); + fwrite(&seed,sizeof(int),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairDSMC::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_global,sizeof(double),1,fp); + fread(&max_cell_size,sizeof(double),1,fp); + fread(&seed,sizeof(int),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&max_cell_size,1,MPI_DOUBLE,0,world); + MPI_Bcast(&seed,1,MPI_INT,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + + // initialize Marsaglia RNG with processor-unique seed + // same seed that pair_style command initially specified + + if (random) delete random; + random = new RanMars(lmp,seed + comm->me); +} + +/*------------------------------------------------------------------------- + rezero and recompute the V_sigma_max values this timestep for use during + the next nrezero timesteps +-------------------------------------------------------------------------*/ + +void PairDSMC::recompute_V_sigma_max(int icell) +{ + int i,j,k; + double Vsigma_max = 0; + + if (number_of_A && number_of_B) { + for (k = 0; k < vsigmamax_samples; k++) { + i = particle_list[itype] + [static_cast(random->uniform()*number_of_A)]; + j = particle_list[jtype] + [static_cast(random->uniform()*number_of_B)]; + if (i == j) continue; + Vsigma_max = MAX(Vsigma_max,V_sigma(i,j)); + } + } + V_sigma_max[itype][jtype] = Vsigma_max; +} + +/*------------------------------------------------------------------------- + VHS model + compute the velocity vector difference between i and j and multiply by + their combined collision cross section, sigma, for neutral-neutral + collisions using the Variable Hard Sphere model +-------------------------------------------------------------------------*/ + +double PairDSMC::V_sigma(int i, int j) +{ + double relative_velocity_sq,relative_velocity,pair_sigma; + double delv[3]; + double *vi = atom->v[i]; + double *vj = atom->v[j]; + + subtract3d(vi,vj,delv); + relative_velocity_sq = dot3d(delv,delv); + relative_velocity = sqrt(relative_velocity_sq); + + // from Bird eq 4.63, and omega=0.67 + // (omega - 0.5) = 0.17 + // 1/GAMMA(2.5 - omega) = 1.06418029298371 + + if (relative_velocity_sq != 0.0) + pair_sigma = sigma[itype][jtype]* + pow(kT_ref/(0.5*reduced_mass*relative_velocity_sq),0.17) * + 1.06418029298371; + else + pair_sigma = 0.0; + + return relative_velocity*pair_sigma; +} + +/*------------------------------------------------------------------------- + generate new velocities for collided particles +-------------------------------------------------------------------------*/ + +void PairDSMC::scatter_random(int i, int j, int icell) +{ + double mag_delv,cos_phi,cos_squared,r,theta; + double delv[3],vcm[3]; + double *vi = atom->v[i]; + double *vj = atom->v[j]; + + subtract3d(vi,vj,delv); + if (itype == jtype) mag_delv = sqrt(dot3d(delv,delv))*0.5; + else mag_delv = sqrt(dot3d(delv,delv)); + + cos_phi = 1.0 - (2.0*random->uniform()); + cos_squared = MIN(1.0,cos_phi*cos_phi); + r = sqrt(1.0 - cos_squared); + delv[0] = cos_phi*mag_delv; + theta = two_pi*random->uniform(); + delv[1] = r*mag_delv*cos(theta); + delv[2] = r*mag_delv*sin(theta); + + if (itype == jtype) { + vcm[0] = (vi[0]+vj[0])*0.5; + vcm[1] = (vi[1]+vj[1])*0.5; + vcm[2] = (vi[2]+vj[2])*0.5; + vi[0] = vcm[0] + delv[0]; + vi[1] = vcm[1] + delv[1]; + vi[2] = vcm[2] + delv[2]; + vj[0] = vcm[0] - delv[0]; + vj[1] = vcm[1] - delv[1]; + vj[2] = vcm[2] - delv[2]; + } else { + vcm[0] = vi[0]*imass_tmass + vj[0]*jmass_tmass; + vcm[1] = vi[1]*imass_tmass + vj[1]*jmass_tmass; + vcm[2] = vi[2]*imass_tmass + vj[2]*jmass_tmass; + vi[0] = vcm[0] + delv[0]*jmass_tmass; + vi[1] = vcm[1] + delv[1]*jmass_tmass; + vi[2] = vcm[2] + delv[2]*jmass_tmass; + vj[0] = vcm[0] - delv[0]*imass_tmass; + vj[1] = vcm[1] - delv[1]*imass_tmass; + vj[2] = vcm[2] - delv[2]*imass_tmass; + } + + total_number_of_collisions++; +} + +/* ---------------------------------------------------------------------- + This method converts the double supplied by the calling function into + an int, which is returned. By adding a random number between 0 and 1 + to the double before converting it to an int, we ensure that, + statistically, we round down with probability identical to the + remainder and up the rest of the time. So even though we're using an + integer, we're statistically matching the exact expression represented + by the double. +------------------------------------------------------------------------- */ + +int PairDSMC::convert_double_to_equivalent_int(double input_double) +{ + if (input_double > INT_MAX) + error->all("Tried to convert a double to int, but input_double > INT_MAX"); + + int output_int = static_cast(input_double + random->uniform()); + return output_int; +} diff -Naur lammps-9Nov09/src/DSMC/pair_dsmc.h lammps-10Nov09/src/DSMC/pair_dsmc.h --- lammps-9Nov09/src/DSMC/pair_dsmc.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-10Nov09/src/DSMC/pair_dsmc.h 2009-11-06 09:48:55.000000000 -0700 @@ -0,0 +1,102 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef PAIR_DSMC_H +#define PAIR_DSMC_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairDSMC : public Pair { + public: + PairDSMC(class LAMMPS *); + virtual ~PairDSMC(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + virtual void write_restart_settings(FILE *); + virtual void read_restart_settings(FILE *); + + private: + double cut_global; + double **cut; + double **sigma; + + double cellx; + double celly; + double cellz; + int ncellsx; + int ncellsy; + int ncellsz; + int total_ncells; + int total_number_of_collisions; + int recompute_vsigmamax_stride; + int vsigmamax_samples; + double T_ref; + double kT_ref; + double two_pi; + double max_cell_size; + + int seed; + int number_of_A; + int number_of_B; + int max_particle_list; + + class RanMars *random; + + int **particle_list; + int **first; + int **number; + + double **V_sigma_max; + + int max_particles; + int *next_particle; + + int itype; + int jtype; + + double imass; + double jmass; + double total_mass; + double reduced_mass; + double imass_tmass; + double jmass_tmass; + double vol; + double weighting; + + void allocate(); + void recompute_V_sigma_max(int); + double V_sigma(int, int); + void scatter_random(int, int, int); + int convert_double_to_equivalent_int(double); + + inline void subtract3d(const double *v1, const double *v2, double *v3) { + v3[0] = v2[0] - v1[0]; + v3[1] = v2[1] - v1[1]; + v3[2] = v2[2] - v1[2]; + } + + inline double dot3d(const double *v1, const double *v2) { + return( v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] ); + } +}; + +} + +#endif diff -Naur lammps-9Nov09/src/DSMC/style_dsmc.h lammps-10Nov09/src/DSMC/style_dsmc.h --- lammps-9Nov09/src/DSMC/style_dsmc.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-10Nov09/src/DSMC/style_dsmc.h 2009-11-06 09:48:55.000000000 -0700 @@ -0,0 +1,20 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PairInclude +#include "pair_dsmc.h" +#endif + +#ifdef PairClass +PairStyle(dsmc,PairDSMC) +#endif diff -Naur lammps-9Nov09/src/Makefile lammps-10Nov09/src/Makefile --- lammps-9Nov09/src/Makefile 2009-11-06 08:38:50.000000000 -0700 +++ lammps-10Nov09/src/Makefile 2009-11-06 09:48:55.000000000 -0700 @@ -13,7 +13,7 @@ # Package variables -PACKAGE = asphere class2 colloid dipole dpd gpu granular \ +PACKAGE = asphere class2 colloid dipole dpd dsmc gpu granular \ kspace manybody meam molecule opt peri poems prd reax xtc PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn \ diff -Naur lammps-9Nov09/src/style.h lammps-10Nov09/src/style.h --- lammps-9Nov09/src/style.h 2009-10-29 16:45:27.000000000 -0600 +++ lammps-10Nov09/src/style.h 2009-11-06 09:48:55.000000000 -0700 @@ -369,26 +369,6 @@ RegionStyle(union,RegUnion) #endif -// style files for standard packages - -#include "style_asphere.h" -#include "style_class2.h" -#include "style_colloid.h" -#include "style_dipole.h" -#include "style_dpd.h" -#include "style_gpu.h" -#include "style_granular.h" -#include "style_kspace.h" -#include "style_manybody.h" -#include "style_meam.h" -#include "style_molecule.h" -#include "style_opt.h" -#include "style_peri.h" -#include "style_poems.h" -#include "style_prd.h" -#include "style_reax.h" -#include "style_xtc.h" - // package and user add-ons #include "style_packages.h" diff -Naur lammps-9Nov09/src/style_packages.h lammps-10Nov09/src/style_packages.h --- lammps-9Nov09/src/style_packages.h 2009-09-25 09:15:47.000000000 -0600 +++ lammps-10Nov09/src/style_packages.h 2009-11-06 09:48:55.000000000 -0700 @@ -11,13 +11,14 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -// style flies for packages +// style files for standard packages #include "style_asphere.h" #include "style_class2.h" #include "style_colloid.h" #include "style_dipole.h" #include "style_dpd.h" +#include "style_dsmc.h" #include "style_gpu.h" #include "style_granular.h" #include "style_kspace.h" @@ -27,5 +28,6 @@ #include "style_opt.h" #include "style_peri.h" #include "style_poems.h" +#include "style_prd.h" #include "style_reax.h" #include "style_xtc.h"