SIMOLANT 02/2024
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.
- When SIMOLANT starts, the temperature is set high and you see
molecules moving freely in a box. This is a (two-dimensional)
gas.
- Cool the system (using slider T)
to a temperature about T=0.6
(watch T= in the right top panel).
The molecules will slow down, the attraction forces prevail and
cause the molecules to condense to a liquid droplet.
- Hence, cool the system to about T=0.3.
Molecules will arrange spontaneously to a hexagonal lattice.
- Try more from menu Prepare system.
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.txt. You need FLTK V1.3 and gcc/g++.
Aims of SIMOLANT
Elementary and high
school
A number of phenomena are shown using a two-dimensional molecular model of matter:
- Condensation of gas and crystallization of liquid on cooling.
- Melting and evaporation on heating.
- Mixing of fluids and gases.
- Capillary action.
- Crystal defects in motion.
- Gas in a gravitational field.
- Impact of a solid body (crystal) to a wall.
Universities
Basic concepts of statistical thermodynamics and molecular simulations can be elucidated:
- Ergodic and deterministic dynamic systems.
- Nucleation, Ostwald ripening.
- Molecular dynamics at constant energy, temperature (several thermostats), and pressure.
- Monte Carlo at constant energy, temperature, and pressure.
- Dependence if the (extended) integrals of motion (total energy) on timestep and thermostat/barostat time constants.
- Flying icecube problem of the Berendsen thermostat
- Radial distribution function and coordination number.
- Walls and periodic boundary conditions.
- Simulation workouts: verification of the Clausius-Clapeyron equation, critical point derermination (by near-critical isotherms), line tension ("surface tension" in 2D).
Molecular model
Interaction potential
The interaction potential is a function of
interparticle distance r:
u(r) = 4/r^8 − 4/r^4
Term 4/r^8 is positive and therefore repulsive. The
term −4/r^4 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).
This is a two-dimensional analogue of the well-known Lennard-Jones
potential. We do not use any real units – the units are defined
by the above formula (and the Boltzmann constant k=1).
Because of consistent implementation of NPT ensembles, the potential
is truncated and smoothed. The truncated potential equals the original
one for r<C, otherwise
u(r) = A(r² − cutoff²)² for r ∈ (C,cutoff),
u(r) = 0 for r > C,
where constants A and C (C<cutoff) are determined
so that both the potential and forces are continuous (the derivate of
forces has jumps). No cutoff corrections are added. The default cutoff=4.
for cutoff<L/2, the nearest-image algorithm is used, otherwise a sum over all periodic images is calculated.
The critical point is approximately T=0.85, P=0.06, ρ=0.3 for cutoff=4.
We have also two types of walls, attractive and repulsive. The
attractive wall is built of the same particles, 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).
The particle–wall potential is thus
uwall(d) = πρ (5/24 d^6−1/d^2)
The repulsive wall is similar, only the attractive (negative) term is omitted.
The wall potentials are not truncated
How realistic is this all?
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 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 similar to vapor condensation in both 2D
and 3D. However, when some 2D systems are cooled, there is no
interface between frozen and liquid phases (Costerlitz-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...).
(E.g., 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 al low
temperatures.
Menus
File
- Open...
Load previously saved file with a configuration.
- Save as...
Save the configuration and (selected) simulation parameters to
a file. Extension .sim is added if not present.
- Protocol name...
Change the name of the file used for printing measured quantities
(see the [ ▯ record ] button).
Extension .txt is added if not present.
- Quit
Quit SIMOLANT. Also pressing Esc will quit.
Prepare system menu
The simulation start, methods, and parameters of simulation, and
quantities to observe are prepared for you to observe phenomena.
- Gas
A dilute system at temperature T=1
represents a (two dimendional) gas.
Try to cool it to about T=0.6 and
observe the changes.
- + Diffusion
As above, and molecules are colored so that you can observe
diffusion.
- + Gravity
Gas in a gravitational field. Watch the vertical density profile.
- Vapor-liquid equilibrium A "vessel" periodic in the
x-direction is prepared with liquid at bottom. Watch pressure and
density as you cool or heat the system a bit. You may try Metropolis
Monte Carlo.
See also clausius-clapeyron.pdf workout
- Nucleation
Dilute gas is prepared and quickly cooled (quenched). Several
small droplets soon form. They coalesce to bigger droplets. For a
more realistic simulation corresponding better to the condensation
of water vapor in wet air, the Langevin thermostat is used – heat
transfer between air and droplets is modeled by random forces. If
you are patient, you may observe another phenomenen
responsible for growth of large droplets at the expense of small
ones – Ostwald ripening. The trick is that small droplets
evaporate faster than the big ones because of their curved
surface. On the contrary, bigger droplets catch more eagerly
flying molecules and grow faster.
- Liquid
A small droplet is produced at
temperature T=0.6 in the periodic
boundary consitions. It will solidify on cooling and evaporate on
heating up.
With a smaller number of particles and lower density ans small τ, the
flyinf icecube artifact can be observed with the Berendesen
thermostat.
- Capillary
Capillary action under a microscope. The left, right, and bottom
walls are attractive and a small gravity is turned on. You can see
a meniscus. Try to make the left and right walls repulsive and
observe the changes!
- Bubble (cavity)
A cavity (commonly called "bubble") in liquid.
- Hexagonal crystal
A perfect piece of hexagonal crystal
at T=0.2 is prepared.
The number of particles may change to match a regular hexagon.
Molecules are colored according to the number of neighbors.
Hints:
- Melt the crystal (increase T).
- Impact it to the ground (use gravity) – it will
recrystallize because of low T.
- Smash it to the ground at constant energy (switch to
the NVE molecular dynamics and then use gravity). A crystal will
melt and perhaps evaporate; set dt=0.02 (using cmd: input field)
in case of problems with too large forces.)
- + Edge defect
A piece of crystal with an edge dislocation. The defect quickly
travels through the crystal and eventually disappears at the
crystal surface.
- + Vacancy
One missing atom in a perfect crystal. This defect diffuses only
slowly.
- + Interstitial
One additional atom in a perfect crystal. This atom has a high
energy and again diffuses out to the surface.
Method menu
Here the simulation method (MC/MD), ensemble (NVE, NVT, NPT) and
thermostat/barostat flavor are selected.
Boundary conditions menu
- Box (soft walls)
Molecules are contained in a square box bounded by four walls
(attractive or repulsive – see below).
- Slit (x-periodic, y-walls)
The system is periodic in the x-direction: Any molecule
leaving the box to the right will appear from the left.
- Periodic
The system is periodic in both the x- and y-directions. The
topology is the same as of a tire.
Magic numbers of particles for crystals: 15, 56, 209, 780. This means
that a crystal without defects matching well a square periodic box is
possible. NPT Berendsen method is recommended.
Show menu
Select quantities to measure and show are a graph to draw. If button
[ ▯ record ] is pressed, the averages of the selected quantities and
graphs are recorded.
- Minimum
Measuring most quantities is turned off, shown are only the
potential energy and temperature/kinetic energy (if applicable) and
basic method parameters. You should not use this option
with ▯ record because most quantities
are not measure.
- Quantities
As above and additional quantities are measured and shown as text:
Pressure (calculated from the virial of force); this pressure is used in barostats
Diagonal components of the pressure tensor (calculated from the virial of force)
Components of pressure on walls (calculated from forces on walls, if applicable)
Line tension ("surface tension" in 2D).
γ = Ly*(Pyy−Pxx) (makes sense only in the slab geometry).
All the following options include these extended measurements, which
are exported if ▯ record is pressed.
- Energy/enthalpy convergence profile
The running convergence profile of energy (potential + kinetic; for
Nosé-Hoover thermostat and MTK without the extended degrees of
freedom; for the Creutz method with the bag energy) or enthalpy (in
the NPT ensembles) is shown. The graph is scaled to the running
miminum–maximum range (as printed) and the standard deviation of the
shown data is printed.
Hints:
Reset the graph by button [ reset ].
Move slider simulation speed so that stride=1.
Use slider measurement block as needed.
- Temperature
The running convergence profile of temperature (kinetic in MD,
averaged bag in MC).
- Volume convergence profile
The volume convergence profile (NPT only).
- Pressure convergence profile
The convergence profile of virial pressure.
- Integral of motion convergence profile
The conserved energy in MD (NVE: total energy, Nosé-Hoover and MTK:
extended energy).
- Radial distribution function (RDF)
The RDF represents the probability of finding a pair of molecules a
given distance r apart, normalized so
that for interactionless particles (ideal gas) it is 1. On
the r-axis,
rmin denotes the potential minimum.
The RDF curve is scaled from zero to the maximum. To look better,
the curve is smoothed by applying two (for block=1) or one (block>1)
three-point formulas.
The RDF exported (after recording) is not smoothed
- Vertical density profile
The vertical density profile is a horizontally averaged number
density as a function of the y-coordinate of atoms (from
bottom). If the system is y-periodic, the profile (of a
horizontal slab) is centered at every step (i.e., it does not matter
when the slab diffuses up and down). The printed profile is placed
so that the box center is zero (range [−L/2,L/2], the shown profile
is in [0,L]. If the box dimensions fluctuate in the NVT ensembles,
the y-box side is always divided in the same number of
histogram bins (the bin width fluctutes); the final graph uses the
averaged volume (area) for the y-coordinate as well as
scaling. The y-range of the shown graph is dynamically
adjusted; overshooted values are orange. Smoothing applies as for
RDF. Instead of homogeneous pressure, both components Pxx and Pyy of
the pressure tensor are shown.
- Droplet radial density profile
Droplet radial density profile shows the full angle (0–360°)
averaged number density from the box center (if walls are
present), or from the droplet or slab center defined as
(L/2π)*atan2(∑ sin(xi*2π/L),∑ cos(xi*2π/L)).
- Cavity radial density profile
As above, but in periodic boundary conditions with respect to
the cavity center,
(L/2π)*(atan2(∑ sin(xi*2π/L),∑ cos(xi*2π/L))+π).
Help menu
- Help
This help.
- About
Copyright and acknowledgements.
Right Panel
Symbols
Top line in the right panel shows the number of particles and
the MC|MD/ensemble/method.
MC→MD means that MC is
running to remove overlaps and will be switched to MD.
Symbol in
brackets [MC|MD/ensemble/method] denotes that a switch
to MC and back is possible in case of problems.
Variables
N | Number of particles. |
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). |
g | Acceleration of gravity. |
d | Radius of the trial displacement (MC only), in the units of Lh=L/2. |
dV | Max. 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). |
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) |
Pid | Pressure calculated from the ideal gas equation of state, Pid = ρ kT (here, k=1). |
Z | Compressiility factor, Z=P/Pid. |
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). |
γ | 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. |
H | Enthalpy. |
Sliders panel
Sliders can be moved by a mouse or (when the mouse pointer is over)
by arrows (ctrl-arrows for fine moves).
Instead of mouse, Tab, Shift-Tab, Enter, etc., work as usual.
- T Temperature.
In the Metropolis Monte Carlo method it is directly a control
parameter of the method. In molecular dynamics a thermostat is
coupled to the system with a typical (correlation or coupling)
time τ.
- τ Thermostat coupling time
(MD only).
Long time means that the system is well insulated from the
thermostat and the system temperature approaches the predefined
value T slowly. Short time means that heat can flow between
the system and thermostat very fast.
- d Trial displacement radius
(MC only).
Too long displacements change the configuration a lot, however,
they are (at lower temperatures) rarely accepted and the
simulation is inefficient. Too short displacements are almost
always accepted, however, the old and new configurations are
almost identical. See the acceptance ratio. Value d=0 (bottom)
means that the trial displacement radius is dynamically adjusted
to reach acc.r.=1/3; in addition, long moves (new random position)
with probability 2 % are performed.
- g Acceleration of gravity.
Negative values point downwards (as in the Earth gravitational
field), positive upwards. In the periodic boundary conditions the
forces of gravity are inactive.
- ρ Number density ρ=N/Area.
Since N is kept constant, increasing
density actually decreased the box size L, and vice versa. However,
for better view, the box is displayed in the same size and the
particle diameters are changed instead. (Imagine you view the screen
from such a distance that you see the particles as discs of the
original size.) Particle coordinates are rescaled with the box.
Area is defined (approximately) as
follows:
- Area = (L−1.2)² in the box (with walls)
- Area = L (L−1.2) in the slit
- Area = L² in the periodic boundary conditions
- P Pressure.
Parameter of NPT simulations.
- N Number of particles.
Change the number of particles.
Note that the density is kept constant so that the box is
rescaled, too.
If MD is running, it will likely be switched to Monte Carlo.
- decrease N: randomly selected particles are removed
- increase N: particles are inserted at random places
Walls
Clicking the respective buttons will toggle repulsive (red) and attractive
(green) walls.
Expert
Some more advanced functions.
Switches
- molecule size:
Controls the size of molecules.
- Real – The disk diameter equals the collision radius,
i.e., the distance where u(r)=0 (it
is r=1).
- Small – 50 % of the collision radius
- Dot – dot of 5 pixels, useful with Lines or Trace
- Pixel – one pixels
- draw mode:
- Movie – Molecules in motion.
- Traces – Trajectories are visible and gradually erased.
- Lines – Trajectories are visible.
- Nothing – Do not draw molecules nor box. Some speed is gained, especially for small systems.
- color mode:
- Black – Molecules are colored black (and this color is
kept).
- y-split – Top half is colored red, bottom half
blue (and these colors are kept). Useful to study diffusion.
- Neighbors – Molecules are drawn in colors according
to the number of neighbors (molecules lying closer than 1.55:
a tick will appear in the RDF graph to mark this distance).
New colors are calculated every step (unless you select Keep).
- Random – Colorize atoms randomly (and keep these colors).
- Keep – Previously selected colors are kept.
- reset
Restore the default behavior (movie of black
molecules). In addition, red hints are removed. The Energy
convergece profile is reset, too (if shown).
- set MC move
In MC simulations: The
trial move size is adjusted automatically. If this button is
pressed, the obtained d (and in NPT also dV) is set and the MC simulation continues with
this constant d and dV.
- run
Stop the simulation temporarily. All other functions
(parameter change etc.) work as ususal and become active as soon
as the simulation is resumed. Useful to chage parameters if the
current setup keeps launching errors.
Simulation speed and measurement blocking
- simulation speed
This slider controls the speed of the simulation/movie and the stride
of calculations (of quantities).
- Leftmost: every frame is shown and delays are added between frames.
- Left: 1–2 delays are added and moderate delays.
- Center: short delays, every 3rd configuration shown.
- Right: short delays, only every n-th configuration is shown.
- Rightmost: every 15-th configuration is shown.
- measurement block
Quantities of interest (kinetic temperature, pressure, RDF) are
averaged over certain number of simulation steps (sweeps in MC) and
then displayed at once. The block length is in range [1,100], the
default is 10.
Command line options
FLTK options are recognized. Note a space between the option and its (optional) value. Useful ones:
-g WxH+X+Y | geometry (W,H=window width and height, X,Y = position). Do not use less than -g 958x683 |
-title "WINDOW TITLE" | change window title |
-ico | start iconized |
-notooltips | do not show tooltips (help text appearing over buttons and sliders) |
SIMOLANT options are in uppercase. The following number/value is without space! No check is done whether the values given are in a reasonable range!
-Dsteps | debug 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.) |
-Ifile | (file=name.sim) open sim-file on start |
-NN | initial number of particles |
-Pcomma-separated assignments | Commands (assignments) as in the "Expert" mode, comma-separated. Executed after initialization (-I, -N) Example: -Pwall=1,th=1 |
-Sspeed | initial value of the simulation speed slider [default=0.7], affects the timer and (for speed>0) the number of the skipped frames. Normal range = [-1.5,2.5] |
-Ttimer | timer delay in s for speed=0 [default=0.05] |
-Wrho | wall number density [default=0.75] |
file | (file=name.sim) open sim-file on start (when used as a single argument; use -Ifile with other arguments) |
Known bugs and caveats
- There are timing, thread-priority, and/or efficiency problems.
The original dumb X11-only version was faster.
- Simple inefficient O(N²) algorithm (no linked-cell list, no
neighbor list...).
- Both / and \ treated as directory separator (on all platforms).
- Inconsistent command line option format.
- Under hard conditions (fast changes of N, density), an
out-of-box particle may pass unnoticed and cause various
errors.
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
- 02/2024 small bug ~O(dt) fixed in the Bussi thermostat, default timestep changed, default qtau for MTK increased
- 04/2023 all whites ignored in a command (cmd: field), hints removed, bug (possible string overflow) fixed
- 03/2022 cumulative distribution functions added to the additional recorded output
control of these additional outputs moved to a separate widget (not in menu File)
bug in surface tension fixed
default thermostat changed to BUSSI (was BERENDSEN)
- 12/2021 bug in surface tension introduced (returned to original good value 03/2022)
NPT L limited to 1e10 to prevent crash in rare cases
vertical density profile centered
vertical density profile better NPT-compatible (Ly=sqrt(⟨L²⟩) used)
- 09/2021 variables _bc and _th (in Command:) renamed to th and bc
bug fixed (Ekin for MC NPT was left from the previous step and not calculated from the T; thus, Pvir and H were incorrect)
- 10/2020 bugs fixed (wall force, show τ slider)
- 07/2020 barostats, cutoff, many GUI and control changes
- 10/2019 Bussi thermostat (Canonical Sampling through Velocity Rescaling) added
- 05/2019 bug fixed (Creutz's bag), changes to satisfy the paranoic C++ compiler
- 04/2018 bug fixed (command g=), small improvements
- 03/2018 Maxwell-Boltzmann thermostat added, small improvements
- 11/2017 convergence profile of total (kinetic+potential) energy added
- 03/2017 comma accepted in cmd: input, toggle [,] will use decimal comma instead of period in the recorded files
- 11/2016 bug fixed (rescaled window + periodic b.c. not properly shown)
- 10/2015 parameter input from keyboard added
- 03/2015 export of convergence profiles for further analysis added
Supported by
- 2012–2013 TUL Liberec, EU: esf (evropský sociální fond v ČR), MŠMT: OP Vzdělání pro konkurenceschopnost (FLTK version)
- 2001–2022 University of Chemistry and Technology, Prague (DOS+Turbo C/X11+*nix versions+FLTK versions)
- 1985–2001 Institute of Chemical Process Fundamentals (DOS+Pascal version)
Acknowledgements
- Ivo Nezbeda (the first version of SIMOLANT was written in Pascal
for a simulation seminar which was held at the Institute of
Chemical Technology in about 1985)
- Karel Matas (FLTK advise)