SIMOLANT 12/2024

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
• Menus
  • File
  • Prepare system
  • Method
  • Boundary conditions
  • Show
  • Tinker
  • Help
• Right Panel
• Symbols
• Variables
• 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.

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)
  • Dependence of the line tension on temperature and the total line tension (surface enthalpy)
  • Kelvin equation
  • 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. Any value of c≥2.232 selects this potential.

For c<L/2 (L is the square box side), the nearest-image algorithm is used, otherwise a sum over all periodic images is calculated.

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 triuncation), 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>(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 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...). 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.

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). Any c ∈ (1,2) is transformed to c = 2^(1/4) = 1.1892071 and selects this potential.
• Penetrable disks and ideal gas
  The PD (penetrable spheres, here penetrable disks) is a simple model of a polymer coil in a good solvent. It is selected by c=2:
  
      u(r) = a(r²−4)²     for r < 2
      u(r) = 0                for r ≥ 2
  
  Normally, a=1/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 tesselation:
  
      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 selected so that ∫_1^∞ u(r)² dr = 512/1155 is the same as for the 2D Lennard-Jones. This potential is selected by c<−1. Parameters a and b are the nodes (zero values of the potential), it should hold 1 < a < b < |c| Note that no check is performed of this condition.

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) and 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, the extension .txt is changed to .csv.
• Export force field
  

  Prints a file (tab-separated TXT or CSV, according to 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.

Force field

• Lennard-Jones c=4 (recommended)
  Standard 2D Lennard-Jones potential smoothly cut off at c=4. Parameter c cn 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=0.5,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=−0.5,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 (−∞%lt;c−1 selects the double-well potential), it must hold 1 < a < b < |c|.

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

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

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.
  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 and measure

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.

• Minimum
  Measuring most quantities is turned off, shown are only the potential energy and temperature/kinetic energy (if applicable) and basic method parameters. If ▯ record is requested, Quantities (full measurements) are selected.
• 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. This pressure for thermostats conserving the linear momenta in the periodic and slit boundary conditions include the “kinettip 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.
  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). 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) = ∑ [r(0).x−r(t).x]² / 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).
• 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) 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. This is an approximate value which is based on measuring simulation and drawing time and calculating the remaining timeout. If the simulation takes a longer time, the FPS is smaller. Note that the speed is also affected by slider Simulation speed. It is void (maximum speed) if this slider is at the rightmost position. See also option -F.
• 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.
Symbol 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.

Sliders panel

Sliders can be moved by a mouse or (when the mouse pointer is over) finely by arrows.

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

Walls

Expert

• record
  Turning the [▯ record] button on will start recording the measured values (of energies, pressure, etc.).
  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.
• ,
  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 ';'.
• include:
  The following verbose information is generated:
  • Nothing – No additional information; i.e., only statistics of thermodynamic quantities is recorded.
  • Convergence prof. – 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. The format is
      VARIABLE=NUMBER
  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. It is still checked 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). 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.
  g	slider	Gravity.
  L		Box size (density ρ gets recalculated).
  measure	menu	Show and measure, integer from 0 in the order defined in the menu.
  method	menu	Method (MC/MD, thermostat..), integer from 0 in the order defined in the menu.
  N	slider	Number of particles.
  nbr	Atom-atom distance limit to define neighbors (for coloring)
  P	slider	Barostat pressure.
  qtau, qτ		Ratio τ(barostat)/τ(thermostat).
  rho, ρ	slider	Number density, ρ=N/V; L gets rcalculated.
  show		The same as “measure”.
  stride	slider	Show each stride-th configuration.
  T	slider	Thermostat temperature.
  tau, τ	slider	Thermostat time constant (MD).
  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.
  • 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).
  • 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 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

• 12/2024 More potentials added: WCA, double well, penetrable disks/spheres
  also velocities saved to .sim
  MSD added
  tho drops 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)