pylj now works in real units

pylj/md.py

@@ -2,7 +2,7 @@
 from pylj import comp, util
 
 
-def initialise(number_of_particles, temperature, box_length, init_conf, timestep_length=5e-3):
+def initialise(number_of_particles, temperature, box_length, init_conf, timestep_length=1e-14):
     """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.
 
@@ -29,10 +29,11 @@ def initialise(number_of_particles, temperature, box_length, init_conf, timestep
     """
     system = util.System(number_of_particles, temperature, box_length, init_conf=init_conf,
                          timestep_length=timestep_length)
-    v = np.sqrt(2 * system.init_temp)
-    theta = 2 * np.pi * np.random.randn(system.particles.size)
-    system.particles['xvelocity'] = v * np.cos(theta)
-    system.particles['yvelocity'] = v * np.sin(theta)
+    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 = comp.compute_forces(system.particles, system.distances,
                                                                             system.forces, system.box_length)
     return system
@@ -94,9 +95,9 @@ def sample(particles, box_length, initial_particles, system):
         arrays.
     """
     temperature_new = util.calculate_temperature(particles)
+    system.temperature_sample = np.append(system.temperature_sample, temperature_new)
     pressure_new = comp.calculate_pressure(particles, box_length, temperature_new)
     msd_new = util.calculate_msd(particles, initial_particles, box_length)
-    system.temperature_sample = np.append(system.temperature_sample, temperature_new)
     system.pressure_sample = np.append(system.pressure_sample, pressure_new)
     system.force_sample = np.append(system.force_sample, np.sum(system.forces))
     system.msd_sample = np.append(system.msd_sample, msd_new)
@@ -151,4 +152,3 @@ def update_velocities(velocities, accelerations, timestep_length):
     velocities[0] += 0.5 * accelerations[0] * timestep_length
     velocities[1] += 0.5 * accelerations[1] * timestep_length
     return [velocities[0], velocities[1]]
-

pylj/sample.py

@@ -221,7 +221,7 @@ def setup_cellview(ax, system):
     """
     xpos = system.particles['xposition']
     ypos = system.particles['yposition']
-    mk = (1052.2 / (system.box_length - 0.78921) - 1.2174)
+    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])
@@ -253,11 +253,10 @@ def setup_rdfview(ax, system):
         The whole system information.
     """
     ax.plot([0], color='#34a5daff')
-    ax.set_xlim([0, system.box_length/2])
+    ax.set_xlim([0, system.box_length / 2])
     ax.set_yticks([])
     ax.set_ylabel('RDF', fontsize=16)
-    ax.set_xlabel('r/Å', fontsize=16)
-
+    ax.set_xlabel('r/m', fontsize=16)
 
 def setup_diffview(ax):
     """Builds the scattering profile visualisation pane.
@@ -269,8 +268,10 @@ def setup_diffview(ax):
     """
     ax.plot([0], color='#34a5daff')
     ax.set_yticks([])
-    ax.set_ylabel('log(I[q])', fontsize=16)
-    ax.set_xlabel('q/Å$^{-1}$', fontsize=16)
+    ax.set_xscale('log')
+    ax.set_yscale('symlog')
+    ax.set_ylabel('I(q)', fontsize=16)
+    ax.set_xlabel('q/m$^{-1}$', fontsize=16)
 
 
 def setup_pressureview(ax):
@@ -282,7 +283,7 @@ def setup_pressureview(ax):
         The axes position that the pane should be placed in.
     """
     ax.plot([0], color='#34a5daff')
-    ax.set_ylabel('Pressure/Pa', fontsize=16)
+    ax.set_ylabel(r'Pressure/$\times10^6$Pa m$^{-1}$', fontsize=16)
     ax.set_xlabel('Time/s', fontsize=16)
 
 
@@ -330,11 +331,11 @@ def update_rdfview(ax, system, average_rdf, r):
     r: array_like
         The radial distribution functions r for each timestep, to later be averaged.
     """
-    hist, bin_edges = np.histogram(system.distances, bins=np.linspace(0, system.box_length/2 + 0.5, 100))
+    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.5 / 2.) * 0.5)
+                 (bin_edges[:-1] + 0.5e-10 / 2.) * 0.5)
     average_rdf.append(gr)
-    x = bin_edges[:-1] + 0.5 / 2
+    x = bin_edges[:-1] + 0.5e-10 / 2
     r.append(x)
 
     line = ax.lines[0]
@@ -357,7 +358,7 @@ def update_diffview(ax, system, average_diff, q):
     q: array_like
         The scattering profile's q for each timestep, to later be averaged.
     """
-    qw = np.linspace(2 * np.pi /(system.box_length/2)*4, system.box_length/10, 100)
+    qw = np.logspace(np.log10(2 * np.pi /(system.box_length)), 10.47, num=1000, base=10)
     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))
@@ -371,7 +372,7 @@ def update_diffview(ax, system, average_diff, q):
     line1.set_xdata(x2)
     line1.set_ydata(y2)
     ax.set_ylim([0, np.amax(y2) + np.amax(y2) * 0.05])
-    ax.set_xlim(0, np.amax(x2))
+    ax.set_xlim(np.amin(x2), np.amax(x2))
 
 
 def update_forceview(ax, system):
@@ -385,11 +386,11 @@ def update_forceview(ax, system):
         The whole system information.
     """
     line = ax.lines[0]
-    line.set_ydata(system.force_sample * 1e-10)
+    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 * 1e-10)-np.amax(system.force_sample * 1e-10) * 0.05,
-                np.amax(system.force_sample * 1e-10)+np.amax(system.force_sample * 1e-10) * 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_tempview(ax, system):
@@ -421,11 +422,12 @@ def update_pressureview(ax, system):
         The whole system information.
     """
     line = ax.lines[0]
-    line.set_ydata(system.pressure_sample * 1e10)
+    data = system.pressure_sample * 1e6
+    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(system.pressure_sample * 1e10) - np.amax(system.pressure_sample * 1e10) * 0.05,
-                np.amax(system.pressure_sample * 1e10) + np.amax(system.pressure_sample * 1e10) * 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):
@@ -437,7 +439,7 @@ def setup_msdview(ax):
         The axes position that the pane should be placed in.
     """
     ax.plot([0], color='#34a5daff')
-    ax.set_ylabel('MSD/Å$^2$', fontsize=16)
+    ax.set_ylabel('MSD/m$^2$', fontsize=16)
     ax.set_xlabel('Time/s', fontsize=16)
 
 

pylj/util.py

@@ -23,14 +23,19 @@ class System:
     timestep_length: float (optional)
         Length for each Velocity-Verlet integration step, in seconds.
     """
-    def __init__(self, number_of_particles, temperature, box_length, init_conf='square', timestep_length=5e-3):
+    def __init__(self, number_of_particles, temperature, box_length, init_conf='square', timestep_length=1e-14):
         self.number_of_particles = number_of_particles
         self.init_temp = temperature
         if box_length <= 600:
-            self.box_length = box_length
+            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))
+        if box_length >= 4:
+            self.box_length = box_length * 1e-10
+        else:
+            raise AttributeError('With a box length of {} the cell is too big to really hold more than one '
+                                 'particle.'.format(box_length))
         self.timestep_length = timestep_length
         self.particles = None
         if init_conf == 'square':
@@ -118,8 +123,8 @@ def calculate_temperature(particles):
     for i in range(0, particles['xposition'].size):
         v = np.sqrt(particles['xvelocity'][i] * particles['xvelocity'][i] + particles['yvelocity'][i] *
                     particles['yvelocity'][i])
-        k += 0.5 * v * v
-    return k / particles['xposition'].size
+        k += 66.234e-27 * v * v / ( 1.3806e-23 * 2 * particles['xposition'].size)
+    return k
 
 def calculate_msd(particles, initial_particles, box_length):
     """Determines the mean squared displacement of the particles in the system.