Algorithm for anisotropic diffusion in hydrogen-bonded networks.
Edoardo Milotti
Dipartimento di Fisica, Università di Trieste,
and INFN – Sezione di Trieste, Via Valerio, 2, I-34127 Trieste, Italy∗
(Dated: November 4, 2018)
Abstract
In this paper I describe a specialized algorithm for anisotropic diffusion determined by a field
of transition rates. The algorithm can be used to describe some interesting forms of diffusion that
occur in the study of proton motion in a network of hydrogen bonds. The algorithm produces data
that require a nonstandard method of spectral analysis which is also developed here. Finally, I
apply the algorithm to a simple specific example.
PACS numbers: 05.40.-a, 02.70.-c, 83.10.Rs
∗Electronic address: milotti@ts.infn.it
http://arxiv.org/abs/0704.0623v1
mailto:milotti@ts.infn.it
I. INTRODUCTION
Protons migrating in water have an anomalously high mobility [1] and their diffusion is
actually limited by the continuous rearrangement of hydrogen bonds [1, 2]. Indeed protons
migrating in ice move faster than protons in water, as the transition rate from one water
molecule to the next is enhanced by the higher molecular order in ice. Proton mobility
increases whenever water molecules are constrained, as in carbon nanotubes [3, 4]. Local
electric fields also orientate water molecules, and thus should lead to a local increase of proton
mobility, and indeed it is now known that there is a definite water dipole orientational order
in the hydration water close to ionizable residues in hydrated proteins [5, 6, 7, 8]: this is
a collective property, which is somewhat independent of the individual fluctuations of the
water dipoles.
Here I am not concerned with the detailed simulation of proton motion which is the sub-
ject of several specialized papers like [9, 10, 11], but rather I wish to set up the framework
for a simulation of the random walks performed by protons in some interesting context, like
proton migration on the surface of hydrated proteins [13]. The basic idea is that protons
move faster in the network of hydrogen bonds just where there is a higher molecular order,
i.e., the transition rate is higher where there is higher spatial order, and because of the
continuous rearrangement of the water molecules which make up a fluctuating bond struc-
ture, the random walk performed by a single proton can actually be viewed as a walk in
continuous space and continuous time, as long as the time resolution of the process is longer
than the relaxation time of water dipole motion.
Here I take for granted that there is some induced order in the hydrogen bond network,
like the dipole field described in [5, 6], and I introduce a corresponding field of transition
rates γ(r, t), such that the time-dependent probability density ρ(r, t) and the associated
probability ∆p = ρ∆V of finding a random walker (a proton) in the small volume ∆V at
position r and time t, yield the following equation for the decrease of ∆p, due to random
walker escape from the region,
∆p(r, t)
= −γ(r, t)∆p(r, t) (1)
I also assume that γ(r, t) is a continuous, differentiable function.
The situation is illustrated in figure 1, which shows a random subdivision of a plane
region: a set of positions – marked by the large black dots – is associated to small surrounding
regions; the arrows in the figure mark the flow of random walkers in the central region to and
from the bordering regions. If the area (actually, the line length in this 2D representation) of
the interface between the central region and the kj-th region is Ai,kj , and the total interface
area of the kj-th region is Akj , then it is easy to see that the total derivative of ∆pi is
= −γi∆pi +
Ai,kj
γkj∆pkj (2)
and the global flow in this discretized system is described by a system of coupled linear
differential equations.
More importantly, we can define currents for the inflow and outflow of random walkers
from a modified form of Fick’s law
Ji→k = α
γiρi − γkρk
where ∆i,k denotes the distance between the centroids of the bordering i-th and k-th region,
and the parameter α is akin to the diffusion coefficient, but is measured in different units
(it has the dimensions of a surface) [14]. The current in the previous formula is actually a
projection along the direction that connects the centroids of the bordering region and it is
easy to generalize to the continuum case and find the outflowing current
J = −α∇(γρ) (4)
so that finally one finds the following Fokker-Planck equation from the conservation of the
total number of random walkers
= α∇2(γρ) (5)
assuming that α does not depend on position.
In the following sections I describe an algorithm to simulate this kind of diffusive motion:
first I discuss the angular distribution, then confinement to motion on surfaces, and in
section IV I show how to extend the algorithm for asynchronous updates. In section V I
give a recipe to analyze asynchronous data. In section VI I discuss a simple example, and
finally in section VII I give a short summary and outlook for the utilization of the algorithm.
II. ANGULAR DISTRIBUTION
From equation (4) we see that the current actually contains two contributions
J = −α∇(γρ) = −α (ρ∇γ + γ∇ρ)) (6)
however when we consider the problem at hand – namely, the diffusion of protons in the
network of hydrogen bonds, and we remark that we wish to describe the individual proton
motion, then we notice that we are only interested in situations where ∇ρ = 0. In fact, pro-
tons repel other protons that are too close, and obey a sort of effective exclusion principle –
which is actually independent of their fermionic nature – and the position of the individual
proton corresponds to a peak of the instantaneous probability density: therefore the current
defined in (4) has the same direction as ∇γ in all cases of practical interest. This direction
corresponds to the average proton motion, but for a single transition to a nearby site it can
only define the axis of an angular probability distribution. Here I make the simplest possi-
ble choice, namely that the angular probability distribution is a simple dipole distribution
defined by the normalized conditional probability density for the unit vector n
P (n|n0) =
(1 + ∆Pn · n0) (7)
where I0 is the normalization factor (I0 = 2π in the 2D case and I0 = 4π in the 3D case),
and n0 is the unit vector
n0(r) =
so that the decrease of the density ρ due to the flow in the angular range ∆Ω, during a given
time interval ∆t, is
(1 + ∆Pn · n0)∆t∆Ω (9)
The constant inside the parenthesis corresponds to the isotropic loss term, while the other
term is associated to the current (6). From a comparison of the elementary flows of random
walkers in direction n we find
J · n = γρ
n · n0 (10)
so that
∆P = I0α
and the conditional angular probability density is
P (n|n0) =
1 + I0α
n · n0
The conditional angular probability density (12) can be used to generate random walks
discarding the time information. Here I take the following time-independent expression for
the transition rate
γ(r, t) = γ(r) =
|r|2 + Γ2
+B (13)
which has an obvious symmetry center, located in the origin, which corresponds to the po-
sition of the peak value as well. This transition rate is motivated by the considerations put
forward in the introduction: if protons migrate in a hydrogen bonded network with polariza-
tion centers that create partial ice-like order in their neighborhood, then the transition rate
(13) is highest, and saturates, close to the polarization centers, and decays to a constant
value with a 1/r2 behavior for r ≫ Γ (i.e., it has a radial dependence like the potential of
electric dipole fields). Notice also that the anisotropy coefficient is
|r|2 + Γ2
The techniques to generate random angles which are distributed according to the probability
density (12) are reviewed in appendix A, and figures 2-5 show some examples: in these
examples all length and distance units are in arbitrary units. Figure 2 shows random walks
around a single center with transition rate (13): the random walker starts at the origin,
with a fixed step length = 0.005 arbitrary units; the horizontal and vertical scales are also
labeled with the same arbitrary length units; the parameters of the transition rate function
are the same in these simulations, A = 1, B = 0.1, and Γ = 1, while α changes in the three
cases displayed in the figure. Larger values of α correspond to higher anisotropy, and we see
that as the anisotropy grows, the random walk becomes more and more compact.
Figure 3 shows a random walk with two centers at positions r1 = (−1, 0), r2 = (1, 0)
(arbitrary units): the random walker starts at the origin, with a fixed step length = 0.01
arbitrary units; the horizontal and vertical scales are also labeled with the same arbitrary
length units. In this case the transition rate is similar to (13), but with two centers,
γ(r) =
|r− r1|2 + Γ21
|r− r2|2 + Γ22
+B (15)
with A1 = A2 = 1, B = 0.1, and Γ1 = Γ2 = 1, and α = 1. Here the random walker explores
the regions around both centers.
Figure 4 shows a situation which is similar to figure 3, although it is more complex. The
transition rate is once again similar to (13), but now it has ten centers,
γ(r) =
|r− rk|2 + Γ2k
+B (16)
with Ak = 1, B = 0.1, and Γk = 0.1, and α = 0.025; the step length is = 0.01 arbitrary
units. The centers are scattered randomly, with a lower bound on the minimum distance
between them; the figure shows three snapshots at different times in the simulation, as the
random walker starts from the center of the figure, drifts to one of the centers and later
migrates to other neighboring centers.
Finally figure 5 shows a random walk in space about two centers at r1 = (−1, 0, 0),
r2 = (1, 0, 0) (arbitrary units), which is very similar to the random walk in figure 3: the
transition rate is still given by expression (15), with A1 = A2 = 1, B = 0.1, and Γ1 = Γ2 = 1,
and α = 0.5, with a fixed step length = 0.1 arbitrary units. Once again the random walker
explores the regions around both centers.
III. DIFFUSION ON SURFACES
In many cases it is important to confine the motion of the random walkers to some
particular portion of space, for instance in the case of protons on hydrated proteins the
motion is confined to the thin hydration layer. The simulation method outlined in the
previous section can be adapted to provide such a confinement to a surface: in this case one
can define at each step the tangent plane at the position of the random walker and proceed
as in the 2D case. Obviously the gradient of the transition rate γ used in the formulas of the
previous section must be projected on the tangent plane, and moreover the directions must
be generated according to the 2D angular distribution (see appendix B). Figure 6 shows
a random walk on a spherical surface with 10 centers as in the examples in the previous
section: the sphere has radius 1 (arb. units), the transition rate is given by expression (16),
with Ak = 10, B = 1, and Γk = 0.25, and α = 0.2; the step length is = 0.1 arbitrary units.
IV. ASYNCHRONOUS TRANSITIONS
In the previous sections we have discussed the space behavior of the random walks, but
obviously we can use the transition rate function γ(r, t) to describe the time behavior as
well. It is possible to choose a fixed time step ∆t and use the transition rates to compute the
probability of generating a transition in the time interval (t, t+∆t) (synchronous update),
however this is inconvenient if the function γ(r, t) spans a wide range of values, because it
means that the choice ∆t ≪ min(1/γ(r, t)) which is required for an accurate simulation,
produces very long waiting times where the transition rate is very small (and therefore, very
large amounts of sampled data). It is actually much more practical to use the transition
rate to generate directly the transition times of each step, which we assume to be indepen-
dent (asynchronous update) from one another. With this – rather natural – assumption of
independency, it is very simple to generate the transition times, as explained in appendix B,
although this leads to uneven sampling, and requires a specialized form of spectral analysis.
V. FOURIER ANALYSIS OF ASYNCHRONOUS DATA
Using asynchronous sampling times it is not possible to use the standard Fourier or other
similar spectral analysis techniques [15]. However the signal produced by the time-domain
simulation is “exact”, at least in the sense that there are no algorithmic artifacts due to
sampling and it is desirable to extract as much information as possible. To this end, I
notice that any function of the position of the random walkers must be stationary between
successive transitions, and that it is possible to make direct use of the definition of Fourier
transform
F (ω) =
s(t)e−iωtdt (17)
where s(t) is any signal produced in the simulation, which depends on the positions of the
random walkers (e.g., a component of the electric dipole moment if the random walkers are
charged particles). The signal s(t) has the fixed value sn in the time interval (tn, tn + 1),
where tn is the time of the n-th transition, and we find:
F (ω) =
∫ tn+1
e−iωtdt
−iω(tn+1+tn)/2(tn+1 − tn)
sin [ω(tn+1 − tn)/2]
ω(tn+1 − tn)/2
sn sin
ω(tn+1 − tn)
ω(tn+1 + tn)
sn sin
ω(tn+1 − tn)
ω(tn+1 + tn)
Using equation (19) the Fourier transform can be evaluated exactly for all frequencies, and
without aliasing: in practice this is possible, practical, and actually useful only for a small
finite set of frequencies. If we had used a Discrete Fourier Transform (DFT) algorithm
[15] to analyze N real samples, we would have found N/2 independent Fourier coefficients,
and using a Fast Fourier Transform (FFT) algorithm the time complexity of the calculation
would be O(N logN). If we use the algorithm defined by equation (19) to compute M values
(M ≤ N/2) of the Fourier transform, the time complexity is clearly O(NM), so that in a
practical calculation we can only compute a reduced number of Fourier coefficients. However,
I remark that in addition to being exact, the algorithm has another major advantage over
the standard DFT calculations: there is no limitation to the set of frequencies that can be
computed, and in particular one can choose a set of frequencies that is not evenly spaced and
that is denser close to the origin, which is particularly useful in this case since the random
walk – when considered as a noise process – is expected to produce a spectrum with a large
power-law peak at low frequencies.
This kind of analysis is actually limited by the finite time span of the generated signal:
we see from equation (18) that the Fourier transform of the generated signal is a sum
of sinc functions, and therefore it is not useful to represent the transform for frequencies
lower than ωmin = π/T where T is the signal duration and ωmin is the lowest positive
zero of the corresponding sinc function. With this limitation, we can sample the Fourier
transform at frequency values that are evenly spaced on a logarithmic scale and obtain a
better representation of the transform close to the origin than is possible with conventional
methods.
In this approach we evaluate the Fourier transform of the simulated signal in the time
interval (0, tN) and we implicitly assume that the signal vanishes outside this interval: this
is different from the standard (implicit) assumption in standard DFT analysis, where the
observed signal is assumed to repeat periodically outside the observation interval [15]. If we
introduce a rectangular window with a width equal to the observation interval (tN), we see
that the present method returns a Fourier transform that is the convolution of the transform
of the signal with the transform of the rectangular window, which is
exp(−iωt)dt = tN exp(−iωtN/2)
sin(ωtN/2)
(ωtN/2)
As a consequence of the convolution associated to the rectangular window we see that a
constant nonzero level produces a sharp peak centered at zero frequency, with a shape given
by equation (20); this peak corresponds to the standard DC peak in DFT analysis, and has
tails with a 1/ω2 spectral behavior that may mimic the low-frequency behavior of a standard
Debye relaxation with a very small decay rate. The mean level of the simulated signal is
∫ tn+1
sn (tn+1 − tn) (21)
and thus we can correct for the DC peak by subtracting its transform
exp(−iωtN/2) sin(ωtN/2) (22)
from the signal transform.
VI. RANDOM WALK ABOUT A SINGLE CENTER
As an example, I consider here a complete simulation (3D space and time data) for
random walks about a single center (as defined by the transition rate (13) ) located in the
origin. The transition rate function used in this example is given again by expression (13),
with A = 100, B = 10−8, Γ = 0.1, α = 0.001, and with a step length = 0.01. The values
of A and B have been chosen to maximize the range of γ sampled by the random walker,
while still keeping a rather short simulation time. I have generated 500 random walks and
10000 transitions for each walk; the random walker always starts at the origin (where the
center of (13) is located). The results of the simulation are shown in figures 7-11: figure 7
is the superposition of a few walks, and it is qualitatively clear that the density profile is
very similar to that of the standard random walk in the plane. Figure 8 shows instead the x
projection of the position signal vs. time, and this is not very different, e.g., from an electric
dipole component if the random walkers are charged particles. The insets show parts of the
signal with increasing magnification, and the last inset displays clearly the stationary parts of
the signal between successive transitions. Figure 9 shows the (unnormalized) distribution of
time intervals between successive transitions: the figure demonstrates clearly that although
the times have been generated according to the interval distribution in appendix B, the
distribution in figure 9 is not a simple exponential, but rather it contains two different
power-law regions (marked by the dotted lines in this log-log plot), which reflects the way in
which the transition probability function is sampled by the random walks. One well-known
property of ordinary random walks is that their mean square radius is proportional to time,
i.e., 〈r2〉 ∝ t: here we see (figure 10) that this linearity is recovered only asymptotically, as
random walkers explore regions that are far away from the origin. Finally, I have used the
power spectral estimation method of section V to analyze the x position signals (as those
in figure 8): the result is shown in figure 11. Figure 11a is the power spectrum obtained in
a single realization of the random walk, while figure 11b is the average of 400 spectra. In
each part of figure 11, the thin gray line represents an ideal power-law spectral density with
the same slope as the the average of 400 spectra, i.e., a 1/f 2 spectrum, which is usually
expected in these types of processes. In fact the random walks simulated here effectively
sample asymptotically only a rather limited range of transition rates – even though the A/B
ratio is quite high – and this means that the usual superposition argument that leads to
more general power-law spectra [16] does not apply here and it is quite natural to find a
1/f 2 spectrum.
VII. DISCUSSION
Before concluding this paper it is important to note that the correct continuum formu-
lation of the diffusion equation in an inhomogeneous environment has been the subject of
much discussion in the past and is still debated (see [17, 18], and references therein). The
generalization is unclear because the microscopic details seem to matter [17]. Moreover,
there is also some interest towards the diffusion equation in various forms of anomalous
diffusion [19]. Here I wish to stress again that the results presented in this paper are spe-
cialized and are meant to address diffusion in structures like those described in [5, 6], unlike
other approaches described in the existing literature that deal with more general diffusion
problems [20, 21, 22]; still, the diffusion equation (5) is similar to equation (5) in reference
[17], and therefore it is interesting to give one further look at its structure. In section II we
have seen that the current (4) has two components, and the gradient term that generates
the random walks discussed here roughly corresponds to the so-called “spurious” drift term
(using the terminology of reference [17]). The other term in the current, namely −αγ∇ρ
has been neglected because the single random walkers considered here are charged fermions
and obey a sort of effective exclusion principle. Indeed, in a context like that of [5, 6] pro-
tons repel because of their charge, while their spin structure, and therefore also their true
fermionic character – and the Pauli exclusion principle – do not matter much; this effective
exclusion principle has been used in the past for an Ising-like modeling of proton motion,
where the presence or absence of a proton at a given position is treated like a pseudo-spin
variable (see the discussion in the review paper [12]). The situation would be very different
if space could be filled by a cloud of random walkers: in a case such as this – which roughly
corresponds to random walkers that effectively behave as bosons – the neglected term should
be included. Thus the actual importance of the different terms of the diffusion equation (5)
depends on the bosonic or fermionic character of the random walkers.
One prominently missing term in the diffusion equation (5) is the usual drift term asso-
ciated to external fields, however if we look at the structure of the current (3), we see that
we can easily produce a flow unbalance with a space- (and possibly time-) dependent α, so
that we obtain a modified diffusion equation
= ∇ · [α∇(γρ)] = ∇α · ∇(γρ) + α∇2(γρ) (23)
with an additional α-dependent drift term.
A final, important comment, is that the algorithm presented here has a sort of backward
approach with respect to other existing algorithms for random walks and diffusion in inho-
mogeneous environments, as it starts directly from the transition rates, instead of deriving
them from a given diffusion equation (as, e.g., [20]).
To conclude, in this paper I have described a novel algorithm for anisotropic diffusion,
which is continuous both in space and time, and I have discussed its application to a simple
example, in anticipation of further work that shall be carried out in a realistic simulation of
noise in proton conduction in weakly hydrated proteins [23].
APPENDIX A: ANGULAR DISTRIBUTIONS
Here I consider first the planar angular distribution defined by the normalized probability
density
P (ϕ; g) =
(1 + g cosϕ) (A1)
Using the standard inversion method (described in many textbooks, see, e.g., [24]), one finds
that the solution ϕ of the nonlinear equation
(ϕ+ g sinϕ) +
= y (A2)
has the distribution described by (A1) if y is a uniform variate on the (0, 1) interval.
The generation of a random direction in space from the normalized probability density
P (n; g) =
(1 + gn · n0) (A3)
requires two angles, a zenithal angle θ and azimuthal angle ϕ, where n0 defines the zenithal
axis, so that the probability of finding a unit vector n in the interval (θ, θ+dθ), and (ϕ, ϕ+dϕ)
dP (θ, ϕ; g) =
(1 + g cos θ) d cos θ
(1 + gx) dx
i.e., the probability density is the product of two independent densities, one uniform with
respect to ϕ on the (0, 2π) interval, and the other linear with respect to x = cos θ on the
(−1, 1) interval. Using again the inversion method, one finds that
has the required linear distribution if y is a uniform variate on the (0, 1) interval. Since n0 is
not usually parallel to the z direction, two rotations are also required: one first rotates the
reference frame so that n0 is parallel to the z axis, this is followed by the angle generation
step, and finally one must transform back to the original reference frame.
APPENDIX B: TIME DISTRIBUTION
Time transitions are generated according to the exponential distribution, which has the
probability density
= γ exp(−γt) (B1)
where γ is the transition rate. The standard inversion method [24] can be used again to
generate times
t = −
ln y (B2)
that are exponentially distributed if y is a uniform variate on the (0, 1) interval.
ACKNOWLEDGMENTS
I wish to thank warmly Giorgio Careri for his insightful comments and suggestions: he
was the first to pinpoint the problem that spurred the research described here, and this
paper would not have been written without his encouragement.
I also wish to thank Alessio Del Fabbro for his careful reading of the manuscript and for
several interesting discussions.
[1] N. Agmon, Chem. Phys. Lett. 244 (1995) 456.
[2] N. Agmon, J. Chim. Phys. (Paris) 93 (1996) 1714.
[3] D. J. Mann and M. D. Halls, Phys. Rev. Lett. 90 (2003) 195503.
[4] R. Jay Mashi, Sony Joseph, N. R. Aluru, and Eric Jakobsson, Nano Lett. 3 (2003) 589.
[5] J. Higo, M Sasai, H. Shirai, H. Nakamura, and T. Kugimiya, PNAS 98 (2001) 5961.
[6] J. Higo and M. Nakasako, J. Comput. Chem. 23 (2002) 1323.
[7] T. Yokomizo, M. Nakasako, T. Yamazaki, H. Shindo, and J. Higo, Chem. Phys. Lett. 401
(2005) 332.
[8] P. Kumar, G. Franzese, S. V. Buldyrev, and H. E. Stanley, Phys. Rev. E 73 (2006) 041505.
[9] J. Halding and P. S. Lomdahl, Phys. Rev. A 37 (1988) 2608.
[10] I. Chochliouros and J. Pouget, J. Phys.: Condensed Matter 7 (1995) 8741.
[11] U. W. Schmitt and G. A. Voth, J. Chem. Phys. 111 (1999) 9361.
[12] G. Careri, Prog. Biophys. Mol. Biol. 70 (1998) 223.
[13] G. Careri and E. Milotti, Phys. Rev. E 67 (2003) 051923.
[14] The distance ∆i,k in equation (3) in the text is defined exactly only in the continuous limit,
however this is not important because the discrete formulation in section I is used only to
introduce the diffusion equation and is not relevant for the algorithm.
[15] S. M. Kay and S. L. Marple, Proc. IEEE 69 (1981) 1380.
[16] M. B. Weissman, Rev. Mod. Phys. 60 (1988) 537.
[17] N. G. van Kampen, J. Phys. Chem. Solids 49 (1988) 673.
[18] M. Christensen and J. Boiden Pedersen, J. Chem. Phys. 119 (2003) 5171.
[19] S. Sellers and J. A. Barker, Phys. Rev. E 74 (2006) 061103.
[20] M. Christensen, J. Comp. Phys. 201 (2004) 421.
[21] K. Kikuchi, M. Yoshida, T. Maekawa, and H. Watanabe, Chem. Phys. Lett. 185 (1991) 335.
[22] K. Kikuchi, M. Yoshida, T. Maekawa, and H. Watanabe, Chem. Phys. Lett. 196 (1992) 57.
[23] G. Careri and E. Milotti, in preparation.
[24] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery: “Numerical recipes in C:
the art of scientific computing, 2nd edition”, section 7.2 (Cambridge Univ. Press, Cambridge,
1992).
FIG. 1: Here space has been subdivided in elementary regions (the dots represent the centroids of
the elementary regions): each region is characterized by a particular transition rate γi and by an
elementary volume Vi. In this figure the transitions between the central region i and the adjacent
regions (denoted here by an index kj) are marked by the gray arrows.
FIG. 2: Random walks around a single center as explained in section II. The random walker starts
at the origin, with a fixed step length = 0.005 arbitrary units; the horizontal and vertical scales are
also labeled with the same arbitrary length units. The width parameter is Γ = 1 in all simulations,
and the other parameters have fixed values as well, A = 1, B = 0.1, while a. α = 1; b. α =
10; c.
α = 10. Larger values of α correspond to higher anisotropy, and we see here that as the anisotropy
grows, the random walk becomes more and more compact.
FIG. 3: Random walk with two centers, as explained in the text. The walk is superposed on the
contour plot of the transition rate γ(r). The random walker starts at the origin, with a fixed step
length = 0.01 arbitrary units; the horizontal and vertical scales are also labeled with the same
arbitrary length units.
FIG. 4: Random walk with several (10) centers, a. after 10000 steps; b. after 40000 steps; c. after
70000 steps. The walk is superposed on the contour plot of the transition rate γ(r).
FIG. 5: Random walk in space with two centers, located at the positions marked by the gray
arrows.
FIG. 6: Random walk on a spherical surface, with several (10) centers, some of which have not yet
been visited by the random walker.
FIG. 7: Superposition of several random walks about a single center in the origin. Isodensity
contour lines are superposed on the density plot.
FIG. 8: This figure shows the x projection of the position signal (this is not very different from
an electric dipole component if the random walkers are charged particles) vs. time for one of the
random walks in the example of section VI; both position and time are in arbitrary units. The
insets show parts of the signal with increasing magnification, and the last inset displays clearly the
stationary parts of the signal between successive transitions.
FIG. 9: Unnormalized distribution of the time intervals ∆t between transitions in the set of 400
random walks of 10000 steps each described in the text: a. logarithm of the relative frequency
vs. ∆t (both in arbitrary units), which shows that the distribution is not a simple exponential,
but rather contains two different power-law regions (dotted lines); b. log-log plot of the same
distribution, where the two power laws are identified by the nearly straight sections.
FIG. 10: Mean squared distance 〈r2〉 vs. time. In an ordinary random walk the mean squared
distance is a linear function of time: here we see that linearity is recovered only asymptotically, as
the random walkers explore regions that are further away from the origin, where the γ is nearly
constant.
FIG. 11: Spectral density calculated with the method of section V, sampled at logarithmically
spaced frequencies: a. spectrum obtained from the x position signal for a single random walk;
b. average of 400 spectra. Spectral densities and frequencies are in arbitrary units. The thin
gray line represents an ideal 1/f2 spectrum, which is expected for this kind of processes: the
computed spectrum deviates from the ideal spectrum only at very low frequency, because of the
limited observation time. Notice also that there is no upward bend at very high frequency – a hint
of the absence of aliasing. In part a. it is clearly visible that the spectrum is sampled at (500)
logarithmically spaced frequencies, because there is no crowding at the high end of the spectrum,
and no rarefaction at low frequency, unlike spectra obtained with the Fast Fourier Transform or
other similar algorithms.
	Introduction
	Angular distribution
	Diffusion on surfaces
	Asynchronous transitions
	Fourier analysis of asynchronous data
	Random walk about a single center
	Discussion
	Angular distributions
	Time distribution
	Acknowledgments
	References

  In this paper I describe a specialized algorithm for anisotropic diffusion
determined by a field of transition rates. The algorithm can be used to
describe some interesting forms of diffusion that occur in the study of proton
motion in a network of hydrogen bonds. The algorithm produces data that require
a nonstandard method of spectral analysis which is also developed here.
Finally, I apply the algorithm to a simple specific example.

Introduction
	Angular distribution
	Diffusion on surfaces
	Asynchronous transitions
	Fourier analysis of asynchronous data
	Random walk about a single center
	Discussion
	Angular distributions
	Time distribution
	Acknowledgments
	References
