Visualization

Introduction

When you are running a simulation, it is often useful to see what is going on by visualizing particles in a 3D view or by plotting observables over time. That way, you can easily determine things like whether your choice of parameters has led to a stable simulation or whether your system has equilibrated. You may even be able to do your complete data analysis in real time as the simulation progresses.

Thanks to ESPResSo's Python interface, we can make use of standard libraries like OpenGL (for interactive 3D views) and Matplotlib (for line graphs) for this purpose. We will also use NumPy, which both of these libraries depend on, to store data and perform some basic analysis.

Simulation

First, we need to set up a simulation. We will simulate a simple Lennard-Jones liquid. Particles will be placed randomly in the simulation box. We will energy-minimize the system to remove overlaps, and then thermalize the system with Langevin dynamics. We can measure the energy as a function of time using system.analysis.energy().

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 18})
import numpy as np

import sys
import tqdm
import logging
import espressomd
logging.basicConfig(level=logging.INFO, stream=sys.stdout)
np.random.seed(42)
matplotlib_notebook = True  # toggle this off when outside IPython/Jupyter

espressomd.assert_features("WCA")

# interaction parameters (purely repulsive Lennard-Jones)
lj_eps = 1.0
lj_sig = 1.0

# system
system = espressomd.System(box_l=[10, 10, 10])
system.time_step = 0.0001
system.cell_system.skin = 0.4

# particle parameters (dense liquid)
density = 0.7
n_part = int(system.volume() * density)

# integration
int_steps = 500
int_n_times = 100

#############################################################
#  Setup System                                             #
#############################################################

# interaction setup
system.non_bonded_inter[0, 0].wca.set_params(epsilon=lj_eps, sigma=lj_sig)

# particle setup
system.part.add(pos=np.random.random((n_part, 3)) * system.box_l)

#############################################################
#  Energy Minimization                                      #
#############################################################

system.integrator.set_steepest_descent(f_max=0, gamma=1.0, max_displacement=lj_eps * 0.01)
# minimize until energy difference < 5% or energy < 1e-3
relative_energy_change = float('inf')
relative_energy_change_threshold = 0.05
energy_threshold = 1e-3
energy_old = system.analysis.energy()['total']
logging.info(f'Energy minimization starts')
logging.info(f'energy: {energy_old:.2e}')
for i in range(20):
    system.integrator.run(50)
    energy = system.analysis.energy()['total']
    logging.info(f'energy: {energy:.2e}')
    relative_energy_change = (energy_old - energy) / energy_old
    if relative_energy_change < relative_energy_change_threshold or energy < energy_threshold:
        break
    energy_old = energy
else:
    logging.info(f'Energy minimization did not converge in {i + 1} cycles')
system.integrator.set_vv()

#############################################################
# Thermalization                                            #
#############################################################

system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=42)
system.time = 0  # reset system timer
energies = np.zeros((int_n_times, 2))
logging.info(f'Thermalization starts')
for i in tqdm.tqdm(range(int_n_times)):
    system.integrator.run(int_steps)
    energies[i] = (system.time, system.analysis.energy()['total'])

plt.plot(energies[:,0], energies[:,1])
plt.xlabel("Time")
plt.ylabel("Energy")
plt.show()

INFO:root:Energy minimization starts
INFO:root:energy: 2.08e+13
INFO:root:energy: 6.57e+01
INFO:root:energy: 6.10e+00
INFO:root:energy: 3.92e+00
INFO:root:energy: 2.75e+00
INFO:root:energy: 1.70e+00
INFO:root:energy: 9.51e-01
INFO:root:energy: 2.68e-01
INFO:root:energy: 5.78e-02
INFO:root:energy: 3.74e-02
INFO:root:energy: 3.35e-05
INFO:root:Thermalization starts

100%|██████████| 100/100 [00:09<00:00, 10.57it/s]

Live plotting

We will write a main() callback function to store the total energy at each integration run into a NumPy array. We will also create a function to draw a plot after each integration run.

# setup matplotlib canvas
plt.xlabel("Time")
plt.ylabel("Energy")
plot, = plt.plot([0], [0])
if matplotlib_notebook:
    from IPython import display
else:
    plt.show(block=False)

energies = np.zeros((int_n_times, 2))

# setup matplotlib update function
current_time = -1
def update_plot():
    i = current_time
    if i < 3:
        return None
    plot.set_xdata(energies[:i + 1, 0])
    plot.set_ydata(energies[:i + 1, 1])
    plt.xlim(0, energies[i, 0])
    plt.ylim(energies[:i + 1, 1].min(), energies[:i + 1, 1].max())
    # refresh matplotlib GUI
    if matplotlib_notebook:
        display.clear_output(wait=True)
        display.display(plt.gcf())
    else:
        plt.draw()
        plt.pause(0.01)

# define a callback function
def main():
    global current_time
    for i in range(int_n_times):
        system.integrator.run(int_steps)
        energies[i] = (system.time, system.analysis.energy()['total'])
        current_time = i
        update_plot()
    if matplotlib_notebook:
        display.clear_output(wait=True)

system.time = 0  # reset system timer
main()

if not matplotlib_notebook:
    plt.close()

Live visualization and plotting

To interact with a live visualization, we need to move the main integration loop into a secondary thread and run the visualizer in the main thread (note that visualization or plotting cannot be run in secondary threads). First, choose a visualizer:

import espressomd.visualization
import threading

visualizer = espressomd.visualization.openGLLive(system)

Then, re-define the main() function to run the visualizer:

def main():
    global current_time
    for i in range(int_n_times):
        system.integrator.run(int_steps)
        energies[i] = (system.time, system.analysis.energy()['total'])
        current_time = i
        visualizer.update()

system.time = 0  # reset system timer

Next, create a secondary thread for the main() function. However, as we now have multiple threads, and the first thread is already used by the visualizer, we cannot call update_plot() from the main() anymore. The solution is to register the update_plot() function as a callback of the visualizer:

# setup new matplotlib canvas
if matplotlib_notebook:
    plt.xlabel("Time")
    plt.ylabel("Energy")
    plot, = plt.plot([0], [0])
    plt.xlim(left=0)

# execute main() in a secondary thread
t = threading.Thread(target=main)
t.daemon = True
t.start()

# execute the visualizer in the main thread
visualizer.register_callback(update_plot, interval=int_steps // 2)
visualizer.start()

<MagicMock name='mock.openGLLive().start()' id='134590741967552'>

While the simulation is running, you can move and zoom around with your mouse.

Important: closing the visualizer GUI will exit the Python session!

If the trajectory runs too quickly, try decreasing int_steps by a factor 10.

Alternative: live visualization without plotting

If live plotting is not important, the code in the previous section simplifies to:

import espressomd.visualization
import threading
visualizer = espressomd.visualization.openGLLive(system)

def main():
    for i in range(int_n_times):
        system.integrator.run(int_steps)
        energies[i] = (system.time, system.analysis.energy()['total'])
        visualizer.update()

system.time = 0  # reset system timer

# execute main() in a secondary thread
t = threading.Thread(target=main)
t.daemon = True
t.start()
visualizer.start()

Alternative: live visualization only

If recording the energy as a time series is not important, the code in the previous section simplifies even further:

import espressomd.visualization
visualizer = espressomd.visualization.openGLLive(system)
visualizer.run()

Customizing the OpenGL visualizer

Visualization of more advanced features of ESPResSo is also possible (e.g. bonds, constraints, lattice-Boltzmann) with the OpenGL visualizer. There are a number of optional keywords that can be used to specify the appearance of the visualization, they are simply stated after system when creating the visualizer instance. See the following examples:

# Enables particle dragging via mouse:
visualizer = espressomd.visualization.openGLLive(system, drag_enabled=True)

# Use a white background:
visualizer = espressomd.visualization.openGLLive(system, background_color=[1, 1, 1])

# Use red color for all (uncharged) particles
visualizer = espressomd.visualization.openGLLive(system, particle_type_colors=[[1, 0, 0]])

The visualizer options are stored in the dict visualizer.specs. The following snippet prints out the current configuration nicely:

for k in sorted(visualizer.specs.keys(), key=lambda s: s.lower()):
    print(f"{k:30}  {visualizer.specs[k]}")

All keywords are explained in the user guide section on the OpenGL visualizer. Specific visualization examples for ESPResSo can be found in the /samples folder. You may need to recompile ESPResSo with the required features used in the examples.