SIMOLANT 03/2025

Contents

Quick popular start

Molecules attract each other when they are close together (the "atom disks" almost touch). If the molecules are even closer, they begin to repel (the disks overlap). If the atoms are far away, they (almost) do not feel each other.

Installation

Get the respective zip-file (simolant-win32.zip for Windows or simolant-amd64.zip for 64-bit linux) and unzip it to a folder. Run simolant.exe or simolant from this folder (this is necessary for the online help).

On other platforms, download simolant-sources.zip and follow the instructions in README.md. You need FLTK V1.3 or higher (not the V2 fork) and gcc/g++.

Aims of SIMOLANT

Elementary and high school

A number of phenomena are shown using a two-dimensional molecular model of matter:

Universities

Basic concepts of statistical thermodynamics and molecular simulations can be elucidated:

Molecular models

A molecular model aka interaction potential aka force field is the energy of a pair of atoms r apart expressed as a function, u(r). The total energy is the sum of these pair terms over all pairs of atoms. We do not use any real units – the units are defined by the respective formulas and the Boltzmann constant k=1. The potential type and its parameters are printed at the top of the right panel next to the number of particles; e.g., "N=300 LJ c=4"

Lennard-Jones-like interaction potential

The default interaction potential is a 2D analogue of the Lennard-Jones potential:

    u(r) = 4/r⁸ − 4/r

Term 4/r⁸ is positive and therefore repulsive. The term −4/r⁴ is negative and therefore attractive; it has a longer range than the repulsive term. The minimum of the potential is −1 at r=2^(1/4)=1.189 (line marked as rmin in the graph of RDF).

The potential is truncated at c (can be set in window cmd:) and smoothed:

    u(r) = 4/r⁸ − 4/r⁴   for r > c₁,
    u(r) = A(r² − c²)²   for r ∈ (c₁,c),
    u(r) = 0                 for r > c,

where constants A and c₁ (c₁<c) are determined so that both the potential and forces are continuous (the derivate of forces has jumps). No cutoff corrections are added. The recommended default c=4, the minimum is c=2.232.

The critical point for the default c=4 version is approximately T=0.85, P=0.06, ρ=0.3 for c=4. The Boyle temperature is 3.0138.

Walls and gravity

The walls can keep the particles in a box or a slit. There are two types of walls, attractive and repulsive. The attractive wall is built of the same particles as the 2D Lennard-Jones (without truncation), but imagine them ground to tiny pieces and smoothly distributed on one side of the wall with number density ρ=N/L²=0.75 (can be changed by option -Pwall=VALUE or from cmd:). The particle–wall potential is thus

    uwall(d) = πρ (5/24 d⁶ −1/d²)

The repulsive wall is similar, only the attractive (negative) term is omitted. The walls are active in the respective boundary conditions, the sign can be changed in panel Walls.

If horizontal walls are present, gravity can be turn on. The potential per particle is

    ugravity(y) = −gy

Negative g is the usual gravity.
Bug or feature: the virial of the force of gravity is not included in the pressure tensor calculation.

How realistic is this?

The two-dimensional model used by SIMOLANT is easier to show and simulate than a realistic three-dimensional model of, e.g., water. Nevertheless, basic concepts shown are qualitatively the same as in three dimensions. The only notable exception is a possibly different nature (universality class) of the freezing transition. If water freezes, a nucleus of ice appears and then it grows – there is an interface between both phases. This is the first-order phase transition, other example is vapor condensation in both 2D and 3D. However, when some 2D systems are cooled, there is no interface between ordered (frozen) and disordered (liquid) phases (Kosterlitz–Thouless continuum transition). Instead, hexagonally ordered areas grow anywhere in the system. At the freezing point the size of these areas diverges to infinity, i.e., in equilibrium we obtain a single crystal structure spanning the whole system. Other 2D systems freeze as in 3D (i.e., via the first order phase transition), but the nucleation phenomena are not so apparent; this is likely the case of 2D hard disks. The type of freezing of our 2D system is not known.

The model is more appropriate for 2D systems as surfactants floating on a surface of water. Such molecules form 2D gas, 2D liquid, and hexagonally ordered 2D solids (and even more...). The surface pressure can be easily demonstrated in a kitchen by throwing a few matches to a bowl of water and touching the surface by a piece of soap.

The potential decreases to zero more slowly than London forces in 3D; therefore, gas at moderate densities is less ideal than in 3D – the compressibility factor Z deviates from 1, especially at low temperatures.

More interaction potentials

In addition to the standard Lennard-Jones-like potential, several more or less weird potentials are supported:

Pair sum and linked-cell list

Calculating the interaction potential and forces needs to consider all neighbors of a particle within the cutoff distance c. Three algorithms are available.

  1. If the cutoff is longer than half box size, c>L/2, one particle may interact with different periodic replicas of other particles. Thus, the sum over all pairs of particles is amended by the triple sum over all replicas of the central cell that may contain interacting particles. For c>L, a particle even interacts with its own periodic replicas. However, condition c>L/2 is not typical for most examples with many particles, so that this algorithm is rarely activated.
  2. If the cutoff is a little bit shorter than half box size or for a few particles only, the pair sum is calculated as a simple double sum.
  3. If the cutoff is shorter than half box size, which typically happens for many particles, the pair sum is calculated by the linked-cell list method. It is based on domain decomposition of the box into ncell×ncell cells. For each particle in each cell, only the cells within the cutoff (and “ahead” to avoid counting each pair twice) are considered and scanned for possible interactions. Only calculations of forces in molecular dynamics and the box swelling/shrinking in NPT MC are treated by this algorithm. The sums needed in Metropolis MC displacements and calculating RDF and other quantities at the end of a cycle are evaluated by algorithms 1 or 2.
The algorithm choice is controlled by variable ncell available from input field cmd: Since the default ncell=0 may not be optimum, you are encouraged to try ncell>1 by hand. The top right panel shows the timing sampled every 2 s.

Menus

File

Prepare system

The simulation start, methods, and parameters of simulation, and quantities to observe are prepared for you to observe phenomena.

Force field

Method

Here the simulation method (MC/MD), ensemble (NVE, NVT, NPT) and thermostat/barostat flavors are selected.

Boundary conditions

Show

Select quantities to measure and show and a graph to draw. If button [▯ record] is pressed, the averages of the selected quantities and graphs (if selected in widget include:) are recorded.

Tinker

Technical setup affecting performance.

Help

Right Panel

Symbols

Top line in the right panel shows the number of particles and the MC|MD/ensemble/method. Symbol MC→MD means that MC is running to remove overlaps and will be switched to MD.
The text in brackets [MC|MD/ensemble/method] denotes that a switch to MC and back is possible in case of problems.

Variables

These parameters are shown during simulation in the right panel and also printed to the protocol simolant.txt

N Number of particles.
LJ c= Cutoff for the Lennard-Jones-like potential.
WCALJ Repulsive potential (c=1.1892).
PD a=, c= Penetrable disks spheres, a=potential depth, c=disk size (cutoff)
DWa,b,c Double well with parameters a,b,c.
T Thermostat or method temperature.
τ Thermostat correlation (coupling) time (MD only).
L Box size (as derived from density).
ρ Number density, ρ = N/Area (see below for the definition of the Area).
wall Number density of walls.
g Acceleration of gravity.
d Radius of the trial displacement (MC only), in the units of Lh=L/2.
dV Maximal trial relative volume change in MC NPT.
P Pressure (parameter of a NPT simulation)
stride stride MD steps or MC sweeps is performed between drawing the configuration and measuring quantities.
block Number of measurements in one block (by stride MC sweeps or MD steps).
#.###+#.### ms This is technical info. The first number is the timer constant in millisecond determined from menu Tinker → FPS (Frames Per Second), slider [simulation speed], and menu Tinker → Timeout (for the slider at maximum speed). The second number is the additional time needed: If the CPU time of one cycle (of stride steps) is less than 1/FPS, the it should oscillate around zero (delays are inserted to reach frame rate FPS). If the second number is large and positive, the calculation cannot be tan within one frame and the movie slows down.
ncellSee here for details.

The following quantities are shown during simulation in the right panel and also printed to the protocol and convergence profiles.

t Time or the number of MC sweeps from 0 during recording.
Tkin Kinetic temperature (averaged in blocks, MD only).
Tbag Creutz demon bag temperature (averaged, MC only).
dt Integration time step (MD only); "dt⇒" denotes that dt has been determined automatically from T and τ.
acc.r. Acceptance ratio (ratio of accepted moves to all attempted moved, MC only).
acc.r/(V) Acceptance ratio of volume change in NPT simulations
Pvir Pressure (averaged configurational pressure based on the virial theorem incl. the kinetic term).
V Volume (in NPT)
Z Compressibility factor, Z=P/ρkT (where k=1); it holds Z=1 for ideal gas. In CREUTZ MC, T=⟨Tbag⟩. In NVE MD, T=⟨Tkin⟩. In NVT and NPT simulations, the thermostat T is used. In NPT simulations, the barostat P is used and ρ=N/⟨V⟩.
Pxx Diagonal xx-component of the pressure tensor (similarly Pyy).
P(top_wall) Pressure calculated from averaged force on the top wall (similarly other walls).
MSDx Mean square displacement in x-direction, similarly MSDy.
γ Surface tension (applies to the slab geometry).
Etot Total (potential + kinetic) energy of the whole simulation box.
Epot Potential energy of the whole simulation box.
p |averaged momentum of the simulation box| (0 for momentum-preserving MD, order parameter of the Vicsek model).
H Enthalpy.
Econserved Conserved energy (MD: NVE, NVT Nosé–Hoover, MTK).
pmomentum (velocity) per particle, |∑ v_i|/N.

Parameters (panel of sliders)

Sliders can be moved by a mouse or (when the mouse pointer is over) finely by arrows. The response of most sliders is not linear.

Walls + shifts

Expert

Some more advanced functions.

Switches

Simulation speed and measurement blocking

Command line options

FILE.sim

An argument without leading hyphen is considered as a sim-file and loaded.

FLTK options are recognized. Note a space between the option and its (optional) value. Useful ones:

-g WxH+X+Ygeometry (W,H=window width and height, X,Y = position). The default is -g 1192x736. Do not use less than about -g 1150x700; particularly, a smaller height will cause overlaps of text and widgets.
-title "WINDOW TITLE"change window title
-icostart iconized
-notooltipsdo not show tooltips (help text appearing over buttons and sliders)

SIMOLANT options are UPPERCASE. The following number/value is without space! No check is done whether the values given are in a reasonable range!

filesim-file (configuration) to open
-Cmethodthe circle drawing method: 0=fl_pie, 1=fl_circle, 2=composed of fl_line (default)
-Dstepsdebug mode (timing and technical info printed to stderr every MD or MC steps), disabled in Windows (recompile with -mconsole instead of -mwindows)
-Iinit(init=number) initial system (Gas=0, Diffusion=1, etc.)
-NNinitial number of particles (legacy, also -PN=)
-NPROTOCOLprotocol name without extension, cf. menu File
-Pcomma-separated assignmentsCommands (assignments) as in the "Expert" mode, comma-separated, decimal separator=period. Executed after initialization (-I, -N)
Example: -Pwall=1,method=1
-Tfactorthe timeout of the timer is at least factor times the 1/FPS. Cf. menu Tinker. [default=0.2]
-Udelayadditional delay in s, do not use [default=1e-5]
-Zseedrandom number seed [default=0=time]

Known bugs and caveats

Copying and Credits

This application is © Jiří Kolafa. It is available under the GNU General Public License Version 3

The sources of SIMOLANT can be found on GitHub.

History

Supported by

Acknowledgements