Source code for MMTK.Random

# Functions for finding random points and orientations.
#
# Written by: Konrad Hinsen
# 

"""
Random quantities for use in molecular simulations
"""

__docformat__ = 'restructuredtext'

from Scientific.Geometry import Vector
from Scientific.Geometry.Transformation import Rotation
from Scientific import N
from MMTK import ParticleProperties, Units

try:
    numeric = N.package
except AttributeError:
    numeric = "Numeric"

RNG = None
try:
    if numeric == "Numeric":
        import RNG
    elif numeric == "NumPy":
        import numpy.oldnumeric.rng as RNG
except ImportError:
    pass


if RNG is None:

    if numeric == "Numeric":
        from RandomArray import uniform, seed
    elif numeric == "NumPy":
        from numpy.oldnumeric.random_array import uniform, seed
    random = __import__('random')
    seed(1, 1)
    random.seed(1)

    def initializeRandomNumbersFromTime():
        random.seed()
        seed(0, 0)

    def gaussian(mean, std, shape=None):
        if shape is None:
            x = random.normalvariate(0., 1.)
        else:
            x = N.zeros(shape, N.Float)
            xflat = N.ravel(x)
            for i in range(len(xflat)):
                xflat[i] = random.normalvariate(0., 1.)
        return mean + std*x

else:

    _uniform_generator = \
                    RNG.CreateGenerator(-1, RNG.UniformDistribution(0., 1.))
    _gaussian_generator = \
                    RNG.CreateGenerator(-1, RNG.NormalDistribution(0., 1.))

    def initializeRandomNumbersFromTime():
        global _uniform_generator, _gaussian_generator
        _uniform_generator = \
                       RNG.CreateGenerator(0, RNG.UniformDistribution(0., 1.))
        _gaussian_generator = \
                       RNG.CreateGenerator(0, RNG.NormalDistribution(0., 1.))

    def uniform(x1, x2, shape=None):
        if shape is None:
            x = _uniform_generator.ranf()
        else:
            n = N.multiply.reduce(shape)
            x = _uniform_generator.sample(n)
            x.shape = shape
        return x1+(x2-x1)*x

    def gaussian(mean, std, shape=None):
        if shape is None:
            x = _gaussian_generator.ranf()
        else:
            n = N.multiply.reduce(shape)
            x = _gaussian_generator.sample(n)
            x.shape = shape
        return mean+std*x

del numeric

#
# Random point in a rectangular box centered around the origin
#
[docs]def randomPointInBox(a, b = None, c = None): """ :param a: the edge length of a box along the x axis :type a: float :param b: the edge length of a box along the y axis (default: a) :type b: float :param c: the edge length of a box along the z axis (default: a) :type c: float :returns: a vector drawn from a uniform distribution within a rectangular box with edge lengths a, b, c. :rtype: Scientific.Geometry.Vector """ if b is None: b = a if c is None: c = a x = uniform(-0.5*a, 0.5*a) y = uniform(-0.5*b, 0.5*b) z = uniform(-0.5*c, 0.5*c) return Vector(x, y, z) # # Random point in a sphere around the origin. #
[docs]def randomPointInSphere(r): """ :param r: the radius of a sphere :type r: float :returns: a vector drawn from a uniform distribution within a sphere of radius r. :rtype: Scientific.Geometry.Vector """ rsq = r*r while 1: x = N.array([uniform(-r, r), uniform(-r, r), uniform(-r, r)]) if N.dot(x, x) < rsq: break return Vector(x) # # Random direction (unit vector). #
[docs]def randomDirection(): """ :returns: a vector drawn from a uniform distribution on the surface of a unit sphere. :rtype: Scientific.Geometry.Vector """ r = randomPointInSphere(1.) return r.normal()
[docs]def randomDirections(n): """ :param n: the number of direction vectors :returns: a list of n vectors drawn from a uniform distribution on the surface of a unit sphere. If n is negative, returns a deterministic list of not more than -n vectors of unit length (useful for testing purposes). :rtype: list """ if n < 0: vs = [Vector(1., 0., 0.), Vector(0., -1., 0.), Vector(0., 0., 1.), Vector(-1., 1., 0.).normal(), Vector(-1., 0., 1.).normal(), Vector(0., 1., -1.).normal(), Vector(1., -1., 1.).normal()] return vs[:-n] else: return [randomDirection() for i in range(n)] # # Random rotation. #
[docs]def randomRotation(max_angle = N.pi): """ :param max_angle: the upper limit for the rotation angle :type max_angle: float :returns: a random rotation with a uniform axis distribution and angles drawn from a uniform distribution between -max_angle and max_angle. :rtype: Scientific.Geometry.Transformations.Rotation """ return Rotation(randomDirection(), uniform(-max_angle, max_angle)) # # Random velocity (gaussian) #
[docs]def randomVelocity(temperature, mass): """ :param temperature: the temperature defining a Maxwell distribution :type temperature: float :param mass: the mass of a particle :type mass: float :returns: a random velocity vector for a particle of a given mass at a given temperature :rtype: Scientific.Geometry.Vector """ sigma = N.sqrt((temperature*Units.k_B)/(mass*Units.amu)) return Vector(gaussian(0., sigma, (3,))) # # Random ParticleVector (gaussian) #
[docs]def randomParticleVector(universe, width): """ :param universe: a universe :type universe: :class:`~MMTK.Universe.Universe` :param width: the width (standard deviation) of a Gaussian distribution :type width: float :returns: a set of vectors drawn from a Gaussian distribution with a given width centered around zero. :rtype: :class:`~MMTK.ParticleProperties.ParticleVector` """ data = gaussian(0., 0.577350269189*width, (universe.numberOfPoints(), 3)) return ParticleProperties.ParticleVector(universe, data)