SIMOLANT 03/2025

Contents

• Quick popular start
• Installation
• Aims of SIMOLANT
• Elementary and high schools
• Universities
• Molecular models
• Lennard-Jones-like interaction potential
• Walls and gravity
• How realistic is this?
• More interaction potentials
• Pair sum and linked-cell list
• Menus
  • File
  • Prepare system
  • Force field
  • Method
  • Boundary conditions
  • Show
  • Tinker
  • Help
• Right Panel
• Symbols
• Variables
• Parameters (panel of sliders)
• Walls
• Expert
• Switches
• Simulation speed and measurement blocking
• Command line options
• Known bugs and caveats
• Copying and Credits

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.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:

• 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.
• Swarming (flocking).

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 of the (extended) integrals of motion (total energy) on timestep and thermostat/barostat time constants.
• Flying ice cube problem of the Berendsen thermostat.
• Radial distribution function and coordination number.
• Walls and periodic boundary conditions.
• Soft penetrable disks, repulsive WCA-like potential, double well potential.
• Simulation workouts:
  • Isotherms and the critical point
  • Verification of the Clausius-Clapeyron equation
  • Line tension ("surface tension" in 2D), contact angle
  • Dependence of the line tension on temperature and the total line tension (surface enthalpy)
  • Pressure outside a droplet/inside cavity (Kelvin equation)
  • Surface / interfacial tension, contact angle
  • Mean square displacement
  • Diffusivity, Arrhenius plot, activation energy
  • Coalescence of droplets
  • Can a double-minimum potential give the Penrose quasicrystal?
  • Phase transition in the Vicsek model
  • More under preparation...

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

    u_wall(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

    u_gravity(y) = −gy

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:

• WCA interaction potential
  The WCA potential is the Lennard-Jones (here 2D version) potential truncated at the minimum and shifted:
  
      u(r) = 4/r⁸ − 4/r⁴ +1   for r < c
      u(r) = 0                       for r ≥ c
  
  where c = 2^(1/4).
• Penetrable disks and ideal gas
  Penetrable spheres, here penetrable disks, is a simple model of a polymer coil in a good solvent.
  
      u(r) = a(r²−c²)²/c²     for r < c
      u(r) = 0                     for r ≥ c
  
  The default is a=1, c=2. From cmd:, any value can be set:
      a>0     repulsive penetrable spheres (disks),
      a=0     ideal gas (not too much fun with it),
      a<0     attractive penetrable spheres (disks) are thermodynamically singular (the thermodynamic limit diverges).
• Double wall
  The double wall potential with a good choice of parameters might give a quasicrystal of 2D Penrose tessellation:
  
      u(r) = const (1−r²) (a²−r²) (b²−r²) (c²/r²−1)²   for r < c
      u(r) = 0                                                             for r ≥ c
  
  where, for the sake of temperature scale purpose, const is calculated so that ∫_1^∞ u(r)² dr = 512/1155 is the same as for the 2D Lennard-Jones. Parameters a and b are the nodes (zero values of the potential), it should hold 1 < a < b < c. No check is performed of this condition. The default a=1.25, b=1.4, c=2 set from the Force Field menu is likely close to the Penrose structure, use cmd: to change.

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.

• ncell=0 (default) means that the appropriate method is selected automatically.
• ncell=1 selects the direct pair method 2 (unless c>L/2 for which method 1 is used).
• ncell>1 chooses the linked-list cell method 3 with the selected number of cells per box side (unless c>L/2 for which method 1 is used).

Menus

File

• Open...
  Load a previously saved file with a configuration, see below for the format.
• Save as...
  Save simulation parameters (N bc walls T L wall dt d method P dV tau qtau g a b c) and the configuration (atom positions and velocities) to a file. Extension .sim is added if not present. The sim-file is ASCII and human-readable.
• Protocol name...
  Change the name of the files used for printing measured quantities (see the [▯ record] button). The default is simolant. The browser offers text files (extension .txt), which are always written if recording is selected. If CSV (Comma-Separated Values) files are to be recorded, see [▯ CSV] button, the extension .txt is changed to .csv.
• Export force field
  

  Prints a file (tab-separated TXT or CSV, according to the current setup) with the following columns:
    r = the interatomic distance
    ufull(r) = the full potential
    u(r) = the truncated potential (the same as ufull(r) except for 2DLJ)
    ffull(r) = the full force
    f(r) = the truncated force
    -du(r)/dr_numeric = as above, obtained by a numerical derivative

• Quit
  Quit SIMOLANT. Also pressing Esc will quit.

Prepare system

• 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 (use sliders simulation speed and Measurement block for better precision).
• Vapor-liquid equilibrium
  A "vessel" periodic in the x-direction is prepared with liquid at bottom, T=0.65.
. Watch pressure and density as you cool or heat the system a bit. See also clausius-clapeyron.pdf workout
• Horizontal slab
  A "vessel" periodic in the x-direction is prepared with a liquid slab at center, T=0.6.
• 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 droplet
  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, lower density and small τ, the flying icecube artifact can be observed with the Berendesen thermostat.
• Two droplets
  Two droplets are prepared at temperature T=0.6 in the periodic boundary conditions. They approach each other with velocity v=1 (can be changed from cmd:) The simulation is switched to NVE. Based on the value of v, heating or cooling is observed on the graph of Tkin.
• Bubble (cavity)
  A cavity (commonly called "bubble") in liquid.
• 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!
• Periodic liquid
  Liquid at temperature T=0.7 and the periodic boundary conditions is prepared in the NPT ensemble.
• 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. The colors are chosen from the Okabe and Ito colorblind-safe palette.
  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. Use slide simulation speed if this is too fast.
• + Vacancy
  One missing atom in a perfect crystal. This defect slowly diffuses.
• + Interstitial
  One additional atom in a perfect crystal. This atom has a high energy diffuses very fast.
• Periodic crystal
  An approximate periodic crystal at temperature T=0.3 and the periodic boundary conditions is prepared in the NPT ensemble.
  It is not possible to exactly match a hexagonal lattice to a square boundary conditions; however, the following "magic numbers" of particles lead to less and less distorted crystals:
  
      N = (0, 1,) 4, 15, 56, 209, 780, 2911, 10864, 40545, ...
      Recursively: N_n = 4N_n−1−N_n−2 (N₀ = 0, N₁ = 1)
      Explicitly: N_n = [(2+√3)^n−(2−√3)^n] / (2√3)
  
  The maximum error in the angle (π/3 = 60°) of both basis vectors is 1/2N = 28.65°/N.
  The code first rounds the number of particles to the nearest (in the logarithmic sense) magic number and then it constructs the crystal. The minimum is 15, the maximum (with N= from cmd:) is 2911. Other "less magic" solutions (as doubled magic numbers with the crystal rotated by 45°) are not considered.
  For N=15 or 56, low T and appropriate ρ or P, the best crystal forms spontaneously. For higher N this happens only occasionally and (e.g.) several heating/freezing cycles are needed to obtain a crystal without defects.
  Hint: prepare a crystal, type button [N+] to include an interstitial and [N−] to remove one atom. You will observe a vacancy and interstitial annihilating.
• Vicsek active matter
  Gas of penetrable disks with a=0.2, c=4, low density, and the Vicsek algorithm. This is a model of a flock of birds or school of fish.

Force field

• Lennard-Jones c=4 (recommended)
  Standard 2D Lennard-Jones potential smoothly cut off at c=4. Parameter c can be changed from cmd:, any value c≥2.232 selects this type of the potential.
• Repulsive (WCALJ, c=1.1892)
  Repulsive potential cut at 2D Lennard-Jones minimum and shifted. Acronym WCALJ refers to its use in the Weeks–Andersen–Chandler perturbation theory.
• Penetrable disks (a=1,c=2)
  A model of penetrable disks (spheres in 3D), a 2D analogue of a model of polymer coils in a good solvent. Parameter c=2 selects penetrable disks, the potential height a=1/2 can be changed from cmd:.
• Ideal gas (a=0,c=2)
  As above with no interaction (ideal gas).
• Attractive disks (a=−1,c=2)
  As above with attraction. This model does not have a thermodynamic limit (the Gibbs state is singular).
• Double well (a=1.25,b=1.4,c=2)
  A potential with a double well. Parameters a,b,c can be changed from cmd:; c is the cutoff, it must hold 1 < a < b < c.

Method

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

• Automatic (MC→MD/Bussi CSVR)
  Several first steps are performed by the Metropolis Monte Carlo method at constant temperature T (this is marked by symbol MC→MD in the right top corner). After particle overlaps are removed, the simulation switches to molecular dynamics (symbol [MD/Bussi CSVR], see Bussi CSVR thermostat). If a problem is encountered (particle overlap which would lead to integration failure), several Metropolis Monte Carlo steps are automatically performed again.
• Monte Carlo NVE (Creutz)
  The Monte Carlo (MC) algorithm samples the configurations of molecules using random moves. In the constant-energy Creutz version, a demon with a bag full of energy visits molecules and tries to push them randomly a bit (random displacement in a disk of diameter d). If the energy decreases, the move is accepted and the demon stores the released energy in the bag. If the energy increases, the move is performed only if there is enough energy in the demon's bag, otherwise the move is rejected and the simulation continues with the original unchanged configuration. This procedure is repeated many times (one MC sweep is when every particle is tried to move once).
  If there is (on average) a lot of energy in the bag, the configurations of the system have high energy, too, which means that the temperature is high. Mathematically, the demon's energy is distributed according to the Boltzmann distribution, prob(Ebag) = const exp(-Ebag/kT), hence the temperature is given by averaged bag energy as

      Tbag = ⟨Ebag⟩/k.

  In the right panel you will see this temperature averaged in blocks of one or more sweeps. It fluctuates in the course of the simulation.
  By default, the trial displacement is determined automatically to reach the acceptance ratio of 0.3. For productive runs, this should be turned off by button [▯ set MC move].
• Monte Carlo NVT (Metropolis)
  The more useful Metropolis Monte Carlo algorithm works at constant temperature T, which is a parameter of the simulation. Here the demon visits a termostat before every step and gets energy corresponding to its temperature. Mathematically

      Ebag = −kT ln(rnd),

  where rnd is a random number drawn uniformly from interval (0,1). Once again, high-energy configurations are more probable at high temperatures, low-energy configurations at low temperatures. The demon's temperature Tbag converges to the thermostat temperature T and then fluctuates around it (on average, the energy the demon returns to the thermostat equals the energy he gets back).
• Monte Carlo NPT (Metropolis)
  The volume of the system randomly fluctuates in the isothermal-isobaric (NPT) ensemble. Again, the change of volume is either accepted or rejected to maintain constant pressure on average. Note that the simulation box is always drawn at the same size and the sizes of molecules change instead.
  For experts: We use the volume measure μ=1/V, to be compatible with MTK, see J. Janek and J. Kolafa: J. Chem. Phys. 160, 184111 (2024).
  If the barostat pressure is too low, the volume becomes too large. In this case, an error is reported and the simulation switches to NVT. You should change (increase) P from the cmd: input field. Alternatively, stop simulation by the [▯ run] button, then switch to NPT, then use the P-slider, then restart simulation by the [▯ run] button.
  The trial volume change is not available via a slider; use cmd: dV=VALUE or automatic setup by button [▯ set MC move].
• Molecular Dynamics (NVE)
  Molecular dynamics (MD) calculates trajectories of molecules in time by numerical integration of the Newtonian equations of motion. Molecules thus move "as in the real world".
  Mathematically, we use the leap-frog algorithm given by formulas:

      v(t+dt/2) := v(t−dt/2) + dt*f(t)/h,
      r(t+dt) := r(t) + dt*v(t+dt/2),

  where v(t) is the velocity (velocities of all particles) and r(t) position(s) at time t, respectively, and dt is the timestep. Forces f(t) are calculated from positions r(t).
  The total (potential + kinetic) energy is conserved (subject to small method errors). Temperature is calculated from the kinetic temperature by formula:

      Tkin = 2 Ekin / k Nf,

  where the leap-frog kinetic energy is

      Ekin = ∑ m [v(t−dt/2)²+v(t+dt/2)²]/4,

  the sum is over all particles, both the Boltzmann constant k and particle mass m are here equal 1, and Nf is the number of degrees of freedom (Nf=2N in a box; because of momentum conservation in the NVE MD, Berendsen and Nosé-Hoover thermostats, Nf=2N−2 in the periodic boundary conditions and Nf=2N−1 in the slit).
  The timestep dt is determined from temperature and thermostat type and then kept constant. In NVE and when temperature changes abruptly, the integrator may fail. A fixed timestep can be set on start by option -Pdt=VALUE, or from running simolant from field cmd:, variable dt. The automatic timestep is set by dt=0.
  Before the MD simulation starts, new velocities are drawn from the Maxwell–Boltzmann distribution. This applies to all MD methods.
• Molecular Dynamics NVT (Berendsen)
  In order to maintain constant temperature, a thermostat has to be used. The simplest one is the Berendsen (friction) thermostat. It works by rescaling all velocities by a factor so that Tkin approaches slowly the defined temperature T. This thermostat (with long enough τ) does not affect transport phenomena as diffusivity and is good for fast cooling or heating. The coupling of the thermostat with the system is given by the coupling (correlation) time τ. Too low values may cause problems for some systems as, e.g., the infamous "flying ice cube" (create a crystal in periodic boundary conditions, use τ=0.1 and wait...). It does not give the true canonical ensemble (the same distribution of energies as if the system were at contact with a thermostat) which is a disadvantage in some application.
• Molecular Dynamics NVT (Nosé-Hoover)
  This thermostat is based on a dynamic variable (something like a "piston of entropy") interacting with the system. The coupling (correlation) time τ must be carefully set (less than 0.3 ps is recommended): Too long values lead to long oscillations in temperature and poor convergence (the "piston" has large "inertia mass"); too short values require short timesteps (the integration may fail). If everything is set properly, this thermostat is probably the best.
  The method can be implemented in the Verlet-like methods by several algorithms. Here we use the TRVP (Time-Reversible Velocity Predictor) with k=1, see J. Kolafa and M. Lísal: J. Chem. Theory Comput. 7, 3596 (2011).
• Molecular Dynamics NVT (Andersen)
  Constant temperature is maintained by thermalizing the molecules randomly: A molecule forgets its old velocity and is launched at random direction with random velocity corresponding (on average) to temperature T (so-called Maxwell-Boltzmann distribution). Again, τ determines how often (on average) a molecule is visited by a random Maxwell-Boltzmann kicker. This type of thermostat can be interpreted as an action of random impacts of invisible gas molecules in which our molecules are immersed. The algorithm is subject to errors for short τ and long dt.
• Molecular Dynamics NVT (Maxwell-Boltzmann)
  As above except that all molecules are asigned Maxwell-Boltzmann random velocity corresponding (on average) to temperature T at once every τ.
• Molecular Dynamics NVT (Langevin)
  Small random forces ceaselessly disturb all molecules. Since this would heat the system, a friction is added to slow the molecules. If both forces are balanced, the system maintains given temperature T. A simple version of the algorithm is implemented so that the ensemble is not exactly canonical for short τ and long dt.
• Molecular Dynamics NVT (Bussi CSVR)
  Also called Canonical Sampling through Velocity Rescaling. Similar to the Berendsen thermostat except that the rescaling factor is amended by a random noise so that the resulting ensemble is canonical, see Bussi G. and Parrinello M., Comp. Phys. Comm. 179, 26 (2008).
• Molecular Dynamics NPT (Berendsen)
  Similar as the Berendsen thermostat. Pressure is maintained by box rescaling if the instantaneous virial pressure differs from parameter P. The analogous volume coupling (correlation) time τ is calculated from the temperature time constant τ and the bulk modulus calculated from a very approximate equation of state, cf. variable qtau available via command in window cmd:.
  Hint: use Show → Volume convergence profile with Measurement block slided left.
• Molecular Dynamics NPT (Martyna et al.)
  Similar to the Nosé-Hoover thermostat but with a box-size-related degree of freedom added. You may interpret it as a "piston", but attached to all particles so that it works also in the periodic boundary conditions. The time constant of this piston is qtau×τ, τ<0.3 and the default qtau=5 are recommended. The true isothermal-isobaric ensemble is generated; however, it may be difficult to adjust the parameters; e.g., for high-temperature gas or crystal. Also, it may easily crash if the initial configuration is far from equilibrium. Can be used with non-periodic boundary conditions, but the scaling is still uniform in all coordinates. Try to view the volume convergence profile and switch all three available barostats.
  The method is based on paper G. J. Martyna, D. J. Tobias, M. L. Klein: J. Chem. Phys. 101, 4177 (1994), for our implementation see J. Janek, J. Kolafa: J. Chem. Phys. 160, 184111 (2024).
• Vicsek with NVT (Langevin)
  The same as the Langevin thermostat with the velocity of each particle (“bird” or “fish”) attaining the weighted averaged velocity from all particles within cutoff c:
  
      v_i := ∑_j v_j*w_j / ∑_j w_j, where w_j=c²−r_j²
  
  Then, the velocity is renormalized to |v|=√T. A Langevin thermostat step follows.
  Only the periodic boundary conditions are supported (box and slit are changed to periodic).
  A user may change the potential and its parameters; particularly, the cutoff c (the same value for the potential and neighbor definition); it, however, must not be too long. Thermostat parameter τ controls the chaos. The temperature with penetrable disks and small a essentially scales time and does not change trajectories.

Boundary conditions

• 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.
  

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.

• Quantities
  As above and additional quantities are measured and shown as text:
  Pvir = pressure (calculated from the virial of force); this pressure is also used in barostats. This pressure for thermostats conserving the linear momenta in the periodic and slit boundary conditions include the “kinettic pressure correction”. Thus, the ideal gas in the limit of ρ→0 gives exact compresibility factor, Z=1. See the MACSIMUS manual for details and tests. The respective kinetic energy scaling factors are printed as qx,qy.
  Pxx, Pyy = diagonal components of the pressure tensor as calculated from the virial of force. This is the natural pressure in the periodic boundary conditions, but it is available also with walls; both values should match. (NB: the virial component of gravity is not included.)
  P(top wall), P(bottom wall) = pressure on walls as calculated from forces on walls.
  γ = gamma = line tension ("surface tension" in 2D):
      γ = Ly*(Pyy−Pxx)
  It makes sense only in the slab geometry.
  MSDx, MSDy = mean square displacement (MSD) measurement is started when [▯ record] is pressed. The MSDs in x and y directions are calculated separately. The formula is
      MSDx(t) = ∑_i [x_i(0)−x_i(t)]² / N
  and is applied for the last configuration in a block. The algorithm operates with scaled coordinates (thus, NPT is valid) and follows particle paths in the periodic boundary conditions. If the displacement from the previous position is greater than 0.45 L, error message is issued and the MSD calculation is turned off. The remedy is to decrease the block size and/or stride.
  The slope of the linear part of MSD.x(t) and MSD.y(t) is 2D in the periodic boundary conditions, where D is the diffusivity. This calculation is rather imprecise (error ~1/√N); therefore, several shorter independent measurements should be averaged. (Ideally, overlapping blocks are recommended, see J. Picálek, J. Kolafa: J. Mol. Liq. 134, 29–33 (2007)).
  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 mininum–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).
• Pressure convergence profile
  The convergence profile of virial pressure.
• Volume convergence profile
  The volume convergence profile (NPT only).
• Momentum (e.g., for Vicsek)
  The total momentum (or velocity because m=1) of the simulation box per particle.
• 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) using formula:

      (L/2π)*atan2(∑ sin(yi*2π/L),∑ cos(yi*2π/L)).

  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 fluctuates); 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; overshot 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 above for both x and y.
• 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))+π).

Tinker

Technical setup affecting performance.

• FPS=15
  Set the number of frames per second to 15. The real FPS may be changed by slider [simulation speed]. Also, if the simulation takes a longer time, the real FPS is slower.
• FPS=30
  Set the number of frames per second to 30.
• FPS=60
  Set the number of frames per second to 60. This is the default.
• FPS=120
  Set the number of frames per second to 120.
• Timeout 99% of simulation (use if shaky)
  If the simulation takes a longer time than is given by the FPS selected, some timeout is still required for widget control. Here the timeout is set to 99% of the simulation time. Use this option if the response of widgets is slow and shaky. See also option -T.
• Timeout 50% of simulation
  The timeout is set to 50% of the simulation time.
• Timeout 20% of simulation (faster)
  The timeout is set to 20% of the simulation time. Use this to speed up the calculations, perhaps at the expense of a less precise response.
• Timeout 7% of simulation (fastest)
  The timeout is set to 7% of the simulation time. Use this to speed up the calculations, at the expense of a less precise response.

Help

• Help
  This help. File simolant.html must be in the same directory as the exe file.
• About
  Copyright and acknowledgements.

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

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

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.

• 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 ensues fast heat flow.
• 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.
  In addition, button [g=0] sets zero gravity.
• ρ   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 scaled 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
  In addition, buttons [N+] and [N−] allow for adding and removing single particles.

Walls + shifts

• If walls are present: Clicking the respective top/bottom/left/right buttons will toggle repulsive (red) and attractive (green) walls. Button [invert] swaps repulsive ↔ attractive walls.
• For periodic boundary conditions: Clicking the respective shift left/right/up/down buttons will move the configuration. The opposite moves are not the same to allow for a finer move.

Expert

• record
  Turning the [▯ record] button on will start recording the measured values (of energies, pressure, etc.), including the convergence and density profiles if selected in the [include] pull-down menu.
  If the [▯ record] button is turned off, averages of the measured quantities are displayed (with error estimates). You may choose from:
  • clear (data lost) – measured data are cleared, recording is turned off.
  • save (...) and clear (also Enter) – write the shown data and more info to a file (default name = simolant.txt, can be changed from the menu File → Protocol name...). The file is overwritten if exists, however, any further invocation of save and clear within one session of SIMOLANT will append the results to the same file (unless the file is changed).
  • continue (also Esc) – no action (continue with recording).

  Note that averaged data (see the block slider) are recorded and that correlations between data are ignored in calculating the standard errors.

• CSV
  The protocol will also be printed to a CSV file (1 line per [▯ record]); note that dt and d share the same column.
  Both the simolant.txt and simolant.csv contain the same informatiom, but simolant.txt is verbose (human-readable) and simolant.csv contains just columns with a header.
• comma
  If selected, the recorded file will contain decimal commas, not dots. The same applies to the CSV file(s), where also the separator is changed to ';'. The cmd: field input accepts both comma and period. The on-screen output uses periods only.
• include:
  The following verbose information is generated:
  • Nothing – No additional information; i.e., only statistics of thermodynamic quantities is recorded.
  • Convergence profile – Convergence profile = time series of quantites; i.e., a line of selected block-averaged thermodynamic quantities is printed for every block.
  • Density profile – The active density/distribution function including the cumulative function is printed when recording is stopped. The density/distribution function is selected from menu Show.
  • Both – Both convergence profile and the selected density/distribution function.
  If CSV is highlighted, the respective numbered CSV files are printed. If CSV is deselected, the same information is printed to the output file (simolant.txt or as set in the File menu) with columns separated by Tabs.
• cmd:
  Allows for entering data in the form of numbers instead of sliders as well as queries. The format to assign a NUMBER to a VARIABLE is:
      VARIABLE=NUMBER
  To print the value of VARIABLE:
      VARIABLE
  Type Enter after entering the data. For instance:
      rho=0.01
      ρ
  The values can be outside the respective slider range; the slider will then show a minimum or maximum. A check is still made whether the variable entered is within some reasonable range, but not so strictly; the simulation may thus crash in some cases. Both decimal comma and dot are accepted (irrespective of the comma toggle); note that this does not apply for option -P. Available variables (input is not case-sensitive):

  variable	access	explanation
  a	(menu)	Parameter of the double well and penetrable disks potentials.
  b	(menu)	Parameter of the double well potential.
  bc	menu	Boundary conditions (0=box, 1=y-periodic, 2=periodic).
  block	slider	Measurements/block for statistics and graphs.
  circle	option -C	circle drawing method (0=fl_pie, 1=fl_circle, 2=composed of fl_line)
  c	(menu)	Potential cut-off.
  d	slider	MC displacement (if [▯ set MC moves] is off).
  dV		MC change of ln(V) (if [▯ set MC moves] is off).
  dt, h		Time step; 0=automatic.
  ff	menu	Force field (0=2DLJ,1=WCALJ,2=PD,3=DW).
  FPS	menu	Animation speed in Frames Per Second, see also Tinker.
  g	slider	Gravity.
  L		Box size (density ρ gets recalculated).
  measure	menu	REMOVED in V03/2025, use show instead.
  method	menu	Method (MC/MD, thermostat..), integer from 0 in the order defined in the menu.
  N	slider	Number of particles; also option -N.
  nbr		Atom-atom distance limit to define neighbors (for coloring)
  ncell		Number of cells per box side, or 1 for the direct pair sum, or 0 for automatic setup. See here for details.
  P	slider	Barostat pressure.
  qtau, qτ		Ratio τ(barostat)/τ(thermostat).
  rho, ρ	slider	Number density, ρ=N/V; L gets recalculated.
  show	menu	Show quantities/graphs, integer from 0 in the order defined in the menu.
  speed	slider	Show each stride-th configuration.
  T	slider	Thermostat temperature.
  tau, τ	slider	Thermostat time constant (MD).
  trace		Trace length for draw mode Traces and Art color mode (default=64).
  v		Droplet–droplet velocity for Prepare system → Two droplets.
  wall		Wall number density (smoothed to make the material of walls).

  All these variables can also be used in option -P at start.

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 Traces
  • Pixel – dot of one pixel
• draw mode:
  • Movie – Molecules in motion.
  • Traces – Trajectories are visible and gradually erased, the length of the trace can be controlled by variable trace from cmd:.
  • Lines – Trajectories are visible and permanent.
  • 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).
  • One – Only one molecule is black, the remaining ones are orange. Useful to illustrate Brownian motion.
  • 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).
  • Art – Cycle colors on the RGB wheel, best with draw mode Lines. The speed of change can be controlled by variable trace from cmd:
  • Keep – Previously selected colors are kept.
• reset
  Restore the default behavior (movie of black molecules). The convergence profiles are reset, too (if shown).
• set MC moves
  In MC simulations: The trial move size (and in NPT also dV) is adjusted automatically to reach the acceptance ratio of 0.3. Should be turned off before a productive run.
• run
  Stop the simulation temporarily. All other functions (parameter change etc.) work as usual and become active as soon as the simulation is resumed. Useful to change 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). The fourth line in the right panel contains the value of stride and the requested time between frames in ms.
  • Leftmost: very slow, every frame is shown with a low FPS (frames per second).
  • Left: every frame is shown with the selected FPS (see menu Tinker).
  • Center: every 3rd configuration is shown with the selected FPS. This is the default.
  • Right: every tenth configuration is shown with the selected FPS.
  • Rightmost: every tenth configuration is shown as fast as possible (no delays added).
  Note that for many particles, the calculation and showing may be slower than the selected FPS.
• 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. The selected stride (see slide simulation speed) still applied; therefore, stride*block=10*100 means that 1000 MD steps are performed to get one block summary data.

Command line options

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:

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

Known bugs and caveats

• There are still timing, thread-priority, and/or efficiency problems. The original dumb X11-only version was faster.
• 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

• 03/2025
Option nomeasure removed, menu renamed to Show, always Pvir etc. measured at end of sweep due to code and control simplicity, linked-cell list implemented (except Metropolis MC and measurements), some polishing of timing, buttons [N+] [g=0] [N] added, configuration shift, menu items reordered, new items added to Prepare system

• 02/2025
Vicsek model of active matter added
two droplets colored, velocity added (variable v), momentum added (for Vicsek)
colormode Art (rainbow) added, variable trace added (for Traces and Art)
• 01/2025 Bug fixed: After change of N while MC, the number of degrees of freedom was not updated leading to wrong results
  MSD also for nonperiodic b.c.
• 12/2024 More potentials added: WCA, double well, penetrable disks/spheres
  also velocities saved to .sim
  MSD added
  two droplets added
• 11/2024 Cutoff smoothing bug fixed
  number of degrees of freedom for MTK was fixed; now lim_ρ→0Z=1 also for small N
  nonstochastic NVT MD kinetic pressure corection added so that now lim_ρ→0Z=1 also for small N
  outputs cleaned; particularly for Z
• 10/2024 Missing header after change of protocol name fixed
  g↔P swapped (was wrong in .txt protocol)
  NPT volume measure μ=const replaced by μ=1/V in MC NPT, to be compatible with MTK
  
• 09/2024 help window and panel enlarged
  Tinker menu added
  timeouts based on the measured simulation times
  new button [invert] walls
  CSV output added
  Okabe and Ito colorblind-safe palette
  backups on save removed
• 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–now 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)