doc/src/compute_rdf.txt

   1 "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
   2
   3 :link(lws,http://lammps.sandia.gov)
   4 :link(ld,Manual.html)
   5 :link(lc,Section_commands.html#comm)
   6
   7 :line
   8
   9 compute rdf command :h3
  10
  11 [Syntax:]
  12
  13 compute ID group-ID rdf Nbin itype1 jtype1 itype2 jtype2 ... :pre
  14
  15 ID, group-ID are documented in "compute"_compute.html command
  16 rdf = style name of this compute command
  17 Nbin = number of RDF bins
  18 itypeN = central atom type for Nth RDF histogram (see asterisk form below)
  19 jtypeN = distribution atom type for Nth RDF histogram (see asterisk form below) :ul
  20
  21 [Examples:]
  22
  23 compute 1 all rdf 100
  24 compute 1 all rdf 100 1 1
  25 compute 1 all rdf 100 * 3
  26 compute 1 fluid rdf 500 1 1 1 2 2 1 2 2
  27 compute 1 fluid rdf 500 1*3 2 5 *10 :pre
  28
  29 [Description:]
  30
  31 Define a computation that calculates the radial distribution function
  32 (RDF), also called g(r), and the coordination number for a group of
  33 particles.  Both are calculated in histogram form by binning pairwise
  34 distances into {Nbin} bins from 0.0 to the maximum force cutoff
  35 defined by the "pair_style"_pair_style.html command.  The bins are of
  36 uniform size in radial distance.  Thus a single bin encompasses a thin
  37 shell of distances in 3d and a thin ring of distances in 2d.
  38
  39 NOTE: If you have a bonded system, then the settings of
  40 "special_bonds"_special_bonds.html command can remove pairwise
  41 interactions between atoms in the same bond, angle, or dihedral.  This
  42 is the default setting for the "special_bonds"_special_bonds.html
  43 command, and means those pairwise interactions do not appear in the
  44 neighbor list.  Because this fix uses the neighbor list, it also means
  45 those pairs will not be included in the RDF. This does not apply when
  46 using long-range coulomb ({coul/long}, {coul/msm}, {coul/wolf} or
  47 similar.  One way to get around this would be to set special_bond
  48 scaling factors to very tiny numbers that are not exactly zero
  49 (e.g. 1.0e-50). Another workaround is to write a dump file, and use
  50 the "rerun"_rerun.html command to compute the RDF for snapshots in the
  51 dump file.  The rerun script can use a
  52 "special_bonds"_special_bonds.html command that includes all pairs in
  53 the neighbor list.
  54
  55 The {itypeN} and {jtypeN} arguments are optional.  These arguments
  56 must come in pairs.  If no pairs are listed, then a single histogram
  57 is computed for g(r) between all atom types.  If one or more pairs are
  58 listed, then a separate histogram is generated for each
  59 {itype},{jtype} pair.
  60
  61 The {itypeN} and {jtypeN} settings can be specified in one of two
  62 ways.  An explicit numeric value can be used, as in the 4th example
  63 above.  Or a wild-card asterisk can be used to specify a range of atom
  64 types.  This takes the form "*" or "*n" or "n*" or "m*n".  If N = the
  65 number of atom types, then an asterisk with no numeric values means
  66 all types from 1 to N.  A leading asterisk means all types from 1 to n
  67 (inclusive).  A trailing asterisk means all types from n to N
  68 (inclusive).  A middle asterisk means all types from m to n
  69 (inclusive).
  70
  71 If both {itypeN} and {jtypeN} are single values, as in the 4th example
  72 above, this means that a g(r) is computed where atoms of type {itypeN}
  73 are the central atom, and atoms of type {jtypeN} are the distribution
  74 atom.  If either {itypeN} and {jtypeN} represent a range of values via
  75 the wild-card asterisk, as in the 5th example above, this means that a
  76 g(r) is computed where atoms of any of the range of types represented
  77 by {itypeN} are the central atom, and atoms of any of the range of
  78 types represented by {jtypeN} are the distribution atom.
  79
  80 Pairwise distances are generated by looping over a pairwise neighbor
  81 list, just as they would be in a "pair_style"_pair_style.html
  82 computation.  The distance between two atoms I and J is included in a
  83 specific histogram if the following criteria are met:
  84
  85 atoms I,J are both in the specified compute group
  86 the distance between atoms I,J is less than the maximum force cutoff
  87 the type of the I atom matches itypeN (one or a range of types)
  88 the type of the J atom matches jtypeN (one or a range of types) :ul
  89
  90 It is OK if a particular pairwise distance is included in more than
  91 one individual histogram, due to the way the {itypeN} and {jtypeN}
  92 arguments are specified.
  93
  94 The g(r) value for a bin is calculated from the histogram count by
  95 scaling it by the idealized number of how many counts there would be
  96 if atoms of type {jtypeN} were uniformly distributed.  Thus it
  97 involves the count of {itypeN} atoms, the count of {jtypeN} atoms, the
  98 volume of the entire simulation box, and the volume of the bin's thin
  99 shell in 3d (or the area of the bin's thin ring in 2d).
 100
 101 A coordination number coord(r) is also calculated, which is the number
 102 of atoms of type {jtypeN} within the current bin or closer, averaged
 103 over atoms of type {itypeN}.  This is calculated as the area- or
 104 volume-weighted sum of g(r) values over all bins up to and including
 105 the current bin, multiplied by the global average volume density of
 106 atoms of type jtypeN.
 107
 108 The simplest way to output the results of the compute rdf calculation
 109 to a file is to use the "fix ave/time"_fix_ave_time.html command, for
 110 example:
 111
 112 compute myRDF all rdf 50
 113 fix 1 all ave/time 100 1 100 c_myRDF\[*\] file tmp.rdf mode vector :pre
 114
 115 [Output info:]
 116
 117 This compute calculates a global array with the number of rows =
 118 {Nbins}, and the number of columns = 1 + 2*Npairs, where Npairs is the
 119 number of I,J pairings specified.  The first column has the bin
 120 coordinate (center of the bin), Each successive set of 2 columns has
 121 the g(r) and coord(r) values for a specific set of {itypeN} versus
 122 {jtypeN} interactions, as described above.  These values can be used
 123 by any command that uses a global values from a compute as input.  See
 124 "Section 6.15"_Section_howto.html#howto_15 for an overview of
 125 LAMMPS output options.
 126
 127 The array values calculated by this compute are all "intensive".
 128
 129 The first column of array values will be in distance
 130 "units"_units.html.  The g(r) columns of array values are normalized
 131 numbers >= 0.0.  The coordination number columns of array values are
 132 also numbers >= 0.0.
 133
 134 [Restrictions:]
 135
 136 The RDF is not computed for distances longer than the force cutoff,
 137 since processors (in parallel) don't know about atom coordinates for
 138 atoms further away than that distance.  If you want an RDF for larger
 139 distances, you can use the "rerun"_rerun.html command to post-process
 140 a dump file and set the cutoff for the potential to be longer in the
 141 rerun script.  Note that in the rerun context, the force cutoff is
 142 arbitrary, since you aren't running dynamics and thus are not changing
 143 your model.  The definition of g(r) used by LAMMPS is only appropriate
 144 for characterizing atoms that are uniformly distributed throughout the
 145 simulation cell. In such cases, the coordination number is still
 146 correct and meaningful.  As an example, if a large simulation cell
 147 contains only one atom of type {itypeN} and one of {jtypeN}, then g(r)
 148 will register an arbitrarily large spike at whatever distance they
 149 happen to be at, and zero everywhere else.  Coord(r) will show a step
 150 change from zero to one at the location of the spike in g(r).
 151
 152 [Related commands:]
 153
 154 "fix ave/time"_fix_ave_time.html
 155
 156 [Default:] none