From 4455078a1b095a48d952aecd22a644b2885cc5a4 Mon Sep 17 00:00:00 2001 From: arm61 Date: Wed, 19 Dec 2018 14:05:12 +0000 Subject: [PATCH 1/4] add numba jit compliation to md things --- pylj/forcefields.py | 4 +++- pylj/md.py | 2 ++ pylj/pairwise.py | 5 ++++- pylj/util.py | 2 ++ requirements.txt | 1 + setup.py | 4 ++-- 6 files changed, 14 insertions(+), 4 deletions(-) diff --git a/pylj/forcefields.py b/pylj/forcefields.py index e1987b0..d3044ff 100644 --- a/pylj/forcefields.py +++ b/pylj/forcefields.py @@ -1,6 +1,7 @@ import numpy as np +from numba import jit - +@jit def lennard_jones(dr, constants, force=False): r"""Calculate the energy or force for a pair of particles using the Lennard-Jones (A/B variant) forcefield. @@ -34,6 +35,7 @@ def lennard_jones(dr, constants, force=False): np.power(dr, -6)) +@jit def buckingham(dr, constants, force=False): r""" Calculate the energy or force for a pair of particles using the Buckingham forcefield. diff --git a/pylj/md.py b/pylj/md.py index d35dd56..b33c274 100644 --- a/pylj/md.py +++ b/pylj/md.py @@ -1,6 +1,7 @@ import numpy as np from pylj import pairwise as heavy from pylj import forcefields as ff +from numba import jit def initialise(number_of_particles, temperature, box_length, init_conf, @@ -60,6 +61,7 @@ def initialize(number_particles, temperature, box_length, init_conf, return a +@jit def velocity_verlet(particles, timestep_length, box_length, cut_off, constants, forcefield, mass): """Uses the Velocity-Verlet integrator to move forward in time. The diff --git a/pylj/pairwise.py b/pylj/pairwise.py index 3db4d97..dd5f170 100644 --- a/pylj/pairwise.py +++ b/pylj/pairwise.py @@ -5,8 +5,10 @@ except ImportError: print("WARNING, using slow force and energy calculations") from pylj import pairwise as heavy + from numba import jit +@jit def compute_force(particles, box_length, cut_off, constants, forcefield, mass): r"""Calculates the forces and therefore the accelerations on each of the particles in the simulation. @@ -279,6 +281,7 @@ def heat_bath(particles, temperature_sample, bath_temp): return particles +@jit def dist(xposition, yposition, box_length): """Returns the distance array for the set of particles. Parameters @@ -317,7 +320,7 @@ def dist(xposition, yposition, box_length): k += 1 return drr, dxr, dyr - +@jit def pbc_correction(position, cell): """Correct for the periodic boundary condition. Parameters diff --git a/pylj/util.py b/pylj/util.py index 4318359..535a417 100644 --- a/pylj/util.py +++ b/pylj/util.py @@ -2,6 +2,7 @@ import numpy as np import webbrowser from pylj import md, mc +from numba import jit class System: @@ -148,6 +149,7 @@ def compute_energy(self): self.constants, self.forcefield) + @jit def integrate(self, method): """Maps the chosen integration method. Parameters diff --git a/requirements.txt b/requirements.txt index 4bfaa4a..617045e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,5 +2,6 @@ jupyter cython numpy matplotlib==2.2.3 +numba coverage coveralls diff --git a/setup.py b/setup.py index 5be3961..614f5ef 100644 --- a/setup.py +++ b/setup.py @@ -33,8 +33,8 @@ 'author_email': 'arm61@bath.ac.uk', 'packages': packages, 'include_package_data': True, - 'setup_requires': ['jupyter', 'numpy', 'matplotlib', 'cython'], - 'install_requires': ['jupyter', 'numpy', 'matplotlib', 'cython'], + 'setup_requires': ['jupyter', 'numpy', 'matplotlib', 'cython', 'numba'], + 'install_requires': ['jupyter', 'numpy', 'matplotlib', 'cython', 'numba'], 'version': VERSION, 'license': 'MIT', 'long_description': long_description, From baeac165441092560afe92784208cb82f7f89d26 Mon Sep 17 00:00:00 2001 From: arm61 Date: Wed, 19 Dec 2018 14:06:21 +0000 Subject: [PATCH 2/4] black formating --- pylj/forcefields.py | 19 ++-- pylj/mc.py | 43 +++++--- pylj/md.py | 193 ++++++++++++++++++--------------- pylj/pairwise.py | 87 ++++++++------- pylj/sample.py | 137 +++++++++++++---------- pylj/tests/test_comp.py | 2 +- pylj/tests/test_forcefields.py | 8 +- pylj/tests/test_mc.py | 28 ++--- pylj/tests/test_md.py | 103 +++++++++--------- pylj/tests/test_pairwise.py | 62 ++++++----- pylj/tests/test_util.py | 96 ++++++++++------ pylj/util.py | 148 +++++++++++++++---------- 12 files changed, 534 insertions(+), 392 deletions(-) diff --git a/pylj/forcefields.py b/pylj/forcefields.py index d3044ff..43f0abf 100644 --- a/pylj/forcefields.py +++ b/pylj/forcefields.py @@ -1,6 +1,7 @@ import numpy as np from numba import jit + @jit def lennard_jones(dr, constants, force=False): r"""Calculate the energy or force for a pair of particles using the @@ -28,11 +29,11 @@ def lennard_jones(dr, constants, force=False): The potential energy or force between the particles. """ if force: - return 12 * constants[0] * np.power(dr, -13) - (6 * constants[1] * - np.power(dr, -7)) + return 12 * constants[0] * np.power(dr, -13) - ( + 6 * constants[1] * np.power(dr, -7) + ) else: - return constants[0] * np.power(dr, -12) - (constants[1] * - np.power(dr, -6)) + return constants[0] * np.power(dr, -12) - (constants[1] * np.power(dr, -6)) @jit @@ -62,8 +63,10 @@ def buckingham(dr, constants, force=False): the potential energy or force between the particles. """ if force: - return constants[0] * constants[1] * np.exp(-constants[1] * dr) - \ - 6 * constants[2] / np.power(dr, 7) + return constants[0] * constants[1] * np.exp(-constants[1] * dr) - 6 * constants[ + 2 + ] / np.power(dr, 7) else: - return constants[0] * np.exp(-constants[1] * dr) \ - - constants[2] / np.power(dr, 6) + return constants[0] * np.exp(-constants[1] * dr) - constants[2] / np.power( + dr, 6 + ) diff --git a/pylj/mc.py b/pylj/mc.py index 7c44fd4..8dc41cc 100644 --- a/pylj/mc.py +++ b/pylj/mc.py @@ -2,9 +2,15 @@ from pylj import forcefields as ff -def initialise(number_of_particles, temperature, box_length, init_conf, - mass=39.948, constants=[1.363e-134, 9.273e-78], - forcefield=ff.lennard_jones): +def initialise( + number_of_particles, + temperature, + box_length, + init_conf, + mass=39.948, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, +): """Initialise the particle positions (this can be either as a square or random arrangement), velocities (based on the temperature defined, and # calculate the initial forces/accelerations. @@ -34,11 +40,18 @@ def initialise(number_of_particles, temperature, box_length, init_conf, System information. """ from pylj import util - system = util.System(number_of_particles, temperature, box_length, - constants, forcefield, mass, - init_conf=init_conf) - system.particles['xvelocity'] = 0 - system.particles['yvelocity'] = 0 + + system = util.System( + number_of_particles, + temperature, + box_length, + constants, + forcefield, + mass, + init_conf=init_conf, + ) + system.particles["xvelocity"] = 0 + system.particles["yvelocity"] = 0 return system @@ -87,8 +100,10 @@ def select_random_particle(particles): The current position of the chosen particle. """ random_particle = np.random.randint(0, particles.size) - position_store = [particles['xposition'][random_particle], - particles['yposition'][random_particle]] + position_store = [ + particles["xposition"][random_particle], + particles["yposition"][random_particle], + ] return random_particle, position_store @@ -110,8 +125,8 @@ def get_new_particle(particles, random_particle, box_length): Information about the particles, updated to account for the change of selected particle position. """ - particles['xposition'][random_particle] = np.random.uniform(0, box_length) - particles['yposition'][random_particle] = np.random.uniform(0, box_length) + particles["xposition"][random_particle] = np.random.uniform(0, box_length) + particles["yposition"][random_particle] = np.random.uniform(0, box_length) return particles @@ -149,8 +164,8 @@ def reject(position_store, particles, random_particle): Information about the particles, with the particle returned to the original position """ - particles['xposition'][random_particle] = position_store[0] - particles['yposition'][random_particle] = position_store[1] + particles["xposition"][random_particle] = position_store[0] + particles["yposition"][random_particle] = position_store[1] return particles diff --git a/pylj/md.py b/pylj/md.py index b33c274..df1279b 100644 --- a/pylj/md.py +++ b/pylj/md.py @@ -4,9 +4,16 @@ from numba import jit -def initialise(number_of_particles, temperature, box_length, init_conf, - timestep_length=1e-14, mass=39.948, - constants=[1.363e-134, 9.273e-78], forcefield=ff.lennard_jones): +def initialise( + number_of_particles, + temperature, + box_length, + init_conf, + timestep_length=1e-14, + mass=39.948, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, +): """Initialise the particle positions (this can be either as a square or random arrangement), velocities (based on the temperature defined, and calculate the initial forces/accelerations. @@ -38,32 +45,42 @@ def initialise(number_of_particles, temperature, box_length, init_conf, System information. """ from pylj import util - system = util.System(number_of_particles, temperature, box_length, - constants, forcefield, mass, - init_conf=init_conf, - timestep_length=timestep_length) + + system = util.System( + number_of_particles, + temperature, + box_length, + constants, + forcefield, + mass, + init_conf=init_conf, + timestep_length=timestep_length, + ) v = np.random.rand(system.particles.size, 2, 12) - v = np.sum(v, axis=2) - 6. + v = np.sum(v, axis=2) - 6.0 mass_kg = mass * 1.6605e-27 v = v * np.sqrt(1.3806e-23 * system.init_temp / mass_kg) v = v - np.average(v) - system.particles['xvelocity'] = v[:, 0] - system.particles['yvelocity'] = v[:, 1] + system.particles["xvelocity"] = v[:, 0] + system.particles["yvelocity"] = v[:, 1] return system -def initialize(number_particles, temperature, box_length, init_conf, - timestep_length=1e-14): +def initialize( + number_particles, temperature, box_length, init_conf, timestep_length=1e-14 +): """Maps to the md.initialise function to account for US english spelling. """ - a = initialise(number_particles, temperature, box_length, init_conf, - timestep_length) + a = initialise( + number_particles, temperature, box_length, init_conf, timestep_length + ) return a @jit -def velocity_verlet(particles, timestep_length, box_length, - cut_off, constants, forcefield, mass): +def velocity_verlet( + particles, timestep_length, box_length, cut_off, constants, forcefield, mass +): """Uses the Velocity-Verlet integrator to move forward in time. The Updates the particles positions and velocities in terms of the Velocity Verlet algorithm. Also calculates the instanteous temperature, pressure, @@ -81,30 +98,31 @@ def velocity_verlet(particles, timestep_length, box_length, util.particle_dt, array_like: Information about the particles, with new positions and velocities. """ - xposition_store = list(particles['xposition']) - yposition_store = list(particles['yposition']) - pos, prev_pos = update_positions([particles['xposition'], - particles['yposition']], - [particles['xprevious_position'], - particles['yprevious_position']], - [particles['xvelocity'], - particles['yvelocity']], - [particles['xacceleration'], - particles['yacceleration']], - timestep_length, box_length) - [particles['xposition'], particles['yposition']] = pos - [particles['xprevious_position'], particles['yprevious_position']] = pos - xacceleration_store = list(particles['xacceleration']) - yacceleration_store = list(particles['yacceleration']) + xposition_store = list(particles["xposition"]) + yposition_store = list(particles["yposition"]) + pos, prev_pos = update_positions( + [particles["xposition"], particles["yposition"]], + [particles["xprevious_position"], particles["yprevious_position"]], + [particles["xvelocity"], particles["yvelocity"]], + [particles["xacceleration"], particles["yacceleration"]], + timestep_length, + box_length, + ) + [particles["xposition"], particles["yposition"]] = pos + [particles["xprevious_position"], particles["yprevious_position"]] = pos + xacceleration_store = list(particles["xacceleration"]) + yacceleration_store = list(particles["yacceleration"]) particles, distances, forces, energies = heavy.compute_force( - particles, box_length, cut_off, constants, forcefield, mass) - [particles['xvelocity'], particles['yvelocity']] = update_velocities( - [particles['xvelocity'], particles['yvelocity']], + particles, box_length, cut_off, constants, forcefield, mass + ) + [particles["xvelocity"], particles["yvelocity"]] = update_velocities( + [particles["xvelocity"], particles["yvelocity"]], [xacceleration_store, yacceleration_store], - [particles['xacceleration'], particles['yacceleration']], - timestep_length) - particles['xprevious_position'] = xposition_store - particles['yprevious_position'] = yposition_store + [particles["xacceleration"], particles["yacceleration"]], + timestep_length, + ) + particles["xprevious_position"] = xposition_store + particles["yprevious_position"] = yposition_store return particles, distances, forces, energies @@ -128,16 +146,19 @@ def sample(particles, box_length, initial_particles, system): arrays. """ temperature_new = calculate_temperature(particles, system.mass) - system.temperature_sample = np.append(system.temperature_sample, - temperature_new) + system.temperature_sample = np.append(system.temperature_sample, temperature_new) pressure_new = heavy.calculate_pressure( - particles, box_length, temperature_new, system.cut_off, - system.constants, system.forcefield) + particles, + box_length, + temperature_new, + system.cut_off, + system.constants, + system.forcefield, + ) msd_new = calculate_msd(particles, initial_particles, box_length) system.pressure_sample = np.append(system.pressure_sample, pressure_new) system.force_sample = np.append(system.force_sample, np.sum(system.forces)) - system.energy_sample = np.append( - system.energy_sample, np.sum(system.energies)) + system.energy_sample = np.append(system.energy_sample, np.sum(system.energies)) system.msd_sample = np.append(system.msd_sample, msd_new) return system @@ -157,31 +178,32 @@ def calculate_msd(particles, initial_particles, box_length): float: Mean squared deviation for the particles at the given timestep. """ - xpos = np.array(particles['xposition']) - ypos = np.array(particles['yposition']) - dxinst = xpos - particles['xprevious_position'] - dyinst = ypos - particles['yprevious_position'] - for i in range(0, particles['xposition'].size): + xpos = np.array(particles["xposition"]) + ypos = np.array(particles["yposition"]) + dxinst = xpos - particles["xprevious_position"] + dyinst = ypos - particles["yprevious_position"] + for i in range(0, particles["xposition"].size): if np.abs(dxinst[i]) > 0.5 * box_length: if xpos[i] <= 0.5 * box_length: - particles['xpbccount'][i] += 1 + particles["xpbccount"][i] += 1 if xpos[i] > 0.5 * box_length: - particles['xpbccount'][i] -= 1 - xpos[i] += box_length * particles['xpbccount'][i] + particles["xpbccount"][i] -= 1 + xpos[i] += box_length * particles["xpbccount"][i] if np.abs(dyinst[i]) > 0.5 * box_length: if ypos[i] <= 0.5 * box_length: - particles['ypbccount'][i] += 1 + particles["ypbccount"][i] += 1 if ypos[i] > 0.5 * box_length: - particles['ypbccount'][i] -= 1 - ypos[i] += box_length * particles['ypbccount'][i] - dx = xpos - initial_particles['xposition'] - dy = ypos - initial_particles['yposition'] + particles["ypbccount"][i] -= 1 + ypos[i] += box_length * particles["ypbccount"][i] + dx = xpos - initial_particles["xposition"] + dy = ypos - initial_particles["yposition"] dr = np.sqrt(dx * dx + dy * dy) return np.average(dr ** 2) -def update_positions(positions, old_positions, velocities, - accelerations, timestep_length, box_length): +def update_positions( + positions, old_positions, velocities, accelerations, timestep_length, box_length +): """Update the particle positions using the Velocity-Verlet integrator. Parameters ---------- @@ -209,19 +231,20 @@ def update_positions(positions, old_positions, velocities, """ old_positions[0] = np.array(positions[0]) old_positions[1] = np.array(positions[1]) - positions[0] += ( - velocities[0] * timestep_length) + ( - 0.5 * accelerations[0] * timestep_length * timestep_length) - positions[1] += ( - velocities[1] * timestep_length) + ( - 0.5 * accelerations[1] * timestep_length * timestep_length) + positions[0] += (velocities[0] * timestep_length) + ( + 0.5 * accelerations[0] * timestep_length * timestep_length + ) + positions[1] += (velocities[1] * timestep_length) + ( + 0.5 * accelerations[1] * timestep_length * timestep_length + ) positions[0] = positions[0] % box_length positions[1] = positions[1] % box_length return [positions[0], positions[1]], [old_positions[0], old_positions[1]] -def update_velocities(velocities, accelerations_old, accelerations_new, - timestep_length): +def update_velocities( + velocities, accelerations_old, accelerations_new, timestep_length +): """Update the particle velocities using the Velocity-Verlet algoritm. Parameters ---------- @@ -239,10 +262,12 @@ def update_velocities(velocities, accelerations_old, accelerations_new, (2, N) array_like: Updated velocities. """ - velocities[0] += 0.5 * (accelerations_old[0] + - accelerations_new[0]) * timestep_length - velocities[1] += 0.5 * (accelerations_old[1] + - accelerations_new[1]) * timestep_length + velocities[0] += ( + 0.5 * (accelerations_old[0] + accelerations_new[0]) * timestep_length + ) + velocities[1] += ( + 0.5 * (accelerations_old[1] + accelerations_new[1]) * timestep_length + ) return [velocities[0], velocities[1]] @@ -260,8 +285,10 @@ def calculate_temperature(particles, mass): boltzmann_constant = 1.3806e-23 # joules/kelvin atomic_mass_unit = 1.660539e-27 # kilograms mass_kg = mass * atomic_mass_unit # kilograms - v = np.sqrt((particles['xvelocity'] * particles['xvelocity']) + - (particles['yvelocity'] * particles['yvelocity'])) + v = np.sqrt( + (particles["xvelocity"] * particles["xvelocity"]) + + (particles["yvelocity"] * particles["yvelocity"]) + ) k = 0.5 * np.sum(mass_kg * v * v) t = k / (particles.size * boltzmann_constant) return t @@ -297,12 +324,9 @@ def compute_force(particles, box_length, cut_off, constants, forcefield, mass): float, array_like Current energies between pairs of particles in the simulation. """ - part, dist, forces, energies = heavy.compute_force(particles, - box_length, - cut_off, - constants, - forcefield, - mass=mass) + part, dist, forces, energies = heavy.compute_force( + particles, box_length, cut_off, constants, forcefield, mass=mass + ) return part, dist, forces, energies @@ -331,11 +355,9 @@ def compute_energy(particles, box_length, cut_off, constants, forcefield): float, array_like Current energies between pairs of particles in the simulation. """ - dist, energies = heavy.compute_energy(particles, - box_length, - cut_off, - constants, - forcefield) + dist, energies = heavy.compute_energy( + particles, box_length, cut_off, constants, forcefield + ) return dist, energies @@ -359,6 +381,5 @@ def heat_bath(particles, temperature_sample, bath_temperature): util.particle_dt, array_like Information about the particles with new, rescaled velocities. """ - particles = heavy.heat_bath(particles, temperature_sample, - bath_temperature) + particles = heavy.heat_bath(particles, temperature_sample, bath_temperature) return particles diff --git a/pylj/pairwise.py b/pylj/pairwise.py index dd5f170..927945c 100644 --- a/pylj/pairwise.py +++ b/pylj/pairwise.py @@ -1,5 +1,6 @@ from __future__ import division import numpy as np + try: from pylj import comp as heavy except ImportError: @@ -39,24 +40,25 @@ def compute_force(particles, box_length, cut_off, constants, forcefield, mass): float, array_like Current energies between pairs of particles in the simulation. """ - particles['xacceleration'] = np.zeros(particles['xacceleration'].size) - particles['yacceleration'] = np.zeros(particles['yacceleration'].size) - pairs = int((particles['xacceleration'].size - 1) * - particles['xacceleration'].size / 2) + particles["xacceleration"] = np.zeros(particles["xacceleration"].size) + particles["yacceleration"] = np.zeros(particles["yacceleration"].size) + pairs = int( + (particles["xacceleration"].size - 1) * particles["xacceleration"].size / 2 + ) forces = np.zeros(pairs) distances = np.zeros(pairs) energies = np.zeros(pairs) atomic_mass_unit = 1.660539e-27 # kilograms mass_amu = mass # amu mass_kg = mass_amu * atomic_mass_unit # kilograms - distances, dx, dy = heavy.dist(particles['xposition'], - particles['yposition'], box_length) + distances, dx, dy = heavy.dist( + particles["xposition"], particles["yposition"], box_length + ) forces = forcefield(distances, constants, force=True) energies = forcefield(distances, constants, force=False) - forces[np.where(distances > cut_off)] = 0. - energies[np.where(distances > cut_off)] = 0. - particles = update_accelerations(particles, forces, mass_kg, dx, dy, - distances) + forces[np.where(distances > cut_off)] = 0.0 + energies[np.where(distances > cut_off)] = 0.0 + particles = update_accelerations(particles, forces, mass_kg, dx, dy, distances) return particles, distances, forces, energies @@ -97,12 +99,12 @@ def update_accelerations(particles, f, m, dx, dy, dr): Information about the particles with updated accelerations. """ k = 0 - for i in range(0, particles.size-1): + for i in range(0, particles.size - 1): for j in range(i + 1, particles.size): - particles['xacceleration'][i] += second_law(f[k], m, dx[k], dr[k]) - particles['yacceleration'][i] += second_law(f[k], m, dy[k], dr[k]) - particles['xacceleration'][j] -= second_law(f[k], m, dx[k], dr[k]) - particles['yacceleration'][j] -= second_law(f[k], m, dy[k], dr[k]) + particles["xacceleration"][i] += second_law(f[k], m, dx[k], dr[k]) + particles["yacceleration"][i] += second_law(f[k], m, dy[k], dr[k]) + particles["xacceleration"][j] -= second_law(f[k], m, dx[k], dr[k]) + particles["yacceleration"][j] -= second_law(f[k], m, dy[k], dr[k]) k += 1 return particles @@ -127,6 +129,7 @@ def second_law(f, m, d1, d2): """ return (f * d1 / d2) / m + def lennard_jones_energy(A, B, dr): """pairwise.lennard_jones_energy has been deprecated, please use forcefields.lennard_jones instead @@ -147,8 +150,10 @@ def lennard_jones_energy(A, B, dr): float: The potential energy between the two particles. """ - print("pairwise.lennard_jones_energy has been deprecated, please use " - "forcefields.lennard_jones instead") + print( + "pairwise.lennard_jones_energy has been deprecated, please use " + "forcefields.lennard_jones instead" + ) return A * np.power(dr, -12) - B * np.power(dr, -6) @@ -172,8 +177,10 @@ def lennard_jones_force(A, B, dr): float: The force between the two particles. """ - print("pairwise.lennard_jones_energy has been deprecated, please use " - "forcefields.lennard_jones with force=True instead") + print( + "pairwise.lennard_jones_energy has been deprecated, please use " + "forcefields.lennard_jones with force=True instead" + ) return 12 * A * np.power(dr, -13) - 6 * B * np.power(dr, -7) @@ -202,19 +209,22 @@ def compute_energy(particles, box_length, cut_off, constants, forcefield): float, array_like Current energies between pairs of particles in the simulation. """ - pairs = int((particles['xacceleration'].size - 1) * - particles['xacceleration'].size / 2) + pairs = int( + (particles["xacceleration"].size - 1) * particles["xacceleration"].size / 2 + ) distances = np.zeros(pairs) energies = np.zeros(pairs) - distances, dx, dy = heavy.dist(particles['xposition'], - particles['yposition'], box_length) + distances, dx, dy = heavy.dist( + particles["xposition"], particles["yposition"], box_length + ) energies = forcefield(distances, constants, force=False) - energies[np.where(distances > cut_off)] = 0. + energies[np.where(distances > cut_off)] = 0.0 return distances, energies -def calculate_pressure(particles, box_length, temperature, - cut_off, constants, forcefield): +def calculate_pressure( + particles, box_length, temperature, cut_off, constants, forcefield +): r"""Calculates the instantaneous pressure of the simulation cell, found with the following relationship: .. math:: @@ -241,15 +251,19 @@ def calculate_pressure(particles, box_length, temperature, float: Instantaneous pressure of the simulation. """ - distances, dx, dy = heavy.dist(particles['xposition'], - particles['yposition'], box_length) + distances, dx, dy = heavy.dist( + particles["xposition"], particles["yposition"], box_length + ) forces = forcefield(distances, constants, force=True) - forces[np.where(distances > cut_off)] = 0. + forces[np.where(distances > cut_off)] = 0.0 pres = np.sum(forces * distances) boltzmann_constant = 1.3806e-23 # joules / kelvin - pres = 1. / (2 * box_length * box_length) * pres + ( - particles['xposition'].size / (box_length * box_length) * - boltzmann_constant * temperature) + pres = 1.0 / (2 * box_length * box_length) * pres + ( + particles["xposition"].size + / (box_length * box_length) + * boltzmann_constant + * temperature + ) return pres @@ -274,10 +288,8 @@ def heat_bath(particles, temperature_sample, bath_temp): Information about the particles with new, rescaled velocities. """ average_temp = np.average(temperature_sample) - particles['xvelocity'] = particles['xvelocity'] * np.sqrt(bath_temp / - average_temp) - particles['yvelocity'] = particles['yvelocity'] * np.sqrt(bath_temp / - average_temp) + particles["xvelocity"] = particles["xvelocity"] * np.sqrt(bath_temp / average_temp) + particles["yvelocity"] = particles["yvelocity"] * np.sqrt(bath_temp / average_temp) return particles @@ -308,7 +320,7 @@ def dist(xposition, yposition, box_length): dyr = np.zeros(int((xposition.size - 1) * xposition.size / 2)) k = 0 for i in range(0, xposition.size - 1): - for j in range(i+1, xposition.size): + for j in range(i + 1, xposition.size): dx = xposition[i] - xposition[j] dy = yposition[i] - yposition[j] dx = pbc_correction(dx, box_length) @@ -320,6 +332,7 @@ def dist(xposition, yposition, box_length): k += 1 return drr, dxr, dyr + @jit def pbc_correction(position, cell): """Correct for the periodic boundary condition. diff --git a/pylj/sample.py b/pylj/sample.py index 7361052..b7861e5 100644 --- a/pylj/sample.py +++ b/pylj/sample.py @@ -13,14 +13,17 @@ class Scattering(object): # pragma: no cover system: System The whole system information. """ + def __init__(self, system): fig, ax = environment(4) self.average_rdf = [] self.r = [] self.average_diff = [] self.q = [] - self.initial_pos = [system.particles['xposition'], - system.particles['yposition']] + self.initial_pos = [ + system.particles["xposition"], + system.particles["yposition"], + ] setup_cellview(ax[0, 0], system) setup_rdfview(ax[0, 1], system) @@ -75,14 +78,17 @@ class Phase(object): # pragma: no cover system: System The whole system information. """ + def __init__(self, system): fig, ax = environment(4) self.average_rdf = [] self.r = [] self.average_diff = [] self.q = [] - self.initial_pos = [system.particles['xposition'], - system.particles['yposition']] + self.initial_pos = [ + system.particles["xposition"], + system.particles["yposition"], + ] setup_cellview(ax[0, 0], system) setup_rdfview(ax[1, 1], system) @@ -129,6 +135,7 @@ class Interactions(object): # pragma: no cover system: System The whole system information. """ + def __init__(self, system): fig, ax = environment(4) @@ -167,6 +174,7 @@ class JustCell(object): # pragma: no cover system: System The whole system information. """ + def __init__(self, system): fig, ax = environment(1) @@ -205,11 +213,12 @@ class CellPlus(object): # pragma: no cover ylabel: string The label for the y-axis of the custom plot. """ + def __init__(self, system, xlabel, ylabel): fig, ax = environment(2) setup_cellview(ax[0], system) - ax[1].plot([0], color='#34a5daff') + ax[1].plot([0], color="#34a5daff") ax[1].set_ylabel(ylabel, fontsize=16) ax[1].set_xlabel(xlabel, fontsize=16) @@ -249,6 +258,7 @@ class Energy(object): # pragma: no cover system: System The whole system information. """ + def __init__(self, system): fig, ax = environment(2) @@ -283,14 +293,17 @@ class RDF(object): # pragma: no cover system: System The whole system information. """ + def __init__(self, system): fig, ax = environment(2) self.average_rdf = [] self.r = [] self.average_diff = [] self.q = [] - self.initial_pos = [system.particles['xposition'], - system.particles['yposition']] + self.initial_pos = [ + system.particles["xposition"], + system.particles["yposition"], + ] setup_cellview(ax[0], system) setup_rdfview(ax[1], system) @@ -348,8 +361,7 @@ def environment(panes): # pragma: no cover elif panes == 4: fig, ax = plt.subplots(2, 2, figsize=(8, 8)) else: - AttributeError('The only options for the number of panes are 1, 2, or ' - '4') + AttributeError("The only options for the number of panes are 1, 2, or " "4") return fig, ax @@ -363,13 +375,12 @@ def setup_cellview(ax, system): # pragma: no cover system: System The whole system information. """ - xpos = system.particles['xposition'] - ypos = system.particles['yposition'] + xpos = system.particles["xposition"] + ypos = system.particles["yposition"] # These numbers are chosen to scale the size of the particles nicely to # allow the particles to appear to interact appropriately - mk = (6.00555e-8 / (system.box_length - 2.2727e-10) - 1e-10) - ax.plot(xpos, ypos, 'o', markersize=mk, markeredgecolor='black', - color='#34a5daff') + mk = 6.00555e-8 / (system.box_length - 2.2727e-10) - 1e-10 + ax.plot(xpos, ypos, "o", markersize=mk, markeredgecolor="black", color="#34a5daff") ax.set_xlim([0, system.box_length]) ax.set_ylim([0, system.box_length]) ax.set_xticks([]) @@ -384,9 +395,9 @@ def setup_forceview(ax): # pragma: no cover ax: Axes object The axes position that the pane should be placed in. """ - ax.plot([0], color='#34a5daff') - ax.set_ylabel('Force/N', fontsize=16) - ax.set_xlabel('Time/s', fontsize=16) + ax.plot([0], color="#34a5daff") + ax.set_ylabel("Force/N", fontsize=16) + ax.set_xlabel("Time/s", fontsize=16) def setup_energyview(ax): # pragma: no cover @@ -397,9 +408,9 @@ def setup_energyview(ax): # pragma: no cover ax: Axes object The axes position that the pane should be placed in. """ - ax.plot([0], color='#34a5daff') - ax.set_ylabel('Energy/J', fontsize=16) - ax.set_xlabel('Step', fontsize=16) + ax.plot([0], color="#34a5daff") + ax.set_ylabel("Energy/J", fontsize=16) + ax.set_xlabel("Step", fontsize=16) def setup_rdfview(ax, system): # pragma: no cover @@ -412,11 +423,11 @@ def setup_rdfview(ax, system): # pragma: no cover system: System The whole system information. """ - ax.plot([0], color='#34a5daff') + ax.plot([0], color="#34a5daff") ax.set_xlim([0, system.box_length / 2]) ax.set_yticks([]) - ax.set_ylabel('RDF', fontsize=16) - ax.set_xlabel('r/m', fontsize=16) + ax.set_ylabel("RDF", fontsize=16) + ax.set_xlabel("r/m", fontsize=16) def setup_diffview(ax): # pragma: no cover @@ -427,10 +438,10 @@ def setup_diffview(ax): # pragma: no cover ax: Axes object The axes position that the pane should be placed in. """ - ax.plot([0], color='#34a5daff') + ax.plot([0], color="#34a5daff") ax.set_yticks([]) - ax.set_ylabel('I(q)', fontsize=16) - ax.set_xlabel('q/m$^{-1}$', fontsize=16) + ax.set_ylabel("I(q)", fontsize=16) + ax.set_xlabel("q/m$^{-1}$", fontsize=16) def setup_pressureview(ax): # pragma: no cover @@ -441,9 +452,9 @@ def setup_pressureview(ax): # pragma: no cover ax: Axes object The axes position that the pane should be placed in. """ - ax.plot([0], color='#34a5daff') - ax.set_ylabel(r'Pressure/$\times10^6$Pa m$^{-1}$', fontsize=16) - ax.set_xlabel('Time/s', fontsize=16) + ax.plot([0], color="#34a5daff") + ax.set_ylabel(r"Pressure/$\times10^6$Pa m$^{-1}$", fontsize=16) + ax.set_xlabel("Time/s", fontsize=16) def setup_tempview(ax): # pragma: no cover @@ -454,9 +465,9 @@ def setup_tempview(ax): # pragma: no cover ax: Axes object The axes position that the pane should be placed in. """ - ax.plot([0], color='#34a5daff') - ax.set_ylabel('Temperature/K', fontsize=16) - ax.set_xlabel('Time/s', fontsize=16) + ax.plot([0], color="#34a5daff") + ax.set_ylabel("Temperature/K", fontsize=16) + ax.set_xlabel("Time/s", fontsize=16) def update_cellview(ax, system): # pragma: no cover @@ -469,8 +480,8 @@ def update_cellview(ax, system): # pragma: no cover system: System The whole system information. """ - x3 = system.particles['xposition'] - y3 = system.particles['yposition'] + x3 = system.particles["xposition"] + y3 = system.particles["yposition"] line = ax.lines[0] line.set_ydata(y3) line.set_xdata(x3) @@ -493,11 +504,15 @@ def update_rdfview(ax, system, average_rdf, r): # pragma: no cover averaged. """ hist, bin_edges = np.histogram( - system.distances, - bins=np.linspace(0, system.box_length / 2 + 0.5e-10, 100)) - gr = hist / (system.number_of_particles * (system.number_of_particles / - system.box_length ** 2) * - np.pi * (bin_edges[:-1] + 0.5e-10 / 2.) * 0.5) + system.distances, bins=np.linspace(0, system.box_length / 2 + 0.5e-10, 100) + ) + gr = hist / ( + system.number_of_particles + * (system.number_of_particles / system.box_length ** 2) + * np.pi + * (bin_edges[:-1] + 0.5e-10 / 2.0) + * 0.5 + ) average_rdf.append(gr) x = bin_edges[:-1] + 0.5e-10 / 2 r.append(x) @@ -528,8 +543,9 @@ def update_diffview(ax, system, average_diff, q): # pragma: no cover qw = np.linspace(2 * np.pi / system.box_length, 10e10, 1000)[20:] i = np.zeros_like(qw) for j in range(0, len(qw)): - i[j] = np.sum(3.644 * (np.sin(qw[j] * system.distances)) / ( - qw[j] * system.distances)) + i[j] = np.sum( + 3.644 * (np.sin(qw[j] * system.distances)) / (qw[j] * system.distances) + ) if i[j] < 0: i[j] = 0 x2 = qw @@ -557,10 +573,10 @@ def update_forceview(ax, system): # pragma: no cover line.set_ydata(system.force_sample) line.set_xdata(np.arange(0, system.step) * system.timestep_length) ax.set_xlim(0, system.step * system.timestep_length) - ax.set_ylim(np.amin(system.force_sample)-np.amax(system.force_sample) * - 0.05, - np.amax(system.force_sample)+np.amax(system.force_sample) * - 0.05) + ax.set_ylim( + np.amin(system.force_sample) - np.amax(system.force_sample) * 0.05, + np.amax(system.force_sample) + np.amax(system.force_sample) * 0.05, + ) def update_energyview(ax, system): # pragma: no cover @@ -578,14 +594,15 @@ def update_energyview(ax, system): # pragma: no cover y = system.energy_sample + 1.3806e-23 * system.temperature_sample line.set_xdata(np.arange(0, system.step) * system.timestep_length) ax.set_xlim(0, system.step * system.timestep_length) - ax.set_xlabel('Time/s', fontsize=16) + ax.set_xlabel("Time/s", fontsize=16) else: y = system.energy_sample - line.set_xdata(np.arange(0, system.step+1)) + line.set_xdata(np.arange(0, system.step + 1)) ax.set_xlim(0, system.step) line.set_ydata(y) - ax.set_ylim(np.amin(y) - np.abs(np.amax(y)) * 0.05, np.amax( - y) + np.abs(np.amax(y)) * 0.05) + ax.set_ylim( + np.amin(y) - np.abs(np.amax(y)) * 0.05, np.amax(y) + np.abs(np.amax(y)) * 0.05 + ) def update_tempview(ax, system): # pragma: no cover @@ -602,10 +619,10 @@ def update_tempview(ax, system): # pragma: no cover line.set_ydata(system.temperature_sample) line.set_xdata(np.arange(0, system.step) * system.timestep_length) ax.set_xlim(0, system.step * system.timestep_length) - ax.set_ylim(np.amin(system.temperature_sample)-np.amax( - system.temperature_sample) * 0.05, - np.amax(system.temperature_sample)+np.amax( - system.temperature_sample) * 0.05) + ax.set_ylim( + np.amin(system.temperature_sample) - np.amax(system.temperature_sample) * 0.05, + np.amax(system.temperature_sample) + np.amax(system.temperature_sample) * 0.05, + ) def update_pressureview(ax, system): # pragma: no cover @@ -625,8 +642,9 @@ def update_pressureview(ax, system): # pragma: no cover line.set_ydata(data) line.set_xdata(np.arange(0, system.step) * system.timestep_length) ax.set_xlim(0, system.step * system.timestep_length) - ax.set_ylim(np.amin(data) - np.amax(data) * 0.05, - np.amax(data) + np.amax(data) * 0.05) + ax.set_ylim( + np.amin(data) - np.amax(data) * 0.05, np.amax(data) + np.amax(data) * 0.05 + ) def setup_msdview(ax): # pragma: no cover @@ -637,9 +655,9 @@ def setup_msdview(ax): # pragma: no cover ax: Axes object The axes position that the pane should be placed in. """ - ax.plot([0], color='#34a5daff') - ax.set_ylabel('MSD/m$^2$', fontsize=16) - ax.set_xlabel('Time/s', fontsize=16) + ax.plot([0], color="#34a5daff") + ax.set_ylabel("MSD/m$^2$", fontsize=16) + ax.set_xlabel("Time/s", fontsize=16) def update_msdview(ax, system): # pragma: no cover @@ -657,5 +675,4 @@ def update_msdview(ax, system): # pragma: no cover line.set_ydata(system.msd_sample) line.set_xdata(np.arange(0, system.step) * system.timestep_length) ax.set_xlim(0, system.step * system.timestep_length) - ax.set_ylim(0, np.amax(system.msd_sample)+np.amax(system.msd_sample) * - 0.05) + ax.set_ylim(0, np.amax(system.msd_sample) + np.amax(system.msd_sample) * 0.05) diff --git a/pylj/tests/test_comp.py b/pylj/tests/test_comp.py index 99460dc..9f811bb 100644 --- a/pylj/tests/test_comp.py +++ b/pylj/tests/test_comp.py @@ -8,7 +8,7 @@ class TestPairwise(unittest.TestCase): def test_dist(self): xpos = np.array([0, 1]) ypos = np.array([0, 1]) - box_length = 5. + box_length = 5.0 dr, dx, dy = comp.dist(xpos, ypos, box_length) assert_almost_equal(dr, [np.sqrt(2)]) assert_almost_equal(dx, [-1]) diff --git a/pylj/tests/test_forcefields.py b/pylj/tests/test_forcefields.py index f101fd5..32321b6 100644 --- a/pylj/tests/test_forcefields.py +++ b/pylj/tests/test_forcefields.py @@ -5,17 +5,17 @@ class TestForcefields(unittest.TestCase): def test_lennard_jones_energy(self): - a = forcefields.lennard_jones(2., [1., 1.]) + a = forcefields.lennard_jones(2.0, [1.0, 1.0]) assert_almost_equal(a, -0.015380859) def test_lennard_jones_force(self): - a = forcefields.lennard_jones(2., [1., 1.], force=True) + a = forcefields.lennard_jones(2.0, [1.0, 1.0], force=True) assert_almost_equal(a, -0.045410156) def test_buckingham_energy(self): - a = forcefields.buckingham(2., [1., 1., 1.]) + a = forcefields.buckingham(2.0, [1.0, 1.0, 1.0]) assert_almost_equal(a, 0.1197103832) def test_buckingham_force(self): - a = forcefields.buckingham(2., [1., 1., 1.], force=True) + a = forcefields.buckingham(2.0, [1.0, 1.0, 1.0], force=True) assert_almost_equal(a, 0.08846028324) diff --git a/pylj/tests/test_mc.py b/pylj/tests/test_mc.py index 879127c..922d074 100644 --- a/pylj/tests/test_mc.py +++ b/pylj/tests/test_mc.py @@ -5,50 +5,50 @@ class TestMc(unittest.TestCase): def test_initialise_square(self): - a = mc.initialise(2, 300, 8, 'square') + a = mc.initialise(2, 300, 8, "square") assert_equal(a.number_of_particles, 2) assert_almost_equal(a.box_length, 8e-10) assert_almost_equal(a.init_temp, 300) - assert_almost_equal(a.particles['xposition']*1e10, [2, 2]) - assert_almost_equal(a.particles['yposition']*1e10, [2, 6]) + assert_almost_equal(a.particles["xposition"] * 1e10, [2, 2]) + assert_almost_equal(a.particles["yposition"] * 1e10, [2, 6]) def test_initialize_square(self): - a = mc.initialize(2, 300, 8, 'square') + a = mc.initialize(2, 300, 8, "square") assert_equal(a.number_of_particles, 2) assert_almost_equal(a.box_length, 8e-10) assert_almost_equal(a.init_temp, 300) - assert_almost_equal(a.particles['xposition']*1e10, [2, 2]) - assert_almost_equal(a.particles['yposition']*1e10, [2, 6]) + assert_almost_equal(a.particles["xposition"] * 1e10, [2, 2]) + assert_almost_equal(a.particles["yposition"] * 1e10, [2, 6]) def test_sample(self): - a = mc.initialise(2, 300, 8, 'square') + a = mc.initialise(2, 300, 8, "square") a = mc.sample(300, a) assert_almost_equal(a.energy_sample, [300]) def test_select_random_particle(self): - a = mc.initialise(2, 300, 8, 'square') + a = mc.initialise(2, 300, 8, "square") b, c = mc.select_random_particle(a.particles) self.assertTrue(0 <= b < 2) self.assertTrue(0 <= c[0] <= 8e-10) self.assertTrue(0 <= c[1] <= 8e-10) def test_get_new_particle(self): - a = mc.initialise(2, 300, 8, 'square') + a = mc.initialise(2, 300, 8, "square") b, c = mc.select_random_particle(a.particles) d = mc.get_new_particle(a.particles, b, a.box_length) - self.assertTrue(0 <= d['xposition'][b] <= 8e-10) - self.assertTrue(0 <= d['yposition'][b] <= 8e-10) + self.assertTrue(0 <= d["xposition"][b] <= 8e-10) + self.assertTrue(0 <= d["yposition"][b] <= 8e-10) def test_accept(self): a = mc.accept(300) assert_almost_equal(a, 300) def test_reject(self): - a = mc.initialise(2, 300, 8, 'square') + a = mc.initialise(2, 300, 8, "square") b = [1e-10, 1e-10] c = mc.reject(b, a.particles, 1) - assert_almost_equal(c['xposition'][1]*1e10, 1) - assert_almost_equal(c['yposition'][1]*1e10, 1) + assert_almost_equal(c["xposition"][1] * 1e10, 1) + assert_almost_equal(c["yposition"][1] * 1e10, 1) def test_metropolis_energy_reduce(self): a = mc.metropolis(300, 100, 1) diff --git a/pylj/tests/test_md.py b/pylj/tests/test_md.py index caa5670..663609b 100644 --- a/pylj/tests/test_md.py +++ b/pylj/tests/test_md.py @@ -5,79 +5,78 @@ class TestMd(unittest.TestCase): def test_initialise_square(self): - a = md.initialise(2, 300, 8, 'square') + a = md.initialise(2, 300, 8, "square") assert_equal(a.number_of_particles, 2) assert_almost_equal(a.box_length, 8e-10) assert_almost_equal(a.init_temp, 300) - assert_almost_equal(a.particles['xposition']*1e10, [2, 2]) - assert_almost_equal(a.particles['yposition']*1e10, [2, 6]) + assert_almost_equal(a.particles["xposition"] * 1e10, [2, 2]) + assert_almost_equal(a.particles["yposition"] * 1e10, [2, 6]) def test_velocity_verlet(self): - a = md.initialise(2, 300, 8, 'square') + a = md.initialise(2, 300, 8, "square") a.particles, a.distances, a.forces, a.energies = md.velocity_verlet( - a.particles, 1, a.box_length, a.cut_off, a.constants, a. - forcefield, a.mass) - assert_almost_equal(a.particles['xprevious_position']*1e10, [2, 2]) - assert_almost_equal(a.particles['yprevious_position']*1e10, [2, 6]) + a.particles, 1, a.box_length, a.cut_off, a.constants, a.forcefield, a.mass + ) + assert_almost_equal(a.particles["xprevious_position"] * 1e10, [2, 2]) + assert_almost_equal(a.particles["yprevious_position"] * 1e10, [2, 6]) def test_update_positions(self): - a = md.initialise(2, 300, 8, 'square') - a.particles['xvelocity'] = 1e4 - a.particles['yvelocity'] = 1e4 - a.particles['xacceleration'] = 1e4 - a.particles['yacceleration'] = 1e4 - b, c = md.update_positions([a.particles['xposition'], - a.particles['yposition']], - [a.particles['xprevious_position'], - a.particles['yprevious_position']], - [a.particles['xvelocity'], - a.particles['yvelocity']], - [a.particles['xacceleration'], - a.particles['yacceleration']], - a.timestep_length, a.box_length) - assert_almost_equal(b[0][0]*1e10, 3) - assert_almost_equal(b[1][0]*1e10, 3) - assert_almost_equal(b[0][1]*1e10, 3) - assert_almost_equal(b[1][1]*1e10, 7) + a = md.initialise(2, 300, 8, "square") + a.particles["xvelocity"] = 1e4 + a.particles["yvelocity"] = 1e4 + a.particles["xacceleration"] = 1e4 + a.particles["yacceleration"] = 1e4 + b, c = md.update_positions( + [a.particles["xposition"], a.particles["yposition"]], + [a.particles["xprevious_position"], a.particles["yprevious_position"]], + [a.particles["xvelocity"], a.particles["yvelocity"]], + [a.particles["xacceleration"], a.particles["yacceleration"]], + a.timestep_length, + a.box_length, + ) + assert_almost_equal(b[0][0] * 1e10, 3) + assert_almost_equal(b[1][0] * 1e10, 3) + assert_almost_equal(b[0][1] * 1e10, 3) + assert_almost_equal(b[1][1] * 1e10, 7) def test_update_velocities(self): - a = md.initialise(2, 300, 8, 'square') - a.particles['xvelocity'] = 1e-10 - a.particles['yvelocity'] = 1e-10 - a.particles['xacceleration'] = 1e4 - a.particles['yacceleration'] = 1e4 + a = md.initialise(2, 300, 8, "square") + a.particles["xvelocity"] = 1e-10 + a.particles["yvelocity"] = 1e-10 + a.particles["xacceleration"] = 1e4 + a.particles["yacceleration"] = 1e4 xacceleration_new = 2e4 yacceleration_new = 2e4 - b = md.update_velocities([a.particles['xvelocity'], - a.particles['yvelocity']], - [xacceleration_new, yacceleration_new], - [a.particles['xacceleration'], - a.particles['yacceleration']], - a.timestep_length) - assert_almost_equal(b[0][0]*1e10, 2.5) - assert_almost_equal(b[1][0]*1e10, 2.5) - assert_almost_equal(b[0][1]*1e10, 2.5) - assert_almost_equal(b[1][1]*1e10, 2.5) + b = md.update_velocities( + [a.particles["xvelocity"], a.particles["yvelocity"]], + [xacceleration_new, yacceleration_new], + [a.particles["xacceleration"], a.particles["yacceleration"]], + a.timestep_length, + ) + assert_almost_equal(b[0][0] * 1e10, 2.5) + assert_almost_equal(b[1][0] * 1e10, 2.5) + assert_almost_equal(b[0][1] * 1e10, 2.5) + assert_almost_equal(b[1][1] * 1e10, 2.5) def test_calculate_temperature(self): - a = md.initialise(1, 300, 8, 'square') - a.particles['xvelocity'] = [1e-10] - a.particles['yvelocity'] = [1e-10] - a.particles['xacceleration'] = [1e4] - a.particles['yacceleration'] = [1e4] + a = md.initialise(1, 300, 8, "square") + a.particles["xvelocity"] = [1e-10] + a.particles["yvelocity"] = [1e-10] + a.particles["xacceleration"] = [1e4] + a.particles["yacceleration"] = [1e4] b = md.calculate_temperature(a.particles, mass=39.948) assert_almost_equal(b * 1e23, 4.8048103702737945) def test_calculate_msd(self): - a = md.initialise(2, 300, 8, 'square') - a.particles['xposition'] = [3e-10, 3e-10] - a.particles['yposition'] = [3e-10, 7e-10] + a = md.initialise(2, 300, 8, "square") + a.particles["xposition"] = [3e-10, 3e-10] + a.particles["yposition"] = [3e-10, 7e-10] b = md.calculate_msd(a.particles, a.initial_particles, a.box_length) assert_almost_equal(b, 2e-20) def test_calculate_msd_large(self): - a = md.initialise(2, 300, 8, 'square') - a.particles['xposition'] = [7e-10, 3e-10] - a.particles['yposition'] = [7e-10, 7e-10] + a = md.initialise(2, 300, 8, "square") + a.particles["xposition"] = [7e-10, 3e-10] + a.particles["yposition"] = [7e-10, 7e-10] b = md.calculate_msd(a.particles, a.initial_particles, a.box_length) assert_almost_equal(b, 10e-20) diff --git a/pylj/tests/test_pairwise.py b/pylj/tests/test_pairwise.py index 45d274b..8efbe07 100644 --- a/pylj/tests/test_pairwise.py +++ b/pylj/tests/test_pairwise.py @@ -6,18 +6,16 @@ class TestPairwise(unittest.TestCase): - def test_update_accelerations(self): part_dt = util.particle_dt() particles = np.zeros(2, dtype=part_dt) ones = np.array([1]) dist = np.array([np.sqrt(2)]) - particles = pairwise.update_accelerations(particles, ones, 1, ones, - ones, dist) - assert_almost_equal(particles['xacceleration'][0], 0.707106781) - assert_almost_equal(particles['yacceleration'][0], 0.707106781) - assert_almost_equal(particles['xacceleration'][1], -0.707106781) - assert_almost_equal(particles['yacceleration'][1], -0.707106781) + particles = pairwise.update_accelerations(particles, ones, 1, ones, ones, dist) + assert_almost_equal(particles["xacceleration"][0], 0.707106781) + assert_almost_equal(particles["yacceleration"][0], 0.707106781) + assert_almost_equal(particles["xacceleration"][1], -0.707106781) + assert_almost_equal(particles["yacceleration"][1], -0.707106781) def test_second_law(self): a = pairwise.second_law(1, 1, 1, np.sqrt(2)) @@ -30,38 +28,52 @@ def test_separation(self): def test_compute_forces(self): part_dt = util.particle_dt() particles = np.zeros(2, dtype=part_dt) - particles['xposition'][0] = 1e-10 - particles['xposition'][1] = 5e-10 + particles["xposition"][0] = 1e-10 + particles["xposition"][1] = 5e-10 particles, distances, forces, energies = pairwise.compute_force( - particles, 30, 15, constants=[1.363e-134, 9.273e-78], - forcefield=ff.lennard_jones, mass=39.948) + particles, + 30, + 15, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, + mass=39.948, + ) assert_almost_equal(distances, [4e-10]) assert_almost_equal(energies, [-1.4515047e-21]) assert_almost_equal(forces, [-9.5864009e-12]) - assert_almost_equal(particles['yacceleration'], [0, 0]) - assert_almost_equal(particles['xacceleration'][0]/1e14, 1.4451452) - assert_almost_equal(particles['xacceleration'][1]/1e14, -1.4451452) + assert_almost_equal(particles["yacceleration"], [0, 0]) + assert_almost_equal(particles["xacceleration"][0] / 1e14, 1.4451452) + assert_almost_equal(particles["xacceleration"][1] / 1e14, -1.4451452) def test_compute_energy(self): part_dt = util.particle_dt() particles = np.zeros(2, dtype=part_dt) - particles['xposition'][0] = 1e-10 - particles['xposition'][1] = 5e-10 - d, e = pairwise.compute_energy(particles, 30, 15, - constants=[1.363e-134, 9.273e-78], - forcefield=ff.lennard_jones) + particles["xposition"][0] = 1e-10 + particles["xposition"][1] = 5e-10 + d, e = pairwise.compute_energy( + particles, + 30, + 15, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, + ) assert_almost_equal(d, [4e-10]) assert_almost_equal(e, [-1.4515047e-21]) def test_calculate_pressure(self): part_dt = util.particle_dt() particles = np.zeros(2, dtype=part_dt) - particles['xposition'][0] = 1e-10 - particles['xposition'][1] = 5e-10 - p = pairwise.calculate_pressure(particles, 30, 300, 15, - constants=[1.363e-134, 9.273e-78], - forcefield=ff.lennard_jones) - assert_almost_equal(p*1e24, 7.07368869) + particles["xposition"][0] = 1e-10 + particles["xposition"][1] = 5e-10 + p = pairwise.calculate_pressure( + particles, + 30, + 300, + 15, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, + ) + assert_almost_equal(p * 1e24, 7.07368869) def test_pbc_correction(self): a = pairwise.pbc_correction(1, 10) diff --git a/pylj/tests/test_util.py b/pylj/tests/test_util.py index afbafb9..df8c474 100644 --- a/pylj/tests/test_util.py +++ b/pylj/tests/test_util.py @@ -6,34 +6,45 @@ class TestUtil(unittest.TestCase): def test_system_square(self): - a = util.System(2, 300, 8, mass=39.948, - constants=[1.363e-134, 9.273e-78], - forcefield=ff.lennard_jones) + a = util.System( + 2, + 300, + 8, + mass=39.948, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, + ) assert_equal(a.number_of_particles, 2) assert_equal(a.init_temp, 300) assert_almost_equal(a.box_length * 1e10, 8) assert_almost_equal(a.timestep_length, 1e-14) - assert_almost_equal(a.particles['xposition'] * 1e10, [2, 2]) - assert_almost_equal(a.particles['yposition'] * 1e10, [2, 6]) - assert_almost_equal(a.initial_particles['xposition'] * 1e10, [2, 2]) - assert_almost_equal(a.initial_particles['yposition'] * 1e10, [2, 6]) + assert_almost_equal(a.particles["xposition"] * 1e10, [2, 2]) + assert_almost_equal(a.particles["yposition"] * 1e10, [2, 6]) + assert_almost_equal(a.initial_particles["xposition"] * 1e10, [2, 2]) + assert_almost_equal(a.initial_particles["yposition"] * 1e10, [2, 6]) assert_almost_equal(a.cut_off * 1e10, 4.0) assert_equal(a.distances.size, 1) assert_equal(a.forces.size, 1) assert_equal(a.energies.size, 1) def test_system_random(self): - a = util.System(2, 300, 8, init_conf='random', mass=39.948, - constants=[1.363e-134, 9.273e-78], - forcefield=ff.lennard_jones) + a = util.System( + 2, + 300, + 8, + init_conf="random", + mass=39.948, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, + ) assert_equal(a.number_of_particles, 2) assert_equal(a.init_temp, 300) assert_almost_equal(a.box_length * 1e10, 8) assert_almost_equal(a.timestep_length, 1e-14) - self.assertTrue(0 <= a.particles['xposition'][0] * 1e10 <= 8) - self.assertTrue(0 <= a.particles['yposition'][0] * 1e10 <= 8) - self.assertTrue(0 <= a.particles['xposition'][1] * 1e10 <= 8) - self.assertTrue(0 <= a.particles['yposition'][1] * 1e10 <= 8) + self.assertTrue(0 <= a.particles["xposition"][0] * 1e10 <= 8) + self.assertTrue(0 <= a.particles["yposition"][0] * 1e10 <= 8) + self.assertTrue(0 <= a.particles["xposition"][1] * 1e10 <= 8) + self.assertTrue(0 <= a.particles["yposition"][1] * 1e10 <= 8) assert_almost_equal(a.cut_off * 1e10, 4.0) assert_equal(a.distances.size, 1) assert_equal(a.forces.size, 1) @@ -41,27 +52,48 @@ def test_system_random(self): def test_system_too_big(self): with self.assertRaises(AttributeError) as context: - util.System(2, 300, 1000, mass=39.948, - constants=[1.363e-134, 9.273e-78], - forcefield=ff.lennard_jones) - self.assertTrue('With a box length of 1000 the particles are probably ' - 'too small to be seen in the viewer. Try something ' - '(much) less than 600.' in str(context.exception)) + util.System( + 2, + 300, + 1000, + mass=39.948, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, + ) + self.assertTrue( + "With a box length of 1000 the particles are probably " + "too small to be seen in the viewer. Try something " + "(much) less than 600." in str(context.exception) + ) def test_system_too_small(self): with self.assertRaises(AttributeError) as context: - util.System(2, 300, 2, mass=39.948, - constants=[1.363e-134, 9.273e-78], - forcefield=ff.lennard_jones) - self.assertTrue('With a box length of 2 the cell is too small to ' - 'really hold more than one particle.' in str( - context.exception)) + util.System( + 2, + 300, + 2, + mass=39.948, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, + ) + self.assertTrue( + "With a box length of 2 the cell is too small to " + "really hold more than one particle." in str(context.exception) + ) def test_system_init_conf(self): with self.assertRaises(NotImplementedError) as context: - util.System(2, 300, 100, init_conf='horseradish', mass=39.948, - constants=[1.363e-134, 9.273e-78], - forcefield=ff.lennard_jones) - self.assertTrue('The initial configuration type horseradish is not ' - 'recognised. Available options are: square or ' - 'random' in str(context.exception)) + util.System( + 2, + 300, + 100, + init_conf="horseradish", + mass=39.948, + constants=[1.363e-134, 9.273e-78], + forcefield=ff.lennard_jones, + ) + self.assertTrue( + "The initial configuration type horseradish is not " + "recognised. Available options are: square or " + "random" in str(context.exception) + ) diff --git a/pylj/util.py b/pylj/util.py index 535a417..f4fffd5 100644 --- a/pylj/util.py +++ b/pylj/util.py @@ -35,9 +35,19 @@ class System: forcefield: function (optional) The particular forcefield to be used to find the energy and forces. """ - def __init__(self, number_of_particles, temperature, box_length, - constants, forcefield, mass, init_conf='square', - timestep_length=1e-14, cut_off=15): + + def __init__( + self, + number_of_particles, + temperature, + box_length, + constants, + forcefield, + mass, + init_conf="square", + timestep_length=1e-14, + cut_off=15, + ): self.number_of_particles = number_of_particles self.init_temp = temperature self.constants = constants @@ -46,32 +56,38 @@ def __init__(self, number_of_particles, temperature, box_length, if box_length <= 600: self.box_length = box_length * 1e-10 else: - raise AttributeError('With a box length of {} the particles are ' - 'probably too small to be seen in the ' - 'viewer. Try something (much) less than ' - '600.'.format(box_length)) + raise AttributeError( + "With a box length of {} the particles are " + "probably too small to be seen in the " + "viewer. Try something (much) less than " + "600.".format(box_length) + ) if box_length >= 4: self.box_length = box_length * 1e-10 else: - raise AttributeError('With a box length of {} the cell is too ' - 'small to really hold more than one ' - 'particle.'.format(box_length)) + raise AttributeError( + "With a box length of {} the cell is too " + "small to really hold more than one " + "particle.".format(box_length) + ) self.timestep_length = timestep_length self.particles = None - if init_conf == 'square': + if init_conf == "square": self.square() - elif init_conf == 'random': + elif init_conf == "random": self.random() else: - raise NotImplementedError('The initial configuration type {} is ' - 'not recognised. Available options are: ' - 'square or random'.format(init_conf)) + raise NotImplementedError( + "The initial configuration type {} is " + "not recognised. Available options are: " + "square or random".format(init_conf) + ) if box_length > 30: self.cut_off = cut_off * 1e-10 else: self.cut_off = box_length / 2 * 1e-10 self.step = 0 - self.time = 0. + self.time = 0.0 self.distances = np.zeros(self.number_of_pairs()) self.forces = np.zeros(self.number_of_pairs()) self.energies = np.zeros(self.number_of_pairs()) @@ -93,8 +109,7 @@ def number_of_pairs(self): int: Number of pairwise interactions in the system. """ - return int((self.number_of_particles - 1) * - self.number_of_particles / 2) + return int((self.number_of_particles - 1) * self.number_of_particles / 2) def square(self): """Sets the initial positions of the particles on a square lattice. @@ -107,8 +122,8 @@ def square(self): for i in range(0, m): for j in range(0, m): if n < self.number_of_particles: - self.particles[n]['xposition'] = (i + 0.5) * d - self.particles[n]['yposition'] = (j + 0.5) * d + self.particles[n]["xposition"] = (i + 0.5) * d + self.particles[n]["yposition"] = (j + 0.5) * d n += 1 def random(self): @@ -117,22 +132,21 @@ def random(self): part_dt = particle_dt() self.particles = np.zeros(self.number_of_particles, dtype=part_dt) num_part = self.number_of_particles - self.particles['xposition'] = np.random.uniform(0, self.box_length, - num_part) - self.particles['yposition'] = np.random.uniform(0, self.box_length, - num_part) + self.particles["xposition"] = np.random.uniform(0, self.box_length, num_part) + self.particles["yposition"] = np.random.uniform(0, self.box_length, num_part) def compute_force(self): """Maps to the compute_force function in either the comp (if Cython is installed) or the pairwise module and allows for a cleaner interface. """ - part, dist, forces, energies = md.compute_force(self.particles, - self.box_length, - self.cut_off, - self.constants, - self.forcefield, - self.mass - ) + part, dist, forces, energies = md.compute_force( + self.particles, + self.box_length, + self.cut_off, + self.constants, + self.forcefield, + self.mass, + ) self.particles = part self.distances = dist self.forces = forces @@ -143,11 +157,13 @@ def compute_energy(self): is installed) or the pairwise module and allows for a cleaner interface. """ - self.distances, self.energies = md.compute_energy(self.particles, - self.box_length, - self.cut_off, - self.constants, - self.forcefield) + self.distances, self.energies = md.compute_energy( + self.particles, + self.box_length, + self.cut_off, + self.constants, + self.forcefield, + ) @jit def integrate(self, method): @@ -158,14 +174,19 @@ def integrate(self, method): The integration method to be used, e.g. md.velocity_verlet. """ self.particles, self.distances, self.forces, self.energies = method( - self.particles, self.timestep_length, self.box_length, - self.cut_off, self.constants, self.forcefield, self.mass) + self.particles, + self.timestep_length, + self.box_length, + self.cut_off, + self.constants, + self.forcefield, + self.mass, + ) def md_sample(self): """Maps to the md.sample function. """ - self = md.sample(self.particles, self.box_length, - self.initial_particles, self) + self = md.sample(self.particles, self.box_length, self.initial_particles, self) def heat_bath(self, bath_temperature): """Maps to the heat_bath function in either the comp (if Cython is @@ -175,8 +196,9 @@ def heat_bath(self, bath_temperature): target_temperature: float The target temperature for the simulation. """ - self.particles = md.heat_bath(self.particles, self.temperature_sample, - bath_temperature) + self.particles = md.heat_bath( + self.particles, self.temperature_sample, bath_temperature + ) def mc_sample(self): """Maps to the mc.sample function. @@ -191,14 +213,15 @@ def select_random_particle(self): """Maps to the mc.select_random_particle function. """ self.random_particle, self.position_store = mc.select_random_particle( - self.particles) + self.particles + ) def new_random_position(self): """Maps to the mc.get_new_particle function. """ - self.particles = mc.get_new_particle(self.particles, - self.random_particle, - self.box_length) + self.particles = mc.get_new_particle( + self.particles, self.random_particle, self.box_length + ) def accept(self): """Maps to the mc.accept function. @@ -208,15 +231,15 @@ def accept(self): def reject(self): """Maps to the mc.reject function. """ - self.particles = mc.reject(self.position_store, self.particles, - self.random_particle) + self.particles = mc.reject( + self.position_store, self.particles, self.random_particle + ) def __cite__(): # pragma: no cover """This function will launch the website for the JOSE publication on pylj.""" - webbrowser.open( - 'http://jose.theoj.org/papers/58daa1a1a564dc8e0f99ffcdae20eb1d') + webbrowser.open("http://jose.theoj.org/papers/58daa1a1a564dc8e0f99ffcdae20eb1d") def __version__(): # pragma: no cover @@ -224,7 +247,7 @@ def __version__(): # pragma: no cover major = 1 minor = 2 micro = 0 - print('pylj-{:d}.{:d}.{:d}'.format(major, minor, micro)) + print("pylj-{:d}.{:d}.{:d}".format(major, minor, micro)) def particle_dt(): @@ -236,11 +259,18 @@ def particle_dt(): - xforce and yforce - energy """ - return np.dtype([('xposition', np.float64), ('yposition', np.float64), - ('xvelocity', np.float64), ('yvelocity', np.float64), - ('xacceleration', np.float64), - ('yacceleration', np.float64), - ('xprevious_position', np.float64), - ('yprevious_position', np.float64), - ('energy', np.float64), ('xpbccount', int), - ('ypbccount', int)]) + return np.dtype( + [ + ("xposition", np.float64), + ("yposition", np.float64), + ("xvelocity", np.float64), + ("yvelocity", np.float64), + ("xacceleration", np.float64), + ("yacceleration", np.float64), + ("xprevious_position", np.float64), + ("yprevious_position", np.float64), + ("energy", np.float64), + ("xpbccount", int), + ("ypbccount", int), + ] + ) From 0e927dbfcbffa4f2774742b068e32d5795703005 Mon Sep 17 00:00:00 2001 From: arm61 Date: Wed, 19 Dec 2018 14:07:09 +0000 Subject: [PATCH 3/4] version number bump --- docs/source/conf.py | 2 +- pylj/util.py | 2 +- setup.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 09a6c38..f7e99a2 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -60,7 +60,7 @@ # built documents. # # The short X.Y version. -version = '1.2.0' +version = '1.2.1' # The full version, including alpha/beta/rc tags. release = version diff --git a/pylj/util.py b/pylj/util.py index f4fffd5..707d60f 100644 --- a/pylj/util.py +++ b/pylj/util.py @@ -246,7 +246,7 @@ def __version__(): # pragma: no cover """This will print the number of the pylj version currently in use.""" major = 1 minor = 2 - micro = 0 + micro = 1 print("pylj-{:d}.{:d}.{:d}".format(major, minor, micro)) diff --git a/setup.py b/setup.py index 614f5ef..5eb87a5 100644 --- a/setup.py +++ b/setup.py @@ -17,7 +17,7 @@ # versioning MAJOR = 1 MINOR = 2 -MICRO = 0 +MICRO = 1 ISRELEASED = True VERSION = '%d.%d.%d' % (MAJOR, MINOR, MICRO) From db3d102ac81a4f921586af8472f5741fda1fa699 Mon Sep 17 00:00:00 2001 From: arm61 Date: Wed, 19 Dec 2018 14:17:41 +0000 Subject: [PATCH 4/4] move jit import out of exception --- pylj/pairwise.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pylj/pairwise.py b/pylj/pairwise.py index 927945c..26783af 100644 --- a/pylj/pairwise.py +++ b/pylj/pairwise.py @@ -1,12 +1,11 @@ from __future__ import division import numpy as np - +from numba import jit try: from pylj import comp as heavy except ImportError: print("WARNING, using slow force and energy calculations") from pylj import pairwise as heavy - from numba import jit @jit