Some STP water for you thirsty gamer boys
Let's consider the following simulation, 100 water molecules at STP:
senpai --in examples/water.mol --out render.xyz --copy 100
The --copy
flag tells SENPAI to insert 100
copies of the substrate described in examples/water.mol
into the universe. But what is the process behind this?
SENPAI's way of loading the simulation
universe_init
is responsible for initialising the universe before the simulation. It itself calls universe_load
, which loads atoms from the input file into a special array, universe->ref_atom
. This array has universe->ref_atom_nb
elements. It contains each atom from the input file, and their information.
universe_init
will then call universe_populate
, a function that places copies of universe->ref_atom
at random locations within the universe.
The universe, as a set of matter, is nothing but an array, universe->atom
, of size (args->copies)*(universe->ref_atom_nb)
. The universe->atom
array therefore is a concatenation of several modified universe->ref_atom
arrays. The modification I was referrencing is the introduction of an offset, identical between all copies of the universe->ref_atom
array.
TL;DR: SENPAI places copies of the loaded substrate at random locations within space.
The issue here is "at random". SENPAI doesn't care if atoms are located at the same physical coordinates. As such, another step is required before starting the simulation: potential reduction.
SENPAI will apply algebraic transformations to the coordinates of each atom within space, so as to lower the potential energy. This process can be tuned through headers/config.h
:
/* PRE-SIMULATION POTENTIAL ENERGY REDUCTION
*
* Before starting a simulation, SENPAI will use a two-stage algorithm to reduce
* the potential energy of the system.
*
* STAGE 1: COARSE (brute force)
* The first stage consists of iterating through the particles and relocating
* them to random offsets, discarding the relocation should the total potential
* increase. The magnitude of the relocation offset gets lowered after a set
* number of attempts. The first stage ends when the potential reaches a
* defined multiple of the target potential.
* UNIVERSE_REDUCEPOT_COARSE_STEP_MAGNITUDE: start magnitude of the relocation
* UNIVERSE_REDUCEPOT_COARSE_MAX_ATTEMPTS: Max attempts before lowering the
* magnitude.
* UNIVERSE_REDUCEPOT_COARSE_MAGNITUDE_MULTIPLIER: Reduce the magnitude by
* multiplying it.
* UNIVERSE_REDUCEPOT_COARSE_THRESHOLD: Target (multiple of the pre-reduction
potential)
*
* STAGE 2: FINE (fine)
* The second stage consists of tuning the coordinates of each atom so as to
* lower the total potential energy. Gradient descent is used here, since the
* analytical gradient is easily computable from the force (F = -nabla*U).
* UNIVERSE_REDUCEPOT_FINE_MAX_STEP: Maximum step an atom can take
* UNIVERSE_REDUCEPOT_FINE_TIMESTEP: Used to compute the step (dx=(dt^2)/2)
*/
#define UNIVERSE_REDUCEPOT_COARSE_STEP_MAGNITUDE (double)1E-9
#define UNIVERSE_REDUCEPOT_COARSE_MAX_ATTEMPTS (size_t)1E2
#define UNIVERSE_REDUCEPOT_COARSE_MAGNITUDE_MULTIPLIER (double)1E-1
#define UNIVERSE_REDUCEPOT_COARSE_THRESHOLD (double)2E0
#define UNIVERSE_REDUCEPOT_FINE_MAX_STEP (double)1E-10
#define UNIVERSE_REDUCEPOT_FINE_TIMESTEP (double)1E-15
It is mathematically impossible to find the absolute minima of the potential function. As such, SENPAI has to rely on a randomization stage AND a gradient descent stage to properly reduce the potential.
However, this method is extremely computationally intensive. So much that I considered moving it out of SENPAI into another project that would spit out a .mol
file that SENPAI could load and not care about. But for now, it shall stay here.
The issue being described here is the following: this method is inefficient, and sometimes SENPAI will spend more time reducing the potential than actually simulating the universe. This has to change, any ideas?