Skip to content

GitLab

Commit 3e974dbd authored by eric pellegrini's avatar eric pellegrini
Browse files

removed oldnumeric deps from Random,Trajectory and Universe

parent 315fb52b
Hide whitespace changes
Inline Side-by-side
MMTK/ChemicalObjects.py
View file @ 3e974dbd
......@@ -10,12 +10,14 @@ Atoms, groups, molecules and similar objects
__docformat__ = 'restructuredtext'
import copy
import numpy
from MMTK import Bonds, Collections, Database, \
Units, Utility, Visualization
from Scientific.Geometry import Vector
from Scientific.Geometry import Objects3D
from Scientific import N
import copy
#
# The base class for all chemical objects.
......@@ -589,8 +591,8 @@ class Atom(ChemicalObject):
except AttributeError:
return None
else:
if N.logical_or.reduce(
N.greater(self.array[self.index, :],
if numpy.logical_or.reduce(
numpy.greater(self.array[self.index, :],
Utility.undefined_limit)):
return None
else:
......@@ -875,7 +877,7 @@ class Molecule(CompositeChemicalObject, ChemicalObject):
# default C-H bond length and X-C-H angle
_ch_bond = 1.09*Units.Ang
_hch_angle = N.arccos(-1./3.)*Units.rad
_hch_angle = numpy.arccos(-1./3.)*Units.rad
_nh_bond = 1.03*Units.Ang
_hnh_angle = 120.*Units.deg
_oh_bond = 0.95*Units.Ang
......@@ -1025,7 +1027,7 @@ class Molecule(CompositeChemicalObject, ChemicalObject):
def _tetrahedralH(self, atom, known, unknown, bond):
r = atom.position()
n = (known[0].position()-r).normal()
cone = Objects3D.Cone(r, n, N.arccos(-1./3.))
cone = Objects3D.Cone(r, n, numpy.arccos(-1./3.))
sphere = Objects3D.Sphere(r, bond)
circle = sphere.intersectWith(cone)
others = filter(lambda a: a.symbol != 'H', known[0].bondedTo())
......
MMTK/Random.py
View file @ 3e974dbd
......@@ -9,83 +9,13 @@ Random quantities for use in molecular simulations
__docformat__ = 'restructuredtext'
import numpy
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
#TODO(Remi): Test all the functions with the one using oldnumeric by providing a give seed
#
# Random point in a rectangular box centered around the origin
......@@ -104,9 +34,9 @@ def randomPointInBox(a, b = None, c = None):
"""
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)
x = numpy.random.uniform(-0.5*a, 0.5*a)
y = numpy.random.uniform(-0.5*b, 0.5*b)
z = numpy.random.uniform(-0.5*c, 0.5*c)
return Vector(x, y, z)
#
......@@ -122,8 +52,8 @@ def randomPointInSphere(r):
"""
rsq = r*r
while 1:
x = N.array([uniform(-r, r), uniform(-r, r), uniform(-r, r)])
if N.dot(x, x) < rsq: break
x = numpy.random.uniform(-r, r,3)
if numpy.dot(x, x) < rsq: break
return Vector(x)
#
......@@ -153,12 +83,12 @@ def randomDirections(n):
Vector(0., 1., -1.).normal(), Vector(1., -1., 1.).normal()]
return vs[:-n]
else:
return [randomDirection() for i in range(n)]
return [randomDirection() for _ in range(n)]
#
# Random rotation.
#
def randomRotation(max_angle = N.pi):
def randomRotation(max_angle = numpy.pi):
"""
:param max_angle: the upper limit for the rotation angle
:type max_angle: float
......@@ -167,7 +97,7 @@ def randomRotation(max_angle = N.pi):
-max_angle and max_angle.
:rtype: Scientific.Geometry.Transformations.Rotation
"""
return Rotation(randomDirection(), uniform(-max_angle, max_angle))
return Rotation(randomDirection(), numpy.random.uniform(-max_angle, max_angle))
#
# Random velocity (gaussian)
......@@ -182,8 +112,8 @@ def randomVelocity(temperature, 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,)))
sigma = numpy.sqrt((temperature*Units.k_B)/(mass*Units.amu))
return Vector(numpy.random.normal(0., sigma, (3,)))
#
# Random ParticleVector (gaussian)
......@@ -198,5 +128,5 @@ def randomParticleVector(universe, width):
with a given width centered around zero.
:rtype: :class:`~MMTK.ParticleProperties.ParticleVector`
"""
data = gaussian(0., 0.577350269189*width, (universe.numberOfPoints(), 3))
data = numpy.random.normal(0., 0.577350269189*width, (universe.numberOfPoints(), 3))
return ParticleProperties.ParticleVector(universe, data)
MMTK/Trajectory.py
View file @ 3e974dbd
......@@ -9,10 +9,11 @@ Trajectory files and their contents
__docformat__ = 'restructuredtext'
import numpy
from MMTK import Collections, Units, Universe, Utility, \
ParticleProperties, Visualization
from Scientific.Geometry import Vector
from Scientific import N
import copy, os, sys
# Report error if the netCDF module is not available.
......@@ -152,7 +153,7 @@ class Trajectory(object):
else:
if mode == 'r':
raise ValueError("can't read trajectory for a non-universe")
index_map = N.array([a.index for a in object.atomList()])
index_map = numpy.array([a.index for a in object.atomList()])
inverse_map = universe.numberOfPoints()*[None]
for i in range(len(index_map)):
inverse_map[index_map[i]] = i
......@@ -211,7 +212,7 @@ class Trajectory(object):
def __getitem__(self, item):
if not isinstance(item, int):
return SubTrajectory(self, N.arange(len(self)))[item]
return SubTrajectory(self, numpy.arange(len(self)))[item]
if item < 0:
item += len(self)
if item >= len(self):
......@@ -240,7 +241,7 @@ class Trajectory(object):
if data.has_key('configuration'):
box = data.get('box_size', None)
if box is not None:
box = box.astype(N.Float)
box = box.astype(numpy.float)
conf = data['configuration']
data['configuration'] = \
ParticleProperties.Configuration(conf.universe, conf.array, box)
......@@ -258,13 +259,13 @@ class Trajectory(object):
return TrajectoryVariable(self.universe, self, name)
elif 'box_size_length' in var.dimensions:
if 'minor_step_number' in var.dimensions:
bs = N.transpose(var[:], [0, 2, 1])
bs = N.reshape(bs, (bs.shape[0]*bs.shape[1], bs.shape[2]))
bs = numpy.transpose(var[:], [0, 2, 1])
bs = numpy.reshape(bs, (bs.shape[0]*bs.shape[1], bs.shape[2]))
return bs[:len(self)]
else:
return var[:]
else:
return N.ravel(N.array(var))[:len(self)]
return numpy.ravel(numpy.array(var))[:len(self)]
def defaultStep(self):
try:
......@@ -483,7 +484,7 @@ class TrajectoryVariable(object):
def __getitem__(self, item):
if not isinstance(item, int):
return SubVariable(self, N.arange(len(self)))[item]
return SubVariable(self, numpy.arange(len(self)))[item]
item = int(item) # gets rid of numpy.intXX objects
if item < 0:
item = item + len(self.trajectory)
......@@ -494,9 +495,9 @@ class TrajectoryVariable(object):
box = None
elif len(self.box_size.shape) == 3:
bs = self.trajectory.block_size
box = self.box_size[item/bs, :, item%bs].astype(N.Float)
box = self.box_size[item/bs, :, item%bs].astype(numpy.float)
else:
box = self.box_size[item].astype(N.Float)
box = self.box_size[item].astype(numpy.float)
array = ParticleProperties.Configuration(self.universe,
self.trajectory.trajectory.readParticleVector(self.name, item),
box)
......@@ -615,10 +616,10 @@ class TrajectorySet(object):
def __getitem__(self, item):
if not isinstance(item, int):
return SubTrajectory(self, N.arange(len(self)))[item]
return SubTrajectory(self, numpy.arange(len(self)))[item]
if item >= len(self):
raise IndexError
tindex = N.add.reduce(N.greater_equal(item, self.nsteps))-1
tindex = numpy.add.reduce(numpy.greater_equal(item, self.nsteps))-1
return self.trajectories[tindex][item-self.nsteps[tindex]]
def __getslice__(self, first, last):
......@@ -634,8 +635,8 @@ class TrajectorySet(object):
data = []
for t in self.trajectories:
var = t.trajectory.file.variables[name]
data.append(N.ravel(N.array(var))[:len(t)])
return N.concatenate(data)
data.append(numpy.ravel(numpy.array(var))[:len(t)])
return numpy.concatenate(data)
def readParticleTrajectory(self, atom, first=0, last=None, skip=1,
variable = "configuration"):
......@@ -666,35 +667,35 @@ class TrajectorySet(object):
and self.cell_parameters[0] is not None:
jump = pt.array[0]-total.array[-1]
mult = -(jump/self.cell_parameters[i-1]).astype('i')
if len(N.nonzero(mult)) > 0:
if len(numpy.nonzero(mult)) > 0:
t._boxTransformation(pt.array, pt.array, 1)
N.add(pt.array, mult[N.NewAxis, : ],
numpy.add(pt.array, mult[numpy.newaxis, : ],
pt.array)
t._boxTransformation(pt.array, pt.array, 0)
jump = pt.array[0] - total.array[-1]
mask = N.less(jump,
mask = numpy.less(jump,
-0.5*self.cell_parameters[i-1])- \
N.greater(jump,
numpy.greater(jump,
0.5*self.cell_parameters[i-1])
if len(N.nonzero(mask)) > 0:
if len(numpy.nonzero(mask)) > 0:
t._boxTransformation(pt.array, pt.array, 1)
N.add(pt.array, mask[N.NewAxis, :],
numpy.add(pt.array, mask[numpy.newaxis, :],
pt.array)
t._boxTransformation(pt.array, pt.array, 0)
elif variable == "box_coordinates" \
and self.cell_parameters[0] is not None:
jump = pt.array[0]-total.array[-1]
mult = -jump.astype('i')
if len(N.nonzero(mult)) > 0:
N.add(pt.array, mult[N.NewAxis, : ],
if len(numpy.nonzero(mult)) > 0:
numpy.add(pt.array, mult[numpy.newaxis, : ],
pt.array)
jump = pt.array[0] - total.array[-1]
mask = N.less(jump, -0.5)- \
N.greater(jump, 0.5)
if len(N.nonzero(mask)) > 0:
N.add(pt.array, mask[N.NewAxis, :],
mask = numpy.less(jump, -0.5)- \
numpy.greater(jump, 0.5)
if len(numpy.nonzero(mask)) > 0:
numpy.add(pt.array, mask[numpy.newaxis, :],
pt.array)
total.array = N.concatenate((total.array, pt.array))
total.array = numpy.concatenate((total.array, pt.array))
else:
self.steps_read.append(0)
return total
......@@ -741,10 +742,10 @@ class TrajectorySetVariable(TrajectoryVariable):
def __getitem__(self, item):
if not isinstance(item, int):
return SubVariable(self, N.arange(len(self)))[item]
return SubVariable(self, numpy.arange(len(self)))[item]
if item >= len(self.trajectory_set):
raise IndexError
tindex = N.add.reduce(N.greater_equal(item,
tindex = numpy.add.reduce(numpy.greater_equal(item,
self.trajectory_set.nsteps))-1
step = item-self.trajectory_set.nsteps[tindex]
t = self.trajectory_set.trajectories[tindex]
......@@ -840,7 +841,7 @@ class ParticleTrajectory(object):
:param vector: the vector to be added
:type vector: Scientific.Geometry.Vector
"""
N.add(self.array, vector.array[N.NewAxis, :], self.array)
numpy.add(self.array, vector.array[numpy.newaxis, :], self.array)
#
# Rigid-body trajectory
......@@ -876,21 +877,21 @@ class RigidBodyTrajectory(object):
ref_cms = object.centerOfMass(reference)
atoms = object.atomList()
possq = N.zeros((steps,), N.Float)
cross = N.zeros((steps, 3, 3), N.Float)
rcms = N.zeros((steps, 3), N.Float)
possq = numpy.zeros((steps,), numpy.float)
cross = numpy.zeros((steps, 3, 3), numpy.float)
rcms = numpy.zeros((steps, 3), numpy.float)
# cms of the CONTIGUOUS object made of CONTINUOUS atom trajectories
for a in atoms:
r = trajectory.readParticleTrajectory(a, first, last, skip,
"box_coordinates").array
w = a._mass/mass
N.add(rcms, w*r, rcms)
numpy.add(rcms, w*r, rcms)
if offset is not None:
N.add(rcms, w*offset[a].array, rcms)
numpy.add(rcms, w*offset[a].array, rcms)
# relative coords of the CONTIGUOUS reference
r_ref = N.zeros((len(atoms), 3), N.Float)
r_ref = numpy.zeros((len(atoms), 3), numpy.float)
for a in range(len(atoms)):
r_ref[a] = atoms[a].position(reference).array - ref_cms.array
......@@ -901,18 +902,18 @@ class RigidBodyTrajectory(object):
"box_coordinates").array
r = r - rcms # (a-b)**2 != a**2 - b**2
if offset is not None:
N.add(r, offset[atoms[a]].array,r)
numpy.add(r, offset[atoms[a]].array,r)
trajectory._boxTransformation(r, r)
w = atoms[a]._mass/mass
N.add(possq, w*N.add.reduce(r*r, -1), possq)
N.add(possq, w*N.add.reduce(r_ref[a]*r_ref[a],-1),
numpy.add(possq, w*numpy.add.reduce(r*r, -1), possq)
numpy.add(possq, w*numpy.add.reduce(r_ref[a]*r_ref[a],-1),
possq)
N.add(cross, w*r[:,:,N.NewAxis]*r_ref[N.NewAxis,
numpy.add(cross, w*r[:,:,numpy.newaxis]*r_ref[numpy.newaxis,
a,:],cross)
self.trajectory._boxTransformation(rcms, rcms)
# filling matrix M (formula no 40)
k = N.zeros((steps, 4, 4), N.Float)
k = numpy.zeros((steps, 4, 4), numpy.float)
k[:, 0, 0] = -cross[:, 0, 0]-cross[:, 1, 1]-cross[:, 2, 2]
k[:, 0, 1] = cross[:, 1, 2]-cross[:, 2, 1]
k[:, 0, 2] = cross[:, 2, 0]-cross[:, 0, 2]
......@@ -927,21 +928,21 @@ class RigidBodyTrajectory(object):
for i in range(1, 4):
for j in range(i):
k[:, i, j] = k[:, j, i]
N.multiply(k, 2., k)
numpy.multiply(k, 2., k)
for i in range(4):
N.add(k[:,i,i], possq, k[:,i,i])
numpy.add(k[:,i,i], possq, k[:,i,i])
del possq
quaternions = N.zeros((steps, 4), N.Float)
fit = N.zeros((steps,), N.Float)
quaternions = numpy.zeros((steps, 4), numpy.float)
fit = numpy.zeros((steps,), numpy.float)
from Scientific.LA import eigenvectors
for i in range(steps):
e, v = eigenvectors(k[i])
j = N.argmin(e)
j = numpy.argmin(e)
if e[j] < 0.:
fit[i] = 0.
else:
fit[i] = N.sqrt(e[j])
fit[i] = numpy.sqrt(e[j])
if v[j,0] < 0.: quaternions[i] = -v[j] # eliminate jumps
else: quaternions[i] = v[j]
self.fit = fit
......
MMTK/Units.py
View file @ 3e974dbd
......@@ -56,7 +56,7 @@ The constants defined in this module are:
"""
from Scientific import N as Numeric
import numpy
# Prefixes
......@@ -80,7 +80,7 @@ peta = 1.0e15
# internal unit: radian
rad = 1.
deg = (Numeric.pi/180.)*rad
deg = (numpy.pi/180.)*rad
# Length units
# internal unit: nm
......@@ -173,20 +173,19 @@ atm = 101325*Pa # atmosphere
# Constants
h = 6.626176e-34*J*s
hbar = h/(2.*Numeric.pi)
hbar = h/(2.*numpy.pi)
k_B = 1.3806513e-23*J/K
eps0 = 1./(4.e-7*Numeric.pi)*A**2*m/J/c**2
eps0 = 1./(4.e-7*numpy.pi)*A**2*m/J/c**2
me = 0.51099906*mega*eV/c**2
electrostatic_energy = 1/(4.*Numeric.pi*eps0)
electrostatic_energy = 1/(4.*numpy.pi*eps0)
# CHARMM time unit
akma_time = Numeric.sqrt(Ang**2/(kcal/mol))
akma_time = numpy.sqrt(Ang**2/(kcal/mol))
# "Atomic" units
Bohr = 4.*Numeric.pi*eps0*hbar**2/(me*e**2)
Bohr = 4.*numpy.pi*eps0*hbar**2/(me*e**2)
Hartree = hbar**2/(me*Bohr**2)
# Remove symbol "Numeric"
del Numeric
del numpy
MMTK/Universe.py
View file @ 3e974dbd
......@@ -12,11 +12,12 @@ Universes
__docformat__ = 'restructuredtext'
import numpy
from MMTK import Bonds, ChemicalObjects, Collections, Environment, \
Random, Utility, ParticleProperties, Visualization
from Scientific.Geometry import Transformation
from Scientific.Geometry import Vector, isVector
from Scientific import N
import copy
try:
......@@ -416,7 +417,7 @@ class Universe(Collections.GroupOfAtoms, Visualization.Viewable):
"""
if self._configuration is None:
np = self.numberOfAtoms()
coordinates = N.zeros((np, 3), N.Float)
coordinates = numpy.zeros((np, 3), numpy.float)
index_map = {}
redef = []
for a in self.atomList():
......@@ -518,7 +519,7 @@ class Universe(Collections.GroupOfAtoms, Visualization.Viewable):
if block:
self.releaseWriteStateLock()
def getParticleScalar(self, name, datatype = N.Float):
def getParticleScalar(self, name, datatype = numpy.float):
"""
:param name: the name of an atom attribute
:type name: str
......@@ -528,7 +529,7 @@ class Universe(Collections.GroupOfAtoms, Visualization.Viewable):
:rtype: :class:`~MMTK.ParticleProperties.ParticleScalar`
"""
conf = self.configuration()
array = N.zeros((len(conf),), datatype)
array = numpy.zeros((len(conf),), datatype)
for a in self.atomList():
array[a.index] = getattr(a, name)
return ParticleProperties.ParticleScalar(self, array)
......@@ -544,7 +545,7 @@ class Universe(Collections.GroupOfAtoms, Visualization.Viewable):
:rtype: :class:`~MMTK.ParticleProperties.ParticleScalar`
"""
conf = self.configuration()
array = N.zeros((len(conf),), N.Int)
array = numpy.zeros((len(conf),), numpy.int)
for a in self.atomList():
try:
array[a.index] = getattr(a, name)
......@@ -636,7 +637,7 @@ class Universe(Collections.GroupOfAtoms, Visualization.Viewable):
del self._atom_properties['velocity']
fixed = self.getParticleBoolean('fixed')
np = self.numberOfPoints()
velocities = N.zeros((np, 3), N.Float)
velocities = numpy.zeros((np, 3), numpy.float)
for i in xrange(np):
m = masses[i]
if m > 0. and not fixed[i]:
......@@ -660,7 +661,7 @@ class Universe(Collections.GroupOfAtoms, Visualization.Viewable):
:type block: bool
"""
velocities = self.velocities()
factor = N.sqrt(temperature/self.temperature())
factor = numpy.sqrt(temperature/self.temperature())
if block:
self.acquireWriteStateLock()
try:
......@@ -670,7 +671,7 @@ class Universe(Collections.GroupOfAtoms, Visualization.Viewable):
self.releaseWriteStateLock()
def degreesOfFreedom(self):
return GroupOfAtoms.degreesOfFreedom(self) \
return Collections.GroupOfAtoms.degreesOfFreedom(self) \
- self.numberOfDistanceConstraints()
def distanceConstraintList(self):
......@@ -1171,9 +1172,9 @@ class Universe(Collections.GroupOfAtoms, Visualization.Viewable):
def numberOfFixedAtoms(self):
return self.getParticleBoolean('fixed').sumOverParticles()
def degreesOfFreedom(self):
return 3*(self.numberOfAtoms()-self.numberOfFixedAtoms()) \
- self.numberOfDistanceConstraints()
# def degreesOfFreedom(self):
# return 3*(self.numberOfAtoms()-self.numberOfFixedAtoms()) \
# - self.numberOfDistanceConstraints()
def mass(self):
return self.masses().sumOverParticles()
......@@ -1196,13 +1197,13 @@ class Universe(Collections.GroupOfAtoms, Visualization.Viewable):
def translateBy(self, vector):
conf = self.configuration().array
N.add(conf, vector.array[N.NewAxis, :], conf)
numpy.add(conf, vector.array[numpy.newaxis, :], conf)
def applyTransformation(self, t):
conf = self.configuration().array
rot = t.rotation().tensor.array
conf[:] = N.dot(conf, N.transpose(rot))
N.add(conf, t.translation().vector.array[N.NewAxis, :], conf)
conf[:] = numpy.dot(conf, numpy.transpose(rot))
numpy.add(conf, t.translation().vector.array[numpy.newaxis, :], conf)
def writeXML(self, file):
file.write('<?xml version="1.0" encoding="ISO-8859-1" ' +
......@@ -1319,17 +1320,17 @@ class Periodic3DUniverse(Universe):
def cartesianToFractional(self, vector):
r1, r2, r3 = self.reciprocalBasisVectors()
return N.array([r1*vector, r2*vector, r3*vector])
return numpy.array([r1*vector, r2*vector, r3*vector])
def cartesianToFractionalMatrix(self):
return N.array(self.reciprocalBasisVectors())
return numpy.array(self.reciprocalBasisVectors())