removed oldnumeric deps from Deformation,Dynamics,Environment,Field,...

removed oldnumeric deps from Deformation,Dynamics,Environment,Field, FourierBasis,Geometry,InternalCoordinates

MMTK/Deformation.py

@@ -21,6 +21,8 @@ of the techniques can be found in the following articles:
 
 __docformat__ = 'restructuredtext'
 
+import numpy
+
 try:
     from MMTK_forcefield import NonbondedList
     from MMTK_deformation import deformation, reduceDeformation, \
@@ -28,7 +30,6 @@ try:
 except ImportError:
     pass
 from MMTK import ParticleProperties
-from Scientific import N
 
 #
 # Deformation energy evaluations
@@ -42,7 +43,7 @@ class DeformationEvaluationFunction(object):
         self.cutoff = cutoff
         self.factor = factor
 
-        nothing = N.zeros((0,2), N.Int)
+        nothing = numpy.zeros((0,2), numpy.int)
         self.pairs = NonbondedList(nothing, nothing, nothing,
                                    universe._spec, cutoff)
         self.pairs.update(universe.configuration().array)
@@ -227,7 +228,7 @@ class DeformationEnergyFunction(DeformationEvaluationFunction):
         if gradients is not None:
             if ParticleProperties.isParticleProperty(gradients):
                 g = gradients
-            elif isinstance(gradients, N.array_type):
+            elif isinstance(gradients, numpy.ndarray):
                 g = ParticleProperties.ParticleVector(self.universe, gradients)
             elif gradients:
                 g = ParticleProperties.ParticleVector(self.universe)
@@ -331,7 +332,7 @@ class FiniteDeformationEnergyFunction(DeformationEvaluationFunction):
         if gradients is not None:
             if ParticleProperties.isParticleProperty(gradients):
                 g = gradients
-            elif isinstance(gradients, N.array_type):
+            elif isinstance(gradients, numpy.ndarray):
                 g = ParticleProperties.ParticleVector(self.universe, gradients)
             elif gradients:
                 g = ParticleProperties.ParticleVector(self.universe)

MMTK/Dynamics.py

@@ -9,9 +9,10 @@ Molecular Dynamics integrators
 
 __docformat__ = 'restructuredtext'
 
+import numpy
+
 from MMTK import Environment, Features, Trajectory, Units
 import MMTK_dynamics
-from Scientific import N
 
 #
 # Integrator base class
@@ -130,8 +131,8 @@ class VelocityVerletIntegrator(Integrator):
             t_parameters = thermostat.parameters
             t_coordinates = thermostat.coordinates
         else:
-            t_parameters = N.zeros((0,), N.Float)
-            t_coordinates = N.zeros((2,), N.Float)
+            t_parameters = numpy.zeros((0,), numpy.float)
+            t_coordinates = numpy.zeros((2,), numpy.float)
         if Features.AndersenBarostatFeature in used_features:
             if self.universe.cellVolume() is None:
                 raise ValueError("Barostat requires finite volume universe")
@@ -144,8 +145,8 @@ class VelocityVerletIntegrator(Integrator):
             b_parameters = barostat.parameters
             b_coordinates = barostat.coordinates
         else:
-            b_parameters = N.zeros((0,), N.Float)
-            b_coordinates = N.zeros((1,), N.Float)
+            b_parameters = numpy.zeros((0,), numpy.float)
+            b_coordinates = numpy.zeros((1,), numpy.float)
 
         args = (self.universe,
                 configuration.array, velocities.array,
@@ -194,7 +195,7 @@ class VelocityScaler(Trajectory.TrajectoryAction):
         :type skip: int
         """
         Trajectory.TrajectoryAction.__init__(self, first, last, skip)
-        self.parameters = N.array([temperature, window], N.Float)
+        self.parameters = numpy.array([temperature, window], numpy.float)
         self.Cfunction = MMTK_dynamics.scaleVelocities
 
 class Heater(Trajectory.TrajectoryAction):
@@ -228,7 +229,7 @@ class Heater(Trajectory.TrajectoryAction):
         :type skip: int
         """
         Trajectory.TrajectoryAction.__init__(self, first, last, skip)
-        self.parameters = N.array([temp1, temp2, gradient], N.Float)
+        self.parameters = numpy.array([temp1, temp2, gradient], numpy.float)
         self.Cfunction = MMTK_dynamics.heat
 
 class BarostatReset(Trajectory.TrajectoryAction):
@@ -255,7 +256,7 @@ class BarostatReset(Trajectory.TrajectoryAction):
         :type skip: int
         """
         Trajectory.TrajectoryAction.__init__(self, first, last, skip)
-        self.parameters = N.zeros((0,), N.Float)
+        self.parameters = numpy.zeros((0,), numpy.float)
         self.Cfunction = MMTK_dynamics.resetBarostat
 
 class TranslationRemover(Trajectory.TrajectoryAction):
@@ -317,8 +318,8 @@ class RotationRemover(Trajectory.TrajectoryAction):
 #
 def _constraintArrays(universe):
     nc = universe.numberOfDistanceConstraints()
-    constraints = N.zeros((nc, 2), N.Int)
-    const_distances_sq = N.zeros((nc, ), N.Float)
+    constraints = numpy.zeros((nc, 2), numpy.int)
+    const_distances_sq = numpy.zeros((nc, ), numpy.float)
     if nc > 0:
         i = 0
         c_blocks = [0]
@@ -329,9 +330,9 @@ def _constraintArrays(universe):
                 const_distances_sq[i] = c[2]*c[2]
                 i = i + 1
             c_blocks.append(i)
-        c_blocks = N.array(c_blocks)
+        c_blocks = numpy.array(c_blocks)
     else:
-        c_blocks = N.zeros((1,), N.Int)
+        c_blocks = numpy.zeros((1,), numpy.int)
     return constraints, const_distances_sq, c_blocks
 
 #

MMTK/Environment.py

@@ -13,7 +13,7 @@ fields, etc.
 
 __docformat__ = 'restructuredtext'
 
-from Scientific import N
+import numpy
 
 #
 # The environment object base class
@@ -56,8 +56,8 @@ class NoseThermostat(EnvironmentObject):
         :type relaxation_time: float
         """
         self.arguments = (temperature, relaxation_time)
-        self.parameters = N.array([temperature, relaxation_time])
-        self.coordinates = N.array([0., 0.])
+        self.parameters = numpy.array([temperature, relaxation_time])
+        self.coordinates = numpy.array([0., 0.])
 
     def setTemperature(self, temperature):
         self.parameters[0] = temperature
@@ -91,8 +91,8 @@ class AndersenBarostat(EnvironmentObject):
         :type relaxation_time: float
         """
         self.arguments = (pressure, relaxation_time)
-        self.parameters = N.array([pressure, relaxation_time])
-        self.coordinates = N.array([0.])
+        self.parameters = numpy.array([pressure, relaxation_time])
+        self.coordinates = numpy.array([0.])
 
     def setPressure(self, pressure):
         self.parameters[0] = pressure

MMTK/Field.py

@@ -22,10 +22,11 @@ without overloading it.
 
 __docformat__ = 'restructuredtext'
 
+import numpy
+
 from MMTK import Collections, ParticleProperties, Visualization
 from Scientific.Visualization import Color
 from Scientific.Geometry import Vector, TensorAnalysis
-from Scientific import N
 
 
 class AtomicField(object):
@@ -60,7 +61,7 @@ class AtomicField(object):
             color = default_color
             for a in atoms:
                 d = a.position()-center
-                weight = N.exp(-d*d/self.box.partition_size**2)
+                weight = numpy.exp(-d*d/self.box.partition_size**2)
                 v = v + weight*data[a]
                 total_weight = total_weight + weight
                 c = default_color
@@ -71,16 +72,16 @@ class AtomicField(object):
                 color = color + c
             values.append(v/total_weight)
             colors.append(color)
-        min = N.minimum.reduce(points)
-        max = N.maximum.reduce(points)
-        axes = (N.arange(min[0], max[0]+1, self.box.partition_size),
-                N.arange(min[1], max[1]+1, self.box.partition_size),
-                N.arange(min[2], max[2]+1, self.box.partition_size))
-        array = N.zeros(tuple(map(len, axes)) + data.value_rank*(3,), N.Float)
-        inside = N.zeros(tuple(map(len, axes)), N.Float)
+        min = numpy.minimum.reduce(points)
+        max = numpy.maximum.reduce(points)
+        axes = (numpy.arange(min[0], max[0]+1, self.box.partition_size),
+                numpy.arange(min[1], max[1]+1, self.box.partition_size),
+                numpy.arange(min[2], max[2]+1, self.box.partition_size))
+        array = numpy.zeros(tuple(map(len, axes)) + data.value_rank*(3,), numpy.float)
+        inside = numpy.zeros(tuple(map(len, axes)), numpy.float)
         for p, v in zip(points, values):
-            indices = N.floor((p-min)/self.box.partition_size+0.5)
-            indices = tuple(indices.astype(N.Int))
+            indices = numpy.floor((p-min)/self.box.partition_size+0.5)
+            indices = tuple(indices.astype(numpy.int))
             array[indices] = v
             inside[indices] = 1.
         self.field = self.field_class(axes, array, data.zero())

MMTK/FourierBasis.py

@@ -17,9 +17,10 @@ see
   | Proteins 33 (1998): 417-429
 """
 
+import numpy
+
 from MMTK import ParticleProperties
 from Scientific.Geometry import Vector
-from Scientific import N
 
 class FourierBasis(object):
 
@@ -43,7 +44,7 @@ class FourierBasis(object):
         p1, p2 = universe.boundingBox()
         p2 = p2 + Vector(cutoff, cutoff, cutoff)
         l = (p2-p1).array
-        n_max = (0.5*l/cutoff+0.5).astype(N.Int)
+        n_max = (0.5*l/cutoff+0.5).astype(numpy.int)
 
         wave_numbers = [(nx, ny, nz)
                         for nx in range(-n_max[0], n_max[0]+1)
@@ -54,7 +55,7 @@ class FourierBasis(object):
 
         atoms = universe.atomList()
         natoms = len(atoms)
-        basis = N.zeros((3*len(wave_numbers)+3, natoms, 3), N.Float)
+        basis = numpy.zeros((3*len(wave_numbers)+3, natoms, 3), numpy.float)
         cm = universe.centerOfMass()
         i = 0
         for rotation in [Vector(1.,0.,0.), Vector(0.,1.,0.), Vector(0.,0.,1.)]:
@@ -64,7 +65,7 @@ class FourierBasis(object):
             i += i
         conf = universe.configuration().array-p1.array
         for n in wave_numbers:
-            k = 2.*N.pi*N.array(n)/l
+            k = 2.*numpy.pi*numpy.array(n)/l
             w = self._w(conf[:, 0], k[0]) * self._w(conf[:, 1], k[1]) * \
                 self._w(conf[:, 2], k[2])
             basis[i, :, 0] = w
@@ -80,9 +81,9 @@ class FourierBasis(object):
 
     def _w(self, x, k):
         if k < 0:
-            return N.sin(-k*x)
+            return numpy.sin(-k*x)
         else:
-            return N.cos(k*x)
+            return numpy.cos(k*x)
 
     def __len__(self):
         return self.array.shape[0]
@@ -108,7 +109,7 @@ def countBasisVectors(universe, cutoff):
     p1, p2 = universe.boundingBox()
     p2 = p2 + Vector(cutoff, cutoff, cutoff)
     l = (p2-p1).array
-    n_max = (0.5*l/cutoff+0.5).astype(N.Int)
+    n_max = (0.5*l/cutoff+0.5).astype(numpy.int)
     n_wave_numbers = 0
     for nx in range(-n_max[0], n_max[0]+1):
         for ny in range(-n_max[1], n_max[1]+1):

MMTK/Geometry.py

@@ -13,8 +13,9 @@ and lattice objects, which define a regular arrangements of points.
 
 __docformat__ = 'restructuredtext'
 
+import numpy
+
 from Scientific.Geometry import Vector
-from Scientific import N
 
 
 # Error type
@@ -100,8 +101,8 @@ class Box(GeometricalObject3D):
         :param corner2: the diagonally opposite corner
         :type corner2: Scientific.Geometry.Vector
         """
-        c1 = N.minimum(corner1.array, corner2.array)
-        c2 = N.maximum(corner1.array, corner2.array)
+        c1 = numpy.minimum(corner1.array, corner2.array)
+        c2 = numpy.maximum(corner1.array, corner2.array)
         self.corners = c1, c2
 
     def __repr__(self):
@@ -112,18 +113,18 @@ class Box(GeometricalObject3D):
 
     def volume(self):
         c1, c2 = self.corners
-        return N.multiply.reduce(c2-c1)
+        return numpy.multiply.reduce(c2-c1)
 
     def hasPoint(self, point):
         c1, c2 = self.corners
-        min1 = N.minimum.reduce(N.fabs(point.array-c1))
-        min2 = N.minimum.reduce(N.fabs(point.array-c2))
+        min1 = numpy.minimum.reduce(numpy.fabs(point.array-c1))
+        min2 = numpy.minimum.reduce(numpy.fabs(point.array-c2))
         return min1 < eps or min2 < eps
 
     def enclosesPoint(self, point):
         c1, c2 = self.corners
-        out1 = N.logical_or.reduce(N.less(point.array-c1, 0))
-        out2 = N.logical_or.reduce(N.less_equal(c2-point.array, 0))
+        out1 = numpy.logical_or.reduce(numpy.less(point.array-c1, 0))
+        out2 = numpy.logical_or.reduce(numpy.less_equal(c2-point.array, 0))
         return not (out1 or out2)
 
     def cornerPoints(self):
@@ -161,10 +162,10 @@ class Sphere(GeometricalObject3D):
     __str__ = __repr__
 
     def volume(self):
-	return (4.*N.pi/3.) * self.radius**3
+	return (4.*numpy.pi/3.) * self.radius**3
 
     def hasPoint(self, point):
-        return N.fabs((point-self.center).length()-self.radius) < eps
+        return numpy.fabs((point-self.center).length()-self.radius) < eps
 
     def enclosesPoint(self, point):
         return (point - self.center).length() < self.radius
@@ -193,7 +194,7 @@ class Cylinder(GeometricalObject3D):
         self.height = (center2-center1).length()
 
     def volume(self):
-        return N.pi*self.radius*self.radius*self.height
+        return numpy.pi*self.radius*self.radius*self.height
 
     def __repr__(self):
         return 'Cylinder(' + `self.center1` + ', ' + `self.center2` + \
@@ -205,7 +206,7 @@ class Cylinder(GeometricalObject3D):
         pt = center_line.projectionOf(point)
         if pt is None:
             return 0
-        return N.fabs((point - pt).length() - self.radius) < eps
+        return numpy.fabs((point - pt).length() - self.radius) < eps
 
     def enclosesPoint(self, point):
         center_line = LineSegment(self.center1, self.center2)
@@ -385,7 +386,7 @@ class Line(GeometricalObject3D):
     __str__ = __repr__
 
     def volume(self):
-	return 0.
+        return 0.
 
 
 class LineSegment(Line):
@@ -414,7 +415,7 @@ class LineSegment(Line):
             return pt
         return None
 
-    def isWithin(point):
+    def isWithin(self, point):
         v1 = point - self.point
         v2 = point - self.point2
         if abs(v1 * v2) < eps:
@@ -434,9 +435,9 @@ def _addIntersectFunction(f, class1, class2):
 # Box with box
 
 def _intersectBoxBox(box1, box2):
-    c1 = N.maximum(box1.corners[0], box2.corners[0])
-    c2 = N.minimum(box1.corners[1], box2.corners[1])
-    if N.logical_or.reduce(N.greater_equal(c1, c2)):
+    c1 = numpy.maximum(box1.corners[0], box2.corners[0])
+    c2 = numpy.minimum(box1.corners[1], box2.corners[1])
+    if numpy.logical_or.reduce(numpy.greater_equal(c1, c2)):
         return None
     return Box(Vector(c1), Vector(c2))
 
@@ -453,7 +454,7 @@ def _intersectSphereSphere(sphere1, sphere2):
                              max(sphere1.radius, sphere2.radius):
 	return None
     x = 0.5*(d**2 + sphere1.radius**2 - sphere2.radius**2)/d
-    h = N.sqrt(sphere1.radius**2-x**2)
+    h = numpy.sqrt(sphere1.radius**2-x**2)
     normal = r1r2.normal()
     return Circle(sphere1.center + x*normal, normal, h)
 
@@ -463,9 +464,9 @@ _addIntersectFunction(_intersectSphereSphere, Sphere, Sphere)
 
 def _intersectSphereCone(sphere, cone):
     if sphere.center != cone.center:
-	raise GeomError("Not yet implemented")
-    from_center = sphere.radius*N.cos(cone.angle)
-    radius = sphere.radius*N.sin(cone.angle)
+	    raise GeomError("Not yet implemented")
+    from_center = sphere.radius*numpy.cos(cone.angle)
+    radius = sphere.radius*numpy.sin(cone.angle)
     return Circle(cone.center+from_center*cone.axis, cone.axis, radius)
 
 _addIntersectFunction(_intersectSphereCone, Sphere, Cone)
@@ -501,8 +502,8 @@ def _intersectCirclePlane(circle, plane):
 	if x > circle.radius:
 	    return None
 	else:
-	    angle = N.arccos(x/circle.radius)
-	    along_line = N.sin(angle)*circle.radius
+	    angle = numpy.arccos(x/circle.radius)
+	    along_line = numpy.sin(angle)*circle.radius
 	    normal = circle.normal.cross(line.direction)
 	    if line.distanceFrom(circle.center+normal) > x:
 		normal = -normal
@@ -515,8 +516,8 @@ _addIntersectFunction(_intersectCirclePlane, Circle, Plane)
 # Rotation
 #
 def rotateDirection(vector, axis, angle):
-    s = N.sin(angle)
-    c = N.cos(angle)
+    s = numpy.sin(angle)
+    c = numpy.cos(angle)
     c1 = 1-c
     try:
         axis = axis.direction
@@ -755,23 +756,23 @@ def superpositionFit(confs):
               and the RMS distance after the optimal superposition
     """
     w_sum = 0.
-    wr_sum = N.zeros((3,), N.Float)
+    wr_sum = numpy.zeros((3,), numpy.float)
     for w, r_ref, r in confs:
         w_sum += w
         wr_sum += w*r_ref.array
     ref_cms = wr_sum/w_sum
-    pos = N.zeros((3,), N.Float)
+    pos = numpy.zeros((3,), numpy.float)
     possq = 0.
-    cross = N.zeros((3, 3), N.Float)
+    cross = numpy.zeros((3, 3), numpy.float)
     for w, r_ref, r in confs:
         w = w/w_sum
         r_ref = r_ref.array-ref_cms
         r = r.array
         pos = pos + w*r
-        possq = possq + w*N.add.reduce(r*r) \
-                      + w*N.add.reduce(r_ref*r_ref)
-        cross = cross + w*r[:, N.NewAxis]*r_ref[N.NewAxis, :]
-    k = N.zeros((4, 4), N.Float)
+        possq = possq + w*numpy.add.reduce(r*r) \
+                      + w*numpy.add.reduce(r_ref*r_ref)
+        cross = cross + w*r[:, numpy.newaxis]*r_ref[numpy.newaxis, :]
+    k = numpy.zeros((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]
@@ -787,16 +788,16 @@ def superpositionFit(confs):
             k[i, j] = k[j, i]
     k = 2.*k
     for i in range(4):
-        k[i, i] = k[i, i] + possq - N.add.reduce(pos*pos)
+        k[i, i] = k[i, i] + possq - numpy.add.reduce(pos*pos)
     from Scientific import LA
     e, v = LA.eigenvectors(k)
-    i = N.argmin(e)
+    i = numpy.argmin(e)
     v = v[i]
     if v[0] < 0: v = -v
     if e[i] <= 0.:
         rms = 0.
     else:
-        rms = N.sqrt(e[i])
+        rms = numpy.sqrt(e[i])
     from Scientific.Geometry import Quaternion
     return Quaternion.Quaternion(v), Vector(ref_cms), \
            Vector(pos), rms

MMTK/InternalCoordinates.py

@@ -9,8 +9,9 @@ Manipulation of molecular configurations in terms of internal coordinates
 
 __docformat__ = 'restructuredtext'
 
+import numpy
+
 import MMTK
-from Scientific import N
 
 #
 # The abstract base class
@@ -195,7 +196,7 @@ class BondAngle(InternalCoordinate):
         v1 = self.universe.distanceVector(self.atoms[1], self.atoms[0])
         v2 = self.universe.distanceVector(self.atoms[1], self.atoms[2])
         angle = v1.angle(v2)
-        if N.fabs(angle - N.pi) < 1.e-4:
+        if numpy.fabs(angle - numpy.pi) < 1.e-4:
             raise ValueError("angle too close to pi")
         axis = v1.cross(v2).normal()
         d = angle-value
@@ -279,7 +280,7 @@ class DihedralAngle(InternalCoordinate):
                                         self.atoms[2], self.atoms[3])
         v = self.universe.distanceVector(self.atoms[1], self.atoms[2])
         axis = v.normal()
-        d = N.fmod(angle-value, 2.*N.pi)
+        d = numpy.fmod(angle-value, 2.*numpy.pi)
         cm1, th1 = self.fragment1.centerAndMomentOfInertia()
         r1 = self.universe.distanceVector(self.atoms[1], cm1)
         th1 -= self.fragment1.mass()*((r1*r1)*delta-r1.dyadicProduct(r1))