diff --git a/python_examples/averages_module.py b/python_examples/averages_module.py index 4ac2d39..a11abe2 100644 --- a/python_examples/averages_module.py +++ b/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 diff --git a/python_examples/config_io_module.py b/python_examples/config_io_module.py index 274cf8f..91f13ad 100644 --- a/python_examples/config_io_module.py +++ b/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 diff --git a/python_examples/corfun.py b/python_examples/corfun.py index 8334389..5dcef48 100755 --- a/python_examples/corfun.py +++ b/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 diff --git a/python_examples/diffusion.py b/python_examples/diffusion.py index 69b11a9..6b53ef3 100755 --- a/python_examples/diffusion.py +++ b/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 diff --git a/python_examples/diffusion_test.py b/python_examples/diffusion_test.py index 4a15ef6..cb03dd0 100755 --- a/python_examples/diffusion_test.py +++ b/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.' diff --git a/python_examples/eos_lj_module.py b/python_examples/eos_lj_module.py index 3e33f56..57d5d93 100644 --- a/python_examples/eos_lj_module.py +++ b/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,19 +104,19 @@ 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 ] @@ -124,8 +124,8 @@ def a_res_full ( temp, rho ): 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,19 +157,19 @@ 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 ] @@ -177,8 +177,8 @@ def a_res_cutshift ( temp, rho ): 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 ] diff --git a/python_examples/error_calc.py b/python_examples/error_calc.py index 36dad71..731c83d 100755 --- a/python_examples/error_calc.py +++ b/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 diff --git a/python_examples/ewald.py b/python_examples/ewald.py index 8368b84..b933ad3 100755 --- a/python_examples/ewald.py +++ b/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 diff --git a/python_examples/ewald_module.py b/python_examples/ewald_module.py index 9f975ca..d328591 100644 --- a/python_examples/ewald_module.py +++ b/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) diff --git a/python_examples/fft3dwrap.py b/python_examples/fft3dwrap.py index aa5c547..cea22fe 100755 --- a/python_examples/fft3dwrap.py +++ b/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") diff --git a/python_examples/grint.py b/python_examples/grint.py index 0efba36..5dc3146 100755 --- a/python_examples/grint.py +++ b/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 diff --git a/python_examples/hit_and_miss.py b/python_examples/hit_and_miss.py index b97d42e..2c08cc3 100755 --- a/python_examples/hit_and_miss.py +++ b/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 diff --git a/python_examples/initialize.py b/python_examples/initialize.py index bfeb724..63d0e6b 100755 --- a/python_examples/initialize.py +++ b/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): diff --git a/python_examples/maths_module.py b/python_examples/maths_module.py index ad71d57..f5e1c73 100644 --- a/python_examples/maths_module.py +++ b/python_examples/maths_module.py @@ -42,7 +42,7 @@ def random_vector(): phi = zeta[1] * 2.0*np.pi # Random angle uniformly sampled in range (0,2*pi) - return np.array ( ( s*np.cos(phi), s*np.sin(phi), c ), dtype=np.float_ ) # Random unit vector + return np.array ( ( s*np.cos(phi), s*np.sin(phi), c ), dtype=np.float64 ) # Random unit vector def random_perpendicular_vector ( old ): """Returns a uniformly sampled unit vector perpendicular to the old vector.""" @@ -88,7 +88,7 @@ def random_quaternion(): break f = np.sqrt ( ( 1.0 - norm1 ) / norm2 ) - return np.array ( ( zeta[0], zeta[1], beta[0]*f, beta[1]*f ), dtype=np.float_ ) # Random quaternion + return np.array ( ( zeta[0], zeta[1], beta[0]*f, beta[1]*f ), dtype=np.float64 ) # Random quaternion def random_rotate_quaternion ( angle_max, old ): """Returns a unit quaternion rotated by a maximum angle (in radians) relative to the old quaternion.""" @@ -185,7 +185,7 @@ def rotate_quaternion ( angle, axis, old ): # Standard formula for rotation quaternion, using half angles rot = np.sin(0.5*angle) * axis - rot = np.array([np.cos(0.5*angle),rot[0],rot[1],rot[2]],dtype=np.float_) + rot = np.array([np.cos(0.5*angle),rot[0],rot[1],rot[2]],dtype=np.float64) e = quatmul ( rot, old ) # Apply rotation to old quaternion return e @@ -201,7 +201,7 @@ def quatmul ( a, b ): return np.array ( [ a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3], a[1]*b[0] + a[0]*b[1] - a[3]*b[2] + a[2]*b[3], a[2]*b[0] + a[3]*b[1] + a[0]*b[2] - a[1]*b[3], - a[3]*b[0] - a[2]*b[1] + a[1]*b[2] + a[0]*b[3] ], dtype=np.float_ ) + a[3]*b[0] - a[2]*b[1] + a[1]*b[2] + a[0]*b[3] ], dtype=np.float64 ) def nematic_order ( e ): """Returns a nematic orientational order parameter.""" @@ -242,7 +242,7 @@ def q_to_a ( q ): assert np.isclose(np.sum(q**2),1.0), 'quaternion normalization error {} {} {} {}'.format(*q) # Write out row by row, for clarity - a = np.empty( (3,3), dtype=np.float_ ) + a = np.empty( (3,3), dtype=np.float64 ) a[0,:] = [ q[0]**2+q[1]**2-q[2]**2-q[3]**2, 2*(q[1]*q[2]+q[0]*q[3]), 2*(q[1]*q[3]-q[0]*q[2]) ] a[1,:] = [ 2*(q[1]*q[2]-q[0]*q[3]), q[0]**2-q[1]**2+q[2]**2-q[3]**2, 2*(q[2]*q[3]+q[0]*q[1]) ] a[2,:] = [ 2*(q[1]*q[3]+q[0]*q[2]), 2*(q[2]*q[3]-q[0]*q[1]), q[0]**2-q[1]**2-q[2]**2+q[3]**2 ] diff --git a/python_examples/mc_chain_lj_module.py b/python_examples/mc_chain_lj_module.py index b2ba603..7087472 100644 --- a/python_examples/mc_chain_lj_module.py +++ b/python_examples/mc_chain_lj_module.py @@ -87,8 +87,8 @@ def regrow ( temperature, m_max, k_max, bond, k_spring, r ): w_tol = 1.e-10 # Min weight tolerance n, d = r.shape assert d==3, 'Dimension error in regrow' - r_try = np.empty((k_max,3),dtype=np.float_) - w = np.empty(k_max,dtype=np.float_) + r_try = np.empty((k_max,3),dtype=np.float64) + w = np.empty(k_max,dtype=np.float64) std = np.sqrt(temperature/k_spring) # Spring bond standard deviation d_max = 3.0*std # Impose a limit on variation, say 3*std diff --git a/python_examples/mc_chain_nvt_sw.py b/python_examples/mc_chain_nvt_sw.py index 8886321..45b0b3b 100755 --- a/python_examples/mc_chain_nvt_sw.py +++ b/python_examples/mc_chain_nvt_sw.py @@ -175,7 +175,7 @@ def update_histogram(q,q_min,q_max): print( "{:40}{:15d}".format('Number of pivot tries per step', n_pivot) ) # Calculate energy histogram and Boltzmann exponents, s(q), which are just values of E/kT -e = -np.arange(nq+1,dtype=np.float_) +e = -np.arange(nq+1,dtype=np.float64) e[0]=0.0 s = e /temperature h = np.zeros_like(e) diff --git a/python_examples/mc_chain_sw_module.py b/python_examples/mc_chain_sw_module.py index 93d0161..a872bae 100644 --- a/python_examples/mc_chain_sw_module.py +++ b/python_examples/mc_chain_sw_module.py @@ -63,8 +63,8 @@ def regrow ( s, m_max, k_max, bond, q_range, r, q ): w_tol = 1.e-10 # Min weight tolerance n, d = r.shape assert d==3, 'Dimension error in regrow' - r_try = np.empty((k_max,3),dtype=np.float_) - w = np.empty(k_max,dtype=np.float_) + r_try = np.empty((k_max,3),dtype=np.float64) + w = np.empty(k_max,dtype=np.float64) r_old = np.copy(r) # Store copy of r q_old = q # Store old q @@ -93,7 +93,7 @@ def regrow ( s, m_max, k_max, bond, q_range, r, q ): r0 = np.copy(r[0,:]) r[:n-m,:] = r[:n-m,:] - r0 - w_new = np.float_(1.0) + w_new = np.float64(1.0) for i in range(n-m,n): # Loop to regrow last m atoms, computing new weight for k in range(k_max): # Loop over k_max tries r_try[k,:] = r[i-1,:] + bond * random_vector() # Trial position in random direction from i-1 @@ -121,7 +121,7 @@ def regrow ( s, m_max, k_max, bond, q_range, r, q ): else: # Remove and reconstruct at start r[:,:] = r_old[::-1,:] # Copy and reverse all n atoms - w_old = np.float_(1.0) + w_old = np.float64(1.0) for i in range(n-m,n): # Old position and weight are stored as try 0 @@ -387,7 +387,7 @@ def accept ( s, q_old, q_new, w_old=1.0, w_new=1.0 ): # For Wang-Landau, s(q) is the entropy function, which is updated through the run # In either case, weights (real, but taking positive integer values here) may optionally appear - w_tol = np.float_(0.5) + w_tol = np.float64(0.5) # Check that we are within bounds for the entropy table assert q_new>=0 and q_newi -coltime = np.full ( n, 1.0e9, dtype=np.float_ ) +coltime = np.full ( n, 1.0e9, dtype=np.float64 ) partner = np.full ( n, n-1, dtype=np.int_ ) for i in range(n-1): coltime[i], partner[i] = update ( i, box, r[i:,:], v[i:,:] ) diff --git a/python_examples/md_nvt_lj.py b/python_examples/md_nvt_lj.py index fd06823..a93310c 100755 --- a/python_examples/md_nvt_lj.py +++ b/python_examples/md_nvt_lj.py @@ -243,12 +243,12 @@ def u4_propagator ( t, j_list ): # Initial values of thermal variables g = 3*(n-1) -q = np.full(m,temperature * tau**2,dtype=np.float_) +q = np.full(m,temperature * tau**2,dtype=np.float64) q[0] = g * temperature * tau**2 fmt_string = "{:15.6f}"*m fmt_string = "{:40}"+fmt_string print(fmt_string.format('Thermal inertias Q',*q)) -eta = np.zeros(m,dtype=np.float_) +eta = np.zeros(m,dtype=np.float64) p_eta = np.random.randn(m)*np.sqrt(temperature) p_eta = p_eta * np.sqrt(q) diff --git a/python_examples/md_nvt_poly_lj.py b/python_examples/md_nvt_poly_lj.py index f4f2490..3f5af71 100755 --- a/python_examples/md_nvt_poly_lj.py +++ b/python_examples/md_nvt_poly_lj.py @@ -242,7 +242,7 @@ def calc_variables ( ): r = r - np.rint ( r ) # Periodic boundaries # Calculate all bond vectors -d = np.empty ( (n,na,3), dtype=np.float_ ) +d = np.empty ( (n,na,3), dtype=np.float64 ) norm = np.sqrt ( np.sum(e**2,axis=1) ) # Quaternion norms e = e / norm[:,np.newaxis] # Ensure normalized quaternions for i, ei in enumerate(e): diff --git a/python_examples/md_poly_lj_module.py b/python_examples/md_poly_lj_module.py index 94837f0..f4e16a8 100644 --- a/python_examples/md_poly_lj_module.py +++ b/python_examples/md_poly_lj_module.py @@ -38,7 +38,7 @@ na = 3 db = np.array([[-np.sin(alpha2), 0.0, -np.cos(alpha2)/3.0], [0.0, 0.0, 2.0*np.cos(alpha2)/3.0], - [np.sin(alpha2), 0.0, -np.cos(alpha2)/3.0]], dtype=np.float_) + [np.sin(alpha2), 0.0, -np.cos(alpha2)/3.0]], dtype=np.float64) diameter = 2.0 * np.sqrt ( np.max ( np.sum(db**2,axis=1) ) ) # Molecular diameter # Cutoff distance and force-shift parameters (all private) chosen as per the reference: @@ -50,7 +50,7 @@ lambda1 = 4.0*(7.0*sr_cut6-13.0*sr_cut12) lambda2 = -24.0*(sr_cut6-2.0*sr_cut12)*sr_cut -m = np.array([ 1.0/3.0, 1.0/3.0, 1.0/3.0 ],dtype=np.float_) # Masses add up to 1.0 +m = np.array([ 1.0/3.0, 1.0/3.0, 1.0/3.0 ],dtype=np.float64) # Masses add up to 1.0 # The following section sets the diagonal moments of inertia "realistically" # based on the values of atomic masses and bond vectors (above), with some checking diff --git a/python_examples/mesh.py b/python_examples/mesh.py index bba2d80..d5c8a83 100755 --- a/python_examples/mesh.py +++ b/python_examples/mesh.py @@ -68,7 +68,7 @@ r = np.random.random_sample( (n,3) ) # For illustration we choose +1 and -1 charges, alternately -q = np.empty(n,dtype=np.float_) +q = np.empty(n,dtype=np.float64) q[::2] = 1.0 q[1::2] = -1.0 diff --git a/python_examples/mesh_module.py b/python_examples/mesh_module.py index 67b4cf2..1e45b63 100755 --- a/python_examples/mesh_module.py +++ b/python_examples/mesh_module.py @@ -42,8 +42,8 @@ def mesh_function ( r, q, sc ): h = 1.0 / sc # mesh spacing - rho = np.zeros( (sc,sc,sc), dtype=np.float_ ) # Zero mesh - v = np.empty( (3,3), dtype=np.float_ ) # Array of weights + rho = np.zeros( (sc,sc,sc), dtype=np.float64 ) # Zero mesh + v = np.empty( (3,3), dtype=np.float64 ) # Array of weights for i in range(n): # Loop over charges nr = np.rint(r[i,:]*sc).astype(np.int_) # Nearest mesh point indices diff --git a/python_examples/qmc_pi_lj.py b/python_examples/qmc_pi_lj.py index 0b4bec0..54d4767 100755 --- a/python_examples/qmc_pi_lj.py +++ b/python_examples/qmc_pi_lj.py @@ -197,7 +197,7 @@ def rad_gyr ( r ): print( "{:40}{:15d} ".format('Number of particles', n) ) print( "{:40}{:15.6f}".format('Box length', box) ) print( "{:40}{:15.6f}".format('Density', n/box**3) ) -r = np.empty( (p,n,3), dtype=np.float_ ) +r = np.empty( (p,n,3), dtype=np.float64 ) for k in range(p): n, box, r[k,:,:] = read_cnf_atoms ( cnf_prefix+str(k).zfill(2)+'.'+inp_tag) r = r / box # Convert positions to box units @@ -225,7 +225,7 @@ def rad_gyr ( r ): for i in range(n): # Loop over atoms # Centre of mass move - dc = random_translate_vector ( dc_max/box, np.zeros(3,dtype=np.float_) ) + dc = random_translate_vector ( dc_max/box, np.zeros(3,dtype=np.float64) ) partial_old = PotentialType ( 0.0, 0.0, False ) partial_new = PotentialType ( 0.0, 0.0, False ) for k in range(p): # Loop over ring polymer indices diff --git a/python_examples/qmc_pi_sho.py b/python_examples/qmc_pi_sho.py index 8fc00fd..22d5703 100755 --- a/python_examples/qmc_pi_sho.py +++ b/python_examples/qmc_pi_sho.py @@ -172,7 +172,7 @@ def e_pi_sho ( p, beta ): k_spring = p * temperature**2 np.random.seed() -x = np.zeros(p,dtype=np.float_) # Set up initial positions at origin +x = np.zeros(p,dtype=np.float64) # Set up initial positions at origin # Calculate initial values pot_cl = 0.5 * np.sum ( x**2 ) / p # Classical potential energy diff --git a/python_examples/qmc_walk_sho.py b/python_examples/qmc_walk_sho.py index 421b3bb..26e481d 100755 --- a/python_examples/qmc_walk_sho.py +++ b/python_examples/qmc_walk_sho.py @@ -96,9 +96,9 @@ assert n_target <= n_max, 'n_target exceeds n_max at start' -x = np.empty(n_max,dtype=np.float_) # position of each walker -v = np.empty(n_max,dtype=np.float_) # potential energy of each walker -k = np.empty(n_max,dtype=np.float_) # floating point number of replicas to make +x = np.empty(n_max,dtype=np.float64) # position of each walker +v = np.empty(n_max,dtype=np.float64) # potential energy of each walker +k = np.empty(n_max,dtype=np.float64) # floating point number of replicas to make replica = np.empty(n_max,dtype=np.int_) # integer number of replicas to make alive = np.empty(n_max,dtype=np.bool_) # flag those walkers still alive bin = np.zeros(n_bin,dtype=np.int_) # histogram bins for wavefunction @@ -166,7 +166,7 @@ # Normalize the wavefunction so that the integral over psi**2 = 1.0 dx = edges[1]-edges[0] -psi = bin.astype(np.float_) # un-normalized wave function +psi = bin.astype(np.float64) # un-normalized wave function factor = np.sum ( psi**2 ) * dx # integral, assuming that -x_max .. +x_max catches everything psi = psi / np.sqrt( factor ) # normalizing factor diff --git a/python_examples/sample_mean.py b/python_examples/sample_mean.py index ffd3cd5..5972744 100755 --- a/python_examples/sample_mean.py +++ b/python_examples/sample_mean.py @@ -34,7 +34,7 @@ print('sample_mean') np.random.seed() -r_0 = np.array([1.0,2.0],dtype=np.float_) +r_0 = np.array([1.0,2.0],dtype=np.float64) a_0 = np.prod ( r_0 ) f = 0.0 diff --git a/python_examples/t_tensor.py b/python_examples/t_tensor.py index 9135d26..f4ae037 100755 --- a/python_examples/t_tensor.py +++ b/python_examples/t_tensor.py @@ -93,7 +93,7 @@ def t4_tensor ( r, r5 ): from itertools import product # Define 3x3 unit matrix or Kronecker delta - u = np.zeros((3,3),dtype=np.float_) + u = np.zeros((3,3),dtype=np.float64) u[0,0]=u[1,1]=u[2,2]=1.0 t4 = 105.0 * np.einsum('i,j,k,l->ijkl',r,r,r,r) # Starting point: outer product of four vectors @@ -120,7 +120,7 @@ def t5_tensor ( r, r6 ): from itertools import product # Define 3x3 unit matrix or Kronecker delta - u = np.zeros((3,3),dtype=np.float_) + u = np.zeros((3,3),dtype=np.float64) u[0,0]=u[1,1]=u[2,2]=1.0 t5 = 945.0 * np.einsum('i,j,k,l,m->ijklm',r,r,r,r,r) # Starting point: outer product of five vectors @@ -147,7 +147,7 @@ def t5_tensor ( r, r6 ): def skew ( a ): """Returns contraction of supplied 3x3 matrix with Levi-Civita tensor.""" - b = np.empty(3,dtype=np.float_) + b = np.empty(3,dtype=np.float64) b[0] = a[1,2] - a[2,1] b[1] = a[2,0] - a[0,2] b[2] = a[0,1] - a[1,0] diff --git a/python_examples/test_pot_atom.py b/python_examples/test_pot_atom.py index 18fd7ce..81631b9 100755 --- a/python_examples/test_pot_atom.py +++ b/python_examples/test_pot_atom.py @@ -34,7 +34,7 @@ def random_positions(n): from maths_module import random_vector import sys - r = np.zeros((n,3),dtype=np.float_) + r = np.zeros((n,3),dtype=np.float64) # atom 0 is always at the origin, now place the others randomly for i in range(1,r.shape[0]): for pos_try in range(npos): diff --git a/python_examples/test_pot_bend.py b/python_examples/test_pot_bend.py index e0955e7..f6d5743 100644 --- a/python_examples/test_pot_bend.py +++ b/python_examples/test_pot_bend.py @@ -48,7 +48,7 @@ def force ( r ): # Store C coefficients in a matrix # In the general case we would not need to calculate every pair # and also we would make use of the symmetry cc[a,b]=cc[b,a] - cc = np.zeros((n,n),dtype=np.float_) # Create C array (scalar products) + cc = np.zeros((n,n),dtype=np.float64) # Create C array (scalar products) for a in range(1,n): for b in range(1,n): cc[a,b]=np.dot(d[a,:],d[b,:]) # Compute C array (zero indices not used) diff --git a/python_examples/test_pot_linear.py b/python_examples/test_pot_linear.py index bce42c6..1bf7e7d 100755 --- a/python_examples/test_pot_linear.py +++ b/python_examples/test_pot_linear.py @@ -34,7 +34,7 @@ def random_positions(n): from maths_module import random_vector import sys - r = np.zeros((n,3),dtype=np.float_) + r = np.zeros((n,3),dtype=np.float64) # molecule 0 is always at the origin, now place the others randomly for i in range(1,r.shape[0]): for pos_try in range(npos): @@ -56,7 +56,7 @@ def random_orientations(n): """Returns n random 3-d vectors in a numpy array (n,3).""" import numpy as np from maths_module import random_vector - e = np.empty((n,3),dtype=np.float_) + e = np.empty((n,3),dtype=np.float64) for i in range(e.shape[0]): e[i,:] = random_vector() return e @@ -152,8 +152,8 @@ def random_orientations(n): print( "{:5d}{:>10}{:15.6f}{:15.6f}{:15.4e}".format(i,cf[xyz],f_exact,fnum,f_exact-fnum) ) ct = ['Tx','Ty','Tz'] -axis = np.empty(3,dtype=np.float_) -esave = np.empty(3,dtype=np.float_) +axis = np.empty(3,dtype=np.float64) +esave = np.empty(3,dtype=np.float64) for (i,xyz), t_exact in np.ndenumerate(t): axis[:] = 0.0 diff --git a/python_examples/test_pot_twist.py b/python_examples/test_pot_twist.py index 27d7f98..8089c9b 100644 --- a/python_examples/test_pot_twist.py +++ b/python_examples/test_pot_twist.py @@ -48,7 +48,7 @@ def force ( r ): # Store C coefficients in a matrix # In the general case we would not need to calculate every pair # and also we would make use of the symmetry cc[a,b]=cc[b,a] - cc = np.zeros((n,n),dtype=np.float_) # Create C array (scalar products) + cc = np.zeros((n,n),dtype=np.float64) # Create C array (scalar products) for a in range(1,n): for b in range(1,n): cc[a,b]=np.dot(d[a,:],d[b,:]) # Compute C array (zero indices not used) @@ -56,7 +56,7 @@ def force ( r ): # Store D coefficients in a matrix # In the general case we would not need to calculate every pair # and also we would make use of the symmetry dd[a,b]=dd[b,a] - dd = np.zeros((n,n),dtype=np.float_) # Create D array + dd = np.zeros((n,n),dtype=np.float64) # Create D array for a in range(1,n): for b in range(1,n): dd[a,b]=cc[a,a]*cc[b,b] - cc[a,b]**2 # Compute D array (zero indices not used) diff --git a/python_examples/wl_hist.py b/python_examples/wl_hist.py index 143bcef..f8cca43 100755 --- a/python_examples/wl_hist.py +++ b/python_examples/wl_hist.py @@ -38,7 +38,7 @@ hist_file = input("Enter histogram file name: ") print('Reading histograms from '+hist_file) -e,h,g,s = np.loadtxt(hist_file,dtype=np.float_,unpack=True) +e,h,g,s = np.loadtxt(hist_file,dtype=np.float64,unpack=True) n = input("Enter number of atoms: ") n = int(n)