<jupyter_start><jupyter_text>Automating `DL_MONTE` equation of state calculations with `simtask`  Introduction `simtask` is a Python package for automating complex tasks involving multiple simulations. Here we introduce the package and provide an example illustrating how it could be used. The example involves calculating the energy and density of the Lennard-Jones fluid using grand-canonical Monte Carlo, via the program `DL_MONTE`, for various temperatures *automatically*. The user specifies which temperatures are to be examined, links `simtask` to the relevant `DL_MONTE` executable, and then sets the `simtask` machinery in motion. The result is that plots of energy and density vs. temperature, along with associated statistical uncertainties and taking into account the equilibration time of the simulation, are generated without any input from the user.  How does `simtask` work? We now briefly describe the key features of `simtask` with the aim of elucidating how it works. The `simtask.task` module is the heart of the package. This contains the two key Python classes, `Task` and `TaskInterface`.`Task` - This is the base class corresponding to a simulation 'task', i.e. something you want to do via simulation. E.g. in the example below the 'task' is to evaluate the energy and density via grand-canonical Monte Carlo simulation. `Task` is an abstract class in the sense that it does nothing. It is subclasses of `Task` which correspond to different possible tasks, and contain code for performing these tasks. For instance the `Measurement` class is a subclass of `Task` which calculates one or more observables at a given set of thermodynamic conditions. This is used in our example below. The code for `Task` subclasses for performing simulation tasks is *general*; it is not specific to a simulation program, e.g. `DL_MONTE`. This is achieved by `Task` subclasses drawing only upon functions within the `TaskInterface` class for simulation-program-specific code. These functions constitute common things one may wish to do during a task (e.g. copy input files to a new directory, run a simulation in that directory, extract the energy from the output files of a directory), common things which can be combinedf into complex tasks. To enable `simtask` to be used with a given simulation program, one must write a subclass of `TaskInterface` which overloads the functions in `TaskInterface` with functions which are appropriate for that simulation code. For instance, `DLMonteInterface` is the subclass of `TaskInterface` which contains functions appropriate for the simulation program `DL_MONTE`. To perform a task in conjunction with `DL_MONTE`, the relevant `Task` subclass, e.g. `Measurement`, simply 'links' to the `DLMonteInterface`, which results in `Measurement` calling the `DL_MONTE`-specific functions for the common things during its invocation.This set-up separates the code pertaining to the task one is interested in performing from the code specific to the simulation program one will use to perform that task. `Task` subclasses will in principle work for *any* simulation program which has an existing `TaskInterface`. Conversely, simulation programs which have an existing `TaskInterface` will work for any future `Task` subclasses.  Preliminaries Installing and setting up `simtask``simtask` can be obtained from ...  Once `simtask` has been installed, the Python environment must be configured so that `simtask` can be 'seen' by Python scripts and modules which wish to utilise the package. There are a number of ways in which this can be achieved. One is to manually change the list of directories which Python searches for packages and modules *within* Python. We will employ this approach here.<jupyter_code>import sys
import os

# You will need to set SIMTASK_HOME appropriately for your local system. 
# Note that SIMTASK_HOME should be the directory ABOVE the directory 'simtask' 
# containing the source code, i.e. the files __init__.py, analysis.py, task.py, etc.
# Rather confusingly this directory might also be called 'simtask', as is the case
# here: '/home/tom/Work/Code_workspace/simtask/simtask' contains the source code.
SIMTASK_HOME = "/home/tom/Work/Code_workspace/simtask"

# Set the search path to the SIMTASK_HOME directory
sys.path.append( SIMTASK_HOME )<jupyter_output><empty_output><jupyter_text>The alternative approach is to alter the shell variable `PYTHONPATH`, which contains a list of directories in which to look for Python packages to be imported. To use `simtask`, it is prudent to set `PYTHONPATH` to the directory containing the package, i.e. the directory stored in the `SIMTASK_HOME` variable in the above code snippet. The following shell command (or whatever equivalent is appropriate for your local system) does this:```bash> SIMTASK_HOME="/home/tom/Work/Code_workspace/simtask"bash> export PYTHONPATH=$PYTHONPATH:$SIMTASK_HOME``` Once `simtask` has been obtained and the Python environment has been configured appropriately, simtask can then be imported from within a Python script, and hence is ready to use.<jupyter_code>import simtask<jupyter_output><empty_output><jupyter_text>Dependencies for the example: `DL_MONTE` and `htk` Below we use `simtask` to perform Monte Carlo simulations by interfacing simtask with the general-purpose Monte Carlo program `DL_MONTE`. Therefore `DL_MONTE` is a prerequisite for the examples below. However `simtask` also requires a further Python package, `htk`, which is distributed alongslide `DL_MONTE`, in order to interface with `DL_MONTE`. `htk` contains a Python interface for launching `DL_MONTE` simulations and manipulating `DL_MONTE` input and output files, functionality which is drawn upon by `simtask`. See https://www.ccp5.ac.uk/DL_MONTE and links therein for details regarding how to obtain `DL_MONTE` and `htk`.  Installing and setting up the dependencies `DL_MONTE` can be installed as described in the manual and the website where it is distributed.Similarly to as described above for `simtask`, the Python search directory for modules and packages must be set up so that `htk` can be 'seen'. As we did before for `simtask`, we elect to modify the search path within Python, though of course one could achieve the same result by modifying the shell variable `PYTHONPATH` as described above.<jupyter_code># You will need to set HTK_HOME appropriately for your local system. 
# Note that HTK_HOME should be the directory ABOVE the directory 'htk' 
# containing the package source code, i.e. the files __init__.py, histogram.py, 
# task.py, etc. Rather confusingly this directory might also be called 'htk', 
# as is the case here: '/home/tom/Work/Code_workspace/DL_MONTE_2_dev/dlmonte2/htk/htk' 
# contains the source code.
HTK_HOME = "/home/tom/Work/Code_workspace/DL_MONTE_2_dev/dlmonte2/htk"

# Set the search path to the 'HTK_HOME directory
sys.path.append( HTK_HOME )<jupyter_output><empty_output><jupyter_text>Example 1: Energy and density of Lennard-Jones fluid at a given $\mu$ and $T$ We are now ready to begin the examples. The ultimate aim is to use `simtask` to calculate the energy and density, with associated uncertainties, of the Lennard-Jones fluid for a range of temperatures $T$ at a given chemical potential $\mu$. However we first consider the problem of evaluating these physical quantities at a *single* temperature; in Example 2 later we consider multiple temperatures. More specifically, we will use `simtask.measurement` to perform back-to-back grand-canonical Monte Carlo simulations with `DL_MONTE` until one of two criteria are met: 1) the energy and density (in fact, the mean number of particles in the system) are evaluated to a specified *precision*; 2) a threshold number of back-to-back simulations have been performed. Note that by 'back-to-back' we mean that the next simulation begins in the configuration where the previous simulation ended.To perform a standard grand-canonical `DL_MONTE` simulation three `DL_MONTE` input files are required, `CONTROL`, `CONFIG` and `FIELD`. The relevant files can be found in the directory housing this notebook. A fuller account of the significance of the contents of the input files can be found elsewhere, e.g. the `DL_MONTE` user manual. We discuss the contents of these files only briefly here. Moreover we do not provide here an account of grand-canonical Monte Carlo or Monte Carlo simulation methodology, further details of which can be found elsewhere. We assume the reader has a basic understanding of grand-canonical Monte Carlo and how to run grand-canonical simulations in `DL_MONTE`.  Simulation detailsThe input files CONTROL, CONFIG and FIELD correspond to a grand-canonical Monte Carlo simulation of the Lennard-Jones fluid as follows. The reduced temperature is $T^*=1.55112$ and the thermodynamic activity is $z=0.04$, where $T^*\equiv k_BT/\epsilon$, $z\equiv\exp(\mu/(k_BT))$, and $\epsilon$ denotes the well depth of the Lennard-Jones potential. Note that the $T^*$ considered is above the critical temeprature (which is $T_c^*=1.326$). The system is a cube with sides of length $10\sigma$, where $\sigma$ is the interparticle separation at which the Lennard-Jones is zero. The potential is truncated at a distance of $2.5\sigma$, and no tail corrections are applied. The Monte Carlo move-set is comprised only of insertion and deletion moves; no translation moves were applied. Finally, the length of the simulation is 120,000 moves. (Note however that back-to-back simulations will be performed; hence considerably more moves than 120,000 will be used in total).  Configuring the output of `simtask` We will now describe the `simtask` Python commands required to perform the task at hand. First of all we should import the relevant `simtask` modules. We also import the `logging` module (a standard library module), since `simtask` uses this module to control its output. The aim of `simtask` is to perform a complicated task involving multiple simulations. During such a task, `simtask` outputs information, via the `logging` module, regarding the progress of the task for the benefit of the user. `logging` enables the 'level' of output to be controlled. For debugging one would typically desire more information to be output; for very long tasks one would typically want a more succint description of the progress of the task. This is all covered by the below code snippet.<jupyter_code>import logging

import simtask.dlmonteinterface as interface
import simtask.measurement as measurement            # IS ALL OF THIS REALLY NEEDED?
import simtask.analysis as analysis
import simtask.task as task

# Set up the logger, which determines the nature of information output by the 
# machinery in the 'simtask' package. The line below results in logging 
# information being output to stdout

handler = logging.StreamHandler()

# Below we set up the 'logger' in simtask.measurement, to determine the nature 
# of output when the 'measurement' module is used (see below). Replacing 
# 'logging.INFO' with 'logging.DEBUG' below results in more information being 
# output by the logger. Using 'logging.WARNING' results in less information 
# being output: only 'warnings'

measurement.logger.setLevel(logging.INFO)
measurement.logger.addHandler(handler)<jupyter_output><empty_output><jupyter_text>Linking to the `DL_MONTE` executable We now set up the relevant `TaskInterface` object, which tells the low-level machinery in `simtask` which simulation program will be used to perform the simulations, and how to perform various jobs specific to that program, e.g. extracting the energy from output files created by the program. In this case we use `DL_MONTE` to perform our simulations; thus the `TaskInterface` class we will use is `DLMonteInterface` (`DLMonteInterface` is a subclass of `TaskInterface`; see above). The code below creates a new `DLMonteInterface` object, which is linked to the `DL_MONTE` relevant executable. Later this object will be 'linked' to the relevant `Task` object before the task is performed.<jupyter_code># The line below creates a DL_MONTE-specific interface. Note that the interface 
# must know the location of the DL_MONTE executable - which is specified as the 
# argument to the DLMonteInterface constructor.
# You will have to adapt the location of the executable to suit your local system.

interface = interface.DLMonteInterface("/home/tom/Work/Code_workspace/DL_MONTE_2_dev_2/dlmonte2/bin/DLMONTE-SRL.X")<jupyter_output><empty_output><jupyter_text>Setting up the `Measurement` object We will use the `Measurement` class to automate the task of calculating the energy and number of molecules for various given chemical potentials.  We will set up a `Measurement` object, `m`, and then begin the task with the command ```pythonm.run()```  Task observablesThe `Measurement` constructor requires a list of `Observable` objects (see below), which represent the physical quantities to be tracked during the task. We are interested in the energy and number of particles. Below we create a list `observables` accordingly.<jupyter_code>energy_obs = task.Observable( ("energy",) )
nmol_obs = task.Observable( ("nmol",0) )
observables = [ energy_obs, nmol_obs ]<jupyter_output><empty_output><jupyter_text>Note that the nature of `Observable` objects may vary between simulation codes. For `DL_MONTE` only observables corresponding to variables output periodically in the `YAMLDATA` file are currently supported. For a variable 'foo' specified in the `YAMLDATA` file the corresponding `Observable` object is returned by the command ```pythontask.Observable( ("foo",) )```Note the essential comma after `"foo"`! For a variable in `YAMLDATA` which is an array (e.g., 'nmol'), the observable corresponding to the `n`th element in the array is returned by the command```pythontask.Observable( ("foo",n-1) )```See the code snippet above: `energy_obs` corresponds to the 'energy' variable in `YAMLDATA`, and `nmol_obs` corresponds to the 1st element in the 'nmol' array in `YAMLDATA` (which in fact is the number of molecules belonging to the 1st molecular species).  Termination criteriaThe `Measurement` class will perform back-to-back simulations until one of a number of customisable criteria have been met. The criteria are determined by the arguments passed to the constructor of the `Measurement` object. For now we apply two criteria. The first is that the energy is evaluated to a precision of 0.2 (energy units), i.e. the back-to-back simulations are terminated if the statistical uncertainty in the energy is less than this precision. The second criterion is that the number of back-to-back simulations does not exceed 20. The below code creates a `Measurement` object with the aforementioned criteria. Moreover we set the name of the directory in which the simulations are to be performed as `fixedprecision`.<jupyter_code># A dictionary is used to specify the threshold precisions to apply to the observables. Here we only
# specify a threshold precision for the energy observable: the simulations will terminate when the energy is
# determined to an uncertainty less than 0.2 (energy units). Note that there is no entry in the dictionary
# corresponding to the 'nmol' observable. The significance of this is that there is no threshold precision
# for this observable.

precisions = { energy_obs : 0.2 }

# Set up the relevant Measurement object - which will actually perform the simulations and data analysis.
# Note that all the simulations and output files pertaining to analysis will be created in the directory
# 'fixedprecision' (via the 'outputdir' argument), and that we will impose a threshold precision on the
# energy as implied by the dictionary 'precisions' described above. Furthermore, we specify that no more
# than 20 simulations will ever be performed (via the 'maxsims' argument)

m = measurement.Measurement(interface, observables, maxsims=20, precisions=precisions, outputdir='fixedprecision')<jupyter_output><empty_output><jupyter_text>Running the taskAs mentioned above, the task is invoked via the command```pythonm.run()```Note that when this is called, information is generated and output to stdout at a level set earlier in the discussion pertaining to the `logging` module.<jupyter_code># Remove the old directory 'fixedprecision', if it exists, before running the task:
# the task will not run if the directory already exists; it will not be overwritten.

import shutil

if os.path.exists('fixedprecision'):
    shutil.rmtree('fixedprecision')


# Run the task

m.run()<jupyter_output>Beginning measurement task...


Beginning simulations...


On simulation 1 of max. 20...

Setting up directory 'fixedprecision/sim_1' for simulation...
Running simulation in 'fixedprecision/sim_1'...
Simulation complete

Extracting data for observable 'energy'...
Extraction complete
Analysing time series for observable 'energy'...
Bypassing equilibration test: equilibration assumed from outset
System is deemed equilibrated
  equilibration time = 0
Deducing blocks to use for block averaging...
Results of post-equilibration time series analysis for observable 'energy':
  number of data points = 1200
  statistical inefficiency = 4.11268841227
  correlation time = 2.01515972449
  calculated block size = 21
Calculating block averages...
  number of blocks found = 57
Analysing block averages for observable 'energy' (for mean and standard error)...
Results of analysis of block averages:
  statistical inefficiency of block averages = 1.0
  mean from block averages = -16.2003129454
  std. erro[...]<jupyter_text>Analysing the outputDuring the task the simulations are performed in subdirectories of the directory `fixedprecision` earmarked as the workspace for the task. The first simulation is performed in subdirectory `sim_1`, the second in `sim_2`, etc. The task involved 6 simulations before the task was terminated on account of the threshold energy being reached. Thus within `fixedprecision` there are thus 6 subdirectories, one corresponding to each simulation. Moreover there are two files, `energy_converge.dat` and `nmol_0_converge.dat`. These contain plots of the observables and their uncertainties obtained after each simulation: the first column in each file corresponds to the simulation number, the second to the mean value of the observable, and the third to the uncertainty. Note that data is pooled from all simulations up to that point in order to come up with the mean and uncertainty. Furthermore the mean and uncertainty are calculated only using post-equilibration data: an algorithm is used to detect when the system has equilibrated, and data obtained before that is ignored in calculating the mean and uncertainty. Details of how the mean and uncertainty are calculated, and how equilibration is detected, is beyond the scope of this document; see the source code for details. PLOT THE ENERGY AND NMOL VS simulation number...  Example 2: Energy and density of Lennard-Jones fluid at a given $\mu$ vs. $T$ Example 1 demonstrated how `simtask` could be used to calculate the value and uncertainty of an observable automatically at a given set of thermodynamic conditions: the energy and density was calculated at a specified chemical potential and temperature. We will now go further, automating the calculation of the energy and density at *many* temperatures.The relevant `simtask` class is `MeasurementSweep`, which performs sets of simulations similar to that performed by `Measurement` at a list of different *control parameters*.Similarly to `Measurement`, to invoke the task corresponding to a `MeasurementSweep` object `sweep` the command is```pythonsweep.run()```However, first one must create the object `sweep`, parameterised appropriately for the task at hand.  Setting up the Measurement objectFirst of all we must make a *template* `Measurement` object which will determine the nature of the calculations of energy and density at each temperature: the calculations at each temperature will track the same observables, and have the same termination criteria, as that of the template `Measurement` object. We make a template with the same observables and `TaskInterface` object (pertaining to `DL_MONTE`) as before. However this time we will use slightly different termination criteria. Instead of using the precision of the energy as a termination criterion, we will impose a maximum time of 30s to the task (i.e. if 30s has elapsed, then no further simulations will be performed beyond the current one). Thus the template `Measurement` object is created as follows.<jupyter_code>measurement_template = measurement.Measurement(interface, observables, maxsims=20, maxtime=30)<jupyter_output><empty_output><jupyter_text>Note that a maximum number of simulations of 20 is still a termination criterion, as it was in Example 1.The next step is to create a list of control parameters - in this case temperature - to examine.<jupyter_code>temperatures = [8000,9000,10000,11000,12000,13000,14000,15000,16000,17000]<jupyter_output><empty_output><jupyter_text>Finally we create the `MeasurementSweep` object `sweep` which will actually perform all the simulations and data analysis. The following command achieves this. Note that the list of control parameters is passed to the argument `paramvalues`. Note also that a string `"temperature"` has been passed to the argument `param`, which is the name of the control parameter. The result is that the `DLMonteInterface` object `interface`, which will create and modify `DL_MONTE` input files during the task, 'knows' to treat the temperature as a control parameter during the task. More specifically, in order to modify the input files appropriately for the next control parameter to be considered, it will search for the keyword `"temperature"` in the `DL_MONTE` input file `CONTROL`, and then modify the value of the temperature next to that keyword. Finally, we should mention that the name of the directory in which all simulations are to be performed is set to `fixedtimesweep`.<jupyter_code>sweep = measurement.MeasurementSweep(param="temperature", paramvalues=temperatures,
                                     measurement_template=measurement_template, outputdir="fixedtimesweep")<jupyter_output><empty_output><jupyter_text>Running the taskAs mentioned above, the task is invoked via the command```pythonsweep.run()```Note that when this is called, information is generated and output to stdout at a level set earlier in the discussion pertaining to the `logging` module.Since we are considering 10 temperatures, and have imposed a time limit of 30 s for all simulations at a given temperature, we expect that the whole task should take about 300s = 5min to perform.<jupyter_code># Remove the old directory 'fixedprecision', if it exists, before running the task:
# the task will not run if the directory already exists; it will not be overwritten.

if os.path.exists('fixedtimesweep'):
    shutil.rmtree('fixedtimesweep')

# Run the task
    
sweep.run()<jupyter_output>Beginning measurement sweep task...


Beginning measurement for control parameter value 8000...

Setting up directory 'fixedtimesweep/param_8000' for control parameter value 8000...
Running measurement for control parameter 8000...

Beginning measurement task...


Beginning simulations...


On simulation 1 of max. 20...

Setting up directory 'fixedtimesweep/param_8000/sim_1' for simulation...
Running simulation in 'fixedtimesweep/param_8000/sim_1'...
Simulation complete

Extracting data for observable 'energy'...
Extraction complete
Analysing time series for observable 'energy'...
Bypassing equilibration test: equilibration assumed from outset
System is deemed equilibrated
  equilibration time = 0
Deducing blocks to use for block averaging...
Results of post-equilibration time series analysis for observable 'energy':
  number of data points = 1200
  statistical inefficiency = 248.252656236
  correlation time = 124.125656756
  calculated block size = 1242
Calculating block averages...
[...]
{'fraction_non_alphanumeric': 0.05893117072156374, 'fraction_numerical': 0.010189228529839884, 'mean_word_length': 5.134183673469388, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 49, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '128785', 'license': 'MIT', 'n_tokens_mistral': 6217, 'n_tokens_neox': 5801, 'n_words': 3272, 'path': 'docs/tutorials/dlmonte_gcmc.ipynb', 'repo_name': 'tomlunderwood/simtask', 'size': '24045'}