Merge pull request #23 from Allen-Tildesley/numpy2

Applied all NumPy 2.0 fixes according to Ruff rule NPY201

python_examples/averages_module.py

@@ -103,12 +103,12 @@ def run_begin ( variables ):
 
     # Store method options and add-constants in module variables
     method=np.array([variable.method for variable in variables],dtype=np.int_)
-    add   =np.array([variable.add    for variable in variables],dtype=np.float_)
+    add   =np.array([variable.add    for variable in variables],dtype=np.float64)
 
     # Zero averages and error accumulators at start of run
     run_nrm = 0.0
-    run_avg = np.zeros(n_avg,dtype=np.float_)
-    run_err = np.zeros(n_avg,dtype=np.float_)
+    run_avg = np.zeros(n_avg,dtype=np.float64)
+    run_err = np.zeros(n_avg,dtype=np.float64)
 
     # Write headings
     print()
@@ -127,8 +127,8 @@ def blk_begin():
     global blk_nrm, blk_avg, blk_msd
 
     blk_nrm = 0.0
-    blk_avg = np.zeros(n_avg,dtype=np.float_)
-    blk_msd = np.zeros(n_avg,dtype=np.float_)
+    blk_avg = np.zeros(n_avg,dtype=np.float64)
+    blk_msd = np.zeros(n_avg,dtype=np.float64)
 
 def blk_add ( variables ):
     """Increment block-average variables."""
@@ -138,7 +138,7 @@ def blk_add ( variables ):
 
     assert len(variables)==n_avg, 'Mismatched variable arrays'
 
-    values = np.array([variable.val for variable in variables],dtype=np.float_)
+    values = np.array([variable.val for variable in variables],dtype=np.float64)
     blk_avg = blk_avg + values
     blk_msd = blk_msd + values**2
     blk_nrm = blk_nrm + 1.0

python_examples/config_io_module.py

@@ -38,10 +38,10 @@ def read_cnf_atoms ( filename, with_v=False ):
     rows, cols = rv.shape
     assert rows == n, "{:d}{}{:d}".format(rows,' rows not equal to ',n)
     assert cols >= 3, "{:d}{}".format(cols,' cols less than 3')
-    r = rv[:,0:3].astype(np.float_) # Coordinate array
+    r = rv[:,0:3].astype(np.float64) # Coordinate array
     if with_v:
         assert cols >= 6, "{:d}{}".format(cols,' cols less than 6')
-        v = rv[:,3:6].astype(np.float_) # Velocity array
+        v = rv[:,3:6].astype(np.float64) # Velocity array
         return n, box, r, v
     else:
         return n, box, r
@@ -59,12 +59,12 @@ def read_cnf_mols ( filename, with_v=False, quaternions=False ):
     assert rows == n, "{:d}{}{:d}".format(rows,' rows not equal to ',n)
     cols_re = 7 if quaternions else 6
     assert cols >= cols_re, "{:d}{}{:d}".format(cols,' cols less than ',cols_re)
-    r = revw[:,0:3].astype(np.float_) # Coordinate array
-    e = revw[:,3:cols_re].astype(np.float_) # Orientation array
+    r = revw[:,0:3].astype(np.float64) # Coordinate array
+    e = revw[:,3:cols_re].astype(np.float64) # Orientation array
     if with_v:
         assert cols >= cols_re+6, "{:d}{}{:d}".format(cols,' cols less than',cols_re+6)
-        v = revw[:,cols_re  :cols_re+3].astype(np.float_) # Velocity array
-        w = revw[:,cols_re+3:cols_re+6].astype(np.float_) # Angular velocity/momentum array
+        v = revw[:,cols_re  :cols_re+3].astype(np.float64) # Velocity array
+        w = revw[:,cols_re+3:cols_re+6].astype(np.float64) # Angular velocity/momentum array
         return n, box, r, e, v, w
     else:
         return n, box, r, e

python_examples/corfun.py

@@ -125,7 +125,7 @@ def c_exact(t):
 # Initial values
 vt = 0.0
 s  = 0.0
-v  = np.empty(nstep,dtype=np.float_)
+v  = np.empty(nstep,dtype=np.float64)
 np.random.seed()
 
 for t in range(-nequil, nstep): # Loop over steps including an equilibration period
@@ -144,10 +144,10 @@ def c_exact(t):
 # This is very slow.
 # A much faster direct method uses the NumPy correlate function (see below)
 
-c  = np.zeros(nt+1,dtype=np.float_)
-n  = np.zeros(nt+1,dtype=np.float_)
+c  = np.zeros(nt+1,dtype=np.float64)
+n  = np.zeros(nt+1,dtype=np.float64)
 t0 = np.empty(n0,dtype=np.int_)
-v0 = np.empty(n0,dtype=np.float_)
+v0 = np.empty(n0,dtype=np.float64)
 mk   = -1 # Storage location of time origin
 full = False
 
@@ -179,7 +179,7 @@ def c_exact(t):
 fft_len = 2*nstep # Actual length of FFT data
 
 # Prepare data for FFT
-fft_inp = np.zeros(fft_len,dtype=np.complex_) # Fill input array with zeros
+fft_inp = np.zeros(fft_len,dtype=np.complex128) # Fill input array with zeros
 fft_inp[0:nstep] = v                          # Put data into first part (real only)
 
 fft_out = np.fft.fft(fft_inp) # Forward FFT
@@ -189,7 +189,7 @@ def c_exact(t):
 fft_inp = np.fft.ifft(fft_out) # Backward FFT (the factor of 1/fft_len is built in)
 
 # Normalization factors associated with number of time origins
-n = np.linspace(nstep,nstep-nt,nt+1,dtype=np.float_)
+n = np.linspace(nstep,nstep-nt,nt+1,dtype=np.float64)
 assert np.all(n>0.5), 'Normalization array error' # Should never happen
 c_fft = fft_inp[0:nt+1].real / n
 

python_examples/diffusion.py

@@ -111,13 +111,13 @@ def unfold ( r_old, r ):
 print("{:40}{:15d}  ".format('Number of particles', n)   )
 print("{:40}{:15.6f}".format('Box (in sigma units)',box) )
 
-msd  = np.zeros(nt+1,dtype=np.float_)
-rvcf = np.zeros(nt+1,dtype=np.float_)
-vacf = np.zeros(nt+1,dtype=np.float_)
-norm = np.zeros(nt+1,dtype=np.float_)
+msd  = np.zeros(nt+1,dtype=np.float64)
+rvcf = np.zeros(nt+1,dtype=np.float64)
+vacf = np.zeros(nt+1,dtype=np.float64)
+norm = np.zeros(nt+1,dtype=np.float64)
 t0 = np.empty(n0,dtype=np.int_)
-v0 = np.empty((n0,n,3),dtype=np.float_)
-r0 = np.empty((n0,n,3),dtype=np.float_)
+v0 = np.empty((n0,n,3),dtype=np.float64)
+r0 = np.empty((n0,n,3),dtype=np.float64)
 mk   = -1 # Storage location of time origin
 full = False
 

python_examples/diffusion_test.py

@@ -95,10 +95,10 @@ def o_propagator ( t, v ):
 
 np.random.seed()
 
-r = np.random.rand(n,3).astype(np.float_)  # Random positions
+r = np.random.rand(n,3).astype(np.float64)  # Random positions
 r = r - 0.5                                # Now in range (-1/2,1/2)
 r = r * box                                # Now in range (-box/2,box/2)
-v = np.random.randn(n,3).astype(np.float_) # Random velocities
+v = np.random.randn(n,3).astype(np.float64) # Random velocities
 v = v * np.sqrt(temperature)               # At desired temperature
 
 cnf_prefix = 'cnf.'

python_examples/eos_lj_module.py

@@ -43,7 +43,7 @@ def power ( tau, delta, c ):
     # f[i,:] is differentiated i times with respect to tau, and then multiplied by tau**i
     # f[:,j] is differentiated j times with respect to delta, and then multiplied by delta**j
 
-    f = np.full ( (3,3), c['n'] * (tau**c['t']) * (delta**c['d']), dtype=np.float_ )
+    f = np.full ( (3,3), c['n'] * (tau**c['t']) * (delta**c['d']), dtype=np.float64 )
     f[1,:] = f[1,:] * c['t']
     f[2,:] = f[2,:] * c['t'] * ( c['t'] - 1.0 )
     f[:,1] = f[:,1] * c['d']
@@ -60,7 +60,7 @@ def expon ( tau, delta, c ):
     # f[i,:] is differentiated i times with respect to tau, and then multiplied by tau**i
     # f[:,j] is differentiated j times with respect to delta, and then multiplied by delta**j
 
-    f = np.full ( (3,3), c['n'] * (tau**c['t']) * (delta**c['d']) * np.exp(-delta**c['l']), dtype=np.float_ )
+    f = np.full ( (3,3), c['n'] * (tau**c['t']) * (delta**c['d']) * np.exp(-delta**c['l']), dtype=np.float64 )
     f[1,:] = f[1,:] * c['t']
     f[2,:] = f[2,:] * c['t'] * ( c['t'] - 1.0 )
     f[:,1] = f[:,1] * (c['d']-c['l']*delta**c['l'])
@@ -78,7 +78,7 @@ def gauss ( tau, delta, c ):
     # f[:,j] is differentiated j times with respect to delta, and then multiplied by delta**j
 
     f = np.full ( (3,3), c['n']*(tau**c['t'])*np.exp(-c['beta']*(tau-c['gamma'])**2)
-                      *(delta**c['d'])*np.exp(-c['eta']*(delta-c['epsilon'])**2), dtype=np.float_)
+                      *(delta**c['d'])*np.exp(-c['eta']*(delta-c['epsilon'])**2), dtype=np.float64)
     f[1,:] = f[1,:] * ( c['t'] - 2.0*c['beta']*tau*(tau-c['gamma']) )
     f[2,:] = f[2,:] * ( ( c['t'] - 2.0*c['beta']*tau*(tau-c['gamma']) )**2 - c['t'] - 2*c['beta']*tau**2 )
     f[:,1] = f[:,1] * ( c['d'] - 2.0*c['eta']*delta*(delta-c['epsilon']) )
@@ -104,28 +104,28 @@ def a_res_full ( temp, rho ):
     tau   = temp_crit / temp # Reduced inverse temperature
     delta = rho / rho_crit   # Reduced density
 
-    a = np.zeros ( (3,3), dtype=np.float_ )
+    a = np.zeros ( (3,3), dtype=np.float64 )
 
     # Coefficients taken from Table 2 of 
     # M Thol, G Rutkai, A Koester, R Lustig, R Span, J Vrabec, J Phys Chem Ref Data 45, 023101 (2016)
 
-    cp = np.empty ( 6, dtype = [ ('n',np.float_), ('t',np.float_), ('d',np.float_) ] )
+    cp = np.empty ( 6, dtype = [ ('n',np.float64), ('t',np.float64), ('d',np.float64) ] )
     cp['n'] = [ 0.005208073, 2.186252, -2.161016, 1.452700, -2.041792, 0.18695286 ]
     cp['t'] = [ 1.000, 0.320, 0.505, 0.672, 0.843, 0.898 ]
     cp['d'] = [ 4.0,   1.0,   1.0,   2.0,   2.0,   3.0   ]
     for c in cp:
         a = a + power ( tau, delta, c )
 
-    ce = np.empty ( 6, dtype = [ ('n',np.float_), ('t',np.float_), ('d',np.float_), ('l',np.float_) ] )
+    ce = np.empty ( 6, dtype = [ ('n',np.float64), ('t',np.float64), ('d',np.float64), ('l',np.float64) ] )
     ce['n'] = [ -0.090988445, -0.49745610, 0.10901431, -0.80055922, -0.56883900, -0.62086250 ]
     ce['t'] = [ 1.294, 2.590, 1.786, 2.770, 1.786, 1.205 ]
     ce['d'] = [ 5.0,   2.0,   2.0,   3.0,   1.0,   1.0   ]
     ce['l'] = [ 1.0,   2.0,   1.0,   2.0,   2.0,   1.0   ]
     for c in ce:
         a = a + expon ( tau, delta, c )
 
-    cg = np.empty ( 11, dtype = [ ('n',np.float_), ('t',np.float_), ('d',np.float_), ('eta',np.float_),
-                                  ('beta',np.float_), ('gamma',np.float_), ('epsilon',np.float_)] )
+    cg = np.empty ( 11, dtype = [ ('n',np.float64), ('t',np.float64), ('d',np.float64), ('eta',np.float64),
+                                  ('beta',np.float64), ('gamma',np.float64), ('epsilon',np.float64)] )
     cg['n']       = [ -1.4667177,  1.8914690,  -0.13837010, -0.38696450, 0.12657020, 0.6057810,
                        1.1791890, -0.47732679, -9.9218575,  -0.57479320, 0.0037729230 ]
     cg['t']       = [ 2.830,  2.548, 4.650, 1.385, 1.460, 1.351, 0.660, 1.496,   1.830, 1.616, 4.970 ]
@@ -157,28 +157,28 @@ def a_res_cutshift ( temp, rho ):
     tau   = temp_crit / temp # Reduced inverse temperature
     delta = rho / rho_crit   # Reduced density
 
-    a = np.zeros ( (3,3), dtype=np.float_ )
+    a = np.zeros ( (3,3), dtype=np.float64 )
 
     # Coefficients taken from Table 1 of 
     # M Thol, G Rutkai, R Span, J Vrabec, R Lustig, Int J Thermophys 36, 25 (2015) 
 
-    cp = np.empty ( 6, dtype = [ ('n',np.float_), ('t',np.float_), ('d',np.float_) ] )
+    cp = np.empty ( 6, dtype = [ ('n',np.float64), ('t',np.float64), ('d',np.float64) ] )
     cp['n'] = [ 0.015606084, 1.7917527, -1.9613228, 1.3045604, -1.8117673, 0.15483997 ]
     cp['t'] = [ 1.000, 0.304, 0.583, 0.662, 0.870, 0.870 ]
     cp['d'] = [ 4.0,   1.0,   1.0,   2.0,   2.0,   3.0   ]
     for c in cp:
         a = a + power ( tau, delta, c )
 
-    ce = np.empty ( 6, dtype = [ ('n',np.float_), ('t',np.float_), ('d',np.float_), ('l',np.float_) ] )
+    ce = np.empty ( 6, dtype = [ ('n',np.float64), ('t',np.float64), ('d',np.float64), ('l',np.float64) ] )
     ce['n'] = [ -0.094885204, -0.20092412,  0.11639644, -0.50607364, -0.58422807, -0.47510982 ]
     ce['t'] = [  1.250, 3.000, 1.700, 2.400, 1.960, 1.286 ]
     ce['d'] = [  5.0,   2.0,   2.0,   3.0,   1.0,   1.0   ]
     ce['l'] = [  1.0,   2.0,   1.0,   2.0,   2.0,   1.0   ]
     for c in ce:
         a = a + expon ( tau, delta, c )
 
-    cg = np.empty ( 9, dtype = [ ('n',np.float_), ('t',np.float_), ('d',np.float_), ('eta',np.float_),
-                                 ('beta',np.float_), ('gamma',np.float_), ('epsilon',np.float_)] )
+    cg = np.empty ( 9, dtype = [ ('n',np.float64), ('t',np.float64), ('d',np.float64), ('eta',np.float64),
+                                 ('beta',np.float64), ('gamma',np.float64), ('epsilon',np.float64)] )
     cg['n']       = [  0.0094333106, 0.30444628,  -0.0010820946, -0.099693391, 0.0091193522,
                        0.12970543,   0.023036030, -0.082671073,  -2.2497821 ]
     cg['t']       = [  3.600, 2.080, 5.240, 0.960, 1.360, 1.655, 0.900, 0.860,  3.950 ]

python_examples/error_calc.py

@@ -104,7 +104,7 @@
 
 # Data generation
 np.random.seed()
-a = np.empty(nstep,dtype=np.float_)
+a = np.empty(nstep,dtype=np.float64)
 a_avg = 0.0 # Zero average accumulator
 a_var = 0.0 # Zero mean-squared accumulator
 

python_examples/ewald.py

@@ -73,7 +73,7 @@
 r = r / box        # Work throughout in unit box
 r = r - np.rint(r) # Apply periodic boundaries
 assert n%2 == 0, 'Error, we require even N for charge balance'
-q = np.empty(n,dtype=np.float_)
+q = np.empty(n,dtype=np.float64)
 q[::2] = 1.0
 q[1::2] = -1.0
 print ( "{:40}{:15.6f}".format('Net charge', np.sum(q)) )
@@ -101,13 +101,13 @@
 # Potentials are stored according to squared distance of neighbouring box
 print('Brute force method')
 
-pot_shell = np.zeros(nbox_sq+1, dtype=np.float_) # Zero array of potentials (not all the elements of this array will be filled)
+pot_shell = np.zeros(nbox_sq+1, dtype=np.float64) # Zero array of potentials (not all the elements of this array will be filled)
 
 for xbox, ybox, zbox in product ( range(-nbox,nbox+1), repeat=3 ): # Triple loop over surrounding periodic boxes
     rbox_sq = xbox**2 + ybox**2 + zbox**2
     if rbox_sq > nbox_sq: # Skip if outside maximum sphere of boxes
         continue
-    rbox_vec = np.array([xbox,ybox,zbox], dtype=np.float_)
+    rbox_vec = np.array([xbox,ybox,zbox], dtype=np.float64)
     start_shift = 0 if rbox_sq>0 else 1 # Skip all i==j if in central box
     for shift in range(start_shift,n):
         # Bare Coulomb term, including box vector, no periodic box correction

python_examples/ewald_module.py

@@ -79,7 +79,7 @@ def pot_k_ewald ( nk, r, q, kappa ):
 
         b = 1.0 / 4.0 / kappa / kappa
         k_sq_max = nk**2            # Store this in module data
-        kfac = np.zeros(k_sq_max+1,dtype=np.float_)
+        kfac = np.zeros(k_sq_max+1,dtype=np.float64)
 
         # Lazy triple loop, which over-writes the same values of ksq several times
         # and leaves some locations empty (those which are not the sum of three squared integers)
@@ -97,9 +97,9 @@ def pot_k_ewald ( nk, r, q, kappa ):
     # Double-check value on later calls
     assert k_sq_max == nk**2, 'nk error'
 
-    eikx = np.zeros((n,nk+1),  dtype=np.complex_) # omits negative k indices
-    eiky = np.zeros((n,2*nk+1),dtype=np.complex_) # includes negative k indices
-    eikz = np.zeros((n,2*nk+1),dtype=np.complex_) # includes negative k indices
+    eikx = np.zeros((n,nk+1),  dtype=np.complex128) # omits negative k indices
+    eiky = np.zeros((n,2*nk+1),dtype=np.complex128) # includes negative k indices
+    eikz = np.zeros((n,2*nk+1),dtype=np.complex128) # includes negative k indices
 
     # Calculate kx, ky, kz = 0, 1 explicitly
     eikx[:,   0] = 1.0 + 0.0j
@@ -165,7 +165,7 @@ def pot_k_pm_ewald ( nk, r, q, kappa ):
     # Assign charge density to complex array ready for FFT
     # Assume r in unit box with range (-0.5,0.5)
     # Mesh function expects coordinates in range 0..1 so we add 0.5
-    fft_inp = mesh_function ( r+0.5, q, 2*nk ).astype(np.complex_)
+    fft_inp = mesh_function ( r+0.5, q, 2*nk ).astype(np.complex128)
 
     # Forward FFT incorporating scaling by number of grid points
     fft_out = np.fft.fftn(fft_inp) / (sc**3)

python_examples/fft3dwrap.py

@@ -85,7 +85,7 @@
 # For a function of the form f(x)*f(y)*f(z) such as this, there are faster ways of initializing the array
 # However, we retain the generic triple loop to make it easy to choose a different function
 
-fft_inp = np.empty((sc,sc,sc),dtype=np.complex_)
+fft_inp = np.empty((sc,sc,sc),dtype=np.complex128)
 
 print()
 print("Initial real-space Gaussian")

python_examples/grint.py

@@ -164,9 +164,9 @@ def f2tanh ( z, z_gl, z_lg, width, rho_g, rho_l ):
 r_vals   = np.linspace ( dr/2, (nr-0.5)*dr, nr )
 
 # Zero the accumulator arrays
-dens = np.zeros ( nk, dtype=np.float_ )
-rho1 = np.zeros ( 2*nz+1, dtype=np.float_ )
-rho2 = np.zeros ( ( 2*nz+1, 2*nc+1, nr ), dtype=np.float_ )
+dens = np.zeros ( nk, dtype=np.float64 )
+rho1 = np.zeros ( 2*nz+1, dtype=np.float64 )
+rho2 = np.zeros ( ( 2*nz+1, 2*nc+1, nr ), dtype=np.float64 )
 
 norm = 0
 t    = 0

python_examples/hit_and_miss.py

@@ -34,7 +34,7 @@
 print('hit_and_miss')
 np.random.seed()
 
-r_0 = np.array([1.0,2.0,3.0],dtype=np.float_)
+r_0 = np.array([1.0,2.0,3.0],dtype=np.float64)
 v_0 = np.prod ( r_0 )
 
 tau_hit  = 0

python_examples/initialize.py

@@ -44,17 +44,17 @@ def fcc_positions ( n, box, length, soft, quaternions ):
     assert n==4*nc**3, "{}{:d}{:d}".format('n, nc mismatch ',n,4*nc**3)
     cell = box / nc  # Unit cell
     box2 = box / 2.0 # Half box length
-    r = np.empty((n,3),dtype=np.float_)
-    e = np.empty((n,3),dtype=np.float_)
+    r = np.empty((n,3),dtype=np.float64)
+    e = np.empty((n,3),dtype=np.float64)
 
-    r_fcc = np.array ( [ [0.25,0.25,0.25],[0.25,0.75,0.75],[0.75,0.75,0.25],[0.75,0.25,0.75] ], dtype=np.float_ )
-    e_fcc = np.array ( [ [1.0,1.0,1.0],[1.0,-1.0,-1.0],[-1.0,1.0,-1.0],[-1.0,-1.0,1.0] ],dtype=np.float_)*np.sqrt(1.0/3.0)
+    r_fcc = np.array ( [ [0.25,0.25,0.25],[0.25,0.75,0.75],[0.75,0.75,0.25],[0.75,0.25,0.75] ], dtype=np.float64 )
+    e_fcc = np.array ( [ [1.0,1.0,1.0],[1.0,-1.0,-1.0],[-1.0,1.0,-1.0],[-1.0,-1.0,1.0] ],dtype=np.float64)*np.sqrt(1.0/3.0)
 
     i = 0
 
     for ix, iy, iz in product(range(nc),repeat=3): # triple loop over unit cells
         for a in range(4): # loop over atoms in unit cell
-            r[i,:] = r_fcc[a,:] + np.array ( [ix,iy,iz] ).astype(np.float_) # in range 0..nc
+            r[i,:] = r_fcc[a,:] + np.array ( [ix,iy,iz] ).astype(np.float64) # in range 0..nc
             r[i,:] = r[i,:] * cell                                          # in range 0..box
             r[i,:] = r[i,:] - box2                                          # in range -box2..box2
             e[i,:] = e_fcc[a]
@@ -81,11 +81,11 @@ def ran_positions ( n, box, length, soft, quaternions ):
 
     print('Random positions')
 
-    r = np.empty((n,3),dtype=np.float_)
+    r = np.empty((n,3),dtype=np.float64)
     if quaternions:
-        e = np.empty((n,4),dtype=np.float_)
+        e = np.empty((n,4),dtype=np.float64)
     else:
-        e = np.empty((n,3),dtype=np.float_)
+        e = np.empty((n,3),dtype=np.float64)
 
     for i in range(n):
 
@@ -148,7 +148,7 @@ def ran_velocities ( nn, e, temperature, inertia, quaternions ):
         w = np.random.randn ( n, 3 ) * w_std_dev
     else:
         w_sq_mean = 2.0 * temperature / inertia
-        w = np.empty ( (n,3), dtype=np.float_ )
+        w = np.empty ( (n,3), dtype=np.float64 )
         for i in range(n):
             w[i,:] = random_perpendicular_vector ( e[i,:] ) # Set direction of angular velocity
             w[i,:] = w[i,:] * np.sqrt(np.random.exponential(w_sq_mean))
@@ -166,7 +166,7 @@ def chain_positions ( n, bond, soft ):
 
     print("{:40}{:15.6f}".format('Chain, randomly oriented bonds = ',bond) )
 
-    r = np.empty ( (n,3), dtype=np.float_ )
+    r = np.empty ( (n,3), dtype=np.float64 )
 
     r[0,:] = [0.0,0.0,0.0] # First atom at origin
     r[1,:] = bond*random_vector() # Second atom at random position (bond length away)
@@ -233,7 +233,7 @@ def chain_velocities ( nn, temperature, constraints, r ):
 
     # Compute total angular momentum and moment of inertia tensor
     ang_mom = np.sum ( np.cross ( r, v ), axis=0 )
-    inertia = np.zeros ( (3,3), dtype=np.float_ )
+    inertia = np.zeros ( (3,3), dtype=np.float64 )
     for i in range(n):
         inertia = inertia - np.outer ( r[i,:], r[i,:] )
         for xyz in range(3):