improved testing

benchmark/utilities.py

@@ -4,12 +4,12 @@
 from matplotlib import pyplot as plt
 import numpy as np
 
-fontsize = 45
-fontlegende = 40
+fontsize = 35
+fontlegende = 35
 font = {'family': 'sans', 'color':  'black', 'weight': 'normal', 'size': fontsize}
 plt.rcParams.update({"text.usetex": True, "font.family": "serif", "font.serif": ["Palatino"]})
 
-def complete_panel(ax, xlabel, ylabel, cancel_x=False, cancel_y=False, font=font, fontsize=fontsize, linewidth=2, tickwidth1=2.5, tickwidth2=2, legend=True, ncol=1, locator_x = 2, locator_y = 2, title=None):
+def complete_panel(ax, xlabel, ylabel, gray=[1, 1, 1], cancel_x=False, cancel_y=False, font=font, fontsize=fontsize, linewidth=2, tickwidth1=2.5, tickwidth2=2, legend=True, ncol=1, locator_x = 2, locator_y = 2, title=None):
 
     if xlabel is not None:
         ax.set_xlabel(xlabel, fontdict=font)
@@ -33,29 +33,43 @@ def complete_panel(ax, xlabel, ylabel, cancel_x=False, cancel_y=False, font=font
     ax.yaxis.offsetText.set_fontsize(20)
     ax.minorticks_on()
 
-    ax.tick_params('both', length=10, width=tickwidth1, which='major', direction='in')
-    ax.tick_params('both', length=6, width=tickwidth2, which='minor', direction='in')
+    ax.tick_params('both', length=10, width=tickwidth1, which='major', direction='in', color=gray)
+    ax.tick_params('both', length=6, width=tickwidth2, which='minor', direction='in', color=gray)
     ax.xaxis.set_ticks_position('both')
     ax.yaxis.set_ticks_position('both')
+    ax.xaxis.label.set_color(gray)
+    ax.yaxis.label.set_color(gray)
+
+    ax.spines['left'].set_color(gray)
+    ax.spines['top'].set_color(gray)
+    ax.spines['bottom'].set_color(gray)
+    ax.spines['right'].set_color(gray)
 
-    # border of graph
     ax.spines["top"].set_linewidth(linewidth)
     ax.spines["bottom"].set_linewidth(linewidth)
     ax.spines["left"].set_linewidth(linewidth)
     ax.spines["right"].set_linewidth(linewidth)
 
+    ax.tick_params(axis='y', which='both', colors=gray)
+    ax.tick_params(axis='x', which='both', colors=gray)
+
     minor_locator_x = AutoMinorLocator(locator_x)
     ax.xaxis.set_minor_locator(minor_locator_x)
     minor_locator_y = AutoMinorLocator(locator_y)
     ax.yaxis.set_minor_locator(minor_locator_y)
 
     if legend:
-        ax.legend(frameon=False, fontsize=fontlegende, 
+        leg = ax.legend(frameon=False, fontsize=fontlegende, 
                 loc='best', handletextpad=0.5, ncol=ncol,
                 handlelength = 0.86, borderpad = 0.3, 
                 labelspacing=0.3)
 
+        for text in leg.get_texts():
+            text.set_color(gray)
+
 myblue = np.array([20, 100, 255])/255
-myred = np.array([255, 0, 0])/255
+myred = np.array([255, 20, 20])/255
 mygray = np.array([50, 50, 50])/255
-myorange = np.array([255, 165, 0])/255
+myorange = np.array([255, 165, 0])/255
+lightgray = [0.1, 0.1, 0.1]
+darkgray = [0.9, 0.9, 0.9]

docs/source/chapters/chapter2.rst

@@ -175,11 +175,14 @@ and mass (:math:`m`) as provided by default.
                 print("Unknown variable type", type)
         return variable
 
-Utilities
----------
+Probe temperature
+-----------------
 
 The rescaling of the velocity requires to be able to calculate temperature. This is 
-done using the formula :math:`T = 2 E_\text{kin} / N_\text{dof}`.
+done using the formula :math:`T = 2 E_\text{kin} / N_\text{dof}`. Create the *Utilities*
+class, and two functions: *calculate_kinetic_energy* and *calculate_temperature*. The
+kinetic energy is simply calculated as :math:`E_\text{kin} = \sum_i 1/2 m v_i^2`, where
+:math:`v_i` is the velocity of atom :math:`i`.
 
 .. code-block:: python
 
@@ -198,8 +201,8 @@ done using the formula :math:`T = 2 E_\text{kin} / N_\text{dof}`.
             Ndof = self.dimensions*self.number_atoms-self.dimensions
             self.temperature = 2*self.Ekin/Ndof
 
-Output
-------
+Save atom coordinate to file
+----------------------------
 
 Since one would like to visualize the box that is being create, let us start to fill the 
 *Outputs* class and crate our first function named *write_lammps_data*. *write_lammps_data*
@@ -410,4 +413,45 @@ The currently written code is:
 Tests
 -----
 
-Let us test the *InitializeSimulation* to make sure it does what is expected of it.
+tofix : for the test to work, *Outputs* need to inerit from *InitializeSimulation*... does 
+that make sense? Should I introduce already an empty MD class? Would be cleaner...
+
+Let us test the *InitializeSimulation* to make sure it does what is expected of it.
+Let us general a very large number of particle with a given temperature:
+
+.. code-block:: python
+
+    temperature=500 # Kelvin
+
+    x = Outputs(number_atoms=50000,
+                Lx=12,
+                desired_temperature=temperature,
+                seed=69817,
+                )
+    x.write_lammps_data()
+
+    velocity = x.atoms_velocities*x.reference_distance/x.reference_time
+    norm_velocity = np.sqrt(velocity.T[0]**2 + velocity.T[1]**2 + velocity.T[2]**2)
+
+Let us calculate a histogram from the norm of the velocity:
+
+.. code-block:: python
+
+    proba, vel = np.histogram(norm_velocity, bins=200, range=(0, 0.1))
+    vel = (vel[1:]+vel[:-1])/2
+    proba = proba/np.sum(proba)
+
+Finally let us plot the distribution for 3 different temperatures using 
+matplotlib pyplot:
+
+.. code-block:: python
+
+    plt.plot(vel, proba, 'o', color=color, markersize=15, label=r'T = '+str(temperature)+' K')
+
+.. figure:: ../_static/chapter1/velocity-distributions-dark.png
+    :alt: NVE energy as a function of time
+    :class: only-dark
+
+.. figure:: ../_static/chapter1/velocity-distributions-light.png
+    :alt: NVE energy as a function of time
+    :class: only-light