Merge fc9dbf1 into 1c7c971

pylj.egg-info/SOURCES.txt

@@ -5,6 +5,7 @@ pylj/__init__.py
 pylj/_testutils.py
 pylj/mc.py
 pylj/md.py
+pylj/pairwise.py
 pylj/sample.py
 pylj/util.py
 pylj.egg-info/PKG-INFO

pylj/mc.py

@@ -1,5 +1,4 @@
 import numpy as np
-from pylj import util
 
 def initialise(number_of_particles, temperature, box_length, init_conf):
     """Initialise the particle positions (this can be either as a square or random arrangement), velocities (based on
@@ -23,6 +22,7 @@ def initialise(number_of_particles, temperature, box_length, init_conf):
     System
         System information.
     """
+    from pylj import util
     system = util.System(number_of_particles, temperature, box_length, init_conf=init_conf)
     system.particles['xvelocity'] = 0
     system.particles['yvelocity'] = 0

pylj/md.py

@@ -1,5 +1,9 @@
 import numpy as np
-from pylj import comp, util
+try:
+    from pylj import comp as heavy
+except ImportError:
+    print("WARNING, using slow force and energy calculations")
+    from pylj import pairwise as heavy
 
 
 def initialise(number_of_particles, temperature, box_length, init_conf, timestep_length=1e-14):
@@ -26,14 +30,15 @@ def initialise(number_of_particles, temperature, box_length, init_conf, timestep
     System
         System information.
     """
+    from pylj import util
     system = util.System(number_of_particles, temperature, box_length, init_conf=init_conf,
                          timestep_length=timestep_length)
     v = np.random.rand(system.particles.size, 2) - 0.5
     v2sum = np.average(np.square(v))
     v = (v - np.average(v)) * np.sqrt(2 * system.init_temp / v2sum)
     system.particles['xvelocity'] = v[:, 0]
     system.particles['yvelocity'] = v[:, 1]
-    system.particles, system.distances, system.forces, system.energies = comp.compute_forces(system.particles,
+    system.particles, system.distances, system.forces, system.energies = heavy.compute_forces(system.particles,
                                                                                              system.box_length,
                                                                                              system.cut_off)
     return system
@@ -70,7 +75,7 @@ def velocity_verlet(particles, timestep_length, box_length, cut_off):
                                                                         box_length)
     xacceleration_store = list(particles['xacceleration'])
     yacceleration_store = list(particles['yacceleration'])
-    particles, distances, forces, energies = comp.compute_forces(particles, box_length, cut_off)
+    particles, distances, forces, energies = heavy.compute_forces(particles, box_length, cut_off)
     [particles['xvelocity'], particles['yvelocity']] = update_velocities([particles['xvelocity'],
                                                                           particles['yvelocity']],
                                                                          [xacceleration_store, yacceleration_store],
@@ -101,17 +106,45 @@ def sample(particles, box_length, initial_particles, system):
         Details about the whole system, with the new temperature, pressure, msd, and force appended to the appropriate
         arrays.
     """
-    temperature_new = util.calculate_temperature(particles)
+    temperature_new = calculate_temperature(particles)
     system.temperature_sample = np.append(system.temperature_sample, temperature_new)
-    pressure_new = comp.calculate_pressure(particles, box_length, temperature_new, system.cut_off)
-    msd_new = util.calculate_msd(particles, initial_particles, box_length)
+    pressure_new = heavy.calculate_pressure(particles, box_length, temperature_new, system.cut_off)
+    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.msd_sample = np.append(system.msd_sample, msd_new)
     return system
 
 
+def calculate_msd(particles, initial_particles, box_length):
+    """Determines the mean squared displacement of the particles in the system.
+
+    Parameters
+    ----------
+    particles: util.particle_dt, array_like
+        Information about the particles.
+    initial_particles: util.particle_dt, array_like
+        Information about the initial state of the particles.
+    box_length: float
+        Size of the cell vector.
+
+    Returns
+    -------
+    float:
+        Mean squared deviation for the particles at the given timestep.
+    """
+    dx = particles['xposition'] - initial_particles['xposition']
+    dy = particles['yposition'] - initial_particles['yposition']
+    for i in range(0, particles['xposition'].size):
+        if np.abs(dx[i]) > 0.5 * box_length:
+            dx[i] *= 1 - box_length / np.abs(dx[i])
+        if np.abs(dy[i]) > 0.5 * box_length:
+            dy[i] *= 1 - box_length / np.abs(dy[i])
+    dr = np.sqrt(dx * dx + dy * dy)
+    return np.average(dr ** 2)
+
+
 def update_positions(positions, velocities, accelerations, timestep_length, box_length):
     """Update the particle positions using the Velocity-Verlet integrator.
 
@@ -162,3 +195,29 @@ def update_velocities(velocities, accelerations_old, accelerations_new, timestep
     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]]
+
+
+def calculate_temperature(particles):
+    """Determine the instantaneous temperature of the system.
+
+    Parameters
+    ----------
+    particles: util.particle_dt, array_like
+        Information about the particles.
+
+    Returns
+    -------
+    float:
+        Calculated instantaneous simulation temperature.
+    """
+    k = 0
+    for i in range(0, particles['xposition'].size):
+        v = np.sqrt(particles['xvelocity'][i] * particles['xvelocity'][i] + particles['yvelocity'][i] *
+                    particles['yvelocity'][i])
+        boltzmann_constant = 1.3806e-23  # joules/kelvin
+        atomic_mass_unit = 1.660539e-27 # kilograms
+        mass_of_argon_amu = 39.948 # amu
+        mass_of_argon = mass_of_argon_amu * atomic_mass_unit # kilograms
+        k += mass_of_argon * v * v / (boltzmann_constant * 2 * particles['xposition'].size)
+    return k
+

pylj/pairwise.py

@@ -0,0 +1,189 @@
+import numpy as np
+
+def compute_forces(particles, box_length, cut_off):
+    """Calculates the forces and therefore the accelerations on each of the particles in the simulation. This uses a
+    12-6 Lennard-Jones potential model for Argon with values:
+
+    - A = 1.363e-134 J m :math:`^{12}`
+    - B = 9.273e-78 J m :math:`^6`
+
+    Parameters
+    ----------
+    particles: util.particle_dt, array_like
+        Information about the particles.
+    box_length: float
+        Length of a single dimension of the simulation square, in Angstrom.
+    cut_off: float
+        The distance greater than which the forces between particles is taken as zero.
+
+    Returns
+    -------
+    util.particle_dt, array_like
+        Information about particles, with updated accelerations and forces.
+    float, array_like
+        Current distances between pairs of particles in the simulation.
+    float, array_like
+        Current forces between pairs of particles in the simulation.
+    float, array_like
+        Current energies between pairs of particles in the simulation.
+    """
+    particles['xaccleration'] = np.zeros(particles['xaccleration'].size)
+    particles['yaccleration'] = np.zeros(particles['yaccleration'].size)
+    pairs = int((particles['xaccleration'].size - 1) * particles['xaccleration'].size / 2)
+    forces = np.zeros(pairs)
+    distances = np.zeros(pairs)
+    energies = np.zeros(pairs)
+    k = 0
+    A = 1.363e-134 # joules / metre ^ {12}
+    B = 9.273e-78 # joules / meter ^ {6}
+    atomic_mass_unit = 1.660539e-27 # kilograms
+    mass_of_argon_amu = 39.948 # amu
+    mass_of_argon = mass_of_argon_amu * atomic_mass_unit # kilograms
+    for i in range(0, particles['xposition'].size-1):
+        for j in range(i, particles['xposition'].size):
+            dx = particles['xposition'][i] - particles['xposition'][j]
+            dy = particles['yposition'][i] - particles['yposition'][j]
+            if np.abs(dx) > 0.5 * box_length:
+                dx *= 1 - box_length / np.abs(dx)
+            if np.abs(dy) > 0.5 * box_length:
+                dy *= 1 - box_length / np.abs(dy)
+            dr = np.sqrt(dx * dx + dy * dy)
+            distances[k] = dr
+            if dr <= cut_off:
+                inv_dr_1 = 1.0 / dr
+                inv_dr_6 = np.power(inv_dr_1, 6)
+                f = (12 * A * (inv_dr_1 * inv_dr_6 * inv_dr_6) - 6 * B * (inv_dr_1 * inv_dr_6))
+                forces[k] = f
+                e = (A * (inv_dr_6 * inv_dr_6) - B * inv_dr_6)
+                energies[k] = e
+                particles['xacceleration'][i] += (f * dx / dr) / mass_of_argon
+                particles['yacceleration'][i] += (f * dy / dr) / mass_of_argon
+                particles['xacceleration'][j] -= (f * dx / dr) / mass_of_argon
+                particles['yacceleration'][j] -= (f * dy / dr) / mass_of_argon
+            else:
+                forces[k] = 0.
+                energies[k] = 0.
+            k += 1
+    return particles, distances, forces, energies
+
+def compute_energy(particles, box_length, cut_off):
+    """Calculates the total energy of the simulation. This uses a
+    12-6 Lennard-Jones potential model for Argon with values:
+
+    - A = 1.363e-134 J m :math:`^{12}`
+    - B = 9.273e-78 J m :math:`^6`
+
+    Parameters
+    ----------
+    particles: util.particle_dt, array_like
+        Information about the particles.
+    box_length: float
+        Length of a single dimension of the simulation square, in Angstrom.
+    cut_off: float
+        The distance greater than which the energies between particles is taken as zero.
+
+    Returns
+    -------
+    util.particle_dt, array_like
+        Information about particles, with updated accelerations and forces.
+    float, array_like
+        Current distances between pairs of particles in the simulation.
+    float, array_like
+        Current energies between pairs of particles in the simulation.
+    """
+    pairs = int((particles['xaccleration'].size - 1) * particles['xaccleration'].size / 2)
+    distances = np.zeros(pairs)
+    energies = np.zeros(pairs)
+    k = 0
+    A = 1.363e-134  # joules / metre ^ {12}
+    B = 9.273e-78  # joules / meter ^ {6}
+    for i in range(0, particles['xpositions']-1):
+        for j in range(i, particles['xpositions']):
+            dx = particles['xposition'][i] - particles['xposition'][j]
+            dy = particles['yposition'][i] - particles['yposition'][j]
+            if np.abs(dx) > 0.5 * box_length:
+                dx *= 1 - box_length / np.abs(dx)
+            if np.abs(dy) > 0.5 * box_length:
+                dy *= 1 - box_length / np.abs(dy)
+            dr = np.sqrt(dx * dx + dy * dy)
+            distances[k] = dr
+            if dr <= cut_off:
+                inv_dr_1 = 1.0 / dr
+                inv_dr_6 = np.power(inv_dr_1, 6)
+                e = (A * (inv_dr_6 * inv_dr_6) - B * inv_dr_6)
+                energies[k] = e
+            else:
+                energies[k] = 0.
+            k += 1
+    return particles, distances, energies
+
+def calculate_pressure(particles, box_length, temperature, cut_off):
+    r"""Calculates the instantaneous pressure of the simulation cell, found with the following relationship:
+
+    .. math::
+        p = \langle \rho k_b T \rangle + \bigg\langle \frac{1}{3V}\sum_{i}\sum_{j<i} \mathbf{r}_{ij}\mathbf{f}_{ij} \bigg\rangle
+
+    Parameters
+    ----------
+    particles: util.particle_dt, array_like
+        Information about the particles.
+    box_length: float
+        Length of a single dimension of the simulation square, in Angstrom.
+    temperature: float
+        Instantaneous temperature of the simulation.
+    cut_off: float
+        The distance greater than which the forces between particles is taken as zero.
+
+    Returns
+    -------
+    float:
+        Instantaneous pressure of the simulation.
+    """
+    pres = 0.
+    A = 1.363e-134  # joules / metre ^ {12}
+    B = 9.273e-78  # joules / meter ^ {6}
+    for i in range(0, particles['xpositions'] - 1):
+        for j in range(i, particles['xpositions']):
+            dx = particles['xposition'][i] - particles['xposition'][j]
+            dy = particles['yposition'][i] - particles['yposition'][j]
+            if np.abs(dx) > 0.5 * box_length:
+                dx *= 1 - box_length / np.abs(dx)
+            if np.abs(dy) > 0.5 * box_length:
+                dy *= 1 - box_length / np.abs(dy)
+            dr = np.sqrt(dx * dx + dy * dy)
+            distances[k] = dr
+            if dr <= cut_off:
+                inv_dr_1 = 1.0 / dr
+                inv_dr_6 = np.power(inv_dr_1, 6)
+                f = (12 * A * (inv_dr_1 * inv_dr_6 * inv_dr_6) - 6 * B * (inv_dr_1 * inv_dr_6))
+                pres += f * dr
+    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)
+    return pres
+
+def heat_bath(particles, temperature_sample, bath_temp):
+    r"""Rescales the velocities of the particles in the system to control the temperature of the simulation. Thereby
+    allowing for an NVT ensemble. The velocities are rescaled according the following relationship,
+
+    .. math::
+        v_{\text{new}} = v_{\text{old}} \times \sqrt{\frac{T_{\text{desired}}}{\bar{T}}}
+
+    Parameters
+    ----------
+    particles: util.particle_dt, array_like
+        Information about the particles.
+    temperature_sample: float, array_like
+        The temperature at each timestep in the simulation.
+    bath_temp: float
+        The desired temperature of the simulation.
+
+    Returns
+    -------
+    util.particle_dt, array_like
+        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)
+    return particles

pylj/tests/test_md.py

@@ -51,4 +51,27 @@ def test_update_velocities(self):
         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]
+        b = md.calculate_temperature(a.particles)
+        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]
+        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]
+        b = md.calculate_msd(a.particles, a.initial_particles, a.box_length)
+        assert_almost_equal(b, 10e-20)
 

pylj/tests/test_util.py

@@ -55,27 +55,4 @@ def test_pbc_correction(self):
         a = util.pbc_correction(1, 10)
         assert_almost_equal(a, 1)
         b = util.pbc_correction(11, 10)
-        assert_almost_equal(b, 1)
-
-    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]
-        b = util.calculate_temperature(a.particles)
-        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]
-        b = util.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]
-        b = util.calculate_msd(a.particles, a.initial_particles, a.box_length)
-        assert_almost_equal(b, 10e-20)
+        assert_almost_equal(b, 1)